A first-principles informed approach to the calculation of Lorenz numbers for complex thermoelectric materials is presented and discussed. Example calculations illustrate the importance of using accurate band structures and energy-dependent scattering times. Results obtained by assuming that the scattering rate follows the density-of-states show that in the non-degenerate limit, Lorenz numbers below the commonly assumed lower limit of 2(kB/q)2 can occur. The physical cause of low Lorenz numbers is explained by the shape of the transport distribution. The numerical and physical issues that need to be addressed in order to produce accurate calculations of the Lorenz number are identified. The results of this study provide a general method that should contribute to the interpretation of measurements of total thermal conductivity and to the search for materials with low Lorenz numbers, which may provide improved thermoelectric figures of merit, zT.

Knowledge of the Lorenz number, L, is essential for interpreting measurements of total thermal conductivity, and it has been recently noted that identifying materials with low Lorenz numbers may provide a path to higher thermoelectric figures of merit.1,2 Measurements of L can be done,3–6 but because such measurements are involved, they are not routinely performed. It has been recently pointed out that for parabolic energy bands with acoustic deformation potential (ADP) scattering, there is a direct relation between L and the routinely measured Seebeck coefficient, S.7 The L vs. S characteristic is independent of the value of the effective mass and independent of temperature when bipolar conduction is negligible. This approach provides a convenient way to determine L from a measured S, but the assumption of a simple band structure and scattering process raises concerns about its generality. Calculations of L can be done (e.g., Ref. 8), but they typically assume simplified band structures and scattering processes. Numerical solutions to the Boltzmann Transport Equation in the Relaxation Time Approximation9,10 should provide the most accurate predictions of L when informed by rigorous treatments of electron scattering.9,11–14 In this paper, we present a first-principles informed approach for calculating L for complex thermoelectric materials. As an example of the approach, we present calculations for Bi2Te3 that assume the scattering rate follows the density-of-states (DOS). These calculations show that for non-degenerate conditions, L can be lower than the often-assumed lower limit of 2(kB/q)2, and they contribute to the understanding of the physical mechanisms behind low Lorenz numbers. Under degenerate conditions, L can deviate from the Sommerfeld limit of (π2/3)(kB/q)2, which often describes metals and heavily doped semiconductors.15 The addition of bipolar conduction further complicates the calculation and can result in unusually large Lorenz numbers. Finally, we discuss uncertainties in the calculation of L in order to guide future work directed at the accurate calculation of L for complex thermoelectric materials.

Our approach is to calculate the L vs. S characteristic for complex thermoelectric materials and to compare the results to analytical calculations for simple, parabolic energy bands. It is convenient to use the L vs. S characteristic because both L and S depend on the energy dependence of the scattering time, but not on its magnitude. The standard expressions for the thermoelectric transport coefficients that result from a relaxation time approximation solution to the Boltzmann transport equation are

σ=+σ(E)dE,
(1a)
S=1qT+(EEF)σ(E)dE+σ(E)dE,
(1b)
κ0=1q2T+(EEF)2σ(E)dE,
(1c)
κe=κ0TσS2,
(1d)
LκeσT,
(1e)

where S is the Seebeck coefficient, κo is the short-circuit electronic thermal conductivity, κe is the open-circuit electronic thermal conductivity, and the differential conductivity is

σ(E)=q2Ξ(E)(f0E),
(1f)

and the diagonal component of the transport function is10,16

Ξxx(E)kυx2τm(E)δ(EEk),
(1g)

where an energy-dependent scattering time has been assumed. As shown by Mahan and Soho, all thermoelectric transport coefficients are determined by the transport function,10 so the Lorenz number, L, which is a ratio of transport coefficients as given by (1e), is also determined by the transport function.

We will take (1e) as the definition of the Lorenz number, which is different from the Wiedemann-Franz Law, which assigns a specific value, L=L0=(π2/3)(kB/q)2 to the Lorenz number.17 It is common, however, to refer to κe/σ=LT as the “Wiedemann-Franz Law” and apply it to semiconductors for which L is rarely L0, as well as to metals (see, for example, Ref. 18). As Mahan and Bartkowiak state, the Wiedemann-Franz Law should be regarded as a rule of thumb.18 Violations occur, especially for electronic structures and scattering processes that produce energetically sharp transport functions. One example is the one-dimensional organic crystals considered in Ref. 2. Our goal in this paper is to present a method to compute the Lorenz number as defined by (1e) for thermoelectric materials using the best available information about band structures and scattering times.

In this work, we solve Eq. (1) but find it convenient to write the transport function in the Landauer form as19 (see also the supplementary material)

Ξxx(E)=2h(M(E)/A)λ(E),
(2)

