Heat transport properties of PbTe_1−xSe_x alloys using equivariant graph neural network interatomic potential

Abstract

The suppression of heat transport in disordered crystals arises from a competition between mass fluctuations and bond disorder, but their relative contributions remain difficult to disentangle. We address this challenge using a machine-learned interatomic potential trained on ab initio data across the PbTe_1−xSe_x alloy series. Molecular dynamics simulations with the trained machine-learned interatomic potential reproduce experimental lattice thermal conductivities and density-functional theory phonon dispersions while enabling frequency-resolved analysis of heat transport from 300 to 800 K. We find a narrow window near 1.7–2.0 THz dominates heat transport across all compositions, where heat is primarily carried by mixed longitudinal acoustic and optical modes. Alloying dramatically reduces spectral diffusivity near 2 THz leading to the deterioration of thermal conductivity. For the parent compounds both the spectral diffusivity and overall thermal conductivity decrease at elevated temperatures. However, the thermal degradation is weaker for the mid-range composition (x = 0.5) due to the greater thermal occupation of vibrational modes and increased heat capacity. Alchemical simulations show that force–constant disorder, not mass contrast, plays the dominant role in this suppression. These results highlight the microscopic mechanisms underlying thermal transport breakdown in alloys and demonstrate how machine-learned interatomic potentials now offer a tractable path toward predictive, physics-rich thermal transport modeling in complex disordered solids.

---

New concepts

This work introduces a new approach for understanding thermal transport in disordered crystalline materials by combining machine-learned interatomic potentials with large-scale molecular dynamics simulations. Our method enables frequency-resolved analysis of heat transport with accuracy comparable to ab initio techniques, but at significantly reduced computational cost. We identify a narrow spectral range where mixed acoustic and optical vibrational modes carry most of the heat, and show how alloying and elevated temperatures disrupt this transport channel. Additionally, we find that force–constant disorder rather than mass disorder plays the dominant role in reducing thermal conductivity. This approach offers new insight into the microscopic origins of thermal conductivity degradation in complex, anharmonic materials and provides a predictive framework for modeling thermal transport in disordered solids. More broadly, it highlights the growing potential of machine-learned interatomic potentials to deliver detailed, physics-informed understanding of heat flow in materials where conventional methods are either insufficient or computationally prohibitive.

Introduction

Thermal conductivity is an important property for heat management of electronic devices and the performance of thermoelectric devices.^1–3 Anharmonic crystals, which deviate significantly from the ideal harmonic approximation of atomic vibrations, exhibit complex phonon–phonon interactions that significantly influence heat transport.⁴ Classical interatomic potentials often do not accurately capture the anharmonic effects, especially at elevated temperatures.⁵ Better descriptions of the complex atomic vibrations are required to more accurately calculate the thermal conductivity of anharmonic crystals. In this work, machine learned interatomic potentials are applied to thermal transport in thermoelectric alloys (PbTe_1−xSe_x) to both test the accuracy of the potential and to reveal new material design principles.

Machine learning interatomic potentials have emerged as powerful tools to learn complex potential energy surfaces from ab initio calculations. Recently, the neural network models NequIP⁶ and Allegro⁷ were introduced with built-in Euclidean group E(3) symmetry. These models train interatomic potentials with the symmetry of atomic systems built-in – the potential energy and other scalar quantities are invariant to symmetry group operations but vector quantities such as atomic forces undergo equivariant transformations. The models train more efficiently and outperform rotationally invariant graph neural network interatomic potentials, shallow neural networks, and kernel-based approaches.⁶ Previously, NequIP and Allegro have been used to investigate molecular crystals,⁸ surface adsorption,⁹ and the thermal conductivity of the α- and β-GeTe phase change materials.¹⁰

Rock salt PbTe and PbSe are strongly anharmonic compounds that remain the foundation of some of the best known thermoelectric materials.^11,12 In the absence of alloying or nanostructures, PbTe and PbSe have room temperature lattice thermal conductivities of ∼2 W m⁻¹ K⁻¹ which decline to ∼1 W m⁻¹ K⁻¹ at 700 K.¹³ Beyond their soft bonding and heavy masses, these low thermal conductivities emerge from large anharmonicity and soft transverse optical branches. The latter results in strong acoustic and optical phonon interaction and low thermal conductivity.^14–17

The lattice thermal conductivity of PbTe has been studied with a wide variety of computational methods. Classical interatomic potentials combined with equilibrium molecular dynamics and Green–Kubo linear response formulation have been used to study the lattice thermal conductivity of both pristine and defected PbTe.¹⁸ A major limitation of the used classical interatomic potentials is that they do not capture important features of the PbTe phonon dispersions such as soft longitudinal optical mode near the Γ-point.¹⁹ Several approaches based on first-principles have been used to study the lattice thermal conductivity of PbTe: density functional theory (DFT) combined with anharmonic lattice dynamics and Boltzmann transport theory,²⁰ DFT combined with the temperature-dependent effective potential technique (TDEP),²¹ and a phonon quasiparticle approach combining first-principles molecular dynamics and lattice dynamics.²² The role of quartic anharmonicity for the lattice thermal conductivity of PbTe has also been studied with DFT by combining the quasiharmonic approximation with temperature-induced anharmonic phonon renormalization.²³ More recently, machine-learning interatomic potentials trained with DFT data have been used to study the lattice thermal conductivity of PbTe. Fan et al. used DFT-trained neuroevolution potentials and homogeneous nonequilibrium molecular dynamics to study pristine PbTe,²⁴ while Qin et al. applied a DFT-trained deep neural network potential to study defected PbTe.²⁵ DFT, neuroevolution potentials, and homogeneous nonequilibrium molecular dynamics have also been used to investigate the lattice dynamics and thermal transport of PbTe under high pressure.²⁶ The Green–Kubo formulation, on the other hand, does not impose a temperature gradient that can introduce finite-size effects and boundary artifacts. Furthermore, all phonon modes, phonon anharmonicity, and finite temperature effects are inherently included, making it well-suited for systems with strong anharmonic effects such as PbTe_1−xSe_x alloys.

