Theoretical modeling of phonon transport process in strongly anharmonic materials at a finite temperature needs to accurately capture the effects of lattice anharmonicity. The anharmonicity of potential energy surface would result in not only strong phonon scatterings but also shifts of phonon frequencies and eigenvectors. In this work, we evaluated the roles of anharmonicity-renormalized phonon eigenvectors in predicting phonon transport properties of anharmonic crystals at high temperatures using molecular dynamics-based normal mode analysis (NMA) methods in both time domain and frequency domain. Using PbTe as a model of strongly anharmonic crystal, we analyzed the numerical challenges to extract phonon lifetimes using NMA methods when phonon eigenvectors deviate from their harmonic values at high temperatures. To solve these issues, we proposed and verified a better fitting strategy, Sum-up Spectrum Fitting Method (SSFM) than the original frequency-domain NMA method. SSFM is to project the total spectrum energy density data of all phonon modes onto an inaccurate (harmonic or quasi-harmonic) eigenvector base and then manually sum up the peaks that belong to the same phonon mode (at the same frequency). The SSFM relaxes the requirement for accurate temperature-dependent eigenvectors, making it robust for analyzing strongly anharmonic crystals at high temperatures.

## I. INTRODUCTION

An improved understanding of the thermal transport process at high temperatures will accelerate materials innovation for a wide range of applications, such as thermal barrier coating, nuclear cladding, and thermoelectrics.^{1–5} Theoretical modeling of the phonon transport process at a finite temperature needs to accurately capture the effects of lattice anharmonicity, especially for strong anharmonic materials.^{6–9} As the temperature increases, the vibrational phase space becomes significantly different from that at zero temperature, invalidating the local harmonic approximation (HA). At finite and high temperatures, the anharmonicity of potential energy surface would result in not only phonon scatterings but also frequency shift of the phonons.^{10} While the leading-order effect of frequency shift can be described by the quasi-harmonic approximation (QHA), the phonon–phonon interaction, namely the self-energy of phonons, could also induce an additional frequency shift, which is a result of higher-order force constants. Such a frequency shift is known as phonon renormalization, particularly important for strongly anharmonic materials such as Zr,^{11} SnSe,^{12,13} SnS,^{12} and PbTe.^{14,15}

To capture these effects and accurately predict thermal properties in strongly anharmonic materials at high temperatures, there are recent developments on first-principles methods to predict phonon properties.^{16–18} With accurate interatomic force constants (IFCs) from density function theory (DFT), thermal properties predicted by the anharmonic lattice dynamics-Boltzmann transport equation (ALD-BTE) approach achieves an excellent agreement with the experimental data for both phonon frequency and lifetime at low and medium temperatures or in weakly anharmonic materials.^{19–24} However, at high temperatures or in strongly anharmonic materials, higher-order phonon scattering processes become important.^{16,17,25} The temperature-dependent effective potential (TDEP) was recently developed to iteratively refine the temperature-dependent IFCs based on atomic configurations and the corresponding potential energy surface at finite temperatures.^{26–28} Another method is the self-consistent phonon approach where a direct computation of Green's function and free energy of phonons at finite temperatures is performed to obtain the renormalized phonon frequencies.^{29,30}

Different from first-principles methods, classical molecular dynamics (MD) can intrinsically consider all orders of anharmonicity and model the thermal transport processes in complicated systems even beyond the perturbative assumption when the temperature is well above the respective Debye temperature.^{31–33} Because MD is in the real space and time domains, normal mode analysis (NMA) methods are typically used to analyze the phonon frequency and lifetime.^{34,35} In the MD-based NMA methods, the time-dependent atomic displacements extracted from the MD simulation are projected to the corresponding normal mode coordinates using the phonon frequencies and eigenvectors computed from lattice dynamics. The NMA can be done in the time-domain (i.e., TD-NMA)^{35,36} by computing the mode-wise correlation functions or in the frequency domain by analyzing the spectral energy density (SED).^{37–39} Previous studies using NMA methods have shown that both TD-NMA and SED methods can extract phonon frequency and lifetime, such as in Argon,^{34,35} Si,^{34,36} CNT,^{34} polymer crystals,^{40} and two-dimensional materials.^{41,42} However, few previous works have been done to quantitatively understand the effect of strong anharmonicity at high temperatures using MD-based methods. For example, McGaughey *et al.* compared the phonon properties of Si and Argon at low temperatures calculated from TD-NMA and SED, where the anharmonicity is not strong in the system.^{34,35} Feng *et al.*^{39} showed from analytical derivations that the eigenvectors are unnecessary in the SED method, without computing and comparing the phonon properties explicitly. Puligheddu *et al.*^{31} compared the phonon property predictions from TD-NMA and ALD-BTE approaches but did not compare the TD-NMA and SED methods. Their comparison of methods in PbTe is only up to 300 K where the higher-order phonon interactions begin to play a role.