where M(E)/A is the number of conduction channels per unit cross-sectional area vs. energy. We compute M(E)/A from a density functional Theory (DFT)-generated band structure using the open source tool, LanTraP 2.0.20 The energy-dependent mean-free-path for backscattering is also needed; it is defined as19 

λ(E)2υx2(E)|υx(E)|τm(E),
(3)

where τm(E) is the momentum relaxation time, and the quantity, υx2/|υx|, is an angle-averaged velocity and is computed as a function of energy from the DFT-generated band structure. For ADP scattering in the elastic limit, the scattering rate is isotropic, equal to the momentum relaxation rate, and proportional to the density of states

1τ(E)=1τm(E)=KelphDOS(E).
(4)

More generally, it is often found that the electron-phonon scattering rate follows the DOS in non-polar semiconductors9,21 and even in polar semiconductors when the electron-phonon coupling is strongly screened (see Sec. IV). The density-of-states can be computed directly from the numerical band structure. The electron-phonon coupling parameter, Kelph, is proportional to the deformation potential squared, but its value does not need to be specified because it appears in both the numerator and denominator of the expressions and, therefore, drops out when computing both L and S (which is one reason for focusing on the L vs. S characteristic). When there are multiple scattering mechanisms, the relative strength of each mechanism must be known.

To illustrate the principles involved, we will begin in Sec. III with a parabolic conduction band for which Eq. (1) can be evaluated analytically when power law scattering for the energy dependent mean-free-path for backscattering is assumed

λ(E)=λ0(EEC)r/(kBT).
(5)

In this expression, r is a characteristic exponent determined by the electron scattering processes. The numerical calculations that will be presented in Sec. III evaluate λ(E) directly from (3) and (4) and do not assume the power law form of (5).

Ignoring bipolar conduction and assuming parabolic energy bands, Eq. (1) can be evaluated to find22 

S=(kBq)((r+2)Fr+1(ηF)Fr(ηF)ηF)=S(kB/q),
(6a)
L=(kBq)2Γ(r+3)Γ(r+2)((r+3)Fr+2(ηF)Fr(ηF)(r+2)(Fr+1(ηF)Fr(ηF))2)=L(kB/q)2,
(6b)

where

ηF=(EFEC)/kBT,
(7)

is the reduced Fermi level. In (6a) and (6b), we defined the dimensionless Seebeck coefficient and Lorenz number, S and L, respectively. In these expressions

Fj(ηF)=1Γ(j+1)0ηjdη1+eηηF,
(8)

is the Fermi-Dirac integral of order j as defined by Blakemore.23 Note that the value of the effective mass and the magnitude of the mean-free-path, λ0, do not appear in these expressions. The L vs. S characteristic computed from Eq. (6) depends only on the scattering exponent, r, and is independent of temperature when bipolar conduction is negligible.

In this section, we compute the L vs. S characteristic in three different ways. The first is an analytical calculation that assumes parabolic bands and power law scattering. The second is a full, numerical treatment of silicon (where the bands are approximately parabolic), and the third is a numerical treatment of Bi2Te3, which has a complex band structure. Comparing these three cases sheds light on how the transport distribution determines L.

Figure 1 shows the computed L vs. S characteristic for various values of r. Acoustic deformation potential (ADP) scattering (commonly thought to be the dominant scattering mechanism in thermoelectric materials24) corresponds to r=0. For degenerate semiconductors, ηF0, S0. In the degenerate limit, Fig. 1 shows that the normalized Lorenz numbers, L=L/(kB/q)2, approach the Sommerfeld limit of π2/3.15 For non-degenerate semiconductors, |S| is large, and Fig. 1 shows that for r=0, L approaches the expected limit of 2. For r=2, which corresponds to ionized impurity scattering, L4 in the non-degenerate limit. For r<0, Fig. 1 shows that L<2 in the non-degenerate limit.

FIG. 1.

The L vs. S characteristic for parabolic energy bands showing several different values of the power law scattering exponent, r. For r=0, which corresponds to ADP scattering in 3D, the results are identical to those in Ref. 7.

FIG. 1.

The L vs. S characteristic for parabolic energy bands showing several different values of the power law scattering exponent, r. For r=0, which corresponds to ADP scattering in 3D, the results are identical to those in Ref. 7.

Close modal

As given by Eq. (1), the thermoelectric transport coefficients are integrals over energy of various powers of (EEF). The weighting factor is the transport function, Ξ(E), times the Fermi window, (f0/E). In a Landauer picture, it can be shown that the transport distribution is proportional to the number of channels for conduction, M(E), times the mean-free-path for backscattering, λ(E).19 For three-dimensional electrons in parabolic bands, one can show that M(E) increases linearly with energy, so for power law scattering