Experimentally, alloying between PbTe and PbSe is known to lower the lattice thermal conductivity by a factor of about 2 and leads to improved thermoelectric performance.^27,28 However, modeling this phenomena remains challenging. Previous works that modeled thermal transport in PbTe_1−xSe_x alloys used either a virtual crystal approach²⁹ or empirical mixing of the interatomic force constants.³⁰ While the lattice thermal conductivity's composition dependence is captured,^29–31 such methods rely on empirical models of the alloys and could not employ interatomic potentials with the accuracy of ab initio calculations.

In this paper, we investigate the thermal properties of PbTe_1−xSe_x alloys from 300 to 800 K using machine-learned interatomic potentials trained on density functional theory (DFT) calculations. We employ the Green–Kubo relation³² to study the compositional dependence of the lattice vibrations on the thermal conductivity. The Green–Kubo relation neither uses a virtual crystal approximation nor does it rely on approximations for phonon scattering and dispersion relations such as the Boltzmann Transport Equation (BTE). Instead, the lattice thermal conductivity is calculated using the fluctuation–dissipation theorem. The direct use of molecular dynamics simulations implies that all phonon modes, phonon anharmonicity, and finite temperature effects are inherently included.³³ Our machine learned interatomic potential was trained on alloys and pure compounds and is general to the composition and configurational space, i.e. decorations. It reproduces the phonon properties from DFT calculations at 0 K and the ab initio molecular dynamics at elevated temperatures. By studying the frequency-dependent heat transport across both alloy composition and temperature, a new view of how thermal transport occurs in PbTe-based alloys emerges.

Methods

The heat transport properties of PbTe_1−xSe_x alloys were calculated using machine learned interatomic potentials. The training set data was acquired from ab initio molecular dynamics (AIMD) simulations of alloys having diverse composition, temperature, and configuration, i.e., decoration or ordering. The machine learned interatomic potential was used to calculate the thermal transport and phonon properties. Each step is described in detail below.

Ab initio molecular dynamics

The reference data was generated from AIMD simulations on PbTe_1−xSe_x using CP2K.³⁴ To generate the initial structures, random configurations, i.e., decorations, of PbTe_1−xSe_x were generated for a 3 × 3 × 3 supercell of the 8-atom unit cell (Fm3̄m space group) using the program supercell.³⁵ In addition to the two pristine compositions, ten decorations were selected for each composition (x = 0.25, 0.5, and 0.75). The sampled decorations were verified to be diverse by checking the partial radial distribution function and potential energy predicted by M3Gnet.³⁶

The atomic positions and lattice vectors of the selected structures were relaxed to minimum energies using DFT and periodic boundary conditions. The lattice constants are shown in Table S1 (ESI†). The DFT calculations utilized the PBEsol functional,³⁷ which has been devised for obtaining better equilibrium properties of bulk solids and their surfaces than the local density approximation (LDA) or the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE).³⁸ PBEsol has been shown to more accurately predict the lattice parameter and phonon dispersion relations of PbTe and PbSe than PBE.³⁹ The calculations use the DZVP-basis sets and Goedecker–Teter–Hutter (GTH) pseudopotentials^40,41 on a 2 × 2 × 2 k-mesh. An example alloy is shown in Fig. 1.

Fig. 1 Example of a PbTe_1−xSe_x alloy with x = 0.25 and lattice vectors 3a = 3b = 3c. The Pb atoms are black, Te atoms are orange, and Se atoms are gold.

Next, AIMD was performed using the relaxed structures as the initial coordinates in the canonical ensemble with time-step of 1 fs and Nosé–Hoover thermostat with a time constant of 200 fs. The temperatures were 250, 300, 350, 500, 800, 1600, and 3000 K. The simulations were about 4 ps for 250 to 800 K. Higher temperatures were simulated for only about 500 fs. Additional AIMD simulations at 400, 600, and 700 K were performed to generate the test set. To minimize unnecessary computations, only the pristine structures and two decorations of each composition were performed at temperature other than 300 K.

The interatomic potential was trained using frames selected from the AIMD trajectories. The model had 5000 training frames, 500 validation frames, and 500 test frames. The random sampling of the frames was balanced across all the compositions (x = 0, 0.25, 0.5, 0.75, and 1). 1700 frames were selected from each alloy composition and 450 frames each were chosen from PbTe and PbSe. The selected frames were recalculated with stricter convergence criteria to obtain more accurate energies, forces, and stress tensors. Example input scripts and the training set are provided in the ESI.†

Interatomic potentials

The interatomic potentials were implemented and trained using Allegro.⁷ The potential energy is partitioned into atomic contributions:and the atomic forces and virial stresses are the derivatives of the energy with respect to the atomic positions and the cell dimensions, respectively, using automatic differentiation.⁷

The model was trained using a cutoff radius of 10 Å and a polynomial envelope function with p = 6. The model had a maximum rotation order, l_max, of 2 and o3 full type E(3)-symmetry for the equivariant features. The network contained three tensor product layers with 64 Bessels per basis. The two-body embedding multi-layer perceptron had 32 tensor features with four hidden layers with dimensions [128, 256, 512, 1024] and SiLU non-linearities.⁴² The latent multi-layer perceptron had three hidden layers of dimensions [1024, 1024, 1024] and SiLU non-linearities. The embedding multi-layer perceptron was built from a single linear layer. The final multi-layer perceptron (from Allegro latent space to edge energies) was a single layer with 128 nodes without nonlinearity. The training used the Adam optimizer with an initial learning rate of 0.001 which decayed by a factor of 0.5 if the validation loss did not improve for 3 epochs. The numerical precision was float64 and batch size was 2. The training minimizes the loss function:

where loss coefficients are λ_e = 0.01, λ_f = 1, λ_σ = 100.

Thermal conductivity

The lattice thermal conductivity along lattice direction α is calculated using the Green–Kubo method^2,32 aswhere k_B is the Boltzmann constant, T is temperature, and V is the volume. The term 〈J_α(0)·J_α(t)〉 is the heat flux autocorrelation function (HFACF). The heat flux iswhere v_i is the velocity of atom i. The first two terms on the right hand side are formed from the potential and kinetic energies per atom, E^k_i. The third term is the atomic stress, Ξ_i, which we implemented such thatThe term is obtained by automatic differentiation and was also used for the calculation of the virial tensor.^43–46

