Lattice thermal transport in group II-alloyed PbTe

Thermoelectric (TE) devices, which are capable of converting waste heat into electricity, are ideal alternative renewable energy technologies to overcome limited fossil fuel resources and environmental challenges.¹ Thermoelectric energy conversion efficiency is characterized by the dimensionless figure of merit ZT = S²σT/(κ_e + κ_l), where S, σ, and T are the Seebeck coefficient, electrical conductivity, and temperature, and κ_e and κ_l are the electronic and lattice thermal conductivities, respectively. ZT is generally optimized by maximizing the thermoelectric power factor S²σ and minimizing κ_l.

PbTe-based TE materials are among the highest performing, partly because of strong inherent phonon anharmonicity, leading to low κ_l.² To further enhance ZT, simultaneously reducing κ_l and improving the power factor can be achieved by introducing anionic and cationic dopants. Recent experiments show that p-type PbTe alloyed with group II elements (Mg, Ca, Sr, and Ba) achieved ZT above 1.5.^3–7 It is found in experiments and theoretical calculations that, with group II dopants, electronic properties are improved owing to convergence of multiple valence bands and bandgap widening.^4,5 With respect to lattice heat transport, significant reduction of κ_l is observed, which was attributed to the all-scale hierarchical architecture-induced phonon scattering due to the presence of solid-solution point defects, nanoscale precipitates, and grain boundaries.^2,3 However, a fundamental understanding of the roles of various phonon scattering mechanisms in reducing κ_l is still lacking. Previous first-principles investigations mainly focused on pristine PbTe and its anionic alloys (PbSe_xTe_1–x),^8,9 showing the large intrinsic anharmonicity and the importance of optical phonons which provide strong scattering channels for acoustic phonons. Recently, temperature-induced phonon renormalization was further investigated using quasiharmonic approximation and temperature-dependent effective potential techniques.^10,11 However, lattice thermal transport for group II alloyed-PbTe is not well understood, with existing theoretical modeling limited to the phenomenological Debye-Callaway model.^5,12,13 Generally, various scattering mechanisms are present in the alloyed phases including intrinsic anharmonic scattering and extrinsic scattering from precipitates, strains, and point defects.¹³ It is fundamentally important to understand the role of various phonon scattering mechanisms in order to explore the possibility to further reduce κ_l. Motivated by this point, we use first-principles methods and the recently-developed compressive sensing lattice dynamics (CSLD)¹⁴ approach to investigate the lattice heat transport of PbTe, MTe, and Pb_0.94M_0.06Te (M = Mg, Ca, Sr, and Ba), particularly focusing on the effects of strain and alloying on κ_l.

The lattice thermal conductivity tensor (⁠$κlαβ$⁠) was calculated by summing over contributions from phonon modes in the first Brillouin zone under the relaxation time approximation¹⁵ 

where N, Ω, and f₀ are the number of phonon wave vectors, the volume of the primitive cell, and the Bose-Einstein distribution function, respectively. For each mode λ, $ωλ, vλα$⁠, and τ_λ represent the phonon frequency, the velocity along the α direction, and the mode lifetime, respectively. The lifetime of each mode was calculated using Fermi's golden rule by treating 3rd-order interatomic force constants (IFCs) as perturbations to harmonic phonons,¹⁵ and the linearized Boltzmann transport equation was solved in an iterative manner to take into account non-equilibrium phonon distributions.^16–20 However, the above formula cannot be directly applied to group II-alloyed PbTe, which are essentially disordered solids and intrinsically difficult to model from first principles. Explicitly treating the alloying effects including both IFC disorder and mass disorder at the first-principles level requires simulating large supercells with significantly reduced symmetry, making it computationally intractable with our currently available resources. To mitigate this limitation, we adopted an effective approach, the virtual crystal approximation (VCA),²¹ where the interatomic interactions of alloys are obtained by averaging those from pristine phases according to a given composition. Additional effects of alloying on the phonon spectra and thermal conductivity were treated by taking into account extrinsic phonon scattering due to mass disorder,^22,23 treated as an external elastic scattering term in BTE. The associated phonon mode scattering rates are given in Tamura's theory by²³ 

