Computational prediction of lattice thermal conductivity -- a comparison of molecular dynamics and Boltzmann transport approaches

TL;DR

This paper compares molecular dynamics and Boltzmann transport approaches for predicting lattice thermal conductivity, highlighting their similarities, differences, and the effects of phonon statistics and scattering mechanisms.

Contribution

It provides a detailed comparison of MD and BTE methods, including phonon scattering, renormalization, and quantum vs. classical statistics, revealing the origins of their agreement.

Findings

MD and BTE results show apparent agreement due to error cancellations.

Phonon statistics significantly affect scattering rates and thermal conductivity.

Equipartition governs phonon statistics in MD simulations.

Abstract

The predictive modeling of lattice thermal conductivity is of fundamental importance for the understanding and design of materials for a wide range of applications. Two major approaches, namely molecular dynamics (MD) simulations and calculations solving approximately the Boltzmann transport equation (BTE), have been developed to compute the lattice thermal conductivity. We present a detailed direct comparison of these two approaches, using as prototypical cases MgO and PbTe. The comparison, carried out using empirical potentials, takes into account the effects of fourth order phonon scattering, temperature-dependent phonon frequencies (phonon renormalization), and investigates the effects of quantum vs. classical statistics. We clarify that equipartition, as opposed to Maxwell Boltzmann, govern the statistics of phonons in MD simulations. We find that lattice thermal conductivity values from MD and BTE show an apparent, satisfactory agreement; however such an agreement is the result of error cancellations. We also show that the primary effect of statistics on thermal conductivity is via the scattering rate dependence on phonon populations.

Tables

Table 1. Table 1: Computed thermal conductivity of MgO and PbTe using different expressions for lifetimes ( τ 𝜏 \tau ) and heat capacities ( C 𝐶 C ). All thermal conductivity values are in W/mK . The lifetimes τ M ​ D subscript 𝜏 𝑀 𝐷 \tau_{MD} ’s are obtained from MD simulations according to Equation 6 , τ B ​ E subscript 𝜏 𝐵 𝐸 \tau_{BE} ’s are obtained from BTE according to phonons in a Bose Einstein distribution, whereas τ E ​ Q subscript 𝜏 𝐸 𝑄 \tau_{EQ} ’s are obtained from BTE according to phonon population obeying equipartition. All BTE calculations include three- and four-phonon processes. The heat capacities C E ​ Q subscript 𝐶 𝐸 𝑄 C_{EQ} and C B ​ E subscript 𝐶 𝐵 𝐸 C_{BE} correspond to those obtained from equipartition and Bose-Einstein statitics, respectively, as shown in Fig. 3 . Finally, the thermal conductivity obtained from Green-Kubo (GK) is listed. The different approximations are labeled A1 to A7 in the first row.

		A1	A2	A3	A4	A5	A6	A7
Compound	Temperature	${\tau }_{MD}$ $C_{EQ}$	${\tau }_{MD}$ $C_{BE}$	${\tau }_{BE}$ $C_{EQ}$	${\tau }_{BE}$ $C_{BE}$	${\tau }_{EQ}$ $C_{EQ}$	${\tau }_{EQ}$ $C_{BE}$	GK
MgO	500K	29.4	27.0	35.5	32.8	27.1	25.2	32.6 (2.1)
MgO	750K	16.3	15.6	20.9	20.2	17.4	16.8	19.6 (0.5)
MgO	1000K	10.4	10.2	14.7	14.4	12.6	12.4	14.1 (0.7)
PbTe	100K	11.2	10.7	13.0	12.4	10.2	9.8	10.4 (0.8)
PbTe	150K	7.0	6.8	7.9	7.7	6.6	6.5	6.5 (0.6)
PbTe	300K	3.1	3.1	3.4	3.4	3.0	3.0	3.2 (0.2)

Equations

$$\kappa =\frac{1}{3V{k}_B{T}^2}{\int }_0^{\infty }⟨J(0)\cdot J(t)⟩\text{d}t$$

$$E_{\lambda }(t)=K_{\lambda }(t)+U_{\lambda }(t)$$

$$K_{\lambda }(t)=\frac{1}{2N}\Vert {\dot{q}}_{\lambda }{\Vert }^2{U}_{\lambda }(t)=\frac{1}{2N}{\omega }_{\lambda }^2\Vert q_{\lambda }{\Vert }^2$$

$${\dot{q}}_{\lambda }=\sum _{l,b}\sqrt{\frac{m_b}{N}}{e}^{i\vec{k}{r}_0(l,b)}{\vec{e}}_{b,\lambda }^{\ast }\cdot v(l,b,t)$$

$$q_{\lambda }=\sum _{l,b}\sqrt{\frac{m_b}{N}}{e}^{i\vec{k}{r}_0(l,b)}{\vec{e}}_{b,\lambda }^{\ast }\cdot u(l,b,t)$$

$${\tau }_{\lambda }={\int }_0^{\infty }\frac{⟨E_{\lambda }(0)E_{\lambda }(t)⟩}{⟨E_{\lambda }(0)E_{\lambda }(0)⟩}\text{d}t$$

$$\frac{⟨E_{\lambda }(0)E_{\lambda }(t)⟩}{⟨E_{\lambda }(0)E_{\lambda }(0)⟩}=e^{-t/{\tau }_{\lambda }}$$

$$H_{harm}(\{q,\dot{q}\})=\sum_{\lambda}\big{[}\frac{1}{2N}\|\dot{q}_{\lambda}\|^{2}+\frac{1}{2N}\omega_{\lambda}^{2}\|q_{\lambda}\|^{2}\big{]}$$

$$⟨E_{\lambda }⟩=k_{B}T⟨n_{\lambda }⟩=\frac{k_{B}T}{\hslash {\omega }_{\lambda }}$$

$$⟨E_{\lambda }⟩=\frac{1}{\Gamma }{\int }_0^{\Gamma }{E}_{\lambda }(t)\text{d}t$$

$$⟨n_{\lambda }⟩=\frac{⟨E_{\lambda }⟩}{\hslash {\omega }_{\lambda }}$$

