Roles of Carrier Doping, Band Gap, and Electron Relaxation Time in the Boltzmann Transport Calculations of a Semiconductor’s Thermoelectric Properties

Abstract

Although there is a growing demand for first-principles predictions of the thermoelectric properties of materials, the contribution of various errors in Boltzmann transport calculations is not negligible. We conducted a typical first-principles calculation and a Boltzmann transport analysis on a typical semiconductor (Si) at various temperatures T while varying the band gap ε_g, electron relaxation time τ_el, and phonon thermal conductivity κ_ph to demonstrate how the calculated thermoelectric properties, which are functions of the carrier doping level, are affected by these parameters. Bipolar conduction drastically decreased zT via a degradation of the Seebeck coefficient S and an increase in the effective Lorenz factor L_eff, indicating the importance of a wide enough ε_g (several multiples of k_BT or higher) for high zT. Thus, the underestimation of ε_g, which frequently happens in first-principles calculations, could induce large errors in calculations for narrow-gap semiconductors. The calculation of the electron thermal conductivity without Peltier thermal conductivity was found to limit the zT of typical semiconductors to below 1. A small value of κ_ph/τ_el, where κ_ph/τ_el is the degree to which a material is a phonon-glass electron-crystal, was necessary to achieve a high zT. Fitting the calculations with experimental thermoelectric properties showed that τ_el can vary by an order of magnitude from 10⁻¹⁵ to 10⁻¹⁴ s, depending on both T and the samples. This indicates that the use of a fixed relaxation time is inappropriate for thermoelectric materials.

1. Introduction

Discoveries of new environmentally-friendly thermoelectric materials for solid state power generation and Peltier cooling devices are highly anticipated.¹⁾ Although thermoelectric materials have been intensively sought after for over half a century, many potentially thermoelectric compounds remain unstudied or understudied. Experiments tend to underestimate the potential of compounds for use as thermoelectric materials, since the optimum doping level and the optimum method for introducing phonon scattering are rarely achieved during the first set of the experiments.

The best-known and the most popular benchmark for the efficiency of a thermoelectric material is the dimensionless thermoelectric figure of merit   

\begin{equation} zT = \frac{S^{2}\sigma T}{\kappa} = \frac{S^{2}\sigma T}{\kappa_{\text{el}} + \kappa_{\text{ph}}}, \end{equation} (1)

where S is the Seebeck coefficient, σ is the electrical conductivity, κ_el is the electron thermal conductivity, κ_ph is the phonon thermal conductivity, and T is the temperature. Thermoelectric conversion efficiencies for power generation and Peltier cooling increase monotonically with zT and reach the Carnot efficiency at zT → ∞. While the criterion for applications has long been zT > 1, a number of bulk materials with zT > 2 have recently been reported. The numerator of eq. (1) is referred to as the power factor P ≡ S²σ, which should be high, especially for Peltier cooling applications.

Equation (1) indicates that a high S, high σ, and low κ are needed to achieve a high zT. Of these parameters, S, σ, and κ_el cannot be controlled independently because of their strong dependence on the carrier doping level n ≡ n_e − n_h, where n_e and n_h are the electron and hole concentrations. In doped semiconductors, a decrease in the carrier density increases S and decreases κ_el but also decreases σ. As a result, zT is low both in low-n and high-n limits, and there is an “optimum n” (usually between 10¹⁹ and 10²¹ cm⁻³) where zT is maximized.¹⁾ Comparisons between different materials usually rely on the maximum reported zT values. However, the maximum zT of a material can only be achieved after developing effective methods to (1) dope carriers up to the optimum n and (2) to decrease κ_ph while avoiding a decrease in the electron relaxation time τ_el. Therefore, researchers have to find the correct impurity elements or good defects to control n and scatter phonons without also scattering electrons. Some materials already have intrinsic phonon scattering mechanisms in their pure state and thus exhibit fairly high zT values without additional processes being required to scatter phonons. The search for such intrinsically low-κ_ph materials is currently being conducted with the help of sophisticated phonon calculations and machine learning.²⁾

Data-driven searches for intrinsically high-zT materials within vast lists of materials are rapidly progressing. Recently, high-throughput first-principles calculations^3–7) have been combined with Boltzmann transport calculation codes such as those in the Boltzmann Transport Properties (BoltzTraP⁸⁾) program package. In the TE Design Lab database^4,5) the electronic properties that are related to the thermoelectric properties of over 2,300 compounds have been calculated theoretically. The database was carefully designed to prevent misleading people into believing that a single representative value for zT can be predicted for each compound. Instead of directly predicting the maximum zT, the database provides the values of the β-factor,⁵⁾ which is closely related to zT via zT = vβ/(uβ + 1). Properties that are intrinsic to the electronic structure are included in β, and the extrinsic variables related to chemical potential and scattering mechanisms are excluded (they are represented by v and u, respectively). Compounds with a larger β value have higher chances of exhibiting a high zT after optimization of n and suppression of the electron scattering.

In such calculations of thermoelectric materials it is relatively easy to reproduce the experimental zT values, because we can set arbitrary values for n, τ_el, and κ_ph if necessary. Observing the trends is also possible, because we can fix these unknown variables as constants, such as the common approximation τ_el = 10⁻¹⁴ s. However, predicting zT is not easy, because all these unknown variables need to be determined accurately.

Since S, σ, and κ_el are strongly dependent on the carrier concentration, we need to determine the chemical potential μ that represents the value of n in the sample. We also need to make sure that the calculated value for ε_g is correct, because determining ε_g from first-principles calculations often yields incorrect values, even though many thermoelectric materials are narrow-gap semiconductors.

The solutions of the Boltzmann transport equations^8,9) yield σ, S, and κ as a function of both μ and T by integrating over energy ε such that   

\begin{equation} \sigma (\varepsilon,T) = \frac{e^{2}}{3}D(\varepsilon)v^{2}(\varepsilon)\tau_{\textit{el}}(\varepsilon,T), \end{equation} (2)

  

\begin{equation} \sigma (\mu,T) = \int \sigma (\varepsilon,T)\left(- \frac{\partial f(\varepsilon,T)}{\partial \varepsilon}\right)d\varepsilon, \end{equation} (3)

  