Ξ(E)M(E)λ(E)(EEC)r+1.
(9)

When r=0 (ADP scattering in parabolic bands), the transport distribution increases linearly with energy and L=2 results in the non-degenerate limit. For r<0, the transport distribution increases sub-linearly with energy, and L<2 results in the non-degenerate limit.

To understand why the Lorenz number is sensitive to the characteristic exponent, r, it is useful to write the transport parameters in terms of averages over energy of moments of (EEF). For S and L, the result is22 

S=(EEFkBT),
(10a)
L=(EEFkBT)2(EEFkBT)2,
(10b)

where the brackets denote an average over energy. Equation (10) shows that the Lorenz number emphasizes the higher energies more that the Seebeck coefficient does. For r>0, higher energies are weighted more than for r=0, so L>2. For r<0, higher energies are weighted less than for r=0, so L<2. When we turn next to the numerical computation of the L vs. S characteristic, we will see that the non-degenerate limits can be explained by examining the shape of the numerically calculated transport distribution.

As a final note, we point out that Eq. (1) assumes three-dimensional electrons, but if we use the corresponding equations for 2D and 1D electrons,22 we find that for ADP scattering in parabolic bands, the L vs. S characteristic is identical in 1D, 2D, and 3D. This can be easily seen from (2). The mean-free-path is proportional to velocity times scattering time, and the scattering time is inversely proportional to the DOS. The number of channels is proportional to velocity times the DOS,19,22 so for parabolic bands, we find

Ξxx(E)(M(E)/A)×λ(E)υ(E)DOS(E)×υ(E)/DOS(E)υ2(E)(EEC).
(11)

Since the transport function is independent of dimension for ADP scattering in parabolic bands, all thermoelectric transport coefficients are independent of dimension.

In this case, we evaluate Eq. (1) assuming a DFT-generated band structure. Figure 2 shows the density-of-states computed from the numerical band structure. The relevant portion where the Fermi window overlaps the density-of-states lies very close to the Fermi level, which is chosen to be at the valence band edge in this case. The corresponding distributions of channels, M(E), energy-dependent mean-free-path, λ(E), and transport distribution, Ξ(E), are also shown in Fig. 2.

FIG. 2.

Key quantities computed from the numerical band structure of Si. (a) The density-of-states (solid line) and its parabolic band fit (dash line). (b) The distribution of channels (solid line), M(E), and its parabolic band fit (dashed line). (c) The energy-dependent mean-free-path, λ(E), and (d) the transport distribution, Ξ(E). The Fermi window, which is shown as a dotted line centered at the Fermi level indicates the relevant range of energies. The Fermi level is set at the valence band edge in these plots. All quantities were extracted from a DFT-generated band structure using Quantum Espresso27 (see supplementary material). The bandgap is 1.1 eV with the midgap at E = 0 eV.

FIG. 2.

Key quantities computed from the numerical band structure of Si. (a) The density-of-states (solid line) and its parabolic band fit (dash line). (b) The distribution of channels (solid line), M(E), and its parabolic band fit (dashed line). (c) The energy-dependent mean-free-path, λ(E), and (d) the transport distribution, Ξ(E). The Fermi window, which is shown as a dotted line centered at the Fermi level indicates the relevant range of energies. The Fermi level is set at the valence band edge in these plots. All quantities were extracted from a DFT-generated band structure using Quantum Espresso27 (see supplementary material). The bandgap is 1.1 eV with the midgap at E = 0 eV.

Close modal

The density-of-states shown in Fig. 2(a) confirms that the conduction and valence bands are nearly parabolic in Si. The onset of a second conduction band about 0.2 eV above the bottom of the conduction band can also be seen. The distribution of channels, M(E), shown in Fig. 2(b), is approximately linear with energy, as expected for parabolic energy bands.19,22 A change in slope for M(E) is seen in the conduction band where the second conduction band begins. The valence band for Si is known to have a complex, warped shape,25 but it has also been shown that it can be described by an equivalent parabolic band,26 which is confirmed by the linear behavior of M(E) in the valence band. Figure 2(c) shows the energy-dependent mean-free-path numerically evaluated according to Eqs. (3) and (4). Note that for a parabolic energy band with ADP scattering, λ(E)=λ0 is independent of energy. The numerically calculated λ(E) is approximately independent of energy. Numerically evaluating λ(E) presents challenges near the band edge because λ(E)υ(E)τm(E), and the velocity approaches zero as the scattering time becomes infinite. As discussed in the supplementary material, the use of a fine grid near the band edge is essential. Finally, Fig. 2(d) shows the numerically computed transport distribution, Ξ(E), according to Eq. (2). The use of a fine energy grid near the band edge produces a smooth Ξ(E), which is then used in (1) to evaluate the thermoelectric transport coefficients.