Even though MD intrinsically considers all orders of anharmonicity, it remains unclear how we should implement the NMA methods considering the higher-order effects of anharmonicity. Phonon eigenvectors are important inputs for the NMA calculations but the renormalization of eigenvectors or the temperature dependence of eigenvectors is not usually taken into consideration in the previous literature. In the TD-NMA method, correct phonon frequencies and eigenvectors are required as inputs. In the SED method, applying the eigenvectors can separate two modes with similar frequencies that are not able to be distinguished. In this aspect, whether harmonic eigenvectors or anharmonic eigenvectors that consider both QHA and phonon renormalization need to be used as inputs for phonon lifetime calculations remains unexplored.

In this work, we evaluated the roles of anharmonicity-renormalized phonon eigenvectors in predicting phonon lifetimes at high temperatures in MD-based TD-NMA and SED methods. We selected a representative strongly anharmonic materials system—lead telluride (PbTe) crystal—and calculated its thermal transport properties well above its Debye temperature ($\theta D=177K$) using the NMA methods.^{6} We tested and verified a better fitting strategy than the original SED method and provided a few recommendations for practically computing phonon transport properties of strongly anharmonic materials using MD-based methods, especially regarding the roles of phonon eigenvectors.

## II. MODELS AND METHODS

### A. Equilibrium molecular dynamics (EMD) simulation

We used the same interatomic potential for all the MD and lattice dynamics simulations so that all the comparisons are fair. For PbTe, we used the potential developed by Qiu *et al.*, which has been demonstrated to reasonably reproduce the vibrational properties and lattice thermal conductivity.^{38,43}

All the MD simulations were performed using the open-source software LAMMPS.^{44} The domain sizes for all the bulk systems are 8 × 8 × 8 unit cells (typically 4096 atoms). Periodic boundary conditions were applied to the simulation domain. The time step was 0.5 fs. Twenty independent cases were simulated with different initial random velocities to reduce the statistical uncertainties.^{45,46} The simulated material was first relaxed in an NPT (constant number of atoms, pressure, and temperature) ensemble with temperatures of 300, 500, and 800 K, and a pressure of 0 MPa for 1 ns, to expand the lattice. Then, both NVT (constant number of atoms, volume, and temperature) and NVE (constant number of atoms, volume, and energy) ensembles were applied for equilibrating the material systems for 1 ns each.^{47} After obtaining relaxed structures, the systems were kept in an NVE ensemble for position and velocity sampling. For each simulation of PbTe (out of 20), the sampling time was 2.5, 2, and 4 ns for Green–Kubo EMD thermal conductivity, TD-NMA lifetime, and SED lifetime calculations, respectively. The sampling interval was 5 fs, which is small enough to accurately obtain the statistical information of atomic vibrations in PbTe. Using the Green–Kubo formula for thermal conductivity, we can compute the thermal conductivity of PbTe at each temperature. The details of computing the Green–Kubo EMD thermal conductivity can be found in our previous publications.^{40,41,48}

### B. Lattice dynamics calculation

For lattice dynamics calculations, the HA uses the lattice constants at 0 K and employs the interatomic potential to compute the interatomic force constants and then write the dynamical matrix.^{15} Phonon frequencies and eigenvectors were solved from the dynamical equation. The HA is inadequate to account for thermal expansion, as the equilibrium distance between atoms in such a model is independent of temperature. The QHA takes into account the temperature-dependent volume and solves the dynamical equation.^{49} The phonon dispersions were calculated under these approximations. Under this framework, we calculated HA and QHA dispersion curves by changing the lattice constants of the unit cell. In HA, the lattice constant of the unit cell at 0 K was obtained by linearly extrapolating the lattice constants at 100, 300, 500, 800, and 1000 K from MD.

The dynamical matrix can also be formed using renormalized force constants calculated by TDEP or SED. Then, the renormalized phonon frequencies and eigenvectors can be computed from the dynamical equation.

### C. NMA methods

To compute phonon lifetime, we used the NMA methods, including the TD-NMA and FD-NMA [also called Spectral Energy Density (SED)], with the general details found in our previous publications.^{40} Here, we briefly introduce these two methods and provide the technical information in this work. In the TD-NMA method, each time-dependent phonon mode is

where