\begin{equation} S(\mu,T) = - \frac{1}{eT\sigma (\mu,T)} \int (\varepsilon - \mu)\sigma (\varepsilon,T)\left(- \frac{\partial f(\varepsilon,T)}{\partial \varepsilon}\right)d\varepsilon, \end{equation} (4)

  

\begin{align} \kappa_{\text{el}}(\mu,T) &= \frac{1}{e^{2}T} \int (\varepsilon - \mu)^{2}\sigma (\varepsilon,T)\left(- \frac{\partial f(\varepsilon,T)}{\partial \varepsilon}\right)d\varepsilon \\ &\quad- S(\mu,T)^{2}\sigma (\mu,T)T. \end{align} (5)

Here, e is the charge of an electron, σ(ε, T) is the transport distribution function (spectral conductivity), D(ε) is the electron density of states, v(ε) is the electron group velocity, and f(ε, T) is the Fermi-Dirac distribution. The actual calculation is conducted in the form of tensors, and the values of these tensors are obtained by taking the average of the diagonal components of the tensors.⁸⁾ Although τ_el(ε, T) is included in σ(μ, T) and σ(ε, T), if τ_el(ε, T) does not depend on ε, then τ_el(ε, T) cancels out during the calculation thus yielding exact values for the S(μ, T) tensor. This cancellation of τ_el(ε, T) does not occur when calculating the σ(μ, T) and κ_el(μ, T) tensors. As a result, we can only obtain the representative values for σ and κ_el, as σ/τ_el and κ_el/τ_el as functions of μ and T.

Note that there are two terms for κ_el in eq. (5). The first term κ₀ represents the normal diffusion of heat. The second term corresponds to the self-Peltier effect: the electric current induced by the thermoelectric voltage carries the heat backwards from the low-T to the high-T end. In this paper, we refer to this second term as the Peltier thermal conductivity κ_P ≡ −S²σT. When calculating κ_el, κ_P is often neglected because it is negligibly small in ordinary metals.⁹⁾ However, this term may not be negligible in thermoelectric materials, because they are designed to maximize P.

By using these calculable parameters σ/τ_el, S, and κ_el/τ_el, the full expression for zT becomes   

\begin{equation} zT = \frac{S^{2}\sigma T}{\kappa_{\text{el}} + \kappa_{\text{ph}}} = \frac{S^{2}(\sigma/\tau_{\text{el}})T}{(\kappa_{\text{el}}/\tau_{\text{el}}) + (\kappa_{\text{ph}}/\tau_{\text{el}})}. \end{equation} (6)

In this expression, we can see that an unknown term (κ_ph/τ_el) still remains. While there have been intensive studies on κ_ph, few studies have ever been conducted on τ_el. Sophisticated phonon calculations can predict the κ_ph values of pure compounds;¹⁰⁾ however, the κ_ph values in high-zT samples are significantly lower than those of pure materials owing to the extrinsic phonon scattering centers introduced intentionally or unintentionally during the fabrication process. Experimental κ_ph values estimated using the Wiedemann-Franz law κ_el = LσT are also unreliable, because the law is not applicable to semiconductors, whose transport distribution function is not linear as a function of energy near μ.⁹⁾ It is also known that L in semiconductors tends to deviate from the theoretical L by ∼20%¹⁾ or more.¹¹⁾

Here, the electron scattering time τ_el is the interval between electron scattering events. Because of a lack of knowledge regarding the real value of τ_el, a constant relaxation time approximation with a fixed τ_el is used for all T to roughly estimate σ and κ_el. A widely used approximation is τ_el = 1 × 10⁻¹⁴ s.^3–8) Compared with their respective pure compounds, high-zT thermoelectric materials are expected to have a higher concentration of electron scattering centers as they originate from the materials that were introduced to control n and scatter phonons.

Scattering theories predict that τ_el has a T-dependence that depends on the dominant scattering mechanisms. When acoustic phonon scattering is dominant, τ_el is proportional to T^−1/2 in covalent compounds, T⁻¹ in metals, and T^−3/2 in ionic compounds. When ionic impurity scattering is dominant, τ_el is proportional to T^3/2. There are also other scattering centers such as grain boundaries and crystal defects. The overall τ_el can be calculated from the inverse summation of τ_el of every scattering mechanism (Στ_el⁻¹)⁻¹.

There are a few theoretical approaches to calculate the τ_el values of Si^12,13) and PbTe¹⁴⁾ from first-principles. Their results suggest that τ_el lies between 10⁻¹⁵ and 10⁻¹³ s. Additionally, there have been attempts to estimate the τ_el of Mg₂Si_1−xSn_x,^15,16) skutterudites,¹⁷⁾ SnSe,¹⁸⁾ Ga₂Ru,¹⁹⁾ and AgGaTe₂²⁰⁾ by fitting the Boltzmann transport equations with experimental σ and carrier concentration values. Further, the electron mobility R_Hσ = eτ_el/m*, where R_H is the experimental Hall coefficient and m* is the calculated effective mass, can be used to estimate τ_el in experimental samples, if the correct value for m* is used. Chen et al. report a list of the ratio of both the experimental and calculated mobilities of 31 compounds²¹⁾ under the assumption that τ_el = 1 × 10⁻¹⁴ s; the ratio ranged from 0.0005 to 10, corresponding to τ_el values between 7 × 10⁻¹⁸ s and 9 × 10⁻¹⁴ s. Although this mobility calculation has other possible sources of errors such as effective mass inaccuracy and anisotropy, this implies that τ_el can vary by orders of magnitude between different samples.

When we ignore this unknown term κ_ph/τ_el, the expression becomes z_elT: the theoretical upper limit for zT in the limit of κ_ph → 0 is given by   

\begin{equation} z_{\text{el}}T(\mu,T) \equiv \frac{S^{2}\sigma}{\kappa_{\text{el}}} = \frac{S^{2}(\sigma/\tau_{\text{el}})T}{\kappa_{\text{el}}/\tau_{\text{el}}} = \frac{S(\mu,T)^{2}}{L_{\text{eff}}(\mu,T)}. \end{equation} (7)

The effective Lorenz factor   

\begin{equation} L_{\text{eff}}(\mu,T) \equiv \frac{\kappa_{\text{el}}(\mu,T)}{\sigma (\mu,T)T} = \frac{\kappa_{\text{el}}(\mu,T)/\tau_{\text{el}}}{(\sigma (\mu,T)/\tau_{\text{el}})T} \end{equation} (8)