where e^λ is the polarization vector of phonon mode λ and x_s(a) and M_s(a) are the concentration and atomic mass of the s-th species on the atomic site indexed by a. According to a study investigating the validity of VCA performed by Larkin and McGaughey,²⁴ Tamura's theory could break down for phonon modes with diffusivities (D) predicted to be below the high-scattering limit $DHS=13vsl$⁠, where v_s and l are the sound velocity and interatomic distance, respectively. To examine the validity of applying Tamura's theory to our systems, we performed mode-by-mode analysis for Pb_0.94Mg_0.06Te, which displays the strongest scattering rates due to large mass contrast. We found that very few modes exceed the high-scattering limit and replacing the diffusivities of these modes with D_HS leads to a change in κ_l of less than 1%. Therefore, we conclude that Tamura's theory is likely valid here despite the neglect of higher-order terms. The ShengBTE package²⁵ was used to perform the iterative calculation of phonon lifetimes and κ_l with 24 × 24 × 24 q-point meshes (see supplementary material for detailed discussion of IFC fitting and computational parameters).

Owing to the general low solubility of M in the PbTe matrix, MTe is experimentally found to precipitate out of the matrix, forming both coherent and incoherent nanoscale precipitates.^3–6,26 Since κ_l reduction is suggested to originate from these precipitates, it is useful to compare the κ_l of pristine PbTe and MTe. To shed light on this point, phonon dispersion and κ_l were computed and analyzed. The phonon dispersions shown in Fig. 1(a) share similar features and exhibit enhanced phonon softening with increasing atomic mass from Mg to Ba. A noticeable difference is observed in MgTe, which shows extra softening of acoustic modes near the X and L points in the first Brillouin zone. Compared to the phonon dispersion of PbTe in Fig. 2(a), one important feature of MTe is the absence of low-lying transverse optical (TO) phonons, which are found to be essential to scatter acoustic phonons,⁸ thus making their κ_l higher than that of PbTe. The computed values of κ_l for MgTe, CaTe, SrTe, BaTe [Fig. 1(b)], and PbTe [Fig. 2(b)] are 3.0, 8.5, 10.5, 10.2, and 3.3 W/mK at room temperature using the fully relaxed lattice parameter at 0 K. The value for PbTe is higher than previously reported computed values (1.9 and 2.1 W/mK).^8,10 The discrepancy may be attributed to two factors: (1) the κ_l of PbTe is extremely sensitive to the lattice parameter (we will further illustrate this point by considering the strain effect below) and therefore can be heavily influenced by the adopted xc functional; (2) κ_l is potentially underestimated under the single mode relaxation time approximation (SMRTA) adopted by the previous studies, which, according to our calculations, is found to give a value of κ_l lower by 0.6 W/mK compared to the more accurate iterative solution of the BTE. Consistent with the difference found in phonon dispersion between MTe and PbTe, κ_l is higher for CaTe, SrTe, and BaTe than for PbTe at 300 K. However, MgTe, despite having the lightest element, exhibits the lowest κ_l among MTe, with a value comparable to that of PbTe. Since the two main factors that determine κ_l are lifetimes and group velocities, we show the computed mode-dependent lifetimes in Fig. 1(c) and group velocities in Fig. S5 in the supplementary material. We see that although MgTe has overall higher group velocities, its much shorter lifetimes lead to significantly reduced κ_l. We find that the κ_l of BaTe computed using lifetimes from MgTe has an even lower value of 1.4 W/mK. To further shed light on the anomalously low κ_l of rock-salt MgTe, we also computed the κ_l of the more stable zinc-blende phase of MgTe and found a value of 14.1 W/mK at 300 K. Detailed analysis shows that although the magnitudes of 3rd-order IFCs are reduced with the structure change from zinc blende to rock salt, the allowed three-phonon process is significantly enhanced (nearly twice in the absorption process) with the more suppressed optical modes in the rock-salt structure, giving rise to anomalously low κ_l. The fact that all MTe compounds have larger or similar κ_l compared with PbTe indicates that extrinsic phonon scattering from interfaces, point defects, grain boundaries, and precipitates should be examined to explain the reduced κ_l of group II-alloyed PbTe.