Here, ** k** is the wavevector, $u\alpha l,b$ is the $\alpha th$ component of the displacement of the $bth$ basis atom in the $lth$ unit cell, $e\u2217$ is the complex conjugate of eigenvector component, $r0l$ is the equilibrium position of the $lth$ unit cell, $\nu $ denotes phonon branches, and

*n*and $NC$ are the numbers of total basis atoms and total unit cells. Different eigenvectors (e.g., eigenvectors from HA and QHA) were used to compare their effects on calculated phonon lifetimes. The spectral phonon relaxation time is obtained by fitting the autocorrelation function of the total energy of each mode.

^{35}

The SED of the phonon modes, i.e., the Fourier transform of the autocorrelation function of each mode has a well-defined spectral shape,

where $\Phi \nu (k,\omega )$ is the Fourier transform of the time derivative of $qk,\nu (t)$. Sun and Allen showed that the approximate expression of this power spectrum has a Lorentzian function shape given by^{34}

where $\Phi k,\nu (\omega )$ is the distribution of SED in the frequency domain at a specific wave vector, *A* is the area of the single peak, *w* is the full width at half maximum, and $\omega c$ is the central frequency (phonon frequency) of the mode. The phonon lifetime can be calculated as

The peak positions map out the phonon dispersion. In the case that two modes have similar frequencies so that their SEDs cannot be distinguished, applying eigenvectors can help separate them into two individual peaks.

The number of wavevectors supported in a SED simulation is limited by the MD domain size so that the dispersion curve is not continuous. For example, if the periodic MD domain contains 8 × 8 × 8 unit cells, we can only resolve 8** k** points (exclude the gamma point) in the [1 0 0] direction with reduced wave vectors of

*** = (**

*k**j*/8,0,0), where

*j*is an integer from 1 to 8. To obtain continuous dispersion curves instead of discrete points, we used the software DynaPhoPy.

^{50}DynaPhoPy is a recently developed program by Carreras

*et al.*to extract the phonon dispersion from the samplings of MD simulations. In DynaPhoPy, phonon frequencies calculated at commensurate

**q**-points by SED corresponding to the chosen supercell size and phonon eigenvectors from the QHA are used to interpolate phonon frequencies at incommensurate

**q**-points.

Specifically for PbTe, the non-analytical corrections are applied to calculate the phonon band structure.^{20,27} This nonanalytic contribution accounts for the splitting between longitudinal optical branch (LO) and transverse optical branch (TO), i.e., the removal of degeneracy between LO and TO phonons at the Brillouin zone center.^{51} The Born effective charges and dielectric tensor were needed for the calculations.^{20} We have verified that the phonon frequencies from DynaPhoPy match well with the SED calculations.

### D. Calculation of BTE thermal conductivity

We used the phonon group velocity from the SED-derived dispersion and the phonon lifetime from NMA methods to predict the thermal conductivity, under the BTE framework with the relaxation time approximation. This thermal conductivity value will be compared with that predicted by the Green–Kubo EMD method (with 14% uncertainty), without the concept of phonons. The uncertainty of the EMD method is calculated when the total simulation time in EMD is 40 ns and the correlation time is 200 ps, according to Wang *et al.*'s work.^{45}

For the ALD-BTE framework, the phonon group velocity and lifetime are calculated from interatomic force constants (IFCs). The system Hamiltonian is cut to the third-order interaction for the calculations without phonon renormalization,

Here, $mi$, $pi$, and $ui$ are the mass, momentum, and displacement from equilibrium of atom *i*, respectively, $\alpha \beta \gamma $ designate the Cartesian components, and $\Phi $ and $\Psi $ are the second- and third-order effective IFCs, respectively. Because the analytical expression of the interatomic potential is known in our work, the IFCs can be directly calculated using the finite difference method. In the IFCs calculations, we used 8 × 8 × 8 supercells to extract second- and third-order force constants using the Phonopy software.^{52} Once the IFCs are known, the phonon frequencies and eigenvectors can be solved from the dynamical equation using the second-order force constants. The phonon lifetime using only the third-order force constants can also be calculated using Fermi's golden rule.^{22} We did not calculate the fourth-order phonon scattering processes. We have used the ShengBTE software package for these calculations with IFCs as inputs, which were solved using a **q**-point mesh of 10 × 10 × 10.^{53}

### E. TDEP method