The L vs. S characteristic for Si is computed by using the M(E) and λ(E) presented in Fig. 2 and calculated from Eq. (1). Three cases are shown in Fig. 3: (i) the conduction band (symbols), (ii) the valence band (solid line), and (iii) the analytically calculated parabolic band reference (dashed line). The numerical results for both the conduction and valence bands of Si are very close to those obtained from the analytical expressions. This might be expected for the conduction band, which consists of six parabolic ellipsoids, but it is not so apparent why the complex valence band behaves like a simple, parabolic band. The reason is provided by Mecholsky et al., who show that the warped valence band is mathematically equivalent to a set of ellipsoidal bands.26 The observation from Fig. 2(d) that the transport distribution is nearly linear with energy over most of the Fermi window, indicates that a parabolic band is a good approximation. The numerical calculations were done at 300 K, but the results are only weakly dependent on temperature (see supplementary material). (For the analytical model, the L vs. S characteristic is independent of temperature.) Note that bipolar effects are not considered here, but are later discussed in Sec. IV.

FIG. 3.

The L vs. S characteristic computed numerically for Si. Symbols: the conduction band with DOS scattering. Solid line: the valence band with DOS scattering. Dashed line: the analytically calculated parabolic band reference.

FIG. 3.

The L vs. S characteristic computed numerically for Si. Symbols: the conduction band with DOS scattering. Solid line: the valence band with DOS scattering. Dashed line: the analytically calculated parabolic band reference.

Close modal

In this case, we evaluate Eq. (1) assuming a DFT-generated band structure for Bi2Te3. Figure 4 shows the density-of-states computed from the numerical band structure along with the corresponding distribution of channels, M(E), energy-dependent mean-free-path, λ(E), and transport function, Ξ(E).

FIG. 4.

Key quantities computed from the numerical band structure of Bi2Te3. (a) The density-of-states. (b) The distribution of channels, M(E). (c) The energy-dependent mean-free-path, λ(E). (d) The transport distribution, Ξ(E). The Fermi window is shown as a dashed line centered at the Fermi level and indicates the relevant range of energies. The Fermi level is set at the valence band edge in these plots. All quantities were extracted from a DFT-generated band structure using Quantum Espresso (see supplementary material). The bandgap is 0.18 eV with the midgap at E = 0 eV.

FIG. 4.

Key quantities computed from the numerical band structure of Bi2Te3. (a) The density-of-states. (b) The distribution of channels, M(E). (c) The energy-dependent mean-free-path, λ(E). (d) The transport distribution, Ξ(E). The Fermi window is shown as a dashed line centered at the Fermi level and indicates the relevant range of energies. The Fermi level is set at the valence band edge in these plots. All quantities were extracted from a DFT-generated band structure using Quantum Espresso (see supplementary material). The bandgap is 0.18 eV with the midgap at E = 0 eV.

Close modal

The non-parabolic energy bands are evident in the density-of-states shown in Fig. 4(a). Nevertheless, the distribution of channels, M(E), shown in Fig. 4(b), is roughly linear with energy over the Fermi window. As shown in Fig. 4(c), the drop in the magnitude of λ(E) for energies away from the band edge is much stronger for the valence band than for the conduction band. Finally, Fig. 4(d) shows the numerically computed transport distribution, Ξ(E), which is proportional to the product of M(E) and λ(E). This figure shows that for the conduction band, the increase in Ξ(E) above the conduction band minimum is slightly sublinear. For the valence band, however, the increase in Ξ(E) for energies below EV is distinctly sublinear. The analytical results presented earlier showed that a sub-linear increase in Ξ(E) with energy leads to L<2 for non-degenerate conditions. Based on the results shown in Fig. 4(d), we might expect L for the conduction band to be slightly less than 2 for non-degenerate conditions, but L for the valence band is expected to be distinctly less than 2 for non-degenerate conditions.

The L vs. S characteristic for Bi2Te3 as computed by using the Ξ(E) shown in Fig. 4 with Eq. (1) is presented in Fig. 5. Three cases are shown: (i) the conduction band (symbols), (ii) the valence band (solid line), and (iii) the analytically calculated parabolic band reference (dashed line). As expected, the numerical result for the conduction band is somewhat less than the parabolic band result. For the valence band, however, L is distinctly less than 2 for non-degenerate conditions. As discussed above, the behavior observed in Fig. 5 can be explained in terms of the transport function.

FIG. 5.

The L vs. S characteristic computed numerically for Bi2Te3. Symbols: the conduction band with DOS scattering. Solid line: the valence band with DOS scattering. Dashed line: the analytically calculated parabolic band reference.

FIG. 5.