is another parameter that can be calculated without knowing τ_el.

In this study, we attempted to demonstrate how these unknown parameters and calculation errors affect the Boltzmann transport calculation of a material’s thermoelectric properties as this calculation is widely used in high-throughput first-principles calculations. We selected Si (Si–Ge alloy) as the test material for the demonstration. It should be noted here that more precise analysis models exist for the theoretical and experimental transport properties of Si–Ge alloys.^11,22–24)

2. Calculation Method

The electronic structure of Si was calculated using WIEN2k code,²⁵⁾ which employs the full-potential linearized augmented plane wave (FLAPW) method. The initial lattice parameters and atomic coordinates were obtained at room temperature. For the exchange-correlation potential, the generalized gradient approximation (GGA) by Perdew, Burke, and Ernzerhof²⁶⁾ was applied. In comparison, the Tran-Blaha modified Becke-Johnson (TB-mBJ)²⁷⁾ potential was applied to the calculation via the GGA to achieve more precise values for the band gaps.²⁸⁾ Spin-polarization and spin-orbit interactions were ignored. The number of k-points was set to 8000 in the first Brillouin zone in the initial calculation, and the mesh was later interpolated by a 5× denser mesh in reciprocal space using code in the BoltzTraP⁸⁾ program package. The muffin-tin radius R_MT was set to 2.10 a.u., while the maximum k was set such that R_MTk_max = 7. The core-valence cutoff energy was set to E_cut = −6.0 Ry. The convergence criteria for energy and charge were 0.0001 Ry and 0.0001 e, respectively. The Boltzmann transport analysis was conducted using BoltzTraP⁸⁾ version 1.5.3; the code was slightly modified²⁹⁾ to include κ_P during the calculation of the transport tensors.

3. Results and Discussion

3.1 Electronic structures

Figure 1 shows the density of states DOS and energy dispersion curves for Si calculated using the GGA²⁶⁾ and TB-mBJ²⁷⁾ potentials. When GGA was used the ε_g was underestimated (around 0.571 eV), which is nearly half the experimental ε_g for Si (1.1701 eV).²⁷⁾ This underestimation is a well-known, common problem of first-principles calculations based on density functional theory (DFT).²⁸⁾ The use of TB-mBJ resulted in a better value of ε_g of ∼1.156 eV. Although the raw calculation result provides the Fermi energy ε_F at the top of the valence band, ε_F at n = 0 (undoped) is shifted to the middle of the band gap during transport calculations by BoltzTraP.²⁹⁾ In reality, ε_F moves up with electron doping and moves down with hole doping.

Fig. 1

(a) Density of states DOS and (b) band structure of Si calculated using the FLAPW method using GGA (black) and TB-mBJ (red) exchange-correlation potentials.

Figure 1(b) shows the band structure of Si along representative reciprocal lattice vectors. The valence band top is composed of one heavy band and one light band, which are degenerate at Γ point, so that the number of hole pockets is two in p-type Si. On the other hand, the conduction band bottom is located along the Γ-X line, so there are six electron pockets in n-type Si.

It can be seen that the band structures obtained using TB-mBJ and GGA are essentially similar, except for the value of ε_g. The concavity of the bands or the effective masses m* were slightly greater with TB-mBJ than with GGA. We used the calculation results based on the experimental lattice parameter³⁰⁾ a = 5.4310 Å at 300 K since there were negligible differences in the calculated value of S. Changes in the lattice volume from −5% to 10% resulted in differences in S within ±0.5 µV/K at n ≈ 10²⁰ cm⁻³ in the GGA calculations. An optimization of the lattice parameters from GGA calculations resulted in an increased lattice parameter a = 5.4797 Å (+2.7% in volume); using this value of a, ε_g increased up to 1.212 eV in the TB-mBJ calculation. Since the agreement with the experimental ε_g was poorer with the optimized a, we employed the results with the experimental a in the following calculations.

3.2 Carrier doping level

In this study, we plotted the thermoelectric properties against the carrier doping level n in units of inverse cubic centimeters by dividing the amount of additional charge by the unit cell volume.³²⁾ This corresponds to the carrier doping level or net carrier concentration, i.e., the difference between the electron and hole concentrations. This n is expected to be constant within a sample, because n is related to the initial dopant concentration in the sample and is not affected by the thermal excitation of electron–hole pairs across the band gap. We plotted the materials’ thermoelectric properties against n on a logarithmic scale. This representation is more informative than using a linear scale, which shrinks the information around the optimum n into a “spiky” region around n = 0.

3.3 Seebeck coefficient

Figure 2(a) shows the values of S of Si at different T plotted against log n. Both p-type and n-type Si showed very similar curves, even though the band structures of the valence band and the conduction band are different. At equal n, the n-type Si showed slightly greater |S| than p-type Si, especially in bipolar states. This is possibly due to the difference in number of carrier pockets: more carriers exist near the band edge in n-type Si, so that |S| of electrons becomes greater than |S| of holes.

Fig. 2

The dependence of the absolute values of calculated Seebeck coefficient |S| on the carrier doping level n in p-type and n-type Si (a) at various T between 300 K and 1300 K with intervals of 100 K, (b) at 900 K using scissors-cut band energies with ε_g values between 0 and 1.17 eV in GGA calculation, with TB-mBJ calculation (ε_g ∼ 1.16 eV) as a reference. (c) The comparison between calculated and experimental values of |S| for Si_0.8Ge_0.2 alloys (reported by Vining et al.²²⁾) plotted against Hall carrier densities for p-type and n-type Si at 300 K. The S–n curves of Si were calculated by using TB-mBJ, and GGA with unshifted ε_g ≈ 0.571 eV. The dark solid circles correspond to the samples, which were used in the calculation of τ_el in the following sections.

At 300 K, the absolute value of |S| decreased monotonically with increasing n from 10¹⁶ to 10²² cm⁻³. Extremely high values of |S| were predicted at low n. When we compare |S| at fixed n, the |S|–n curves shifted upwards with increasing T. Above 700 K, peaks started to appear in the |S|–n curves. The peak position moved to higher n with increasing T.