The TDEP method developed by Hellman *et al.* is implemented for phonon renormalization, which is to obtain the best possible second-order and third-order force constants as a fit to the MD potential energy surface at finite temperatures. The fitting is done by comparing the forces of the model and the MD system at each time step and minimizing the difference. The explicit procedure for extracting these IFCs is described in the work by Hellman *et al.*^{26,27} With the TDEP method, we can renormalize the second- and third-order force constants considering the higher-order anharmonicity induced frequency-shifts of the system. The phonon dispersion, group velocity, eigenvector, and lifetime were calculated with the renormalized IFCs. Similar to the calculations of phonon dispersion curves in DynaPhoPy, the non-analytical corrections were applied in TDEP as well.

## III. RESULTS AND DISCUSSION

### A. Phonon renormalization—Frequency and eigenvector

First, we calculate the shift of phonon frequencies and eigenvectors from their harmonic values due to the lattice anharmonicity at high temperatures, i.e., phonon renormalization. Figure 1(a) shows the phonon dispersion curves of PbTe at 300 and 800 K by QHA, SED, and TDEP methods. We also calculated the phonon dispersion curve of PbTe at 0 K by HA as a comparison, which has the highest phonon frequencies without considering lattice thermal expansion. At 300 K, the dispersion curves calculated from QHA, SED, and TDEP methods agree well with each other. This matching indicates that the dominant temperature effect at 300 K (∼$1.7\theta D$) is the volumetric lattice expansion, where the QHA is sufficient to capture the major effect. Figure 1(b) shows the same comparison at 800 K (∼$4.5\theta D$), which clearly demonstrates that the QHA is not able to capture all the temperature effects on anharmonic phonon frequencies. The higher-order anharmonic phonon interactions must be considered to understand the anharmonic phonon frequency shift, which is more evident in the optical branches. The phonon frequencies calculated by the SED and TDEP methods agree with each other at 300 K but there are some numerical discrepancies at 800 K. Because the continuous SED dispersion curves match well with the SED peaks calculated from MD simulations (black circles), we chose to calculate group velocities from the continuous SED dispersion curves at different temperatures for the BTE thermal conductivity calculations later.

We then evaluate how the strong anharmonicity at high temperatures affects phonon eigenvectors. In contrast with anharmonic phonon frequencies, anharmonic eigenvectors are seldom taken into consideration during the calculations. For PbTe with two basis atoms in the primitive cell, the amplitude of the two components $|eTe|$ and $|ePb|$ can be viewed as the weights of the phonon normal mode projections on the two basis atoms. We then used the eigenvector amplitude ratio $|eTe|/|ePb|$ to evaluate the temperature effect. Figure 2 shows the amplitude ratios of phonon eigenvectors of PbTe with increasing temperatures from 0 to 1000 K. The amplitude ratio does not change with increasing temperatures at the Г point, whereas the ratio does shift at many other wavevectors (away from the Brillouin zone center) with a significant amount. This result is similar to the observation in the work by Feng *et al.*^{39} where they found a consistent difference in harmonic and anharmonic eigenvectors in the longitudinal acoustic mode in PbTe at 800 K.

At the same temperature, the TDEP-resolved anharmonicity-renormalized eigenvectors can capture a larger change of amplitude ratio than those solved from QHA. This difference is negligible at 300 K but becomes noticeable at 800 K (see Fig. S1 in the supplementary material). Intuitively, the differences between the eigenvectors at 0 K using HA, using the QHA under the “frozen phonon” assumption, or renormalized at high temperatures could result in different lifetimes. This paper will focus on resolving the effect of temperature-dependent eigenvectors in extracting phonon lifetimes. We will first show the three numerical challenges related to using correct phonon eigenvectors for normal mode analysis, and then test a new strategy to extract the lifetime with much higher accuracy.

### B. Numerical challenges of normal mode analysis

Next, we will examine how phonon anharmonicity affects the prediction of phonon lifetime using NMA methods. Here we show three aspects that need extra attention to obtain robust results of phonon lifetimes.

First, there could be extra uncertainty if one directly fits the SED without decomposing into individual modes, especially when the frequencies are close. In the frequency domain, each peak in the SED spectrum corresponds to a phonon mode and a Lorentzian function can then be used to fit the peak either manually or using an automatic fitting program recognizing the peaks. Phonon eigenvectors are not needed if the individual peaks can be perfectly separated. However, if two phonon modes have similar phonon frequencies, the peaks in the frequency domain will be very close to each other so that distinguishing one peak from the other is practically difficult. Figure 3 demonstrates the necessity of eigenvectors to separate overlap peaks in SED.^{39} Figure 3(a) shows the overlapping peaks of the transverse and longitudinal optical modes (around 3 THz) in the SED spectrum, which can be separated with TDEP-renormalized eigenvectors, shown in Figs. 3(b) and 3(c). This phenomenon is not a coincidence; it can be obviously observed from the dispersion curve of PbTe that these two phonon modes have very similar phonon frequencies. If these details are ignored, the fitting will then result in an underestimated lifetime because of the artificially increased peak width due to overlapping. By using phonon eigenvectors to project SED into different modes, two overlapped peaks are separated and then each of the peaks can be fitted separately.