The molecular dynamics simulations were performed using the LAMMPS code⁴⁷ in the canonical ensemble with a time-step of 2 fs. The temperature was set with a Nosé–Hoover thermostat with damping equal to 0.2 ps. The structures were equilibrated for 1 ns and then the heat flux was recorded for at least 2 ns. To start, a random configuration of Se and Te atoms was placed on the rocksalt lattice. The structure was then relaxed before conducting the molecular dynamics simulation. To assess the degree of ergodicity and quantify error statistics, the heat flux was obtained for 10 independent molecular dynamics simulations at 300 K and 5 simulations at higher temperatures. The heat flux autocorrelation function was obtained over 50 ps time-windows. The lattice thermal conductivity was averaged along its directional components, i.e., (κ_x + κ_y + κ_z)/3. The spectral thermal conductivity was obtained through a Fourier transform of the heat flux autocorrelation function and its amplitude was scaled uniformly across frequencies to be consistent with the thermal conductivity calculated using Green Kubo. The calculated spectral thermal conductivity of GeTe is shown for comparison in Fig. S1 (ESI†).

Phonon properties

The phonon properties of PbTe and PbSe at 0 K were calculated within the harmonic approximation using Phonopy^48,49 and a supercell containing 216 atoms. The phonon properties at finite temperature were calculated using the fix/phonon package of LAMMPS.⁵⁰ Fix/phonon is based on fluctuation–dissipation theory and acquires the dynamical matrix during the molecular dynamics simulation and captures the finite temperature effects and anharmonicity. The dynamical matrix, D_kα,k′β(q), describes how a displacement of atom k′ in the β direction affects the restoring force on atom k in the α direction at a given wave vector q. The fix/phonon method obtains the force constant coefficients in reciprocal space as^50–53

Φ_kα,k′β(q) = k_BTG_kα,k′β⁻¹(q) (6)

where G are the Green's function coefficients measured from the instantaneous and averaged atomic positions of the kth atom, R and 〈R〉, such that

The dynamical matrix is then obtained by⁵⁰

D_kα,k′β(q) = (m_km_k′)^−1/2Φ_kα,k′β(q) (8)

where m are the atomic masses. The eigenvalues of D are the phonon frequencies at q. The phonon density of states, g, and phonon dispersion were extracted using the package phana.⁵⁰

The spectral heat capacity is

C(ω, T) = g(ω)c_v(ω) (9)

where c_v is the phonon mode heat capacity derived from the Bose–Einstein distribution aswith .

The lattice thermal conductivity within the Boltzmann transport formalism in the relaxation-time approximation is:

The thermal transport can be viewed in terms of the spectral diffusivity:

D(ω) = v²(ω)τ(ω) (12)

where the term v²τ can be calculated asThe phonon group velocity is v and τ is the phonon relaxation time. The phonon properties were calculated on a 3 × 3 × 3 supercell of a 64-atom unit cell. The finite temperature phonon dispersion of pristine PbTe and PbSe was obtained using the primitive cell.

Synthesis and measurement of PbTe_1−xSe_x

Un-doped PbTe_1−xSe_x alloys are known to have a carrier concentration far lower than optimal for thermoelectric performance; as such, the experimental literature is quite scarce for un-doped compositions. Herein, a series of bulk polycrystalline samples was prepared specifically for comparing with the above computation predictions.

Synthesis. High-purity elements Te (shots, 5N, 99.999%) and Se (shots, Alfa Aesar, 99.999%) were used without further purification. Pb (granules, Alfa Aesar, 99.99%) was acid-washed in a 1% HCl solution before use. Stoichiometric amounts of elements were weighed into glass ampoules (PbTe_1−xSe_x, x = 0, 0.10, 0.25, 0.50, 0.75, 0.90, 1.0). The ampoules were evacuated to 10⁻³ Torr and subsequently flame sealed. In a standing furnace, the ampoules were ramped over 18 hours to 1100 °C, intermittently agitated at 1100 °C for 8 hours, and then air quenched. The resulting ingots were ball milled for 3 × 15 minutes, and the powders were re-sealed in glass ampoules. The powders were annealed at 600 °C for 4 days and then further ground to <106 μm with a mortar and pestle. The powders were loaded into a graphite die and hot pressed under vacuum at 600 °C with a uniaxial pressure of 40 MPa. The resulting pellets were roughly 12.7 mm in diameter and 1.5 mm tall. The pellets were polished before measurement achieving near parallelism with less than ±0.005 μm thickness variation. All pellets reached >97% of their theoretical density.

Characterization. The thermal diffusivities (D) were measured on the Netzsch LFA 467 HyperFlash. The diffusivities were converted into thermal conductivities (κ) using the geometrically determined density (ρ) and the Dulong–Petit heat capacity (C_p), as shown in the following equation:

κ = DC_pρ (14)

The Lorenz number (L) was calculated over the full temperature range of PbTe using the following approximation:

where L is in units of 10⁻⁸ W Ω K⁻².⁵⁴ To calculate the lattice thermal conductivity (κ_l), the electronic thermal conductivity (κ_ele = LσT) was subtracted from κ. Due to the extremely small variance in L across the alloy series, the L estimated for PbTe was used to calculate κ_ele for all sample stoichiometries. Maintaining a consistent L also allowed for reasonable κ_l calculations for the bipolar PbTe_0.90Se_0.10 sample. Information on electronic transport properties and crystal structure is available in the ESI.†

Results and discussion

In this section, the trained interatomic potential is first validated against DFT results and experimental measurements of the stoichiometric compounds (i.e., PbTe and PbSe). With this validation, the heat transport properties of PbTe_1−xSe_x alloys are then investigated at temperatures ranging from 300 to 800 K.

Interatomic potential validation

The development of the Allegro interatomic potential was an iterative process to ensure that the final potential was robust. Attention was paid to both to accuracy (i.e., forces and energies across a range of temperatures) and stability in time (i.e., >50 ps). The use of AIMD simulations from 250–3000 K as training and test data was vital for both broadly sampling the configuration space and learning the minimum accurately.