These behaviors in S are in agreement with bipolar conduction occurring; this conduction is caused by the thermal excitation of electron–hole pairs across the band gap. The drastic decrease in |S| is induced by the increase in the total carrier concentration and by the cancellation of S by carriers of the opposite sign. This effect is significant at low n, because thermal excitation adds equal concentrations of electrons and holes, which dominate over the carrier concentration induced by dopants.

The degradation of |S| by bipolar conduction was more significant when ε_g was smaller. Figure 2(b) shows the calculated S at 900 K for artificially set ε_g values (from the experimental value 1.17 eV down to 0 eV) using the ‘scissors-cut’ operation implemented in BoltzTraP. Once ε_g was set to 1.17 eV in GGA, it roughly matched the value of S calculated using TB-mBJ. However, when ε_g was set to smaller values, |S| was reduced on the left side of the plot, indicating that the decrease in |S| at low n is influenced by bipolar conduction. When ε_g was decreased down to 0.75 eV (10 k_BT), it started to influence S at n ≈ 10¹⁹ cm⁻³, which is the optimum range of n for thermoelectric materials. When ε_g was decreased down to 0.25 eV (3.3 k_BT), the influence of bipolar conduction on S was as high as n ≈ 10²⁰ cm⁻³. When ε_g was set to 0, the peak of |S| was around 10²¹ cm⁻³; however, the peak value of |S| was very small.

These results tell us two things about ε_g. First, we should select a parent compound for thermoelectric materials that has a wide enough ε_g. This concept is known as the 10 k_BT-rule.³³⁾ At least for Si, bipolar conduction lowers the |S| of the optimum range of n when the ε_g is less than 10 k_BT. Second, when the calculated value of ε_g is small, the predictions of the thermoelectric properties will contain a large uncertainty. This highlights the difficulty of using first-principles calculations to determine a material’s thermoelectric properties. Although many thermoelectric materials are narrow-gap semiconductors, their band gaps are often underestimated in DFT calculations.

Figure 2(c) shows the agreement of the experimental S and calculated S(n) at 300 K. The figure shows the calculated S for T = 300 K in comparison with the experimental S of p-type and n-type Si_0.8Ge_0.2 alloys.²²⁾ Although there was a high agreement between calculated and experimental |S|, there were slight differences in |S|, especially for n-type Si_0.8Ge_0.2. One of the possible origins for this is that electronic structure of Si_0.8Ge_0.2 is different from the one for pure Si, especially in the conduction band.

3.4 σ/τ_el, κ_el/τ_el, and the effective Lorenz factor

Figure 3(a) and (b) show the n-dependence of the calculated value of σ/τ_el of Si. Almost no difference in σ/τ_el was observed between p-type and n-type Si. At 300–600 K, we can see that σ/τ_el is almost proportional to n. The quantity σ/τ_el was independent of T, as expected from the Drude model σ/τ_el = n_totale²/m*, where n_total = n_e + n_h $ \simeq $ n under monopolar conduction. Above 700 K, the σ/τ_el–n curves showed a bend at low n. This is because of the increase in n_total caused by bipolar conduction. It should be noted that the thermal excitation of carriers does not change n.

Fig. 3

Carrier density dependences of calculated values of (a) σ/τ_el in p-type Si, (b) σ/τ_el in n-type Si, (c) κ_el/τ_el in p-type Si, (d) κ_el/τ_el in n-type Si, (e) L_eff in p-type Si, and (f) L_eff in n-type Si between 300 and 1300 K in intervals of 100 K.

Figure 3(c) and (d) show the n-dependence of κ_el/τ_el calculated using two different expressions for κ_el. At the temperature where the bipolar conduction is negligible (300–600 K), κ_el/τ_el increased monotonically with n, just as σ/τ_el does (see Fig. 3(a) and (b)). This similarity can be explained via the Wiedemann-Franz law (κ_el = L₀σT).

Large discrepancies in the calculated results were observed when we excluded or included the Peltier term κ_P = −S²σT. When κ_el/τ_el was calculated according to κ_el/τ_el = κ₀ + κ_P, the value of κ_el/τ_el was lower than when κ_el/τ_el was calculated according to κ_el/τ_el = κ₀. This difference was the largest around n ≈ 10²⁰ cm⁻³, where P = −κ_P/T tended to reach its maximum. This shows that the approximation κ_el = κ₀ is unsuitable, especially for thermoelectric materials that are designed to achieve a maximum P. In the metallic limit (n > 10²¹ cm⁻³) and in the bipolar limit (n < 10¹⁸ cm⁻³ for T = 1300 K), the value of κ_el/τ_el calculated with κ_P was almost equal that calculated without κ_P. The bipolar effect increases κ_el not only by increasing n_total, but also by increasing the energy per carrier particle by ε_g. Thermally excited carriers possess additional energy equal to ε_g from the high-T end to the low-T end of the sample. Therefore, zT is suppressed under bipolar conduction both because of the cancellation of S and because of the increase in κ_el.

By using both the calculated κ_el/τ_el and σ/τ_el values, we calculated the effective Lorenz factor L_eff = κ_el/(σT) = (κ_el/τ_el)(σ/τ_el)⁻¹T⁻¹ and compared the results to those of the theoretical Lorenz factor L₀ in the free-electron model   

\begin{equation} L_{0} = \frac{\pi^{2}k_{\text{B}}{}^{2}}{3e^{2}} \approx 2.45 \times 10^{-8}\,\text{V$^{2}$\,K$^{-2}$}. \end{equation} (9)

Figure 3(e) and (f) show the values of L_eff calculated using two different expressions for κ_el. When we calculated κ_el using κ₀, L_eff varied greatly and differed significantly from L₀, except for in the metallic region where n > 10²¹ cm⁻³. However, when we included κ_P, L_eff was almost constant over a much wider range of n, except in the region influenced by bipolar conduction. When bipolar conduction was taking place, the increase in κ_el/τ_el resulted in an exponential increase in L_eff. The L_eff in the semiconducting region was calculated to be around 1.67 × 10⁻⁸ V² K⁻² for p-type Si and 1.86 × 10⁻⁸ V² K⁻² for n-type Si, which were 32% and 24% lower than L₀. Such decreases in L_eff have been theoretically predicted for Si by Bera et al.¹¹⁾ They reported that the L_eff of Si is dependent on calculation methods, and it may be further decreased by using exact solutions for the Boltzmann transport equations beyond the relaxation time approximation.¹¹⁾