$$\frac{⟨K_{\lambda }(0)K_{\lambda }(t)⟩}{⟨K_{\lambda }(0)K_{\lambda }(0)⟩}=e^{-t/{\tau }_{\lambda }}{\cos }^2{\omega }_{\lambda }(T)$$

$$\frac{\partial n_{\lambda }}{\partial t}(\mathbf{r})=\frac{\partial n_{\lambda }}{\partial t}(\mathbf{r}{)}_{\text{diff}}+\frac{\partial n_{\lambda }}{\partial t}(\mathbf{r}{)}_{\text{scatt}}+\frac{\partial n_{\lambda }}{\partial t}(\mathbf{r}{)}_{\text{ext}}=0.$$

$$-\frac{\partial n_{\lambda }}{\partial t}(\mathbf{r}{)}_{\text{scatt}}=\frac{n_{\lambda }-n_{\lambda }^0}{{\tau }_{\lambda }},$$

$$n_{\lambda }-n_{\lambda }^0={\mathbf{v}}_{\lambda }\frac{\partial n_{\lambda }}{\partial T}\nabla T{\tau }_{\lambda },$$

$$J^{\alpha }=-\sum _{\beta }{\kappa }_l^{\alpha \beta }(\nabla T{)}^{\beta },$$

$$J^{\alpha }=\sum _{\lambda }\int n_{\lambda }\hslash {\omega }_{\lambda }{\mathbf{v}}_{\lambda }^{\alpha }\frac{d\mathbf{k}}{(2\pi {)}^3},$$

$${\kappa }_l^{\alpha \beta }=\frac{1}{NV{k}_{\text{B}}{T}^2}\sum _{\lambda }{n}_{\lambda }^0(n_{\lambda }^0+1)(\hslash {\omega }_{\lambda }{)}^2{v}_{\lambda }^{\alpha }{v}_{\lambda }^{\beta }{\tau }_{\lambda },$$

$$\begin{array}{rl}{\tau }_{3,\lambda }^{-1}& =\sum _{{\lambda }_1{\lambda }_2}\left\{\frac{1}{2}\left(1+n_{{\lambda }_1}^0+n_{{\lambda }_2}^0\right){\zeta }_{-}+\left(n_{{\lambda }_1}^0-n_{{\lambda }_2}^0\right){\zeta }_{+}\right\},\end{array}$$

$$\begin{array}{rl}{\tau }_{4,\lambda }^{-1}& =\sum _{{\lambda }_1{\lambda }_2{\lambda }_3}\left\{\frac{1}{6}\frac{n_{{\lambda }_1}^0{n}_{{\lambda }_2}^0{n}_{{\lambda }_3}^0}{n_{\lambda }^0}{\zeta }_{--}+\frac{1}{2}\frac{\left(1+n_{{\lambda }_1}^0\right)n_{{\lambda }_2}^0{n}_{{\lambda }_3}^0}{n_{\lambda }^0}{\zeta }_{+-}+\frac{1}{2}\frac{\left(1+n_{{\lambda }_1}^0\right)\left(1+n_{{\lambda }_2}^0\right)n_{{\lambda }_3}^0}{n_{\lambda }^0}{\zeta }_{++}\right\},\end{array}$$

$$\begin{array}{rl}{\zeta }_{\pm }& =\frac{\pi \hslash }{4N}|V_{\pm }^{(3)}{|}^2{\Delta }_{\pm }\frac{\delta ({\omega }_{\lambda }\pm {\omega }_{{\lambda }_1}-{\omega }_{{\lambda }_2})}{{\omega }_{\lambda }{\omega }_{{\lambda }_1}{\omega }_{{\lambda }_2}},\end{array}$$

$$\begin{array}{rl}{\zeta }_{\pm \pm }& =\frac{\pi {\hslash }^2}{8N^2}|V_{\pm \pm }^{(4)}{|}^2{\Delta }_{\pm \pm }\frac{\delta ({\omega }_{\lambda }\pm {\omega }_{{\lambda }_1}\pm {\omega }_{{\lambda }_2}-{\omega }_{{\lambda }_3})}{{\omega }_{\lambda }{\omega }_{{\lambda }_1}{\omega }_{{\lambda }_2}{\omega }_{{\lambda }_3}},\end{array}$$

$$\begin{array}{rl}{V}_{\pm }^{(3)}& =\sum _{b,l_1{b}_1,l_2{b}_2}\sum _{\alpha {\alpha }_1{\alpha }_2}{\Phi }_{0b,l_1{b}_1,l_2{b}_2}^{\alpha {\alpha }_1{\alpha }_2}\frac{e_{\alpha b}^{\lambda }{e}_{{\alpha }_1{b}_1}^{\pm {\lambda }_1}{e}_{{\alpha }_2{b}_2}^{-{\lambda }_2}}{\sqrt{m_b{m}_{b_1}{m}_{b_2}}}{e}^{\pm i{\mathbf{k}}_1\cdot {\mathbf{r}}_{l_1}}{e}^{-i{\mathbf{k}}_2\cdot {\mathbf{r}}_{l_2}},\end{array}$$

$$\begin{array}{rl}{V}_{\pm \pm }^{(4)}& =\sum _{b,l_1{b}_1,l_2{b}_2,l_3{b}_3}\sum _{\alpha {\alpha }_1{\alpha }_2{\alpha }_3}{\Phi }_{0b,l_1{b}_1,l_2{b}_2,l_3{b}_3}^{\alpha {\alpha }_1{\alpha }_2{\alpha }_3}\frac{e_{\alpha b}^{\lambda }{e}_{{\alpha }_1{b}_1}^{\pm {\lambda }_1}{e}_{{\alpha }_2{b}_2}^{\pm {\lambda }_2}{e}_{{\alpha }_3{b}_3}^{-{\lambda }_3}}{\sqrt{m_b{m}_{b_1}{m}_{b_2}{m}_{b_3}}}{e}^{\pm i{\mathbf{k}}_1\cdot {\mathbf{r}}_{l_1}}{e}^{\pm i{\mathbf{k}}_2\cdot {\mathbf{r}}_{l_2}}{e}^{-i{\mathbf{k}}_3\cdot {\mathbf{r}}_{l_3}}.\end{array}$$

