Origin of anomalous anharmonic lattice dynamics of lead telluride

The origin of the anomalous anharmonic lattice dynamics of lead telluride is investigated using molecular dynamics simulations with interatomic force constants (IFCs) up to quartic terms obtained from first principles. The calculations reproduce the peak asymmetry of the radial distribution functions and the double peaks of transverse optical phonon previously observed with neutron diffraction and scattering experiments. They are identified to be due to the extremely large nearest-neighbor cubic IFCs in the [100] direction. The outstanding strength of the nearest-neighbor cubic IFCs relative to the longer-range ones explains the reason why the distortion in the radial distribution function is local.

💡 Research Summary

This paper investigates the origin of the anomalous anharmonic lattice dynamics observed in lead telluride (PbTe), a material renowned for its intrinsically low lattice thermal conductivity and high thermoelectric performance. Experimental studies using neutron diffraction, inelastic neutron scattering (INS), and extended X‑ray absorption fine structure (EXAFS) have reported two striking phenomena: (i) a pronounced asymmetry in the nearest‑neighbor Pb–Te pair distribution function (PDF) that grows with temperature, and (ii) a double‑peak structure in the transverse‑optical (TO) phonon mode near the Brillouin‑zone centre, accompanied by an avoided crossing with the longitudinal acoustic (LA) branch. However, the microscopic cause of these features has remained debated.

To address this, the authors first compute interatomic force constants (IFCs) up to fourth order from first‑principles density‑functional theory (DFT) using the real‑space displacement method. Harmonic (second‑order) IFCs are retained out to the sixth neighbour, while cubic (third‑order) and quartic (fourth‑order) terms are limited to the first neighbour. These IFCs provide a non‑empirical, highly accurate force field that can be directly inserted into classical molecular‑dynamics (MD) simulations.

MD simulations are performed on two supercell sizes: a 4 × 4 × 4 cell (512 atoms) for the radial distribution function and a larger 10 × 10 × 10 cell (8 000 atoms) for the dynamical structure factor S(Q, ω). For each temperature (50, 100, 150, 200, 250 K) ten independent trajectories are generated, equilibrated, and then run in the microcanonical ensemble for 600 ps (small cell) or 2 ns (large cell) with a 2 fs timestep. Ensemble averaging reduces statistical noise.

The computed PDFs reveal that the first‑neighbor Pb–Te peak broadens and becomes markedly non‑Gaussian as temperature rises. Quantitative analysis using skewness and kurtosis shows a monotonic increase in both metrics for the Pb–Te pair, while Pb–Pb and Te–Te pairs remain essentially Gaussian. This indicates that the anharmonic distortion is highly localized to the nearest‑neighbor bond.

The dynamical structure factor S(Q, ω) is evaluated for wavevectors along the