To assess the performance of the machine-learned interatomic potential, we first evaluated its accuracy on the training set. The predicted energies, forces, and stresses closely match the PBEsol reference values, indicating an excellent fit. Model generalization was tested on independent configurations drawn from ab initio molecular dynamics (AIMD) simulations at 400, 600, and 700 K, as described previously. On this test set, the mean absolute errors (MAEs) were 0.62 meV per atom for energies, 6.65 meV Å⁻¹ for forces, and 0.03 meV Å⁻³ for stresses; the corresponding root mean squared errors (RMSEs) were 0.62 meV per atom, 8.38 meV Å⁻¹, and 0.05 meV Å⁻³, respectively. Data leakage is expected to be minimal due to the diversity and size of the training and validation datasets, which include thousands of decorrelated frames.

The radial distribution functions obtained with molecular dynamics simulations using the trained interatomic potential were consistent with those from AIMD simulations at the same temperature. The radial distribution functions of PbTe and the PbTe_0.25Se_0.75 alloy at 300 K and the PbTe_0.5Se_0.5 alloy at 300 K, 500 K, and 800 K are shown in Fig. S2 (ESI†). The interatomic potential successfully captures the interatomic separation and peak shape across different compositions.

The phonon dispersions of PbTe and PbSe obtained with the trained interatomic potential are in excellent agreement with PBEsol as shown in Fig. 2. At the Γ-point, PbTe has frequencies of 0 and 1.189 THz for the trained interatomic potential and 0 and 1.160 THz for PBEsol. PbSe has frequencies of 0 and 1.608 THz for the trained interatomic potential and 0 and 1.604 THz for PBEsol. The optical modes in PbSe are higher in frequency than PbTe due to the lighter atomic mass of Se and can be seen in the atom-projected phonon density of states. Furthermore, features such as the intersection of the longitudinal acoustic and transverse optical phonon modes in the Γ–X segment are consistent. We note that deviations occurred for smaller cutoff radii underlining the importance of long range interactions (Fig. S3, ESI†). Interatomic potentials trained with fewer tensors or Bessel functions had less reliable prediction of heat transport properties.

Fig. 2 Phonon dispersion of (A) PbTe and (B) PbSe at 0 K calculated with PBEsol and Allegro. The phonon density of states projected on the Pb and Te or Se atoms calculated with Allegro is also shown. The unit cells have Fm3̄m space group symmetry.

Transport in PbTe and PbSe

The long timescale simulations enabled by the Allegro interatomic potential allowed us to compute thermal transport coefficients using the Green–Kubo formalism. Before examining the PbTe_1−xSe_x alloys, we first consider stoichiometric PbTe and PbSe. As shown in Fig. 3, the normalized heat flux autocorrelation function (HFACF) for PbTe at 300 K exhibits pronounced oscillations at short times that decay to zero at longer correlation times. Within approximately 30 ps, the HFACF falls below 1% of its initial value. Including the atomic stress term in the heat flux expression (eqn (4)) was essential for accurately capturing these correlations.

Fig. 3 (A) Normalized heat flux autocorrelation function at 300 K and (B) lattice thermal conductivity of PbTe and PbSe from 300–800 K calculated using the Green–Kubo method compared with experimental results. The error bars are the standard error.

To reduce statistical noise, the heat flux was sampled over multiple short windows of 50 ps. This window length was chosen to capture the short-time behavior of the HFACF while ensuring its decay to zero, enabling convergence of the Green–Kubo integral and physical validity of the resulting transport coefficients. To further improve statistical reliability, we averaged the HFACF over several independent simulations with different initial conditions, thereby approximating a true ensemble average. This combined approach mitigates long-term drift and suppresses non-equilibrium artifacts.

The lattice thermal conductivity of PbTe and PbSe was calculated from 300–800 K using the Green–Kubo method. This method inherently includes anharmonic effects and finite temperature effects. The thermal conductivity was verified to be isotropic by inspection in each direction (Fig. S4, ESI†). The lattice thermal conductivity broadly decayed due to increased phonon–phonon interactions with increasing temperature.

Excellent summaries of other κ_L calculations and measurements of PbTe are found in ref. 22 and 25. However, prior experimental measurements either involved rather outdated instrumentation from the 1960s or involved doped samples. As such, we instead synthesized undoped samples across the alloy space.

The resulting polycrystalline ingots had densities >97% of the theoretical density and did not show evidence of secondary phases in X-ray diffraction (Fig. S5, ESI†). Resistivity values were consistently >10 mΩ cm at 300 K (Fig. S6, ESI†); as such the electronic contribution to thermal conductivity is less than 0.01 W m⁻¹ K⁻¹ (Fig. S7, ESI†). At high temperature, the small band gaps of these chalcogenides led to significant carrier excitation, with the Hall carrier concentration rising from ∼5 × 10⁻¹⁷ cm⁻³ at 300 K to ∼1 × 10⁻¹⁹ cm⁻³ at 800 K (Fig. S8, ESI†).

Laser flash diffusivity measurements on bulk polycrystalline PbTe and PbSe pellets indicate that thermal conductivity is primarily due to lattice vibrations. These experimental results show reasonably good agreement with the lattice thermal conductivity obtained from Allegro-based Green–Kubo calculations (Fig. 3B). Due to the polycrystallinity and defects within the experimental sample, the calculated thermal conductivities are slightly higher than the experimental values at room temperature. The difference becomes smaller at high temperature as phonon–phonon scattering dominates.

Parent spectral properties

Beyond predicting the aggregate lattice thermal conductivity, the Green–Kubo method allows for direct calculation of the frequency-dependent contributions, i.e., the spectral thermal conductivity. From the heat flux autocorrelation function, the spectral thermal conductivity is obtained from a Fourier transform. In PbTe, it is dominated by a sharp peak at around 1.7 THz (Fig. 4A). This spectral window consists of mixed optical modes and longitudinal acoustic modes with high band velocity (cf. Fig. 2A). PbSe shows very similar behavior, with a spectral thermal conductivity centered at approximately 2 THz. This peak shift is consistent with the stiffening found in the dispersion upon replacing tellurium with selenium (Fig. 2).

Fig. 4 (A) Spectral thermal conductivity of PbTe and PbSe at 300 K (B) spectral heat capacity and (C) spectral diffusivity, D, of PbTe at various temperatures. The plots of D are offset by increments of 1 × 10⁻⁵ m² s⁻¹.