$$n_{\lambda }^0=\frac{1}{e^{\hslash {\omega }_{\lambda }/k_{\text{B}}T}-1},$$

$$C_{\lambda }=k_{\text{B}}{\left(\frac{\hslash {\omega }_{\lambda }}{k_{\text{B}}T}\right)}^2\frac{e^{\hslash {\omega }_{\lambda }/k_{\text{B}}T}}{{\left(e^{\hslash {\omega }_{\lambda }/k_{\text{B}}T}-1\right)}^2}.$$

$$\begin{array}{rl}{\mathbf{\Phi }}^{\text{CS}}& ={\arg \min }_{\mathbf{\Phi }}\Vert \mathbf{\Phi }{\Vert }_1+\frac{\mu }{2}\Vert \mathbf{F}-\mathbb{A}\mathbf{\Phi }{\Vert }_2^2,\end{array}$$

Full text

I Introduction

Thermal transport properties are important for myriad materials, including thermoelectrics, (opto)electronics, photovoltaic and photoelectrochemical cells, and batteries. To facilitate fundamental understanding and materials design, the accurate prediction of thermal transport coefficients is critical. A multitude of methodologies have been introduced to facilitate such prediction. Of particular significance are two categories of methods, based either on considering atomic motion (molecular dynamics, MD) or collective vibrational excitations, i.e. phonons in solids (anharmonic lattice dynamics, ALD).

Over the years, several MD-based approaches have been developed to compute the thermal conductivity in bulk crystalline and amorphous materials, nanostructures, and fluids, with the central premise of tracking heat flow through atomic motion. In the Green-Kubo (GK) method, based on the fluctuation-dissipation theorem, the thermal conductivity is calculated from the fluctuations of the heat correlation function during an equilibrium simulation, in the NVE ensemble. Non Equilibrium (NE) methods compute the thermal conductivity in the steady state from the response of the system to a perturbation. They can be classified according to the nature of the perturbation: in the older direct method it is a temperature gradient causing a heat flux Evans and Morriss (2007); in reverse NEMD, also known as the Muller-Plathe approach,the constant heat flux acts as the perturbation Muller-Plathe (1997). The ”approach to equilibrium” methods compose the third class of MD based method. In these approaches, the thermal conductivity is computed from the time response of the system to an instantaneous perturbation. This perturbation can be in the form of a square Lampin et al. (2013); Melis et al. (2014) or sinusoidal temperature profile Puligheddu, Gygi, and Galli (2017). Meanwhile, phonon-based ALD approaches for calculating lattice thermal conductivity, applicable to crystalline systems, have been developed based on the Boltzmann transport equation (BTE), with the central assumption that collective vibrations can be thought of as quasi-particles termed phonons which have explicit group velocities and scatter with each other Ziman (1960); Srivastava (1990); Peierls (1996). Within these approaches, the lattice thermal conductivity is calculated by solving the BTE, under the relaxation time approximation, with relaxation times obtained using perturbation theory and considering anharmonic phonon-phonon interactions, such as three- and four-phonon scattering processes Maradudin and Fein (1962); Debernardi, Baroni, and Molinari (1995); Deinzer, Birner, and Strauch (2003); Broido, Ward, and Mingo (2005); Feng and Ruan (2016). While MD and ALD are fundamentally distinct and practically disparate approaches, they have nonetheless both been shown to give reasonably accurate predictions for similar systems Turney et al. (2009); He et al. (2012); Feng and Ruan (2016).

Typically, GK calculations are performed using classical MD where the potential energy is given by an interatomic potential as a function of the positions of the atoms. Recently, Marcolongo et. al. Marcolongo, Umari, and Baroni (2015) and Carbogno et al. Carbogno, Ramprasad, and Scheffler (2017) showed that the Green-Kubo formalism can be applied to ab initio MD, although the method of Ref. Carbogno, Ramprasad, and Scheffler (2017) is limited to solids and contains approximations on the partition of slow and fast vibrational modes. We note that the Green-Kubo formalism does not rely on the definition of phonons and is generally applicable: indeed it can describe conductive heat transport in crystalline or amorphous materials, materials with defects, nanostructures, and fluids. The GK method implicitly accounts for all anharmonic terms in the potential energy. In BTE language, this would translate to including all phonons processes in the perturbative expansion (instead of only up to the third or fourth order), thus allowing for greater accuracy at high temperatures. Nonetheless, the GK method suffers from three deficiencies: (i) being based on classical MD, it is not known how to include quantum effects such as zero point energy or quantization of the phonon energy levels, leading to lower accuracy at low temperatures. This problem is of a theoretical nature and, to our knowledge, has not been solved, although semi-empirical corrections have been proposed Wang, Chan, and Ho (1990); Lee et al. (1991); (ii) long simulation times (order of ns) are needed for convergence, especially for systems with high values of thermal conductivity, mostly due to noise in the long time tail of the heat flux autocorrelation function. Acceleration is possible based on the methods proposed in Ref. Ercole, Marcolongo, and Baroni (2017) or Carbogno, Ramprasad, and Scheffler (2017), although the latter one is limited to crystals and based on phonon theory; and (iii) the value of the thermal conductivity calculated in a GK simulation is a function of the number of atoms included in the MD cell, converging in the infinite limit to the bulk value. Currently, no theory or model exists to predict or extrapolate to the thermodynamic limit, leading to the necessity of simulating at various cell sizes with increasing numbers of atoms until satisfactory convergence is obtained. We refer to He et al. (2012) for a detailed study.

