Predicting phonon dispersion relations and lifetimes from the spectral energy density

We derive and validate a technique for predicting phonon dispersion relations and lifetimes from the atomic velocities in a crystal using the spectral energy density. This procedure, applied here to carbon nanotubes, incorporates the full anharmonicity of the atomic interactions into the lifetime and frequency predictions. It can also account for nonperiodic interactions between phonons and nonbonded molecules near the solid surface. We validate the technique using phonon properties obtained from anharmonic lattice dynamics calculations and thermal conductivities obtained from nonequilibrium molecular-dynamics simulation.

Predicting and describing the behavior of semiconductors and dielectric materials requires knowledge of their phonon dispersion relations and lifetimes.At low temperatures, the phonon properties of a bulk crystal can be evaluated using harmonic and third-order anharmonic lattice dynamics calculations. 1,2In carbon nanotubes ͑CNTs͒, silicon nanowires, and other nanoscale materials in realistic environments, however, the phonons will scatter via interactions with nonbonded molecules at surfaces. 3,4Since these interactions are nonperiodic in time and space, their effects on the phonon properties cannot be directly incorporated into a lattice dynamics calculation.Challenges in predicting the phonon properties also arise with increasing temperature, as the magnitude of the atomic displacements grows and the influence of higher-order phonon scattering events becomes increasingly important. 5Since third-order anharmonic lattice dynamics calculations only account for three-phonon scattering events, the phonon lifetimes they predict at elevated temperatures are overestimated. 1 Herein, we derive and validate a straightforward technique for predicting phonon dispersion relations and lifetimes directly from the velocities of the atoms in a crystal using the phonon spectral energy density.This technique, which extends formulations developed by others, [6][7][8][9] naturally includes the full temperature-dependent anharmonicity of the atomic interactions and can be applied to both bulk crystals and crystals interacting with nonbonded molecules.We apply this technique to classical molecular-dynamics ͑MD͒ simulations of empty and water-filled ͑8,8͒ CNTs and use the predicted phonon properties to calculate the mode-by-mode contributions to the total thermal conductivity.Next, we compare the fully anharmonic phonon properties to those obtained from third-order anharmonic lattice dynamics calculations.Our findings indicate that, even at room temperature, neglecting higher-order scattering events in a CNT leads to an overestimation of the acoustic phonon lifetimes and the CNT length required to obtain fully diffusive phonon transport.We then use the spectral energy density to identify how interactions with confined water molecules shift the phonon frequencies, lower the phonon lifetimes, and reduce the CNT thermal conductivity.
The empty and water-filled CNT thermal conductivities predicted using the phonon properties obtained from the spectral energy density are in excellent agreement with what we predicted using a direct application of the Fourier law in a MD simulation. 3The number of atoms and simulation runtime required to predict the spectral energy density, however, are both at least one order of magnitude smaller than required by MD-based techniques that predict only the thermal conductivity.Thus, although applied here to classical MD simulation, the techniques outlined in this work can be used to predict fully anharmonic phonon properties using atomic velocities obtained from ab initio MD simulations driven by density-functional theory calculations.
The spectral energy density is derived by projecting the positions of the atoms in a crystal onto the normal modes of vibration, q͑ ; t͒, where is the wave vector, denotes the polarization branch, and t is time.The contribution to the normal-mode amplitude stemming from u ␣ ͑ n x,y,z b ; t͒, the displacement in the ␣ direction of atom b ͑with mass m b ͒ inside unit cell n x,y,z , is 10

͑1͒
where e is the corresponding component of the polarization vector, r͑ n x,y,z 0 ͒ is the equilibrium position of each unit cell, and N T ͑=N x N y N z ͒ is the total number of unit cells.The total normal-mode amplitude is qͩ and the average kinetic energy of normal mode ͑ ͒ is

͑2͒
where q ‫ء˙‬ is the complex conjugate of q ˙.The kinetic energy can be transformed from the time domain to the frequency domain by Parseval's theorem 11 to give

͑3͒
where is angular frequency.Substituting the time derivative of Eq. ͑1͒ into Eq.͑3͒ and averaging over N T gives the spectral energy density,