In a previous study using a classical interatomic potential for PbTe, the maximum of the spectral κ was reported to occur between 1 and 1.5 THz (at 300 K).¹⁸ There were also significant contributions to κ between 0 and 1 THz and minor contributions between 3 and 4 THz (optical modes). A more recent study using neuroevolution potentials trained with the DFT-PBE method reported a relatively broad peak in spectral κ at around 1 THz, with smaller contributions between 0 and 0.5 THz and 1.5 and 2.5 THz.²⁶ Compared to the DFT-PBEsol-based interatomic potential used here, the DFT-PBE-trained neuroevolution potential shows somewhat larger differences in the phonon dispersion relations predicted by DFT or the interatomic potential.

To resolve the impact of temperature on the spectral diffusivity, D, we begin by calculating the spectral heat capacity from the phonon density of states using eqn (9) (Fig. 4B). At high temperature, the thermal energy populates the higher frequency modes and the molar heat capacities approach the Dulong–Petit value. At 800 K the molar heat capacity of PbTe is 24.85 J mol atom⁻¹ K⁻¹.

Next, the spectral diffusivity of PbTe was calculated from the spectral heat capacity and thermal conductivity using eqn (13). The diffusivity has a sharp peak which lies on the high frequency shoulder of the heat capacity (Fig. 4C). The available thermal energy is not effectively transferred at other frequencies. To consider the implications of these results, we review the classical approximations for the frequency dependence of phonon scattering within a power law approximation (τ ∝ ωⁿ): both phonon-boundary (n: 0 to −1) and phonon–phonon (n: −2) exhibit comparatively weak scattering at high frequencies to phonon-point defect (n: −4).⁵⁵ As such, engineering point defects in PbTe and PbSe should be quite effective at lowering lattice thermal conductivity. Indeed, a prior computational study by Qin et al. considering the effect of point defects with deep neural network potentials predicted a decrease in κ_L up to 50%.²⁵

Fig. 4C shows how the diffusivity decays with increasing temperature in PbTe. The peak in the spectral diffusivity decays by approximately one order of magnitude when the temperature is increased from 300 to 800 K. Viewing these results on logarithmic scale (Fig. S9, ESI†), we find that this order of magnitude reduction is consistent across all frequencies. The monotonic shift in the spectral diffusivity towards higher frequencies with temperature is unexpected; we hypothesize that this arises from shifting phonon branches with increasing temperature. Such shifts are also visible in the evolution of the heat capacity (Fig. 4B).

The spectral diffusivity of PbSe has a narrow peak similar to PbTe (Fig. S10, ESI†). Its narrow spectral window occurs at the corresponding region of mixed optical modes and longitudinal acoustic with high band velocity at 2.0 THz (Fig. 2B). The phonon dispersions were also calculated at 300 K, and as expected the band velocity is high for both PbTe and PbSe at these frequencies (Fig. S11, ESI†). These band dispersions include finite temperature and anharmonicity.

Transport in alloys

Next, the heat transport properties of PbTe_1−xSe_x alloys were investigated using the trained interatomic potential. The compositional dependence of the lattice thermal conductivity, κ_L, of the alloys at 300 K is shown in Fig. 5. Alloying causes a sharp decline until x = 0.5. This qualitative decline is expected from classical transport models such as the Klemens model.⁵⁶ The trained interatomic potential obtained thermal conductivities for compositions outside the training set. For example, a 10% substitution decreases the lattice thermal conductivity by about a third. Considering the decay of the heat flux autocorrelation function, the alloys typically exhibited a decay within 2 ps, reflecting the increased scattering pathways of the alloys (Fig. S12, ESI†).

Fig. 5 Calculated and experimental lattice thermal conductivity of the PbTe_1−xSe_x alloy at 300 K. The results of this work are compared with experimental PbTe_1−xSe_x⁵⁹ and Pb_0.96Sb_0.02Te_1−xSe_x alloys.²⁷ The error bars are the standard error.

Our calculated results have remarkably good agreement with both our own and previously published experimental measurements. The experimental lattice thermal conductivity of the Te-rich alloys are slightly smaller than the calculated values. While the difference is slightly smaller for a larger supercell size (Fig. S13, ESI†), it is more likely that the polycrystallinity and defects lower the thermal conductivity in the experimental sample as stated above. These effects would be more pronounced in the pristine structures. It should also be noted that experiments have variability,^28,57,58 and we also show results from PbTe_1−xSe_x standard solid solutions⁵⁹ and more recent works containing known Na or Sb impurities.²⁷ The impact on alloying on κ_L was less pronounced in previous computational works,^29,30,60 which may be related to our interatomic potential being trained directly on alloys, our inherent inclusion of anharmonic effects, or the exchange–correlation functional.

The reduction in lattice thermal conductivity from alloying is, as shown above, quite strong. Broadly, the reduction arises from the breakdown in periodicity of both the nuclear masses and spring constants. To this end, the relative contributions were probed with alchemical calculations consisting of “virtual” structures with the spring constants of a x = 0.5 alloy but having average mass on the anionic sites. The thermal conductivity of the virtual alloy without mass disorder was nearly as low as the x = 0.5 alloy (0.61 W m⁻¹ K⁻¹ and 0.46 W m⁻¹ K⁻¹, respectively, and using a cell containing 216 atoms). This indicates that the spring constant disorder plays a larger role in reducing the conductivity of the alloys than mass disorder.

Alloy spectral properties

The alchemical approach above highlights the important role that disordered bonding plays in reducing lattice thermal conductivity. To further resolve the underlying dynamics, we consider the spectral thermal conductivity and constituent properties therein.

To begin, we consider the vibrational density of states across the alloy series (Fig. 6). PbTe has a broad peak at 2.6 THz and two peaks at lower frequency (1.1 and 1.6 THz). These latter two peaks are associated with the transverse acoustic modes near U–W–X and L (Fig. 2). In contrast, the energy difference between these peaks collapses for PbSe, as can be seen by the nearly flat L–W–X transverse acoustic branches. Across the series, the alloys' transverse acoustic DOS show a smooth transition from the two end members. In addition, the higher frequency peaks blueshift and more peaks are resolved for the Se-rich alloys. Such shifts are expected for alloys containing the lighter Se and the corresponding stiffening of the bonds.