The BTE-ALD approach is more advantageous compared with MD-based method in the following respect: (i) the quantum (Bose-Einstein) statistics is utilized, thus (presumably) attaining better accuracy at low temperatures; (ii) first-principles calculations may be used to obtain third and fourth order force constants in a relatively straightforward manner for a variety of systems, to accurately model heat transport in bulk crystals with relatively low computational costs, and computation for larger systems in some cases can be parallelized Debernardi, Baroni, and Molinari (1995); Broido et al. (2007); Li et al. (2014); Zhou et al. (2014); and (iii) convergence issues with respect to system size and sampling are less severe than in MD. However, BTE-ALD suffers from several deficiencies: (i) the approach depends on the assumption that lattice vibrations can be treated as quasiparticles, i.e., both the anharmonic frequency shift and broadening (real and imaginary part of the phonon self-energy) are relatively small compared to the phonon frequency, and fails when quasiparticles are not well-defined, such as in the localization limit in systems with strong intrinsic disorder Allen et al. (1999); Seyf et al. (2017); (ii) the approach is based on perturbation theory, which is accurate only in the perturbative limit and may fail in severely anharmonic systems or at high temperatures since the scattering processes are limited to third or at most fourth order; and (iii) it is difficult to use the phonon BTE method to model heat transport in non-crystalline systems including defects and nanostructures, and the approach cannot be applied to fluids or even solids when diffusion occurs.

In an earlier study Turney et al. (2009), Turney et al. presented and compared methods in several categories: (1) quasiharmonic and anharmonic lattice dynamics calculations, (2) a combination of quasiharmonic lattice dynamics calculations and molecular dynamics simulations, and (3) Green-Kubo and direct molecular dynamics, to assess their validity. They pointed out that the lattice dynamics calculations tend to underestimate lattice thermal conductivity at above half of the Debye temperature. However, their lattice dynamics calculations excluded higher-than-third-order anharmonic phonon-phonon interactions. In addition, the impact of different statistics (quantum vs. classical) on calculated thermal conductivity was not clarified. In a follow-on study, the same authors Turney, McGaughey, and Amon (2009) assessed different corrections to MD simulations to account for quantum statistics, using silicon and the Stillinger Weber potential. They found that these corrections failed at low temperatures due to the classical distribution of phonon modes. He et. al.He et al. (2012) compared the GK-MD and BTE approaches for the computation of thermal conductivity of Si, Ge, and Si-Ge alloys, and found that the results are consistent for the pure compounds, though alloy systems prove to be problematic for BTE. In this study, some investigations into the effects of quantum statistics were carried out via the use of classical statistics in BTE. However, phonon occupation numbers were assumed to obey the Maxwell-Boltzmann distribution, which we will show below to be incorrect. The effects of temperature-induced anharmonic phonon renormalization or higher-than-third-order interactions were again not considered. Since several of the previous studies were carried out on Si or related systems, it would also be of interest to test the results on more diverse types of systems.

In the current work, we aim to provide a controlled, comprehensive, and systematic comparison of MD vs. ALD-BTE based approaches for the prediction of lattice thermal conductivity. In order to accomplish a meaningful comparison, we select two representative systems: a small-gap semiconductor with low lattice thermal conductivity (PbTe), and an insulator with higher lattice thermal conductivity (MgO). These two materials have the further advantage of having simple structures. For each material, we use MD and BTE-ALD, with the same interatomic potential, to evaluate the temperature-dependent values of the lattice thermal conductivity well below, near, and well above their respective Debye temperatures. We further investigate the source of the differences between the two approaches by deriving lifetimes from MD simulations, enumerating heat capacities and phonon occupancies under different treatment of statistics, and using these quantities in the BTE expression. The goal is to determine the underlying physical reasons for divergences and convergences between the two approaches.

II Methods

II.1 Model Systems and Interatomic Potentials

To ensure that any differences in computed thermal conductivity arise solely from the treatment of heat transport (i.e. atoms vs. phonons), we used the same interatomic potential for both MD and BTE calculations. For MgO, we used the potential described in Ref. Shukla et al. (2008); for PbTe, we used the potential from Ref. Qiu et al. (2012). Both potentials are of the Buckingham-Coulomb type. The MgO potential was shown Shukla et al. (2008) to predict experimental lattice constant, thermal expansion and thermal conductivity reasonably well in the 300K-1500K range. The PbTe potential developed by Qiu et al.Qiu et al. (2012) was shown to reasonably reproduce the mechanical and phonon properties of PbTe bulk crystal, as well as lattice thermal conductivity. We studied MgO at 500K, 750K, and 1000K and PbTe at 100K, 150K, and 300K. The experimental Debye temperatures of MgO and PbTe are 743K Beg (1976) and 177K Houston, Strakna, and Belson (1968), respectively. To account for thermal expansion, the lattice parameters computed from NPT MD simulations at each temperature were used.

II.2 Calculation of Thermal Conductivity using Molecular Dynamics

In the GK formalism, the thermal conductivity is given by

[TABLE]

where the brackets $\langle\cdot\rangle$ indicate the thermodynamic average and $J$ is the heat flux due to spontaneous fluctuations.

All classical simulations were performed using LAMMPS Plimpton (1995). The MgO samples contain 32768 (500-750K) or 4096 (1000K) atoms, with the smaller number of atoms for higher temperature due to the shorter mean free path at higher temperature. The PbTe samples contain 8192 atoms. For all systems, Nose-Hoover NVT equilibration runs of 20 ps for MgO and 100 ps for PbTe were followed by NVE simulations of 3 ns to obtain the lattice thermal conductivity ($\kappa$). For each temperature and material, 4-9 GK MD runs were performed, to give a total of 12-27 ns of statistics. The time step used for MgO is 1 fs whereas that for PbTe is 0.5 fs. The MgO 500K and 1000K data are from the Supporting Information in Ref.Puligheddu, Gygi, and Galli (2017).

II.2.1 Calculation of Phonon Lifetimes

Before starting our simulations, we compute the phonon frequencies and eigenvectors for the target temperature using the phonon renormalization (PRN) scheme described in Ref. Xia (2018). Once the phonon frequencies $\omega_{\lambda}$ and eigenvectors $e_{\lambda}$ are known, the energy of each phonon mode can be computed during a MD simulation as a function of time. The energy of a phonon mode $\lambda$ is calculated as

[TABLE]

where $K$ and $V$ are, respectively, the kinetic and potential components of the energy of the phonon mode. The kinetic component is calculated from the projection $\dot{q}$ of the phonon mode on the velocities $v$ of the atoms. The potential component is calculated by projecting onto the displacement $u$ of the atoms from their equilibrium positions $r_{0}$, i.e.

[TABLE]