The L vs. S characteristic computed numerically for Bi2Te3. Symbols: the conduction band with DOS scattering. Solid line: the valence band with DOS scattering. Dashed line: the analytically calculated parabolic band reference.

Close modal

Four issues will be discussed in this section. The first is the sensitivity of the results to the band structure. The second issue is the sensitivity of the results to scattering processes, and the third issue concerns the influence of bipolar effects, which have been neglected in the analyses discussed above. Finally, numerical issues are briefly discussed here and at more length in the supplementary material.

The differences between the analytically and numerically calculated L vs. S characteristics are the result of both band structure and scattering. In this section, the band structure and scattering effects are examined for Si and Bi2Te3. Figure 6(a) shows the L vs. S characteristic for hole-only conduction in the valence bands of Si and compares different treatments of scattering with the analytical results that assume parabolic bands and a constant MFP. When the MFP is computed by assuming that the scattering rate follows the density of states, the small difference between the analytical result and the numerical result is mainly due to the carrier velocity averages as given by (3). The use of a constant MFP approximation (CMFPA) also produces results that are close to the analytical results. The warped valence bands produce a linear M(E), but the energy-dependent angle-averaged velocity differs from that of a parabolic band. It may be surprising that the numerical results agree so closely with the analytical results because the valence bands of Si are complicated, but as discussed by Mecholsky et al., the warped valence bands are mathematically equivalent to a set of ellipsoidal bands.26 Finally, we consider a commonly used approximation when calculating thermoelectric properties—the constant relaxation time approximation (CRTA). This assumption simplifies the analysis but is non-physical. Figure 6(a) shows that the results for the CRTA are much different from the others. The use of a constant relaxation time is non-physical—even for parabolic bands.

FIG. 6.

The L vs. S characteristic for (a) silicon and (b) Bi2Te3 with constant mean free path approximation (“CMFPA,” crosses), constant relaxation time approximation (“CRTA,” solid circles), and scattering rate that is proportional to the DOS (“DOS,” solid). In both cases, a reference calculation with a single parabolic band and constant mean free path (“Analytical,” dashed) is shown. These calculations treat the valence band alone and do not consider bipolar conduction.

FIG. 6.

The L vs. S characteristic for (a) silicon and (b) Bi2Te3 with constant mean free path approximation (“CMFPA,” crosses), constant relaxation time approximation (“CRTA,” solid circles), and scattering rate that is proportional to the DOS (“DOS,” solid). In both cases, a reference calculation with a single parabolic band and constant mean free path (“Analytical,” dashed) is shown. These calculations treat the valence band alone and do not consider bipolar conduction.

Close modal

Next, we turn to Fig. 6(b), the case of Bi2Te3, which has a band structure that is more complicated and significantly non-parabolic as compared to silicon. In this case, the results for a numerical band structure with a constant MFP or a constant scattering time are close to the analytical, parabolic band model. Significantly lower Lorenz numbers are, however, obtained when DOS scattering is assumed. DOS scattering is thought to be the most physical of the three scattering models assumed.

Although the technique described in this paper can make use of arbitrary numerical tables of band structure and energy-dependent scattering times, we have illustrated it by assuming that the scattering rate follows the DOS. This has been found to be a good approximation for non-polar semiconductors (e.g., Refs. 21 and 9) but its validity for complex thermoelectric materials is less clear. Techniques to rigorously evaluate electron-phonon scattering rates in complex materials are available11–14 and can address this question. Calculations for p-type SnSe are shown in Fig. 7 below. (For a description of the methods used, see supplementary material.) Each point in the plot represents the scattering rate for a particular k-state with all relevant acoustic and optical phonon modes included. The polar phonon interactions are screened according to the Debye length with the Fermi level being near the valence band edge. The results show that the scattering rate is k-dependent with the variation with k increasing as the energy increases. A comparison of the scattering rate with the momentum relaxation rate shows that the two are approximately the same. Without screened polar interactions, the scattering rate is significantly larger than the momentum relaxation rate. The solid line is the computed DOS for SnSe. These results show that even for a polar material (with a relatively high carrier concentration; 4 × 1019 cm−3 in this case) with a complex band structure, the assumption that the momentum relaxation time follows the DOS is reasonable and much better than the often-assumed constant scattering time. In general, however, a detailed examination of material-specific scattering processes will be needed to provide accurate scattering times for the approach described in this paper.

FIG. 7.

Computed electron-phonon scattering rate for p-type SnSe. (See supplementary material for a description of the techniques used.) Each point represents the momentum scattering rate for a particular k-state. The solid line is the computed DOS for SnSe. The Fermi level is 64.6 meV below the top of the valence band, and corresponds to a hole concentration of 4 × 1019 cm−3.