Second, it is challenging to obtain the exact temperature-dependent eigenvectors that perfectly separate the SED peaks without any residual side peak. Since the MD simulation time is finite, it is difficult to ensure ergodicity in the accessible vibration phase space, and the eigenvectors obtained by TDEP could intrinsically contain small uncertainties. Using these inexact eigenvectors will cause numerically artificial modes, which are essentially leaked from other modes in this decomposition process. Figure 3 shows an example of decomposing the SED using the TDEP eigenvectors. Although a majority of the SED is successfully separated, a small amount of SED that should be in the optical transverse mode is partially projected onto an acoustic transverse residual mode and vice versa. In Fig. 9 in Appendix A, we show an example of the SED spectrum for residual peaks at different wavevectors on the same branch. We find that the modes further from the center of the Brillouin zone generate residual peaks with larger magnitudes and leak more SED information from the main peaks. This is reasonable because the effect of anharmonicity is stronger away from the zone center (see Figs. 1 and 2). Unfortunately, in our calculations, none of the eigenvectors from HA, QHA, and TDEP can perfectly separate branches without generating residual peaks.

Third, linewidths of long-lived low-frequency eigenmodes could be very small and the uncertainty of lifetimes could be large due to the scarce of data points around the peaks (see Figs. 10 and 11 in Appendixes B and C using Si as an example). Although using a longer time of trajectory sampling helps to increase the frequency resolution, fitting the SED peaks with the Lorentzian function would still result in uncertainties in phonon lifetime or linewidths when the sharp peak is accompanied by a residual side peak.

Figure 4 shows an example of the aforementioned numerical challenges in NMA methods and their consequences in the calculation of phonon lifetime. Figures 4(a) and 4(b) show that residual peaks appear in the mode-dependent SED spectrum when using TDEP-resolved eigenvectors to separate the SED spectrum. If one brutally deals with this issue in the SED method by ignoring the residual side peaks, the inconsistency of phonon lifetimes will be observed between the TD-NMA and SED as shown in Figs. 4(c)–4(e). Such inconsistency could be attributed to that (1) the major SED peak might not be the “main” peak shape since a small part of the spectral energy could be leaked into other eigenmodes as well due to the numerical inaccuracies of eigenvectors, and (2) it is difficult to judge the convergence of integration when computing lifetime in TD-NMA since the autocorrelation function carried a small amount of information of other modes. Due to these two limitations of SED and TD-NMA, we found that the lifetimes of optical modes are severely overestimated and the lifetimes of acoustic modes are underestimated using TD-NMA, while the lifetimes are more scattered in SED as shown in Fig. 4(e). Although the TD-NMA and SED methods are mathematically equivalent,^{34,35} this fitting example demonstrates that extracting lifetime using TD-NMA or SED method is sensitive to the eigenvectors, but obtaining a perfect eigenvector using a finite amount of MD trajectory is challenging. In Sec. III C, we resolve this challenge based on the idea that the total SED should be identical regardless of the set of orthonormal vectors for mode decomposition.^{39}

### C. The “Sum-up Spectrum Fitting Method” (SSFM)

This section proposes a practical strategy to resolve the numerical challenges of obtaining the temperature-dependent eigenvectors and lifetimes using MD, which is called the “Sum-up Spectrum Fitting Method” (SSFM). This idea is based on the fact that the total SED is invariant, whatever the choice of orthonormal vectors for decomposition.^{39} Figure 5 demonstrates the SSFM process. We can first apply eigenvectors to separate overlapped peaks and attribute these peaks to different phonon modes of different branches, four out of six phonon modes with their main peaks shown in Fig. 5(a) and their corresponding residual peaks shown in Fig. 5(b). As expected, the residual peaks sit on other existing phonon modes. Based on the idea of SED invariance of decomposition, we can restore the perfectly decomposed peak shape by summing over the main SED peaks associated with the eigenvectors and the residual side peaks leaked from other modes but sitting at the same frequencies.