3.5 z_elT: the electronic upper limit for zT

The n-dependence of z_elT is thought to be material-dependent and a function of both n and T; it is believed that this dependence can be determined purely from its electronic structure. However, as presented in Fig. 4(a) and (b), the n-dependence of z_elT calculated using different expressions for κ_el exhibited very different behaviors. When we used κ_el = κ₀ + κ_P, the z_elT of Si could significantly exceed 1 at low n. The simple expression z_elT = S²/L₀ also yielded similar z_elT curves to the ones obtained with κ_P, although some discrepancies were observed in the low-n region where bipolar conduction was dominant.

Fig. 4

Dependence of z_elT on the carrier doping level of (a) p-type and (b) n-type Si calculated using three different expressions for κ_el between temperatures of 300 and 1300 K in intervals of 100 K.

However, when we used the approximation κ_el = κ₀, a strange upper limit of z_elT appeared at z_elT = 1, indicating that there was no chance of achieving zT > 1. In our calculations this limit was observed for almost any semiconductors with steep band edges. This is consistent with past theoretical calculations that included κ_P^9,34–37) or used the Wiedemann-Franz law¹⁸⁾ to successfully achieve z_elT > 1, while studies that only used κ₀^7,38–41) only found values of z_elT below 1.

The origin of this false “wall of z_elT = 1” can be explained as follows. In the low-n limit of non-degenerate semiconductors, μ (the center of the Fermi-Dirac distribution) is located nearly at the middle of the band gap, where no hole/electron states exist. In normal semiconductors the DOS at the band edges are very steep. As a result, the hole/electron states are concentrated near the band edge energy ε₀, which is either the top of the valence band (ε_V) for holes or the bottom of the conduction band (ε_C) for electrons. At low T, only a few carriers can be thermally excited up to these states, and the difference in ε of each carrier from ε₀ is negligibly small, compared to the large ε₀-μ. Then we can treat ε-μ as a constant (ε₀-μ) for all carriers, which can be taken out of the integrals of eqs. (3)–(5).   

\begin{align} S(\mu,T) \xrightarrow{| \varepsilon_{0} - \mu | \gg \varDelta \varepsilon} &- \frac{1}{eT\sigma(\mu,T)}(\varepsilon_{0} - \mu)\sigma(\varepsilon_{0},T) \\ &\quad \rightarrow - \frac{1}{eT}(\varepsilon_{0} - \mu) \end{align} (10)

  

\begin{equation} S(\mu,T)^{2}\sigma (\mu,T)T \xrightarrow{| \varepsilon_{0} - \mu | \gg \varDelta \varepsilon}\frac{1}{e^{2}T}(\varepsilon_{0} - \mu)^{2}\sigma(\varepsilon_{0},T) \end{equation} (11)

  

\begin{equation} \kappa_{0}(\mu,T) \xrightarrow{| \varepsilon_{0} - \mu | \gg \varDelta \varepsilon}\frac{1}{e^{2}T}(\varepsilon_{0} - \mu)^{2}\sigma(\varepsilon_{0},T) \end{equation} (12)

  

\begin{equation} \kappa_{\text{P}}(\mu,T) \xrightarrow{| \varepsilon_{0} - \mu | \gg \varDelta \varepsilon} - \frac{1}{e^{2}T}(\varepsilon_{0} - \mu)^{2}\sigma(\varepsilon_{0},T) \end{equation} (13)

  