Fig. 6 (A) Phonon density of states of PbTe_1−xSe_x alloys at 300 K. Atom-projected density of states of the alloys with (B) PbTe_0.75Se_0.25 and (C) PbTe_0.25Se_0.75 normalized by the number of atoms of each element.

The atom-projected vibrational density of states of the PbTe_0.75Se_0.25 and PbTe_0.25Se_0.75 alloys are shown in Fig. 6B and C. Note, each density of state is normalized by the number of atoms of each element. The acoustic, low frequency, modes are dominated by Pb contributions and optical modes at high frequency modes are primarily associated with the lighter Te and Se atoms.

Typically, we would view lattice thermal conductivity within the Boltzmann transport formalism in the relaxation-time approximation (eqn (11)). However, the phonon modes of the alloys are smeared and cannot be neatly broken down into their individual velocity components. For example, the phonon dispersion of the PbTe_0.75Se_0.25 alloy at 300 K shown in Fig. S14 (ESI†) has no distinct bands above 0.5 THz. Thus, the reduction of thermal transport is better viewed in terms of the spectral diffusivity (eqn (12)). We thus seek to decompose the spectral thermal conductivity of the alloys into the heat capacity and diffusivity components.

At temperatures well above the Debye temperature, the equipartition theorem implies that the spectral heat capacity is proportional to the phonon density of states. As such, the spectral heat capacity of these alloys (Fig. S15, ESI†) follows a trend similar to the DOS shown in Fig. 6.

Armed with the alloy spectral heat capacities, the spectral thermal diffusivity of alloys at finite temperatures can be examined in detail (Fig. 7). In alloys, the most efficient frequency of heat transport shifts from about 1.7 to 2.0 THz. Despite the massive reduction in the lattice thermal conductivity in alloys seen in Fig. 5, the lattice thermal conductivity still fundamentally arises from the same frequency window at top of the longitudinal branch region. However, there is an order of magnitude loss of diffusivity (Fig. 7B).

Fig. 7 (A) Spectral diffusivity and (B) its maximum for PbTe_1−xSe_x alloys. The plots are truncated when phonon density of states falls below 0.01 states per ų. The plots of D are offset by increments of 2 × 10⁻⁵ m² s⁻¹.

Next, we can examine the thermal transport behavior of the alloys at high temperature. As shown in Fig. 8A, the thermal conductivity of the x = 0.5 alloy only decays at high temperatures, highlighting the glass-like transport in these alloys. This is consistent with the values measured experimentally. Resolving the spectral properties for the PbTe_0.5Se_0.5 alloy shows an increase in the heat capacity from 2.2 to 3.2 THz and a small decrease in the spectral thermal diffusivity at elevated temperature (Fig. 8B and C). The thermal diffusivity of the alloy is already an order of magnitude lower than pristine PbTe at 300 K and the increased thermal occupation is able to limit the deterioration of the thermal conductivity. The maximum diffusivity and the heat capacity at that frequency of the alloy and PbTe are shown in Fig. S16 (ESI†).

Fig. 8 (A) Calculated and experimentally measured lattice thermal conductivity of PbTe_0.5Se_0.5 alloys from 300 to 800 K. (B) Calculated spectral heat capacity and (C) spectral diffusivity of the alloy. The diffusivity plots are truncated when phonon density of states falls below 0.01 states per ų. The plots of D are offset by increments of 1 × 10⁻⁶ m² s⁻¹.

Conclusions

The phonon properties and heat transport of PbTe_1−xSe_x alloys were investigated using an equivariant graph neural network interatomic potential. By training directly on ab initio calculations of the alloys, the interatomic potential achieves both compositional and configurational generality. It shows strong agreement with DFT calculations and experimental measurements, confirming its accuracy and reliability. The long timescale simulations enabled by the machine-learned interatomic potential allowed for the calculation of thermal transport coefficients using the Green–Kubo formalism. The Green–Kubo formulism intrinsically incorporates phonon–phonon interactions, anharmonicity, and finite-temperature behavior without relying on quasiharmonic approximations or empirical corrections. The predicted reduction in lattice thermal conductivity for PbTe_1−xSe_x alloys is consistent with experimental measurements and previous computational results.^29,30,60,61 Additionally, our approach using machine-learned interatomic potential facilitates spectral analysis of thermal conductivity and heat capacity across the 300–800 K temperature range. We identify a narrow spectral window near 1.7–2.0 THz, dominated by mixed longitudinal acoustic and optical modes, as the primary contributor to heat transport across all compositions. Alloying strongly suppresses the diffusivity both in this spectral window and throughout the spectrum resulting in the deterioration of thermal conductivity. At elevated temperatures, both the spectral diffusivity and overall thermal conductivity decrease across the parent compounds. However, in mid-range compositions (x = 0.5), the reduction in spectral diffusivity occurs alongside an increase in heat capacity driven by greater thermal occupation of vibrational modes. This compensating effect contributes to a comparatively weaker suppression of thermal conductivity, indicating a reduced sensitivity to thermal degradation. Using alchemical simulations, we further demonstrate that force–constant disorder, rather than mass contrast, is the dominant mechanism behind this suppression. These findings elucidate the microscopic origins of the thermal transport degradation in PbTe_1−xSe_x alloys. The findings demonstrate how to utilize machine-learned interatomic potentials for predictive, physics-rich thermal transport modeling in both medium and high-temperature thermoelectrics as well as other complex, disordered solids.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data supporting this article have been included as part of the ESI.†

Acknowledgements

We acknowledge CSC – IT Center for Science (Finland) for funding (Finland-Colorado joint collaboration project CNTR0018608) and access to the LUMI supercomputer, owned by the EuroHPC Joint Undertaking, hosted by CSC and the LUMI consortium. E. S. T., T. B., C. G., and A. N. acknowledge the support of the National Science Foundation under award OAC 2118201.

References