The SSFM provides a “post-processing” approach to correctly separate overlapping peaks without residual peaks, which makes the fitting process more robust. This strategy will relax the requirement for accurate anharmonicity-renormalized phonon eigenvectors. Also, the SSFM is helpful to relax the fitting requirement of sharp peaks on phonon modes with high lifetimes in SED. Figure 6(a) compares phonon lifetimes of PbTe obtained by the SSFM with these by the original SED fitting procedure at 300 K. Comparing to fitting two degenerate transverse branches separately and ignoring the residual peaks in the original SED fitting procedure, in SSFM, we sum up the spectrum of these two branches, fit the combined spectral energy distribution, and assume these two branches have the same phonon lifetimes. As a statistical fitting method, SSFM allows information of residual peaks to be used in the fitting and allows averaging the lifetimes from degenerate branches, which reduces the uncertainty of the numerical fitting and provides more accurate phonon lifetimes. Figure 6(b) compares the lifetimes computed by the SSFM with the reported values from the SED method by Qiu *et al.*^{38} and the ALD method by Puligheddu *et al.*^{31} Using the same interatomic potential to describe PbTe, our lifetime result is consistent with both literature data.^{31,38} The literature data by Puligheddu *et al.*^{31} explicitly consider both the third-order and fourth-order phonon scattering processes using the ALD approach. To match the anharmonic frequency, we used the frequency shift data in their original publications to re-plot their lifetime data. The original data in Qiu's *et al*.’s paper^{38} use the SED method but do not consider the anharmonic phonon frequency shift at 300 K. For a better visual comparison, we also re-plot their lifetime as a function of frequency with frequency shift directly captured by the SED method from 0 to 300 K. The magnitudes of our phonon lifetimes match well with the predictions from Qiu *et al.* and acoustic part of the data from Puligheddu *et al*. This comparison confirms the accuracy of the phonon lifetimes predicted by our SED-SSFM. Moreover, we calculated thermal conductivity using the group velocities calculated from the continuous SED dispersion curve and lifetimes predicted from the SSFM and compared them to the thermal conductivity calculated by the EMD method. The EMD-predicted thermal conductivity is 2.54 W/mK; the thermal conductivity prediction from the original SED method is 2.85 W/mK, whereas the one from the SSFM method is 2.58 W/mK. The prediction from the SSFM method agrees excellently well with the EMD-predicted thermal conductivity.

### D. Roles of phonon eigenvectors

We can now examine the roles of phonon eigenvectors in the NMA methods. In contrast with the SED method, TD-NMA intrinsically has a higher requirement for accurate eigenvectors. In the time domain, applying the correct eigenvectors is supposed to be the only approach to perfectly separate one phonon mode from the rest of them. Figure 7(a) shows the comparison of phonon lifetimes predicted using different eigenvectors in the TD-NMA methods. Because phonon lifetimes by TD-NMA with HA eigenvectors are calculated using HA phonon frequencies, they have higher frequencies than other results. Phonon lifetimes calculated from TDEP are slightly higher than those from HA and QHA eigenvectors. It seems that the phonon lifetimes calculated from TD-NMA have large variations, especially when compared with lifetimes of the same mode from the SED method. As analyzed in Sec. III C, the underestimation of lifetimes for acoustic modes and the overestimation of lifetimes for optical modes are due to the inaccurate anharmonic eigenvectors so that residual modes are counted in the TD-NMA fitting.

Figure 7(b) shows the role of eigenvectors in frequency-domain NMA methods in predicting phonon lifetimes at 800 K in PbTe. Without any eigenvectors, the phonon lifetimes at the low frequencies are consistent with using other eigenvectors but the lifetimes are lower at the high frequencies. The lower lifetimes are all due to the overlapping of the SED peaks as shown in Fig. 3(a). The lifetimes calculated by the original fitting method scatter for the acoustic modes, which is attributed to the inaccurate phonon eigenvectors and the ignorance of residual peaks. Using the SSFM, the phonon lifetimes are consistent with each other using eigenvectors from HA, QHA, and TDEP. This is because eigenvectors from HA, QHA, and TDEP can successfully separate transverse and longitudinal branches in this material. Eigenvectors from HA are good enough for SED fittings with the SSFM, which significantly relaxes the requirement for accurate anharmonic eigenvectors at high temperatures in PbTe. If there is a phase transition between 0 K and the target temperature for other highly anharmonic materials, we believe that a set of QHA eigenvectors at the target temperature should work well for predicting phonon lifetimes in materials at the high-temperature phase using the SSFM.