\begin{equation} \therefore z_{\text{e}}T = \frac{S(\mu,T)^{2}\sigma (\mu,T)T}{\kappa_{0}(\mu,T) + \kappa_{\text{P}}(\mu,T)} \xrightarrow{| \varepsilon_{0} - \mu | \gg \varDelta \varepsilon}\left[ \begin{array}{l} 1\ (\text{$\kappa_{0}$ only})\\ \infty\ (\kappa_{0} + \kappa_{\text{P}}) \end{array} \right. \end{equation} (14)

When we include κ_P, the denominator κ₀ + κ_P becomes 0 and z_elT diverges to infinity in the low-n, low-T limit. However, when we neglect κ_P, z_elT does not exceed 1. This implies that zT > 1 can only be achieved when a large power factor (S²σ) cancels out a part of κ₀ via the self-Peltier effect. Mahan and Sofo have proposed that when the transport distribution function $\sigma (\varepsilon )( - \frac{\partial f(\varepsilon ,T)}{\partial \varepsilon })$ is a Dirac delta function and μ is located away from ε, then κ_P completely cancels out κ₀ to yield an extremely high value of zT.³⁴⁾ We consider that ordinary semiconductors in the low-n limit may have transport distribution functions that are like delta-functions because of their steep band edges, thus fulfilling this condition for high zT.

3.6 Effects of κ_ph and τ_el on zT

Although knowing z_elT is useful for predicting the upper limit for zT, we need to estimate the actual values of zT by including the contribution of the ignored term κ_ph/τ_el. A theoretical estimation of this parameter is difficult, since both κ_ph and τ_el are expected to strongly depend on the material and its microstructure since phonons and electrons scatter at impurities, grain boundaries, and crystal defects.

As shown in Fig. 5, the introduction of κ_ph/τ_el results in a large reduction of zT at low n. The strangely high z_elT at low n only occurs when neglecting κ_ph, as it is the dominant term of κ in the low-n region. Therefore, the value of z_elT should not be used to estimate zT in the low-n region. However, at high n the value of z_elT is more related to zT, because then zT is less influenced by introducing κ_ph/τ_el as κ_el is dominant. This unknown term κ_ph/τ_el can be understood to represent the material’s degree of phonon-glass electron-crystal, a concept that was first proposed by Slack.⁴²⁾

Fig. 5

Dependence of zT on the carrier doping level for different values of κ_ph/τ_el in [W m⁻¹ K⁻¹ s⁻¹] for (a) p-type Si at 300 K, (b) n-type Si at 300 K, (c) p-type Si at 900 K, and (d) n-type Si at 900 K. The data for κ_ph/τ_el = 0 are presented as z_elT in Fig. 4.

3.7 Estimation of κ_ph and τ_el from the experimental values

It is difficult to calculate κ_ph and τ_el from first-principles calculations. Even if the values of κ_ph and τ_el of pure crystals were calculated, these values will be decreased in the actual samples owing to extrinsic electron and phonon scattering mechanisms. In this study, we therefore analyzed whether κ_ph and τ_el should instead be obtained by comparing the calculation results with the experimental values.

By combining the calculation results with the reported experimental values for S_exp(T), σ_exp(T), and κ_exp(T) we first tried to find the dataset of values of ε_F that best agrees with the experimental data. This method is useful when the carrier density is not available, although it requires that the calculation reproduces the actual electronic structure. The Boltzmann transport calculation can provide thousands of datasets for different values of (ε_F, T). The parameter n_calc, which is the doped carrier concentration, is determined by ε_F. The other parameters S_calc, (σ/τ_el)_calc, and (κ_el/τ_el)_calc are all obtained as functions of (ε_F, T). We searched for the ε_F that best reproduced S_exp and n_exp. Instead of using the value of n_exp from Hall measurements, the corresponding dataset for T was found by searching for a dataset with S_calc = S_exp, because in general n_exp at higher T was not reported. Two or more datasets were found to have S_calc = S_exp when the S–n curve peaked; we selected the dataset with the most reasonable value for n.

Figure 6(a) shows the T-dependence of ε_V-μ for p-type and ε_C-μ for n-type semiconductors that was estimated by comparing the value of S calculated using TB-mBJ to the values of S_exp of three samples of Si_0.8Ge_0.2 alloys that were reported by Vining et al.²²⁾ and Dismukes et al.⁴¹⁾ Here, the DOS of Si_0.8Ge_0.2 and Si was assumed to be the same. Vining’s p-type (No. 75) and n-type (No. 93) samples were obtained via ball-milling of the alloys followed by hot-pressing. A Dismukes’ n-type sample (No. 1834) was obtained using the zone-levelling crystal growth technique so that the sample would have a higher crystallinity and larger grains compared with the Vining’s samples. Since ε_V-μ for the p-type semiconductor yielded almost 0 meV and ε_C-μ for the n-type yielded about 6 meV, we can see that μ is located almost at the band edge, as expected for semiconductors near the degenerate regions. The n estimated from S in Fig. 6(b) was on the order of 10²⁰ cm⁻³, which is consistent with the experimental values. The cause of the slight increase with T thus appears to be either calculation errors or the thermal excitation of localized carriers.

Fig. 6

Temperature dependence of (a) the chemical potential, where ε_V is the top of the valence band and ε_C is the bottom of the conduction band; (b) the carrier doping level; (c) the effective relaxation time; (d) the inverse of the effective relaxation time, with linear fits of the data points; and (e) the phonon thermal conductivity in comparison with κ_total–L₀σT (displayed as pale lines and symbols) that was estimated using experimental values of S, σ, and κ of 0.08% boron-doped Si_0.8Ge_0.2 alloys reported in the literature. The ball-milled, pressure-sintered samples prepared by Vining et al.²²⁾ is shown as BM, and the zone-levelled, high-crystallinity sample by Dismukes is shown as ZL.

The effective relaxation time $\tau _{\text{el}}^{\text{eff}}$ was obtained as   

\begin{equation} \tau_{\text{el}}^{\text{eff}} = \frac{\sigma_{\text{exp}}}{(\sigma/\tau_{\text{el}})_{\text{calc}}} \end{equation} (15)

by using experimental values of the electrical conductivity σ_exp. Although $\tau _{\text{el}}^{\text{eff}}$ differs from the value of τ_el of the ideal crystals,^12,13) $\tau _{\text{el}}^{\text{eff}}$ and $\tau _{\text{el}}^{\text{ideal}}$ should satisfy Matthiesen’s law in the Boltzmann transport equations,   

\begin{equation} \frac{1}{\tau_{\text{el}}^{\text{eff}}} = \frac{1}{\tau_{\text{el}}^{\text{ideal}}} + \frac{1}{\tau_{\text{el}}^{\text{ext}}}, \end{equation} (16)

where $\tau _{\text{el}}^{\text{ext}}$ is the extrinsic relaxation time, which is determined by impurity doping and the microstructure.

Figure 6(c) shows the T-dependence of $\tau _{\text{el}}^{\text{eff}}$ estimated for the three samples. All samples exhibited a value of $\tau _{\text{el}}^{\text{eff}}$ on the order of 10⁻¹⁵ s, while a monotonic decrease of $\tau _{\text{el}}^{\text{eff}}$ was observed with increasing T. We also estimated the value of τ_el of the other thermoelectric materials by using the values of S_exp, σ_exp, and κ_exp at their best values of zT. Although the calculations were not reliable enough, the values of τ_el were estimated to be between 10⁻¹⁵ and 10⁻¹⁴ s; single crystals⁴³⁾ tended to exhibit longer τ_el values than polycrystalline samples.⁴⁴⁾

This analysis of $\tau _{\text{el}}^{\text{eff}}$ enabled us to determine the dominant scattering mechanism in the samples. Figure 6(d) shows that in the two n-doped samples $(\tau _{\text{el}}^{\text{eff}})^{ - 1}$ and T had a linear relationship to each other. Therefore, the dominant scattering mechanism at high T in these samples was expected to be scattering caused by an acoustic distortion potential (in metals), where τ_el is proportional to T⁻¹. The scattering mechanism in the p-type sample appears to be more independent of T.

In Fig. 6(e), the effective phonon thermal conductivity was calculated according to   

\begin{equation} \kappa_{\text{ph}}^{\text{eff}} = \kappa_{\text{exp}} - \tau_{\text{el}}^{\text{eff}}(\kappa_{\text{el}}/\tau_{\text{el}})_{\text{calc}}, \end{equation} (17)

using the experimental thermal conductivity κ_exp. Since the value of L_eff of the compound was smaller than of L₀, the value of $\kappa _{\text{ph}}^{\text{eff}}$ calculated by this method was higher than that using the free-electron model $\kappa _{\text{ph}}^{\text{FEM}} = \kappa _{\text{exp}} - L_{0}\sigma T$. The Dismukes’ sample with its large grains exhibited a higher $\kappa _{\text{ph}}^{\text{eff}}$ than the Vining’s sample with its small grains. These results imply that κ_ph is strongly dependent on the density of the grain boundaries, which act as phonon scattering centers.