Aside from nanoscale precipitates, the size and strain fluctuations caused by point defects in the solid solution could also play an important role in determining the overall κ_l in alloyed PbTe. The fully relaxed lattice parameters of MgTe, CaTe, SrTe, and BaTe are 5.971, 6.400, 6.723, and 7.093 Å, respectively. When PbTe (fully relaxed lattice parameter of 6.441 Å) forms solid solutions with 6% (12%) MTe, its lattice parameters thus could vary roughly from −0.5% (−1.0%) to 0.5% (1.0%). It was shown in quasiharmonic approximation calculations¹⁰ that κ_l is very sensitive to the PbTe lattice parameter. However, the effects of static strain, specifically caused by alloying with group-II elements, on the κ_l of PbTe were not taken into account in previous studies,^5,12,13 motivating us to investigate its effect on phonon dispersion and κ_l. Figure 2(a) shows the phonon dispersions of PbTe calculated with −1%, 0%, and 1% lattice expansion. One striking feature is the presence of strong softening of optical phonon modes with lattice expansion. Particularly, TO modes near the Γ point fall deep into the acoustic region when the lattice is slightly expanded. The increased overlap between TO and acoustic modes increases the scattering phase space,¹¹ giving rise to further reduction in κ_l. The corresponding κ_l plotted in Fig. 2(b) shows that κ_l can be increased from 3.3 to 3.8 (4.8) W/mK with 0.5 (1.0)% lattice contraction and decreased to 2.5 (1.8) W/mK with 0.5 (1.0)% lattice expansion. Since the relaxed lattice parameter of 6.441 Å is smaller than the experimental value of 6.462 Å at 300 K,²⁷ we computed κ_l using the experimental lattice parameter and found a value of 2.7 W/mK (2.1 W/mK under SMRTA) at 300 K. Our results compare very well with Ref. 10 which used the same xc functional, and our higher value is from solving BTE in an iterative (more accurate) manner. Compared to experiments, our result is close to the measurement performed on single-crystal PbTe of 2.4 W/mK.²⁸ The slight overestimation can be attributed to (1) overestimation of the energy of low-lying TO optical phonon modes (see Fig. S2 in the supplementary material) and (2) neglect of higher-than-third-order phonon-phonon interactions. We note that experimental measurements on polycrystalline samples yield a yet lower κ_l of 2.0 W/mK²⁹ due to additional phonon scattering from defects and grain boundaries. Despite the overestimation presented in pristine PbTe, our result is close to a κ_l of 2.8 W/mK in Na-doped PbTe, which allows us to further investigate the additional alloying effects from group II elements. Our results also reveal that the strain effect alone can significantly influence the intrinsic anharmonicity of the PbTe matrix and different effects should be expected for different group II elements. By examining the phonon lifetime in Fig. 2(c), we find that lattice parameter increase tends to enhance the scattering rate, revealing that lifetime reduction plays a major role in reducing κ_l with lattice expansion. Note that here we only consider the static strain effect and ignore phonon scattering due to spatial strain fluctuations.