FIG. 7.

Computed electron-phonon scattering rate for p-type SnSe. (See supplementary material for a description of the techniques used.) Each point represents the momentum scattering rate for a particular k-state. The solid line is the computed DOS for SnSe. The Fermi level is 64.6 meV below the top of the valence band, and corresponds to a hole concentration of 4 × 1019 cm−3.

Close modal

The study so far has considered L vs. S when only one type of carrier is present. This is acceptable at 300 K for Si, which has a bandgap of 1.1 eV, but for many common TE materials, bipolar effects can become important at the temperatures of interest (e.g., Bi2Te3 at room temperature with EG0.10.2eV). Computing an L vs. S characteristic in the presence of bipolar conduction adds uncertainties because in addition to knowledge of the majority carrier transport properties and bandgap, minority carrier properties are also needed (e.g., the MFP and effective mass in a parabolic band model or the minority carrier band structure and scattering rate in a full, numerical model). We can illustrate the considerations by looking at Bi2Te3 for which considerable experimental data is available. Using available data for S vs. σ for n-Bi2Te3 and p-Bi2Te3,28 we perform full numerical simulations assuming DOS scattering and adjust the electron-phonon coupling term in (4) to match the bipolar experimental data. As shown in the supplementary material, an excellent fit is obtained if a bandgap of EG=0.18 eV is assumed. Figure 8, shows the resulting L vs. S characteristic for Bi2Te3 at 300 K calculated with and without bipolar effects included.

FIG. 8.

The L vs. S characteristic for Bi2Te3 with and without bipolar conduction included (“Bipolar,” symbols) and without bipolar conduction (“Holes,” solid line). DOS scattering is assumed in both cases and the electron-phonon coupling constants are evaluated as discussed in the supplementary material.

FIG. 8.

The L vs. S characteristic for Bi2Te3 with and without bipolar conduction included (“Bipolar,” symbols) and without bipolar conduction (“Holes,” solid line). DOS scattering is assumed in both cases and the electron-phonon coupling constants are evaluated as discussed in the supplementary material.

Close modal

To understand the shape of the bipolar curve in Fig. 8, consider the movement of the Fermi level beginning from inside the valence band and moving toward the middle of the band gap. Deep inside the valence bands, bipolar effects due to conduction band electrons are negligible. In this case, L approaches the degenerate limit of π2/3 and S is small. This is the left part of the plot where S is small and both bipolar and unipolar L vs. S characteristics converge. As the Fermi level moves above the valence maximum and toward midgap, S increases and L decreases, until bipolar effects become significant, L begins to increase, and the bipolar characteristic diverges from the unipolar characteristic. For S200μV/K, we would deduce L1.75 from the unipolar characteristics and L2.1 from the bipolar curve. Unless carefully modeled, bipolar effects will introduce a significant uncertainty in the estimation of the Lorenz number from a measured Seebeck coefficient. An accurately computed L vs. S. characteristic could provide a good way to deduce the Lorenz number from the measured Seebeck coefficient, but the mapping of L to S is greatly complicated when bipolar effects are present. Accurate deduction of the Lorenz number would require careful modeling of bipolar effects or their suppression.

As illustrated by the examples presented, accurate calculations of L vs. S require a careful treatment of the band structure and energy-dependent scattering time. Numerical issues are also important. A good resolution of the band structure and scattering rates near the band edges is essential. Even with a highly refined k-mesh, however, obtaining accurate results for quantities near the band edge can still be challenging. This can be seen in Fig. 2(c), the MFP calculated numerically for silicon. The MFP shown in Fig. 2(c) is very close to the band edge where the scattering rate, which is proportional to the DOS, rapidly approaches zero. The MFP is proportional to velocity times the scattering time (the inverse of the scattering rate). Near the band edge, the velocity approaches zero and the scattering time approaches infinity, which presents numerical challenges. The results displayed in Fig. 2(c) used a refined k-mesh produced with an iterative technique that is described in the supplementary material. The numerical fluctuations in the MFP path are however, rather small, and what matters finally is the transport function, which is the product of MFP and number of conducting channels. As shown in Fig. 2(d), the transport function is smooth near the band edge. Finally, we note that when an exactly parabolic band structure is used in the numerical procedure, the results are in close agreement with the analytical results, which indicates that the iterative k-grid refinement procedure described in the supplementary material is effective.

In this paper, we presented a technique to compute the Lorenz number for complex thermoelectric materials. The method makes use of DFT-calculated band structures and assumes acoustic deformation potential scattering (or more generally, a scattering rate that follows the density-of-states). The accuracy of the numerical method was established by comparing the computed L vs. S characteristics for a material (Si) with nearly parabolic energy bands to analytical results that assume parabolic bands and a constant mean-free-path. We also showed that when the transport distribution increases less slowly than linearly with energy away from the band edge, then non-degenerate Lorenz numbers less than 2(kB/q)2 result.