͑4͒
The polarization vectors do not appear in Eq. ͑4͒ because they are orthonormal. 10The spectral energy density is the average kinetic energy per unit cell as a function of wave vector and frequency.Although similar expressions for calculating phonon properties have appeared in other reports, 8,9 Eq. ͑4͒ and its derivation are formally motivated by lattice dynamics.For the one-dimensional CNT systems considered here, as shown in Fig. 1, the Brillouin zone is one dimensional.The nondegenerate wave vectors are thus limited to z =2n z / a z N z , where a z is the lattice constant, N z is the total number of unit cells in the z direction, and n z is an integer ranging from −N z + 1 to N z .
In Fig. 2͑a͒ we present the spectral energy density for an empty 12.3-nm-long ͑N z = 50͒ 1.10-nm-diameter CNT.The resolution of the frequency axis is 0.001 THz and the resolution of the wave-number axis, which is based on the number of unit cells, is 0.02 z / ͑2 / a z ͒.The atomic velocities were obtained from NVE ͑constant mass, volume, and en-ergy͒ MD simulation at a temperature of 298 K. Interactions between carbon atoms are modeled using the reactive empirical bond order ͑REBO͒ potential 12 and periodic boundary conditions are imposed in the axial direction.We integrate the equations of motion using the velocity Verlet scheme with a 0.5 fs time step and allow the system to equilibrate for 1 ns prior to data collection.The shading on the plot indicates the magnitude of the spectral energy density for each ͑ z , ͒ combination corresponding to a total integration time 0 of 1 ns.We find that integrating beyond 1 ns, which is approximately five times greater than the longest phonon lifetime, does not change the predicted phonon properties.
In Fig. 3͑a͒, we present the spectral energy density along z = 0 for frequencies below 5 THz.Similar peak-and-valley profiles exist along the entire frequency axis at each allowed wave vector.The locations of the peaks in the spectral energy density are in excellent agreement with the phonon frequencies we obtain from harmonic lattice dynamics calculations, indicating that phonon frequency shifts due to anharmonic effects are negligible.For each phonon mode, the range of frequencies accessed by the carbon atoms is related to the anharmonicity of the interatomic potential and the corresponding rate of multiphonon scattering processes.The shape of this frequency spread for each mode is the Lorentzian function where I is the peak magnitude, c is the frequency at the peak center, and ␥ is the half-width at half-maximum.Equation ͑5͒ can be derived by substituting the time-dependent anharmonic normal-mode coordinates 1 into Eq.͑2͒. 13The phonon lifetime is then defined by =1/ 2␥.We superimpose fits of Eq. ͑5͒ onto the spectral energy density shown in Fig. 3͑a͒ and estimate the uncertainty of each fit to be 10% of the phonon lifetime.
In Fig. 2͑b͒, we present the phonon dispersion curves, which are formed by connecting the peaks in the spectral energy density with continuous and differentiable curves along z .For many-atom unit cells, like the 32-atom CNT unit cell used here, guidance from lattice dynamics calculations is needed to determine how the peaks in the spectral energy density should be joined and to identify degenerate branches.The axial group velocity is obtained from v g,z = ‫ץ‬ / ‫ץ‬ z , which we evaluate numerically using a central difference technique.
The relationship between phonon properties and the axial CNT thermal conductivity k z comes by combining the Boltzmann transport equation under the relaxation-time approximation with the Fourier law to give k z = ͚ ͚ c ph v g,z 2 , where c ph is the phonon specific heat. 1 For the classical system considered here, c ph = k B / V, where k B is the Boltzmann constant and V͑=DbN z a z ͒ is the system volume.Using phonon properties obtained from the spectral energy density, we predict a total CNT thermal conductivity of 393 W / m K.As illustrated in Fig. 2͑b͒, we find that the four acoustic branches are responsible for 173 W / m K ͑44%͒ of the thermal conductivity, with the majority of this contribution coming from acoustic modes with reduced wave vectors less than 0.2 z / ͑2 / a z ͒.The midfrequency optical branches, defined here as the 55 optical branches with ⌫-point ͑ z =0͒ frequencies less than 30 THz, are responsible for 196 W / m K ͑50%͒.The remaining 37 higher-frequency optical branches provide the remaining 24 W / m K ͑6%͒.The nonuniform distribution of the thermal conductivity contributions is primarily related to the disparate phonon mean-free paths ͑=v g,z ͒, which extend up to 780 nm for the acoustic modes, to 395 nm for the midfrequency optical modes, but only to 9 nm for the higher-frequency optical modes.We find that doubling the number of unit cells in the CNT fragment ͑i.e., doubling its length͒ does not change the predicted phonon lifetimes or the thermal conductivity.Since the spectral energy density is obtained directly from the atomic velocities, there is no distinction between normal and umklapp phonon scattering processes.Both effects are fully included in the phonon lifetime calculations.
In an independent investigation, 3 we used a direct application of the Fourier law in a MD simulation to explore the variation in thermal conductivity with length for an equivalent 1.10-nm-diameter CNT also modeled using the REBO potential.We found that the thermal conductivity reached a fully diffusive ͑i.e., length-independent͒ value of 407 W / m K when the CNT length exceeded 700 nm.This thermal conductivity is in excellent agreement with the value of 393 W / m K calculated here using phonon properties obtained from the spectral energy density.The 700 nm sample length required to obtain fully diffusive transport in the direct-method calculations, needed to eliminate changes in the phonon lifetimes related to boundary scattering, is longer than 98% of the phonon mean-free paths predicted from the spectral energy density.The additional contributions to thermal transport from phonons with longer mean-free paths, present in the direct-method MD simulations, are comparable to the uncertainty in our reported thermal conductivity prediction. 3We note that the thermal conductivities of graphitic materials predicted using the REBO potential are lower than those measured from experiments, which range from 2000 to 10 000 W / m K. 14 For modes with frequencies greater than 2 THz, the fully anharmonic lifetimes calculated using the spectral energy density are within 10% of the lifetimes we calculate using third-order anharmonic lattice dynamics calculations ͓both at a temperature of 298 K; see Fig. 3͑a͔͒.We find that increasing the length of the CNT fragment beyond 50 unit cells ͑the MD system size͒ does not change the magnitude of the thirdorder anharmonic lifetimes for frequencies greater than 2 THz.For frequencies less than 2 THz, however, the lifetimes we calculate from third-order anharmonic lattice dynamics calculations exceed the fully anharmonic lifetimes predicted from the spectral energy density by up to two orders of magnitude.Furthermore, even for CNT fragments 2000 unit cells long, some of the low-frequency third-order anharmonic lifetimes are not fully converged.Using lifetimes extrapolated from the third-order anharmonic lattice dynamics calculations in place of the fully anharmonic MD values, we find that the maximum phonon mean-free path increases from 780 nm to over 10 m, and that the CNT thermal conductivity increases by a factor of 2. Similarly long mean-free paths were reported by Mingo and Broido, 15 who used thirdorder anharmonic lattice dynamics calculations to explore the transition to fully diffusive phonon transport with increasing length in smaller-diameter CNTs.A comparison of these two third-order anharmonic lattice dynamics results to the fully anharmonic phonon properties suggests that neglecting higher-order phonon scattering leads to an overestimation of the CNT thermal conductivity and the length required to observe fully diffusive phonon transport.
The spectral energy density can also be used to study phonon transport in systems that cannot be investigated using lattice dynamics calculations.For example, in Fig. 3͑b͒ we present the spectral energy density along z = 0 for frequencies less than 5 THz obtained from MD simulation of a water-filled 1.10-nm-diameter CNT.Interactions between water molecules are modeled using the TIP4P/2005 ͑Ref.16͒ potential and the water-carbon interactions are modeled using a Lennard-Jones potential. 17Comparing the water-filled CNT spectral energy density to the empty CNT spectral energy density ͓Fig.3͑a͔͒ clearly shows the effects of collision linewidth broadening, whereby interactions with the water molecules smear the sub-10-THz phonon modes over a wider range of frequencies, reducing the phonon lifetimes by a factor of 5-10.For phonon modes with frequencies less than 2 THz, the water molecules also introduce a small frequency shift of 0.05-0.1 THz.We find that for frequencies greater than 10 THz ͑the maximum frequency accessed by water molecules inside the CNT͒, the spectral energy density and the phonon lifetimes are unaffected by interactions with the water molecules. 3sing the group velocity and lifetime data obtained from the complete water-filled CNT spectral energy density, we predict a total thermal conductivity of 311 W / m K.This value is in excellent agreement with the thermal conductivity of 300 W / m K we predicted using a direct application of the Fourier law in a MD simulation. 3Since the spectral energy density naturally captures the effects of nonbonded atoms on phonon transport, mimicking such interactions by damping the atomic velocities is not necessary. 18n addition to providing a more complete description of phonon transport, the computational effort required to predict the thermal conductivity using the spectral energy density is orders of magnitude smaller than that required in other MDbased techniques.For example, using a direct application of the Fourier law to predict the thermal conductivity of a 1.10nm-diameter CNT requires a simulation system cell containing over 100 000 atoms and an equilibration and data collection time of 20 ns. 3 Using the spectral energy density technique, we find that an integration time of 1 ns and a system size of 1600 atoms are sufficient to obtain fully converged and system-size-independent phonon properties.*mcgaughey@cmu.edu

FIG. 1 .FIG. 2 .
FIG.1.͑Color online͒ Carbon nanotube fragment used when evaluating the phonon spectral energy density.The dashed lines identify the 32-atom unit cell used in Eq. ͑4͒.The chirality vector for the CNT is ͑8,8͒ and the zero-temperature nearest-neighbor atomic spacing is 1.42 Å.The van der Waals thickness of the CNT, b, is taken to be 3.4 Å ͑Ref.3͒.

FIG. 3 .
FIG.3.͑Color online͒ ͑a͒ Semilogarithmic plot of the empty CNT spectral energy density along z = 0 for frequencies below 5 THz at T = 298 K.The lifetime of each mode is determined by fitting each peak to Eq. ͑5͒.Sample fits, shown as solid red lines, and the corresponding phonon lifetimes are presented for each peak.The normal-mode frequencies calculated from harmonic lattice dynamics calculations are shown as vertical dashed lines and the phonon lifetimes obtained from third-order anharmonic lattice dynamics calculations are given in parentheses.͑b͒ Water-filled CNT spectral energy density along z =0.