Considering the strain effect only, the κ_l of PbTe should increase with MgTe alloying and decrease with SrTe alloying. However, experimentally κ_l decreases with both group II alloys, indicating that additional scattering mechanisms may dominate over the strain effect. Another important factor introduced by alloying elements is mass fluctuations in crystals. By neglecting the change and disorder of IFCs induced by MTe and assuming the random distribution of M atoms over Pb sites, we examine the phonon scattering due to mass disorder, i.e., alloy scattering. To better illustrate the effect of mass disorder scattering, here we use the experimental lattice parameter of PbTe which shows better agreement between theoretical and experimental κ_l. The computed κ_l of PbTe alloyed with 6% MTe is shown in Fig. 3(a). We find that the largest κ_l reduction can be achieved through alloying with Mg which has the largest mass contrast, while Ba with the most similar mass to Pb reduces κ_l much less effectively. This is due to the fact that mass disorder scattering rates are proportional to mass contrast and increase with the increased relative alloying concentration. Further considering the strain effect on top of mass disorder scattering, as shown in Fig. 3(b) and Tables S1 and S2 in the supplementary material, we find that κ_l changes from 1.39, 1.49, 1.84, and 2.30 W/mK to 1.55, 1.50, 1.72, and 1.67 W/mK for Mg, Ca, Sr, and Ba, respectively, where the decrease/increase in κ_l is in line with the corresponding lattice expansion/contraction. The strain effect on Pb_0.94Ba_0.06Te which leads to significant reduction of κ_l indicates its importance when there is large lattice mismatch. Compared to experiments, our results of Pb_0.94Mg_0.06Te agree well with a recent experimental report that room temperature κ_l is significantly reduced from 2.83 to 1.74 W/mK for Pb_0.94Mg_0.06Te.^4,6 For Pb_0.94Ca_0.06Te, the experimental value of 1.33 W/mK also agrees with our theoretical value of 1.50 W/mK at 300 K.³⁰ For Pb_0.94Sr_0.06Te, the computed κ_l of 1.72 W/mK agrees reasonably with an experimentally reported value of about 2.0 W/mK,⁵ where the authors utilized the non-equilibrium processing technique to extend the solubility from less than 1 mol. % to about 5 mol. %. The slight underestimation in theory could be due to the requirement of full solubility in modeling mass disorder scattering. Our results suggest that mass disorder scattering plays the most important role in reducing κ_l for cationic dopants with large mass contrast (such as Mg, Ca, and Sr) near or slightly exceeding the solubility limit and relatively weakens the extra phonon scattering caused by static strain compared to pristine PbTe. Our results also suggest that extrinsic scattering, for example, induced by precipitated nanostructures and nanostructuring which limits the maximum mean free path, is not required to explain the experimental results when alloying group II atoms are mostly in solid solution form. For cationic dopants with reduced mass contrast such as Ba, however, κ_l reduction cannot be primarily attributed to the strain effect or mass disorder scattering. In contrast to Mg, Ca, and Sr, nanostructures are expected to be responsible for κ_l reduction with Ba. This can be inferred from the following: (1) the large difference between the computed κ_l of 1.67 W/mK for Pb_0.94Ba_0.06Te and the experimentally observed κ_l of about 1.20 W/mK in Pb_0.97Ba_0.03Te (Ref. 30) and (2) Ba has extremely low solubility (less than 0.5%) in PbTe.¹³ Comparison between phonon lifetimes in Fig. 3(c) confirms the strong scattering of optical modes through introducing mass disorder, while acoustic modes are less affected, suggesting that κ_l can be further reduced by creating mesoscale grain boundaries without significantly affecting electronic transport. Note that here we only compare to experimental results at 300 K since we only considered the room-temperature lattice parameter. We expect that our conclusion remains qualitatively sound at high temperatures, where lattice expansion and phonon renormalization (hardening of TO modes) could further affect κ_l.

We have applied a first-principles-based compressive sensing lattice dynamics approach for modeling lattice heat transport in PbTe and its solid solutions formed with group II elements. With the computed harmonic and anharmonic interatomic force constants, we modeled the lattice heat transport by iteratively solving the linearized Boltzmann transport equations with phonon lifetimes obtained from perturbation theory. Lattice thermal conductivities were computed and analyzed for PbTe and its alloys within the framework of virtual crystal approximation, with a particular emphasis on the effects of strain and mass disorder. The results show that the current approach is able to explain the experimentally observed reduction of lattice thermal conductivity, thus allowing further improvement of ZT through phonon engineering.

See supplementary material for detailed descriptions of the computational parameters, compressive sensing lattice dynamics method, and validation of computed phonon dispersions and lattice thermal conductivity of pristine PbTe using the experimental lattice constant.

This work was supported by the Midwest Integrated Center for Computational Materials (MICCoM) as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division (No. 5J-30161-0010A). The use of the Center for Nanoscale Materials, an Office of Science user facility, was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility, supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We gratefully acknowledge the computing resources provided on Blues, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.