4. Remaining Sources of Errors

The sources of errors in first-principles calculation of the thermoelectric properties of materials can be classified into three categories: over-simplification of the crystal structure, inaccuracy in the electronic structure calculation, and inaccuracies in the transport calculation.

An over-simplification of the crystal structure can occur in calculations of almost any thermoelectric material since the highest zT values are usually achieved in extremely complex structures.¹⁾ For carrier concentration control and phonon scattering, the crystals are solid solutions of various impurity elements that are randomly distributed throughout the crystal. Defects and disorders can act in a similar way, making it impossible to define unit cells. A number of thermoelectric materials have incommensurate or misfit-layer structures, which are difficult to define by periodic unit cells. The presence of such structures would modify the electronic structure near the band gap so that it is different from that of the parent compound.

The second type of error is related to the electronic structure calculations. Underestimation (and sometimes, overestimation) of ε_g is an unavoidable error in DFT. This error is critical for narrow-gap semiconductors. However, if the bandgap was at least open, the scissors-cut operation of the band energies in BoltzTraP⁸⁾ could be used to improve the accuracy.²¹⁾ The experimental values of the band gaps can be obtained from databases such as citrination.⁴⁵⁾ The use of other high-cost exchange-correlation potentials that reproduce ε_g may also be helpful, although the errors in m* need to be small enough. The introduction of various interactions, such as spin-polarizations for magnetic elements and spin-orbit interactions for heavy elements, is also necessary to reduce errors. Magnetic insulators exhibit a finite ε_g only in spin-polarized calculations, and Mott insulators may require an electron localization potential U to open the band gaps. Spin-orbit interactions should be dominant in many heavy-element thermoelectric materials; these split the bands into two similar bands separated by several hundred millielectronvolts. At the crossing-points of those bands, band reconstruction occurs to generate bands with small m* that are separated by a small gap. In some compounds, the electronic structures are sensitive to the lattice volume and to atomic positions. In such cases a careful selection of the calculation conditions is necessary.

The third type of error is related to the Boltzmann transport calculations. This type of error can be usually be mitigated by increasing the density of the k-point mesh. In BoltzTraP, the use of a coarse k-point mesh results in oscillations and noise in the thermoelectric properties at high n. The introduction of anisotropy would improve the calculations of anisotropic compounds. A more precise treatment of τ_el, which should be dependent on energy and bands, may also improve accuracy of the Boltzmann transport calculations.

5. Conclusion

We investigated how the Boltzmann transport calculation of the thermoelectric properties of a typical semiconductor (Si) from a first-principles calculation is affected by unknown variables such as n, ε_g, and κ_ph/τ_el. We determined that to achieve a high zT the parent compound should have a wide ε_g, low κ_ph, and long τ_el. We demonstrated that calculated thermoelectric properties such as S, σ/τ_el, κ_el/τ_el, z_elT, and zT can vary by orders of magnitude when n is varied. The prediction of an extremely high value for S was not found to be special, as it was achieved in a normal semiconductor with a wide enough ε_g in the low-n limit. These properties were strongly affected by ε_g, even though the estimation of ε_g is difficult in DFT. This error is critical for thermoelectric materials, because narrow-gap semiconductors are favored for ease of carrier doping. We also found that κ_P is required for calculations of κ_el when we need to show that zT can exceed 1. The calculated L_eff was nearly equal to L₀ in the metallic limit of n > 10²¹ cm⁻³; however, it was ∼20% lower than L₀ when n < 10²⁰ cm⁻³ and increased by orders of magnitude when bipolar conduction became dominant. The values of zT were largely dependent on κ_ph and τ_el, although both parameters are strongly influenced by typical engineering processes such as impurity doping and microstructural control.

We attempted to estimate the T-dependence of μ, n, τ_el, and κ_ph by fitting the calculation results with the T-dependences of S_exp, σ_exp, and κ_exp provided in experimental papers on the thermoelectric properties of certain materials. We showed that τ_el was on the order of several 10⁻¹⁵ s and varied with T as expected for the dominant electron scattering mechanism. It was also shown that the values of τ_el differed between samples fabricated in different ways.

To predict semiconductor’s thermoelectric properties correctly, it is necessary to reduce many errors in the values of parameters related to the crystal structures, electronic structure calculations, and Boltzmann transport calculations. We also found that it is necessary to careful check that the calculated results are in agreement with the experimental data for any calculations of the thermoelectric properties of semiconductors from first-principles.

Acknowledgments

This work was supported by KAKENHI grants (grant numbers 23760647 and 16K14379) from the Ministry of Education, Culture, Sports, Science and Technology and by the “Materials research by Information Integration” Initiative (MI²I) project of the Support Program for Starting Up Innovation Hub from Japan Science and Technology Agency (JST). We thank Jim Bailey, PhD, from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

REFERENCES