[TABLE]

[TABLE]

Here, $l$ and $b$ are indices over the primitive cells and the atoms inside a primitive cell, respectively. $N$ is the number of atoms in the system, $m_{b}$ is the mass of atom $b$, and $\omega_{\lambda}$ is the frequency of the phonon mode $\lambda$. The displacement $u$ is given by $u(t)=r(t)-r_{0}(t)$. The vector $\vec{e}_{b,\lambda}$ describes the direction and phase of the displacement of atom $b$ due to the phonon with wavevector $\vec{k}$ and polarization $\lambda$.

The lifetime of each phonon mode is then calculated from the normalized autocorrelation of the energy:

[TABLE]

In order to reduce the noise due to the tail of the autocorrelation, we compute the lifetime from a fit to the equation

[TABLE]

with the lifetime $\tau_{\lambda}$ as a single fitting parameter. Lifetimes are computed as described above during an NVE simulation, once the system is prepared at the target temperature using a Nose-Hoover thermostat. The approach described to obtain lifetimes relies on the definition of phonons in the quasi-harmonic approximation.

II.2.2 Average Phonon Occupation

The classical Hamiltonian of the system can be rewritten, in the harmonic approximation, as a sum of kinetic and potential energy over phonon modes:

[TABLE]

As a consequence of the equipartition theorem, the average kinetic and potential energy for each phonon mode is equal to $k_{B}T/2$. Each phonon mode has an average energy and phonon number given by:

[TABLE]

From the procedure used to calculate the phonon lifetimes, described above, we can obtain the phonon energy averaged over time

[TABLE]

where $\Gamma$ is the total simulation time, and the mode-by-mode average phonon occupation number $\langle n_{\lambda}\rangle$ is

[TABLE]

Our results for phonon numbers for MgO and PbTe at various temperatures are shown in Figure 1, where we compare the distributions obtained in our calculations against those predicted by equipartition, as well as the Maxwell-Boltzmann (MB) and Bose-Einstein (BE) statistics. The computed distributions deviate significantly from those of the MB distribution, illustrating that classical phonons in MD are distributed not via MB but equipartition (EQ). The computed distributions deviate also from BE statistics but the agreement improves with increasing temperature, as expected.

II.2.3 Frequency Shifts

Phonon frequencies are dependent on temperature due to different degrees of atomic displacements from their equilibrium positions as a function of T, amounting to sampling different portions of the anharmonic potential energy surfaces. This change in phonon frequencies as a function of T (which we call shift in phonon frequencies relative to their T=0 value, or phonon renormalization) manifests in MD simulations. In our analysis, we compute the frequency shift from the kinetic part $K_{\lambda}(t)$ of the phonon mode energy. We fit its normalized time autocorrelation function to the equation

[TABLE]

where we added the dependence $(T)$ to the frequency $\omega_{\lambda}$. The frequency shift is then defined as the difference between the frequency obtained during MD simulations and that obtained from diagonalizing the dynamical matrix at T=0. Note that this procedure does not include effects arising from the temperature-dependent change in lattice constant, as we compute both frequencies at the same lattice constant.

II.3 Anharmonic Lattice Dynamics and Boltzmann Transport Equation

II.3.1 Boltzmann Transport Equation

The Boltzmann transport equation (BTE) can be used to describe the time evolution of the positions and momenta of a system of particles, e.g. phonons. According to the BTE, at equilibrium the evolution of the occupation probability $n_{\lambda}$ of a specific phonon mode $\lambda$ due to diffusion, scattering, and the presence of an external heat current must balance:

[TABLE]

Under the relaxation time approximation (RTA), the scattering term can be expressed as

[TABLE]

where $n_{\lambda}^{0}$ and $\tau_{\lambda}$ denote the equilibrium occupation probability and relaxation time, respectively. Under an external temperature gradient $\Delta T$, the deviation of the occupation number from its equilibrium value $\lambda$ is given by

[TABLE]

where $\mathbf{v}_{\lambda}$ is the phonon group velocity vector. In the linear regime, the lattice thermal conductivity tensor is defined by

[TABLE]

where $J^{\alpha}$ is the heat current along $\alpha$ direction and can be obtained as

[TABLE]

where $\mathbf{k}$ denotes the phonon wave vector. The resulting lattice thermal conductivity tensor within the RTA is

[TABLE]

where $k_{\text{B}}$ is the Boltzmann constant, $V$ is the volume of unit cell, $N$ is the total number of phonon wave vectors included in the summation, $\omega_{\lambda}$ and $v_{\lambda}$ are the frequency and group velocity of phonon mode $\lambda$, with Cartesian coordinates indexed by $\alpha$ and $\beta$. Typically, $\omega_{\lambda}$ and $v_{\lambda}$ are extracted from phonon dispersions in the harmonic approximation, assuming small atomic displacements, by computing second-derivatives of the potential energy with respect to atomic displacements. Within the framework of anharmonic lattice dynamics, $\tau_{\lambda}$ is assumed to arise primarily from intrinsic phonon-phonon scattering events Mahan (2000), with the lowest-order contribution being three-phonon processes.

II.3.2 Anharmonic Lattice Dynamics

Most recently, Feng and Ruan Feng and Ruan (2016, 2018a); Feng, Lindsay, and Ruan (2017) performed rigorous calculations of four-phonon scattering rates in the ALD-BTE framework including fourth order IFCs by extending the derivation of Maradudin et al. Maradudin and Fein (1962). Their studies reveal that even for diamond, silicon, and germanium, which are generally considered low-anharmonicity materials, the contribution of four-phonon scattering rates is comparable to that of three-phonon scattering rates near/above the Debye temperature. Following their derivation based on Fermi’s golden rule (FGR) Feng and Ruan (2016), the scattering rates ($\tau_{3,\lambda}^{-1}$ and $\tau_{4,\lambda}^{-1}$) associated with three- and four-phonon processes (see Fig. 2) in the single mode relaxation time approximation (SMRTA) are given by

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

where $V_{\pm}^{(3)}$ and $V_{\pm\pm}^{(4)}$ are in turn given by

[TABLE]

[TABLE]