To further check the accuracy of our phonon lifetime calculations for PbTe, phonon lifetimes by different methods are used to calculate lattice thermal conductivities under the BTE framework. The same group velocities from the continuous SED dispersion curve were used in all comparisons. Figures 7(c) and 7(d) show the comparison of the BTE thermal conductivities using lifetimes by TD-NMA and SED using different eigenvectors. These values are also compared with the EMD-predicted thermal conductivity, which should serve as the baseline. In the TD-NMA method, using the TDEP-renormalized eigenvectors, the calculated phonon lifetimes lead to the closest BTE thermal conductivity to the EMD thermal conductivity. At the same time, we must admit that it is challenging to obtain accurate eigenvectors (to completely avoid residual peaks) for TD-NMA calculations at high temperatures. In the SED methods, although lifetimes by SED without eigenvectors are smaller than the correct values due to the overlapping peaks at high frequencies as shown in Fig. 7(b), they do not significantly change the final BTE thermal conductivities in PbTe. These optical modes with small lifetimes and group velocities have limited effects on the calculations of BTE thermal conductivities. The inaccurate lifetimes from the original SED fitting lead to an overestimation of thermal conductivity. A similar comparison at 300 K is shown in Fig. S2 in the supplementary material. The SED with SSFM is more robust and insensitive to the uncertainties of phonon eigenvectors, compared with the TD-NMA method, for calculating phonon lifetimes and thermal conductivity.

Figure 8 compares the EMD-predicted thermal conductivity of PbTe with reference values and calculations using both MD-based NMA methods and ALD-BTE approaches at different temperatures.^{38} When compared with the literature value, our prediction is within the uncertainty range. Thermal conductivity predicted using the ALD-BTE approach using only the third-order IFCs is higher than the EMD value (see details in Fig. S3 in the supplementary material). This observation was validated by a recent work by Puligheddu *et al.*,^{31} who shows that the ALD-BTE can predict the correct thermal conductivity of PbTe with temperatures higher than 300 K only when the fourth-order phonon scattering processes are included. By using lifetimes from TD-NMA with TDEP-renormalized eigenvectors and SSFM with HA eigenvectors, BTE thermal conductivities are calculated from 300 to 800 K. These lifetimes, considering all phonon–phonon scatterings, provide thermal conductivities within the uncertainty range of EMD thermal conductivities.

## IV. SUMMARY AND OUTLOOK

Based on the results presented above, we have the following practical recommendations for predicting phonon transport properties of strongly anharmonic crystals at high temperatures using MD-based NMA methods:

At sufficiently high temperatures for anharmonic materials, the phonon renormalization must be considered for the shifts of phonon frequency, phonon eigenvector, and phonon group velocity when the QHA is not enough. Either the TDEP or SED method can be employed for this purpose.

For materials with weak anharmonicity, TD-NMA with phonon renormalized eigenvectors and frequencies can provide a reasonably satisfying phonon lifetime result. The TDEP or other similar approaches

^{31}can implement the renormalization scheme, but for strongly anharmonic phonon modes, calculating the accurate renormalized anharmonic eigenvectors and then computing the phonon lifetimes remains a challenge.For materials with strong anharmonicity, the SED method with the proposed SSFM fitting strategy is recommended. SSFM is to project the total phonon spectrum including all phonon modes onto separate modes using a set of inaccurate (harmonic or quasi-harmonic) phonon eigenvectors and then manually sum up the peaks that belong to the same phonon mode (at the same frequency). In agreement with the previous belief,

^{34}eigenvectors are still necessary to ensure the separation of overlapping peaks and the harmonic or quasi-harmonic eigenvectors should be good enough for this purpose. The proposed SSFM fitting strategy deals with the residual peaks generated when the eigenvectors are not accurate enough. In this perspective, even though theoretically equivalent, the SED method is practically more flexible and robust than the TD-NMA method when computing phonon lifetimes of strongly anharmonic materials.

We showed that phonon transport properties can be accurately calculated for strongly anharmonic crystals at high temperatures using MD-based simulations if the interatomic potential is known and accurate. If we combine the methods described in this work with the recent research efforts toward data-driven machine-learning-based MD force fields with the first-principles accuracy,^{54} a better framework considering all orders of anharmonicity can be established to study transport properties in strongly anharmonic materials at high temperatures, especially when the DFT-ALD-BTE approach becomes challenging.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the comparison between QHA and TDEP for solving PbTe eigenvectors at 300 and 800 K, the comparison between different SED methods for cumulative BTE thermal conductivities of PbTe at 300 K, and phonon lifetimes calculated using ALD with third-order force constants.

**ACKNOWLEDGMENTS**