All thermoelectric transport coefficients are determined by the transport function;10 the Lorenz number is determined by the shape, but not the magnitude of the transport distribution. It is sensitive, therefore, to both the band structure and to energy dependent scattering. It is easy to show mathematically from Eq. (10b) that a delta-function transport distribution gives L=0. The physical reason is clear. If there is a single channel, then when we open-circuit it to measure κe, no electrons flow. Since there is no flow of electrons, there can be no flow of heat. If a transport distribution decreases rapidly with energy away from the band edge, then the transport distribution approximates a delta-function at the band edge, and low Lorenz numbers can be expected. More generally, we showed that whenever the transport function increases less than linearly with energy, a non-degenerate Lorenz number less than the parabolic band with the ADP scattering result of L=2 can be expected. Note also that the influence of energetically offset bands on L as discussed by Thesberg et al.29 can also be understood in terms of the transport distribution.

It has been shown that when parabolic bands with acoustic deformation potential scattering is a good approximation, then the L vs. S characteristics provide a convenient way to deduce the Lorenz number from a measured Seebeck coefficient.7 As illustrated in this paper, the L vs. S characteristic is material-specific. We presented a procedure to compute material-specific Lorenz numbers and L vs. S. characteristics using the best available band structures and scattering rates. We illustrated this general procedure with some specific examples.

The ability to compute a material-specific L vs. S characteristic with high confidence would be useful in two ways. First, it could be used to deduce L from the more readily measured S, which would be useful for analyzing measurements of total thermal conductivity in order to deduce the lattice component. Second, it could be used to identify materials that give a low L for a given S. Such materials might be useful for increasing the thermoelectric figure of merit, zT.1,2 To produce such high confidence calculations requires high accuracy band structures (especially near the band edges), and high-accuracy energy-dependent scattering times and mean-free-paths. Decades of work on band structure calculations suggests that the first challenge can be solved, but computing the energy-dependent mean-free-path for arbitrary, complex thermoelectric materials is still not routine. It will be important to test the assumption that acoustic deformation potential scattering is generally the dominant scattering mechanism in thermoelectric materials.24 First-principles scattering calculations may provide the capabilities needed.11–14 

See supplementary material for discussions on the DFT calculation of silicon and Bi2Te3 band structures, the numerical advantage of the Landauer approach and distribution of mode calculation, the numerical sensitivity of Lorenz number calculation, the temperature dependence in L vs. S, the effect of bipolar conduction on L vs. S, and the SnSe scattering rate calculation from first-principles.

This work was partially supported by the Defense Advanced Research Projects Agency (Award No. HR0011-15-2-0037) and JHU/APL MATRIX funding (Award No. HR0011-16-C-0011). J.M. acknowledges support from the NSERC (Discovery Grant RGPIN-2016-04881) and Compute Canada.

1.
R. W.
McKinney
,
P.
Gorai
,
V.
Stevanović
, and
E. S.
Toberer
, “
Search for new thermoelectric materials with low Lorenz number
,”
J. Mater. Chem. A
5
,
17302
17311
(
2017
).
2.
A.
Casian
, “
Violation of the Wiedemann-Franz law in quasi-one-dimensional organic crystals
,”
Phys. Rev. B
81
,
155415
(
2010
).
3.
G. S.
Kumar
,
G.
Prasad
, and
R. O.
Pohl
, “
Experimental determinations of the Lorenz number
,”
J. Mater. Sci.
28
,
4261
4272
(
1993
).
4.
R.
Venkatasubramanian
, “
Lattice thermal conductivity reduction and phonon localization like behavior in superlattice structures
,”
Phys. Rev. B
61
,
3091
3097
(
2000
).
5.
A. F.
May
,
E. S.
Toberer
,
A.
Saramat
, and
G. J.
Snyder
, “
Characterization and analysis of thermoelectric transport in n-type Ba8Ga16−xGe30+x
,”
Phys. Rev. B
80
,
125205
(
2009
).
6.
K. C.
Lukas
,
W. S.
Liu
,
G.
Joshi
,
M.
Zebarjadi
,
M. S.
Dresselhaus
,
Z. F.
Ren
 et al, “