In the above equations, $l$, $b$, and $\alpha$ are indexes of primitive cells, basis atoms, and cartesian coordinates, respectively; $m$ is the atomic mass, and $\mathbf{r}$ is the lattice vector of the primitive cell; $\mathbf{k}$, $n_{\lambda}^{0}$, $\omega_{\lambda}$, and $e^{\lambda}$ are phonon wave vector, equilibrium occupation number, frequency, and eigenvector, respectively; $\Phi_{0b,l_{1}b_{1},l_{2}b_{2}}^{\alpha\alpha_{1}\alpha_{2}}$ and $\Phi_{0b,l_{1}b_{1},l_{2}b_{2},l_{3}b_{3}}^{\alpha\alpha_{1}\alpha_{2}\alpha_{3}}$ are the third and fourth order IFCs, respectively. In the three-phonon processes, $\zeta_{-}$ represents the splitting process ($\lambda\rightarrow\lambda_{1}+\lambda_{2}$) and $\zeta_{+}$ indicates the combination process ($\lambda+\lambda_{1}\rightarrow\lambda_{2}$). Similarly, $\zeta_{--}$, $\zeta_{+-}$ and $\zeta_{++}$ represent the splitting ($\lambda\rightarrow\lambda_{1}+\lambda_{2}+\lambda_{3}$), redistribution ($\lambda+\lambda_{1}\rightarrow\lambda_{2}+\lambda_{3}$), and combination ($\lambda+\lambda_{1}+\lambda_{2}\rightarrow\lambda_{3}$) in the four-phonon processes Feng and Ruan (2016). In both three- and four-phonon processes, momentum conservation is strictly enforced as indicated by $\Delta_{\pm}$ and $\Delta_{\pm\pm}$, and energy conservation is enforced by $\delta$ functions, which are approximated by adaptive and regular Gaussian smearing Li et al. (2012, 2014) in computing $\tau_{3,\lambda}^{-1}$ and $\tau_{4,\lambda}^{-1}$, respectively.

II.3.3 Computational Details for BTE

In the BTE-ALD calculations, for both MgO and PbTe, we used 6$\times$6$\times$6 supercells to extract harmonic and anharmonic interatomic force constants (IFCs) up to the fourth order using compressive sensing lattice dynamics Zhou et al. (2014). There is no explicit cutoff distance enforced on the harmonic IFCs Parlinski, Li, and Kawazoe (1997). To further verify the extracted harmonic IFCs, we compare the calculated phonon dispersions with those independently obtained by Phonopy Togo and Tanaka (2015). The cutoff distance of the third order IFCs is limited to the seventh nearest neighbour shell, which leads to converged lattice thermal conductivity when only three-phonon interactions are accounted for. Considering the short-range nature of 4th-order anharmonicity Zhou et al. (2014) and the associated combinatorial explosion in the number of parameters, the 4th-order IFCs are limited to the second nearest neighbour shell. We also verified that the additional inclusion of the third nearest neighbour shell in the fourth order IFCs leads to negligible changes in both frequency shifts and lattice thermal conductivity. The phonon BTE with renormalized harmonic IFCs and anharmonic IFCs as input were solved using $\mathbf{q}$-point mesh of 16$\times$16$\times$16 and 12$\times$12$\times$12, respectively, which are equivalent to the supercell structures used in MD simulations.

II.3.4 Iterative and non-iterative solutions of BTE

Phonon BTE under RTA can be linearized; when solved non-self-consistently, their solutions are within the so-called single mode relaxation approximation (SMRTA) and when solved iteratively, they are within the relaxation time approximation (IRTA), respectively. Within the SMRTA, scattering rates arise from both normal ($N$) and Umklapp ($U$) processes, which lead to underestimation of $\kappa_{l}$ since $N$ processes do not introduce thermal resistance directly. For systems with significant $N$ processes, an iterative solution to the linearized BTE should be pursued by accounting for nonequilibrium phonon distribution through iteratively refining phonon populations Omini and Sparavigna (1995, 1996); Broido, Ward, and Mingo (2005). In this study, we calculate three-phonon scattering rates by solving the phonon BTE in an iterative manner. Considering the extremely high computational cost of including additional four-phonon scattering in the iterative solver, we treat four-phonon scattering non-iteratively, and combine them with three-phonon scattering rates based on Matthiessen’s rule to give the total scattering rates. We found that this strategy leads to reasonably accurate lattice thermal conductivity for compounds which are dominated by $U$ processes, such as PbTe, as discussed in detail in an earlier study Xia (2018). It is worth noting that fully iterative soluation of BTE including both three- and four-phonon scatterings has been recently developed and implemented by Feng *et al. * Feng and Ruan (2018b) and Ravichandran *et al. * Ravichandran and Broido (2018),respectively.

II.3.5 Temperature-dependent Frequencies: Anharmonic Phonon Renormalization

Several first-principles phonon renormalization (PRN) schemes based either on real or reciprocal space formalisms have been introduced to treat strong anharmonicity effects on phonon frequencies Souvatzis et al. (2009); Errea, Calandra, and Mauri (2014); Hellman, Abrikosov, and Simak (2011); Tadano and Tsuneyuki (2015); van Roekeghem, Carrete, and Mingo (2016). In this work, we utilize our recently introduced real space based PRN scheme to compute temperature-dependent phonon freqeuncies and eigenvectors Xia (2018). The temperature effects are taken into account via constructing temperature-dependent effective harmonic potential coefficients Wallace (1998) which includes the temperature-dependent corrections from higher-order IFCs on top of harmonic IFCs Xia (2018). We refer the reader to Refs. [Xia, 2018; Xia and Chan, 2018] for detailed discussions.

II.4 Statistics

In the expression of the thermal conductivity obtained from the BTE, the phonon statistics enter both the definition of the heat capacity and that of the scattering rates (inverse of relaxation time). According to Bose-Einstein (BE) statistics, the phonon population of mode $\lambda$ in equilibrium is

[TABLE]

which leads to a mode heat capacity of

[TABLE]