The authors would like to thank Xin Qian and Abel Carreras for helpful discussions and thank Olle Hellman for providing access to the TDEP program. Computational resources were provided by the High-Performance Computing Center at North Carolina State University and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (NSF) (Grant No. ACI-1548562). The authors acknowledge financial support from the National Science Foundation (NSF) under Award Nos. CBET 1943813 and DMR 2011978. This work was also supported by the Faculty Research and Professional Development Fund from North Carolina State University.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: THE “RESIDUAL PEAK” ISSUE IN THE SED FITTING OF PbTe – WAVEVECTOR DEPENDENCE

Residual peaks are generated if the eigenvectors are not accurate enough when separating overlapping branches. Figures 9(a) and 9(b) show that the residual peaks appear in the fitting process of PbTe of the optical and acoustic transverse branches at different wave vectors at 800 K. To show the details of the residual peaks, the main peaks around 3 THz for the optical branch and 1 THz for the acoustic branch are not shown in the figure. In the same branch, we picked modes from the first, fifth, and eighth wave vectors [i.e., ** k1**(0.0625, 0, 0),

**(0.25, 0, 0), and**

*k*4**(0.5, 0, 0)] in the Brillouin zone of PbTe to do the SED fitting using eigenvectors. The modes further from the center of the Brillouin zone generate residual peaks with larger magnitudes and leak more SED information from the main peaks. This is reasonable because the anharmonicity effect is stronger away from the zone center.**

*k*8^{39}

### APPENDIX B: APPLICATION OF NMA METHODS TO Si CRYSTAL AT 300 K

We calculated the phonon lifetimes of Si using the Tersoff potential.^{55} The domain size of bulk Si is 8 × 8 × 8 unit cells (typically 4096 atoms). Periodic boundary conditions were applied to the simulation domain. The time step was 0.5 fs. Twenty independent cases were simulated with different initial random velocities to reduce the statistical uncertainties for TD-NMA fitting.^{45,46} Figure 10 shows the phonon lifetimes of Si at 300 K calculated from ALD, TD-NMA, and SED. The ALD-derived lifetimes were from the ShengBTE software by inputting the third-order interatomic force constants (IFCs), where the quasi-harmonic approximation (QHA) was applied to the lattice.^{53} This result can serve as a baseline for comparison because the effects of high-order phonon scattering in Si at 300 K are negligible. Due to the limited 8 × 8 × 8 simulation domain size, the sampling of phonon modes is limited. When comparing the lifetimes predicted from the SED method with the ALD_QHA baseline, we find that most of the lifetimes are matched with only a few exceptions at the very-low-frequency range. We attribute this discrepancy to the resolution in the peak fitting for a very long phonon lifetime. When comparing the lifetimes predicted from the TD-NMA method with the ALD_QHA baseline, we find that the exceptions are at the very high-frequency range. Because the TD-NMA needs accurate phonon eigenvectors and eigenvectors at the zone boundary are sensitive to the atomic environment, we believe this discrepancy results from the statistical errors of the eigenvectors at the zone boundary with high phonon frequencies. Residual peaks will be generated due to the inaccuracy of the eigenvectors.

### APPENDIX C: THE “SHARP PEAK” ISSUE IN THE SED FITTING OF Si

Figure 11 shows the comparison of SED fittings with and without using the Sum-up Spectrum Fitting Method (SSFM) in phonon lifetime calculations and analysis of the artifacts in the SED fitting. Figure 11(a) shows the distribution of SED of Si at 300 K where this specific mode seems to possess a lower lifetime at the low-frequency range in Fig. 10 (dashed circle). When the phonon lifetime is very high, much larger than 100 ps, the corresponding peak width in SED (frequency domain) is very small. We notice that even a slight difference in the fitting could result in a substantial variation in the lifetime value. For example, the red fitting curve is supposed to be correct with the 370 ps lifetime. If the fitting curve is just slightly different from the red one, such as the blue fitting line, the fitted lifetime could be changed from 370 to 13 263 ps. If one of the peak data is missing or an automatic fitting program is applied here, the lifetime would become 170 ps, which is the one shown in Fig. 11. Therefore, to further improve the resolution of the SED on these modes with extremely high lifetimes for a single-peak fitting, a much smaller time step or a significantly longer sampling time should be set in MD simulations, to increase the number of data points in the frequency domain, which will substantially increase the computation cost for all modes. As an alternative option without increasing the computation cost, the SSFM is helpful to relax the fitting requirement of sharp peaks on modes with high lifetimes in SED by using the artificial information in SED with inaccurate eigenvectors. By summing up acoustic and optical transverse branches of Si at 300 K, more data points are available for fitting in Fig. 11(b) compared with Fig. 11(a). The extra data from SED result in a better statistic fitting for phonon lifetime even within this shape peak.