Experimental determination of the Lorenz number in Cu0.01Bi2Te2.7Se0.3 and Bi0.88Sb0.12
,”
Phys. Rev. B
85
,
205410
(
2012
).
7.
H.-S.
Kim
,
Z. M.
Gibbs
,
Y.
Tang
,
H.
Wang
, and
G. J.
Snyder
, “
Characterization of Lorenz number with Seebeck coefficient measurement
,”
APL Mater.
3
,
041506
(
2015
).
8.
E.
Flage-Larsen
and
Ø.
Prytz
, “
The Lorenz function: Its properties at optimum thermoelectric figure-of-merit
,”
Appl. Phys. Lett.
99
,
202108
(
2011
).
9.
E.
Witkoske
,
X. F.
Wang
,
M.
Lundstrom
,
V.
Askarpour
, and
J.
Maassen
, “
Thermoelectric band engineering: The role of carrier scattering
,”
J. Appl. Phys.
122
,
175102
(
2017
).
10.
G. D.
Mahan
and
J. O.
Sofo
, “
The best thermoelectric
,”
Proc. Natl. Acad. Sci. U.S.A.
93
,
7436
7439
(
1996
).
11.
B.
Qiu
,
Z.
Tian
,
A.
Vallabhaneni
,
B.
Liao
,
J. M.
Mendoza
,
O. D.
Restrepo
 et al, “
First-principles simulation of electron mean-free-path spectra and thermoelectric properties in silicon
,”
Europhys. Lett.
109
,
57006
(
2015
).
12.
F.
Giustino
,
M. L.
Cohen
, and
S. G.
Louie
, “
Electron-phonon interaction using Wannier functions
,”
Phys. Rev. B
76
,
165108
(
2007
).
13.
S.
Poncé
,
E. R.
Margine
,
C.
Verdi
, and
F.
Giustino
, “
EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions
,”
Comput. Phys. Commun.
209
,
116
133
(
2016
).
14.
B.
Liao
,
J.
Zhou
,
B.
Qiu
,
M. S.
Dresselhaus
, and
G.
Chen
, “
Ab initio study of electron-phonon interaction in phosphorene
,”
Phys. Rev. B
91
,
235419
(
2015
).
15.
N. W.
Ashcroft
and
N. D.
Mermin
,
Solid-State Physics
(
Saunders
,
Philadelphia
,
1976
).
16.
W. E.
Bies
,
R. J.
Radtke
,
H.
Ehrenreich
, and
E.
Runge
, “
Thermoelectric properties of anisotropic semiconductors
,”
Phys. Rev. B
65
,
085208
(
2002
).
17.
C.
Kittel
,
Introduction to Solid State Physics
(
Wiley
,
1986
).
18.
G. D.
Mahan
and
M.
Bartkowiak
, “
Wiedemann–Franz law at boundaries
,”
Appl. Phys. Lett.
74
,
953
954
(
1999
).
19.
C.
Jeong
,
R.
Kim
,
M.
Luisier
,
S.
Datta
, and
M.
Lundstrom
, “
On Landauer versus Boltzmann and full band versus effective mass evaluation of thermoelectric transport coefficients
,”
J. Appl. Phys.
107
,
023707
(
2010
).
20.
See https://nanohub.org/groups/needs/lantrap for L
anTraP 2.0
, Purdue University, 2017.
21.
M. V.
Fischetti
, “
Monte Carlo simulation of transport in technologically significant semiconductors of the diamond and zinc-blende structures. I. Homogeneous transport
,”
IEEE Trans. Electron Devices
38
,
634
649
(
1991
).
22.
M.
Lundstrom
and
C.
Jeong
,
Near-Equilibrium Transport: Fundamentals and Applications
(
World Scientific Publishing Company
,
2013
).
23.
J. S.
Blakemore
, “
Approximations for Fermi-Dirac integrals, especially the function used to describe electron density in a semiconductor
,”
Solid-State Electron.
25
,
1067
1076
(
1982
).
24.
H. J.
Goldsmid
,
The Physics of Thermoelectric Energy Conversion
(
Morgan & Claypool Publishers
,
2017
).
25.
J.
Singh
,
Physics of Semiconductors and Their Heterostructures
(
McGraw-Hill College
,
1992
).
26.
N. A.
Mecholsky
,
L.
Resca
,
I. L.
Pegg
, and
M.
Fornari
, “
Theory of band warping and its effects on thermoelectronic transport properties
,”
Phys. Rev. B
89
,
155131
(
2014
).
27.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
 et al, “
QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
28.
M. T.
Pettes
,
J.
Maassen
,
I.
Jo
,
M. S.
Lundstrom
, and
L.
Shi
, “
Effects of surface band bending and scattering on thermoelectric transport in suspended bismuth telluride nanoplates
,”
Nano Lett.
13
,
5316
5322
(
2013
).
29.
M.
Thesberg
,
H.
Kosina
, and
N.
Neophytou
, “
On the Lorenz number of multiband materials
,”
Phys. Rev. B
95
,
125206
(
2017
).

Supplementary Material