In the classical (equipartition) limit, each phonon mode has an energy of $k_{\text{B}}T$ and population of $k_{\text{B}}T/\hbar\omega_{\lambda}$, thus giving rise to a temperature-independent mode heat capacity of $k_{\text{B}}$. The classical limit always overestimates the heat capacity (see Fig. 3), particularly at low temperatures. Specifically in the case of MgO, the classical limit leads to overestimates in the heat capacity of 11.2%, 5.3% and 3.1% at 500 K, 750 K and 1000 K, respectively.

As rigorously derived by Feng and Ruan Feng and Ruan (2016), phonon population in the classical limit (equipartition) cannot be directly used to compute three- and four-phonon scattering rates in Eq. 19 and 20, since the resulting relaxation time is not properly defined from the linearized BTE. However, since classical limit and BE statistics do converge to the same limit at high temperature, we utilized the equipartition phonon occupation in Eq. 19 and 20 to compute scattering rates in order to compare with molecular dynamics simulations performed at high T.

III Results and Discussions

III.1 Comparison of temperature-dependent phonon frequencies

It is instructive to compare the phonon frequencies and their changes as a function of temperature obtained using the two different approaches. Fig. S2 in the Supplemental Materials (SM) shows comparisons of phonon dispersions, showing that the fitted potential using CSLD gives similar dispersions to those obtained using finite displacements with the original interatomic potential. In Fig. 4, we show the mode-resolved frequency shifts as calculated in ALD and in MD, relative to the frequencies obtained from direct diagonalization. Fig. S6 in the SM displays the phonon density of states for these three cases. As expected, the frequency shifts increase with temperature. With the exception of high frequency MgO modes, the frequency shifts in MD are found to be smaller than the ones in ALD. The difference in frequency shifts between MD and ALD increases at high frequency, likely because of the different phonon distributions. This difference persists at high temperatures, despite an expected decrease in the difference between the two phonon distributions with increasing temperatures. Moreover, we note that the agreement between ALD and MD approaches is better in MgO, where frequency shifts are large ($\sim$1 THz), than in PbTe, where frequency shifts are on average 10 times smaller ($\sim$0.1 THz). To further confirm the ALD results, we performed additional self-consistent phonon (SCPH) calculations of frequency shifts in reciprocal space, as derived from many-body Green-function theory Mahan (2000); Errea, Rousseau, and Bergara (2011); Tadano and Tsuneyuki (2015); the results are found to agree with those obtained by the real space-based PRN scheme. This indicates that part of the discrepancies between MD and ALD may be numerical in nature.

III.2 Comparison of phonon lifetimes

We compare the phonon mode-resolved scattering rates of MgO and PbTe obtained from MD and ALD at various temperatures in Fig. 5. To perform a valid comparison, we enforce equipartition for phonon population in ALD calculations, thus following the statistics in MD simulations. Overall, good agreement is found between MD and ALD results, despite the rather different formalisms based on atomic motion or phonon quasiparticles. In general, acoustic modes show much smaller scattering rates than those of optical modes, primarily due to the limited scattering phase space of low-frequency phonon modes. Scattering rates of both acoustic and optical modes increase with enhanced phonon populations at higher temperatures. We notice that, for both MgO and PbTe, scattering rates of acoustic modes obtained from the two methods agree well with each other, while those of optical modes are found to be smaller from MD simulations, but with decreased discrepancy between MD and ALD at high temperatures. We note that the larger discrepancy associated with the high-lying optical modes at relatively low temperatures might be due to the break down of relaxation time approximation when equipartition is assumed for phonon population in the linearized BTE, while high temperature tends to reduce such discrepancy, as pointed out by by Feng and Ruan Feng and Ruan (2016). To shed light on the impact of statistics on lifetimes in ALD calculations, we also compare the scattering rates calculated using phonon populations obeying Bose-Einstein statistics and equipartition, respectively, in Fig. S5 (see Supplementary Materials). Consistently, we find that equipartition leads to higher phonon scattering rates, and again, with decreased difference from those obtained by assuming Bose-Einstein statistics at higher temperatures. This discrepancy is deeply rooted in the fact that phonon populations from Bose-Einstein statistics and equipartition are different, particularly for the high-frequency optical phonons. It is the overall increased phonon populations in the case of equipartition that leads to higher phonon scattering rates. Our results suggest that phonon scattering rates calculated using Bose-Einstein statistics compare better with MD simulations that those from equipartition, particularly for the high-lying optical phonon modes. This agreement, however, may be accidental or suggest error cancellations. depend on the populations $n$ of various phonon modes. The effect of these differences is rather difficult to estimate, even limiting ourselves to three phonon processes in the SMRTA, due for example to the form of equation 19, involving both a sum and a difference between different phonon populations.

To reveal the role of the higher-than-third-order phonon scattering processes in determining overall phonon scattering rates, we show the decomposed total scattering rates by separating them into contributions from three- and four-phonon processes for MgO and PbTe in Fig. 6. The MgO results reveal that the contribution of four-phonon scattering is comparable to that of three-phonon scattering at all temperatures studied here. Particularly, four-phonon scattering rates of optical modes with frequency of about 10 THz are even higher than those of three-phonon processes, highlighting the importance of including four-phonon scattering processes to accurately predict the thermal conductivity. It is also evident from the PbTe results that four-phonon scattering rates have a stronger temperature-dependence than three-phonon scattering rates, with increasing relative contributions to total scattering rates at higher temperatures.

III.3 Thermal Conductivity