1. X. Qian, J. Zhou and G. Chen, Phononengineered extreme thermal conductivity materials, Nat. Mater., 2021, 20, 1188–1202.
2. S. Arabha, Z. S. Aghbolagh, K. Ghorbani, S. M. Hatam-Lee and A. Rajabpour, Recent advances in lattice thermal conductivity calculation using machine-learning interatomic potentials, J. Appl. Phys., 2021, 130, 210903.
3. M. Zhang, J. Cai, F. Gao, Z. Zhang, M. Li, Z. Chen, Y. Wang, D. Hu, X. Tan and G. Liu, et al., Improved Thermo electric Performance of p-Type PbTe by Entropy Engineering and Temperature-Dependent Precipitates, ACS Appl. Mater. Interfaces, 2023, 16, 907–914.
4. P. Parajuli, S. Bhattacharya, R. Rao and A. Rao, Phonon anharmonicity in binary chalcogenides for efficient energy harvesting, Mater. Horiz., 2022, 9, 1602–1622.
5. F. Knoop, T. A. Purcell, M. Scheffler and C. Carbogno, Anharmonicity measure for materials, Phys. Rev. Mater., 2020, 4, 083809.
6. S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat. Commun., 2022, 13, 2453.
7. A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth and B. Kozinsky, Learning local equivariant representations for large-scale atomistic dynamics, Nat. Commun., 2023, 14, 579.
8. A. K. Gupta, M. M. Stulajter, Y. Shaidu, J. B. Neaton and W. A. de Jong, Equivariant neural networks utilizing molecular clusters for accurate molecular crystal lattice energy predictions, ACS Omega, 2024, 9, 40269–40282.
9. J. S. Jestila, N. Wu, F. Priante and A. S. Foster, Accelerated lignocellulosic molecule adsorption structure determination, J. Chem. Theory Comput., 2024, 20, 2297–2312.
10. S.-H. Lee, J. Li, V. Olevano and B. Sklenard, Equivariant graph neural network interatomic potential for Green-Kubo thermal conductivity in phase change materials, Phys. Rev. Mater., 2024, 8, 033802.
11. A. D. LaLonde, Y. Pei, H. Wang and G. Jeffrey Snyder, Lead telluride alloy thermoelectrics, Mater. Today, 2011, 14, 526–532.
12. Q. Yan and M. G. Kanatzidis, High performance thermoelectrics and challenges for practical devices, Nat. Mater., 2022, 21, 503–513.
13. A. El-Sharkawy, A. Abou El-Azm, M. Kenawy, A. Hillal and H. Abu-Basha, Thermophysical properties of polycrystalline PbS, PbSe, and PbTe in the temperature range 300–700 K, Int. J. Thermophys., 1983, 4, 261–269.
14. O. Delaire, J. Ma, K. Marty, A. F. May, M. A. McGuire, M.-H. Du, D. J. Singh, A. Podlesnyak, G. Ehlers and M. Lumsden, et al., Giant anharmonic phonon scattering in PbTe, Nat. Mater., 2011, 10, 614–619.
15. H. Kim and M. Kaviany, Effect of thermal disorder on high figure of merit in PbTe, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 045213.
16. S. Lee, K. Esfarjani, T. Luo, J. Zhou, Z. Tian and G. Chen, Resonant bonding leads to low lattice thermal conductivity, Nat. Commun., 2014, 5, 3525.
17. J. M. Skelton, S. C. Parker, A. Togo, I. Tanaka and A. Walsh, Thermal physics of the lead chalcogenides PbS, PbSe, and PbTe from first principles, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 205203.
18. B. Qiu, H. Bao, G. Zhang, Y. Wu and X. Ruan, Molecular dynamics simulations of lattice thermal conductivity and spectral phonon mean free path of PbTe: bulk and nanostructures, Comput. Mater. Sci., 2012, 53, 278–285.
19. Y. Chen, X. Ai and C. A. Marianetti, First-Principles Approach to Nonlinear Lattice Dynamics: Anomalous Spectra in PbTe, Phys. Rev. Lett., 2014, 113, 105501.
20. T. Shiga, J. Shiomi, J. Ma, O. Delaire, T. Radzynski, A. Lusakowski, K. Esfarjani and G. Chen, Microscopic mechanism of low thermal conductivity in lead telluride, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 155203.
21. A. H. Romero, E. K. U. Gross, M. J. Verstraete and O. Hellman, Thermal conductivity in PbTe from first principles, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 214310.
22. Y. Lu, T. Sun and D.-B. Zhang, Lattice anharmonicity, phonon dispersion, and thermal conductivity of PbTe studied by the phonon quasiparticle approach, Phys. Rev. B, 2018, 97, 174304.
23. Y. Xia, Revisiting lattice thermal transport in PbTe: the crucial role of quartic anharmonicity, Appl. Phys. Lett., 2018, 113, 073901.
24. Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen and T. Ala-Nissila, Neuroevolution machine learning potentials: combining high accuracy and low cost in atomistic simulations and application to heat transport, Phys. Rev. B, 2021, 104, 104309.
25. M. Qin, X. Zhang, J. Zhu, Y. Yang, Z. Ti, Y. Shen, X. Wang, X. Liu and Y. Zhang, A machine learning methodology to investigate the lattice thermal conductivity of defected PbTe, J. Mater. Chem. A, 2023, 11, 10612–10627.
26. R. Cheng, X. Shen, S. Klotz, Z. Zeng, Z. Li, A. Ivanov, Y. Xiao, L.-D. Zhao, F. Weber and Y. Chen, Lattice dynamics and thermal transport of PbTe under high pressure, Phys. Rev. B, 2023, 108, 104306.
27. P. F. Poudeu, J. D’Angelo, H. Kong, A. Downey, J. L. Short, R. Pcionek, T. P. Hogan, C. Uher and M. G. Kanatzidis, Nanostructures versus Solid Solutions: Low Lattice Thermal Conductivity and Enhanced Thermoelectric Figure of Merit in Pb_9.6Sb_0.2Te_10−xSe_x Bulk Materials, J. Am. Chem. Soc., 2006, 128, 14347–14355.
28. Y. Pei, X. Shi, A. LaLonde, H. Wang, L. Chen and G. J. Snyder, Convergence of electronic bands for high performance bulk thermoelectrics, Nature, 2011, 473, 66–69.
29. Z. Tian, J. Garg, K. Esfarjani, T. Shiga, J. Shiomi and G. Chen, Phonon conduction in PbSe, PbTe, and PbTe_1−xSe_x from first-principles calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 184303.
30. T. Murakami, T. Shiga, T. Hori, K. Esfarjani and J. Shiomi, Importance of local force fields on lattice thermal conductivity reduction in PbTe_1−xSe_x alloys, Europhys. Lett., 2013, 102, 46002.
31. Z. Tian, M. Li, Z. Ren, H. Ma, A. Alatas, S. D. Wilson and J. Li, Inelastic X-ray scattering measurements of phonon dispersion and lifetimes in PbTe_1−xSe_x alloys, J. Phys.: Condens. Matter, 2015, 27, 375403.
32. R. Kubo, The fluctuation-dissipation theorem, Rep. Prog. Phys., 1966, 29, 255.
33. Z. Zhang, Y. Guo, M. Bescond, J. Chen, M. Nomura and S. Volz, Heat Conduction Theory Including Phonon Coherence, Phys. Rev. Lett., 2022, 128, 015901.
34. T. D. Kuhne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schutt and F. Schiffmann, et al., CP2K: an electronic structure and molecular dynamics soft ware package-Quickstep: efficient and accurate electronic structure calculations, J. Chem. Phys., 2020, 152, 194103.
35. K. Okhotnikov, T. Charpentier and S. Cadars, Supercell program: a combinatorial structure-generation approach for the local-level modeling of atomic substitutions and partial occupancies in crystals, J. Cheminf., 2016, 8, 1–15.
36. C. Chen and S. P. Ong, A universal graph deep learning interatomic potential for the periodic table, Nat. Comput. Sci., 2022, 2, 718–728.
37. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Phys. Rev. Lett., 2008, 100, 136406.
38. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
39. J. M. Skelton, D. Tiana, S. C. Parker, A. Togo, I. Tanaka and A. Walsh, Influence of the exchange-correlation functional on the quasi-harmonic lattice dynamics of II VI semiconductors, J. Chem. Phys., 2015, 143, 064710.
40. C. Hartwigsen, S. Goedecker and J. Hutter, Relativistic separable dual-space Gaussian pseudopotentials from H to Rn, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641.
41. J. VandeVondele and J. Hutter, Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases, J. Chem. Phys., 2007, 127, 114105.
42. D. Hendrycks and K. Gimpel, Gaussian Error Linear Units (GELUs), arXiv, 2016, e-prints, arXiv:abs/1606.08415 DOI:10.48550/arXiv.1606.08415.
43. A. P. Thompson, S. J. Plimpton and W. Mattson, General formulation of pressure and stress tensor for arbitrary many body interaction potentials under periodic boundary conditions, J. Chem. Phys., 2009, 131, 154107.
44. Z. Fan, L. F. C. Pereira, H.-Q. Wang, J.-C. Zheng, D. Donadio and A. Harju, Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 094301.
45. D. Surblys, H. Matsubara, G. Kikugawa and T. Ohara, Application of Atomic Stress to Compute Heat Flux via Molecular Dynamics for Systems With Many-Body Interactions, Phys. Rev. E, 2019, 99, 051301.
46. H. Yu, Y. Zhong, L. Hong, C. Xu, W. Ren, X. Gong and H. Xiang, Spin dependent graph neural network potential for magnetic materials, Phys. Rev. B, 2024, 109, 144426.
47. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al., LAMMPS – a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
48. A. Togo, L. Chaput, T. Tadano and I. Tanaka, Implementation strategies in phonopy and phono3py, J. Phys.: Condens. Matter, 2023, 35, 353001.
49. A. Togo, First-principles Phonon Calculations with Phonopy and Phono3py, J. Phys. Soc. Jpn., 2023, 92, 012001.
50. L. T. Kong, Phonon dispersion measured directly from molecular dynamics simulations, Comput. Phys. Commun., 2011, 182, 2201–2207.
51. C. Campana and M. H. Muser, Practical Green's function approach to the simulation of elastic semi-infinite solids, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 075420.
52. L. T. Kong, G. Bartels, C. Campana, C. Denniston and M. H. Muser, Implementation of Green's function molecular dynamics: an extension to LAMMPS, Comput. Phys. Commun., 2009, 180, 1004–1010.
53. L. T. Kong, C. Denniston and M. H. Muser, An improved version of the Green's function molecular dynamics method, Comput. Phys. Commun., 2011, 182, 540–541.
54. H.-S. Kim, Z. M. Gibbs, Y. Tang, H. Wang and G. J. Snyder, Characterization of Lorenz number with Seebeck coefficient measurement, APL Mater., 2015, 3, 041506.
55. E. S. Toberer, A. Zevalkink and G. J. Snyder, Phonon engineering through crystal chemistry, J. Mater. Chem., 2011, 21, 15843–15852.
56. P. Klemens, Thermal resistance due to point defects at high temperatures, Phys. Rev., 1960, 119, 507.
57. A. Ioffe and A. Ioffe, Thermal conductivity of semiconductor solid solutions, Sov. Phys.-Solid State, 1960, 2, 719–728.
58. G. Alekseeva, B. Efimova, L. Ostrovskaya, O. Serebryannikova and M. Tsypin, Thermal conductivity of solid solutions based on lead telluride, Sov. Phys. Semicond., 1971, 4, 1122–1125.
59. V. V. Devyatkova and E. D. Tikhonov, Sov. Phys. Solid State, 1965, 7, 1427–1431.
60. H. Wang, A. D. LaLonde, Y. Pei and G. J. Snyder, The criteria for beneficial disorder in thermoelectric solid solutions, Adv. Funct. Mater., 2013, 23, 1586–1596.
61. J. Wang, H. Wang, G. Snyder, X. Zhang, Z. Ni and Y. Chen, Characteristics of lattice thermal conductivity and carrier mobility of undoped PbSe-PbS solid solutions, J. Phys. D: Appl. Phys., 2013, 46, 405301.

---

Footnotes

† Electronic supplementary information (ESI) available: It contains the training sets, example scripts for DFT, relaxed lattice constants, comparison of the radial distribution functions, thermal conductivity along different directions and for different cell sizes, additional experimental characterization, and supplemental spectral heat transport properties. See DOI: https://doi.org/10.1039/d5mh00934k

‡ Current address: VTT Technical Research Centre of Finland Ltd., VTT, FI-02044 Finland.

---