• 1)   G.J.  Snyder and  E.S.  Toberer: Nat. Mater. 7 (2008) 105.
• 2)   A.  Seko,  A.  Togo,  H.  Hayashi,  K.  Tsuda,  L.  Chaput and  I.  Tanaka: Phys. Rev. Lett. 115 (2015) 205901.
• 3)   G.K.H.  Madsen: J. Am. Chem. Soc. 128 (2006) 12140.
• 4)   P.  Gorai,  D.  Gao,  B.  Ortiz,  S.  Miller,  S.A.  Barnett,  T.  Mason,  Q.  Lv,  V.  Stevanović and  E.S.  Toberer: Comput. Mater. Sci. 112 (2016) 368.
• 5)   J.  Yan,  P.  Gorai,  B.  Ortiz,  S.  Miller,  S.A.  Barnett,  T.  Mason,  V.  Stevanović and  E.S.  Toberer: Energy Environ. Sci. 8 (2015) 983.
• 6)   Q.  Hao,  D.  Xu,  N.  Lu and  H.  Zhao: Phys. Rev. B 93 (2016) 205206.
• 7)  M. Núñez-Valdez, Z. Allahyari and A.R. Oganov: arXiv:1610.07824v1.
• 8)   G.K.H.  Madsen and  D.J.  Singh: Comput. Phys. Commun. 175 (2006) 67.
• 9)   T.  Takeuchi: Mater. Trans. 50 (2009) 2359.
• 10)   J.E.  Turney,  E.S.  Landry,  A.J.H.  McGaughey and  C.H.  Amon: Phys. Rev. B 79 (2009) 064301.
• 11)   C.  Bera,  M.  Soulier,  C.  Navone,  G.  Roux,  J.  Simon,  S.  Volz and  N.  Mingo: J. Appl. Phys. 108 (2010) 124306.
• 12)   B.  Qiu,  Z.  Tian,  A.  Vallabhaneni,  B.  Liao,  J.M.  Mendoza,  O.D.  Restrepo,  X.  Ruan and  G.  Chen: Europhys. Lett. 109 (2015) 57006.
• 13)   M.  Fiorentini and  N.  Bonini: Phys. Rev. B 94 (2016) 085204.
• 14)   S.  Ahmad and  S.D.  Mahanti: Phys. Rev. B 81 (2010) 165203.
• 15)   X.J.  Tan,  W.  Liu,  H.J.  Liu,  J.  Shi,  X.F.  Tang and  C.  Uher: Phys. Rev. B 85 (2012) 205212.
• 16)   H.  Wang,  W.  Chu and  H.  Jin: Comput. Mater. Sci. 60 (2012) 224.
• 17)   J.  Yang,  P.  Qiu,  R.  Liu,  L.  Xi,  S.  Zheng,  W.  Zhang,  L.  Chen,  D.J.  Singh and  J.  Yang: Phys. Rev. B 84 (2011) 235205.
• 18)   S.  Hao,  F.  Shi,  V.P.  Dravid,  M.G.  Kanatzidis and  C.  Wolverton: Chem. Mater. 28 (2016) 3218.
• 19)   Y.  Takagiwa,  K.  Kitahara and  K.  Kimura: J. Appl. Phys. 113 (2013) 023713.
• 20)   W.  Wu,  K.  Wu,  Z.  Ma and  R.  Sa: Chem. Phys. Lett. 537 (2012) 62.
• 21)   W.  Chen,  J.-H.  Pöhls,  G.  Hautier,  D.  Broberg,  S.  Bajaj,  U.  Aydemir,  Z.M.  Gibbs,  H.  Zhu,  M.  Asta,  G.J.  Snyder,  B.  Meredig,  M.A.  White,  K.  Persson and  A.  Jain: J. Mater. Chem. C 4 (2016) 4414, Supplementary Material.
• 22)   C.B.  Vining,  W.  Laskow,  J.O.  Hanson,  R.  Van der Beck and  P.D.  Gorsuch: J. Appl. Phys. 69 (1991) 4333–4340.
• 23)   A.J.  Minnich,  H.  Lee,  X.W.  Wang,  G.  Joshi,  M.S.  Dresselhaus,  Z.F.  Ren,  G.  Chen and  D.  Vashaee: Phys. Rev. B 80 (2009) 155327.
• 24)   G.A.  Slack and  M.A.  Hussain: J. Appl. Phys. 70 (1991) 2694.
• 25)  P. Blaha, K. Schwarz, G.K.H. Madsen, D. Kvasnicka and J. Luitz: WIEN2k. An augmented plane wave plus local orbitals program for calculating crystal properties, Vienna University of Technology, Austria (2001).
• 26)   J.P.  Perdew,  K.  Burke and  M.  Ernzerhof: Phys. Rev. Lett. 77 (1996) 3865.
• 27)   F.  Tran and  P.  Blaha: Phys. Rev. Lett. 102 (2009) 226401.
• 28)   D.J.  Singh: Phys. Rev. B 82 (2010) 205102.
• 29)  After the line “thermal(1:3,1:3)=kappa(1:3,1:3)” in fermiintegrals.f in BoltzTraP,⁸⁾ we added the following code: “DO i=1,3; DO j=1,3; DO ialp=1,3; thermal(i,j)=thermal(i,j)+temp*seebeck(ialp,i)*nu(ialp,j); ENDDO; ENDDO; ENDDO”.
• 30)   G.  Celotti,  D.  Nobili and  P.  Ostoja: J. Mater. Sci. 9 (1974) 821.
• 31)   W.  Bludau,  A.  Onton and  W.  Heinke: J. Appl. Phys. 45 (1974) 1846.
• 32)  In BoltzTraP,⁷⁾ the number of additional charge per unit cell is given as “N” in case.trace. The unit cell volume, which needs to be the reduced cell volume, is written in case.outputtrans.
• 33)   G.D.  Mahan: Solid State Phys. 51 (1998) 81.
• 34)   T.  Fang,  S.  Zheng,  T.  Zhou,  L.  Yan and  P.  Zhang: Phys. Chem. Chem. Phys. 19 (2017) 4411.
• 35)   Y.  Takagiwa,  K.  Kitahara,  Y.  Matsubayashi and  K.  Kimura: J. Appl. Phys. 111 (2012) 123707.
• 36)   G.D.  Mahan and  J.O.  Sofo: Proc. Natl. Acad. Sci. USA 93 (1996) 7436–7439.
• 37)   Z.  Fan,  H.Q.  Wang and  J.C.  Zheng: J. Appl. Phys. 109 (2011) 073713.
• 38)   Y.  Katsura and  H.  Takagi: J. Electron. Mater. 42 (2012) 1365.
• 39)   L.F.  Wan and  S.P.  Beckman: Phys. Chem. Chem. Phys. 16 (2014) 25337.
• 40)   B.-Z.  Sun,  Z.  Ma,  C.  He and  K.  Wu: Phys. Chem. Chem. Phys. 17 (2015) 29844.
• 41)   S.  Kessair,  O.  Arbouche,  K.  Amara,  Y.  Benallou,  Y.  Azzaz,  M.  Zemouli,  M.  Bekki,  M.  Ameri and  B.S.  Bouazza: Indian J. Phys. 90 (2016) 1403.
• 42)  G.A. Slack: CRC Handbook of Thermoelectrics 407 (1995).
• 43)   J.P.  Dismukes,  L.  Ekstrom,  E.F.  Steigmeier,  I.  Kudman and  D.S.  Beers: J. Appl. Phys. 10 (1976) 2899.
• 44)   Y.  Katsura: Mater. Trans. 57 (2016) 1035.
• 45)   J.  O’Mara,  B.  Meredig and  K.  Michel: JOM 68 (2016) 2031. URL: https://citrination.com/.