Here we discuss our finding, summarized in Table 1, in which we use lifetimes obtained from MD simulations ($\tau_{MD}$), as well as from three- and four-phonon scattering ($\tau_{BE}$, in the BTE expression for thermal conductivity (Eq. 18). We also use the classical equipartition ($C_{EQ}$) and Bose-Einstein ($C_{BE}$) expressions for heat capacity in Eq. 18. All the above results are listed in comparison to the Green-Kubo (GK) thermal conductivity, and labeled approximations 1-7 or A1-A7 in Table 1. For MgO, the thermal conductivity values obtained from BTE using three- and four-phonon scattering processes, and Bose-Einstein statistics (A4 in Table 1), agrees strikingly well (to within 3%) with the Green-Kubo results (A7). For PbTe, however, the agreement is only satisfactory (within 16%). This agreement is better than previous comparisons (e.g. Ref. Turney et al. (2009)), likely due to the inclusion of fourth-order phonon scattering in the BTE-ALD treatment. The agreement, however, may be due to a number of error cancellations, as explained below.

One significant difference between the MD and BTE approaches is the inclusion of phonon scattering processes to all orders in MD, as opposed to the order-by-order expansion in BTE. Because the rate of scattering processes at different orders are summed, i.e. via Matthiessen’s rule, the inclusion of higher order should reduce the lifetimes and therefore the thermal conductivity. Indeed, substituting lifetimes from MD into the BTE expression with the BE expression for heat capacity (A2) results in lower thermal conductivity values compared to using BTE lifetimes (A4), due to the inclusion of higher-order scattering processes in MD. Because of the computation of the lifetimes using energy autocorrelation functions, the dissipative vs non-dissipative effects of umklapp vs normal processes, respectively, are also preserved.

The effects of statistics on the heat capacity and the corresponding effects on thermal conductivity can be made apparent by considering the changes in predicted values as the heat capacity is evaluated using equipartition (A1, A3, A5) or Bose-Einstein statistics (A2, A4, A6), when the treatment of the lifetimes is held constant. Since classical heat capacities are consistently higher, and especially so at low temperatures, than Bose-Einstein heat capacities, the thermal conductivity values are also higher in A1, A3, and A5 than in A2, A4, and A6, respectively. The difference decreases with increasing temperature as expected. Moreover, we find that using BE heat capacities with MD lifetimes (A2) worsens rather than improves the agreement with the BTE-ALD results (A4). In agreement with Ref. Turney, McGaughey, and Amon (2009), BE treatment of heat capacities together with MD lifetimes is ruled out as a possible quantum correction for classical MD simulations of thermal conductivity.

As mentioned above, phonon occupations also directly affect scattering rates calculated from three- and four-phonon scattering processes. This effect is apparent from a comparison between A4 and A6. It is seen that changing the occupation numbers to those corresponding to equipartition (A6) rather than the Bose-Einstein distribution (A4), without changing the heat capacity, changes the thermal conductivity significantly (up to 24%). Note also, that the effect of statistics on scattering rates is significantly larger than the effects of statistics on heat capacities. For example, the differences between A5 and A6 is smaller than those between A4 and A6.

Interestingly, substituting MD lifetimes as well as heat capacities evaluated with equipartition into the BTE expression (A1) results in thermal conductivity values which are significantly different from those obtained from the Green Kubo expression (A7). The difference is likely due to the single-mode approximation used in the BTE approach. This highlights the fact that apart from treatments of lifetimes and heat capacities, the summation according to phonon modes in and of itself introduces a difference between BTE and MD methods for computing thermal conductivity.

IV Conclusion

In conclusion, we have presented a detailed comparison between molecular dynamics (MD) and Boltzmann Transport Equation using anharmonic lattice dynamics (BTE-ALD) approaches for the computation of thermal conductivity. Improvements in BTE-ALD such as the inclusion of fourth-order phonon scattering processes, and treatment of temperature-dependent phonon frequency shifts (phonon renormalization) were included. Issues regarding the proper treatment of statistics in MD simulations, namely classical equipartition rather than Maxwell Boltzmann, were addressed. Thermal conductivity values were found to agree well between BTE-ALD and Green-Kubo (GK) MD, but a detailed analysis showed that such agreements are introduced by cancellations of different errors. By substituting lifetimes derived from energy-energy autocorrelation function from MD simulations into the BTE expression, we determined the effects of higher-than-4th-order phonon processes, the effects of statistics via the lifetimes and heat capacities, and the effects of single-mode relaxation time approximation itself, to the calculation of thermal conductivity. Significantly, we find that the effects of statistics on thermal conductivity is primarily due to effects on scattering rates. We also find that a full substitution of MD lifetimes and heat capacities in the BTE expression fails to reproduce GK results, indicating a significant effect of single-mode relaxation time treatment on thermal conductivity. One significant lesson from this comparison is that for integrated properties such as transport coefficients, apparent agreement may mask fundamental physical differences, thus caution is advised in the interpretation of the results.

Acknowledgment

This work was supported by the Midwest Integrated Center for Computational Materials (MICCoM) as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division (5J-30161-0010A). Use of the Center for Nanoscale Materials, an Office of Science user facility, was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

I Interatomic force constants from Compressive sensing lattice dynamics

The essential ingredients for the calculation of $\tau_{\lambda}^{-1}$ are harmonic and anharmonic interatomic force constants (IFCs). We utilized our recently-developed compressive sensing lattice dynamics (CSLD) Zhou et al. (2014) approach to extract both harmonic and anharmonic IFCs, which are obtained from a convex optimization problem that minimizes a weighted sum of the $\ell_{1}$ norm of the IFCs $\mathbf{\Phi}$ and the root-mean-square fitting error

[TABLE]

where $\mathbf{F}$ is a vector composed of atomic forces calculated by DFT, and $\mathbb{A}$ is a matrix formed by the products the atomic displacements. The parameter $\mu$ is a tuning parameter used to control the relative weights of the force fitting error versus the sparsity of $\mathbf{\Phi}$ components. Since (1) small $\mu$ leads to sparse solution of $\mathbf{\Phi}$ at the expense of the accuracy of the fitted $\mathbf{\Phi}$ (“underfitting”) and (2) large $\mu$ give a dense solution of $\mathbf{\Phi}$ but with poor predictive power (“overfitting”), an optimal $\mu$ is determined by monitoring the predictive relative error for a leave-out subset of the training data not used in fitting. To reduce the parameter space of $\mathbf{\Phi}$, linear constraints from crystal symmetry and translational symmetry are explored and strictly enforced by finding the null space representation of $\mathbf{\Phi}$. More details on numerical issues and the symmetrization of IFCs can be found in Ref. [Zhou et al., 2014; Xia and Chan, ; Nelson et al., 2013a, b; Tadano and Tsuneyuki, 2015]. The dependence of force prediction error on the order of included interactions is given in S1

II Phonon Dispersion and Density of States