In-plane phonon transport in thin films

The in-plane phonon thermal conductivities of argon and silicon thin films are predicted from the Boltzmann transport equation under the relaxation time approximation. We model the thin films using bulk phonon properties obtained from harmonic and anharmonic lattice dynamics calculations. The input required for the lattice dynamics calculations is obtained from interatomic potentials: Lennard-Jones for argon and Stillinger–Weber for silicon. The effect of the boundaries is included by considering only phonons with wavelengths that fit within the film and adjusting the relaxation times to account for mode-dependent, diffuse boundary scattering. Our model does not rely on the isotropic approximation or any fitting parameters. For argon films thicker than 4.3 nm and silicon films thicker than 17.4 nm, the use of bulk phonon properties is found to be appropriate and the predicted reduction in the in-plane thermal conductivity is in good agreement with results obtained from molecular dynamics simulation and experiment. We include the effects of boundary scattering without employing the Matthiessen rule. We find that the Matthiessen rule yields thermal conductivity predictions that are at most 12% lower than our more accurate results. Our results show that the average of the bulk phonon mean free path is an inadequate metric to use when modeling the thermal conductivity reduction in thin films. © 2010 American Institute of Physics. doi:10.1063/1.3296394


I. INTRODUCTION
The thermal transport properties of micro-and nanostructures can be significantly different than those of the corresponding bulk material.][10] This reduction can be problematic in devices where efficient heat removal is key to reliability ͑e.g., lasers͒, but beneficial when a low thermal conductivity is desired, such as in a thermoelectric energy conversion device.
The thermal conductivity reduction in a thin film is caused by phonon scattering at the film boundaries.4][15][16] These models, however, obscure the underlying phonon physics because they rely on fitting parameters and assume isotropic phonon properties.Studies of thin films using molecular dynamics ͑MD͒ simulation, 11,17 while confirming the thermal conductivity reduction, have not been able to elucidate the underlying phonon-level mechanisms.While MD simulations can be used to predict phonon transport properties in bulk materials, [18][19][20][21][22] it is unclear how the available techniques can be applied to thin films.Additionally, using MD simulation to access film thicknesses large enough to recover the bulk thermal conductivity is a daunting computational task.
Here, we present a procedure for predicting the in-plane phonon thermal conductivity of thin films using lattice dynamics ͑LD͒ calculations. 19Our procedure uses the properties of all the phonons in the Brillouin zone and does not rely upon fitting parameters.In Sec.II, we describe the theory behind LD calculations and present our methodology for modeling phonon transport in thin films.We then apply this procedure to Lennard-Jones 23 ͑LJ͒ argon films with thicknesses ranging from 4.3 to 540 nm and Stillinger-Weber 24 ͑SW͒ silicon films with thicknesses ranging from 17.4 to 71 200 nm in Sec.III.The films are oriented with the crossplane direction parallel to the ͓001͔ crystallographic direction.We predict the in-plane phonon thermal conductivity of LJ argon thin films at a temperature of 20 K and SW silicon thin films at temperatures of 400 K, where predictions from classical MD simulations are available for comparison.SW silicon thin films are also modeled at a temperature of 300 K, for which experimental data are available.From the analysis, we find that the average of the bulk phonon mean free path ͑MFP͒ cannot be used in simple models to predict the thermal conductivity reduction in thin films, as suggested by others. 25,26Mode-dependent phonon properties, as used in this work, are required.phonon thermal conductivity ͑i.e., the thermal conductivity in the x-or y-direction͒, k, can be predicted from which is derived by solving the Boltzmann transport equation under the relaxation time approximation and applying the Fourier law. 19,27The phonon specific heat, c ph , x-component of group velocity, v g,x , and relaxation time, , all depend on the phonon mode, denoted by the wave vector, , and dispersion branch, .These properties can be obtained from harmonic and anharmonic LD calculations. 19armonic LD is a method for determining the natural frequencies of vibration ͑i.e., the phonon frequencies͒, , of a crystal. 28In harmonic LD, the atomic interactions are described using the second derivatives of the total system potential energy with respect to the equilibrium atomic positions.The vibrations take the form of traveling waves with coordinates given by where Q 0,H ͑ ͒ is a constant and t is time.In an infinite crystal lattice, the wave vector can take on any value in the Brillouin zone.In our calculations, we consider a finite number of wave vectors corresponding to a finite crystal lattice with periodic boundary conditions.In this case, the total number of wave vectors, N, is equal to the total number of unit cells and the number of dispersion branches is equal to three times the number of atoms in the unit cell.
In anharmonic LD, the third-and fourth-order derivatives of the potential energy are applied as a perturbation to the harmonic phonon modes.The phonons can then be described as traveling waves that decay and/or grow with time, where ⌬ and ⌫ are the shift in the phonon frequency and the phonon-phonon scattering rate that result from anharmonic effects.The LD techniques assume that the atomic displacements from their equilibrium positions are much smaller than the interatomic spacing.These methods thus lose accuracy at high temperatures ͑above half of the Debye temperature͒, where the atomic displacements become large and the fifthand higher-order derivatives cannot be neglected. 19 compute the phonon properties needed to evaluate Eq. ͑1͒ from the frequencies and scattering rates found from the LD calculations.The x-component of group velocity is where A = + ⌬ and x is the x component of the wave vector.The phonon-phonon relaxation time is related to the phonon-phonon scattering rate through The phonon specific heat is where V is the system volume, k B is the Boltzmann constant, and ϵ ប A ͑ ͒ / k B T for the Planck constant divided by 2, ប, and temperature, T. Additional details on the theory and implementation of the LD methods are discussed elsewhere. 19,29,30

B. Implementation
Because a thin film has a finite thickness, LD calculations will only predict the existence of standing waves in the cross-plane ͑z͒ direction.These standing waves do not transport energy ͑i.e., they are non-propagating phonons͒.The scenario where no phonons propagate in the cross-plane direction is nonphysical.For very thick films, the phonon properties should approach those in a bulk crystal, where phonons propagate in all directions.Due to this issue, we cannot use LD calculations to directly predict the phonon properties of a thin film.
To predict the thermal conductivity of a thin film, we first use LD calculations to predict bulk phonon properties.We use lattice constants, a, of 5.32 Å for argon 21 and 5.43 Å for silicon. 18For both materials, we use 32ϫ 32ϫ 32 conventional unit cells and periodic boundaries in all three Cartesian directions to evaluate the bulk phonon properties.This Brillouin zone resolution is fine enough that increasing it does not appreciably change the bulk thermal conductivity predictions.
The effects of the boundaries are included by ͑i͒ only considering phonons with wavelengths that fit within the film and ͑ii͒ adjusting the relaxation times to account for boundary scattering.The constraint on the wavelength, which is inversely proportional to the wave vector magnitude, is expressed in terms of the cross-plane component of the wave vector, z .This constraint can be written as where N z is the number of unit cells in the cross-plane direction and l z is an integer whose magnitude is less than N z / 2. The thickness of the film, L z , is equal to aN z .Equation ͑7͒ also gives the valid wave vectors in the bulk system, where N z = 32.In the thin films, N z must be a factor or multiple of N z in the bulk calculation.Otherwise, Eq. ͑7͒ gives wave vectors that are not found in the bulk system.We consider films with N z =2 j , where j is a non-negative integer.Using this model, we thus only need to compute the bulk phonon properties ͑c ph , v g,x , and p-p ͒ once for each combination of material and temperature.
To account for boundary scattering, we set where F is a mode-dependent scaling factor given by Here, p ͑ ͒ is the specularity parameter and ␦ ͑ ͒ ͒͒ , with v g,z being the cross-plane ͑z͒ component of the phonon group velocity.We note that ␦ ͑ ͒ is the ratio of the film thicknesses to the z component of the phonon MFP in bulk.The specularity parameter is the probability that a particular phonon mode will reflect off the boundary.It ranges from zero for completely diffuse scattering to unity for completely specular scattering.In the Appendix, we derive Eqs.͑8͒ and ͑9͒ from the Boltzmann transport equation under the relaxation time approximation.Under the additional assumption that the boundary scattering is independent of the phonon-phonon scattering, the Matthiessen rule can be used to write 31,27 1 where M is the effective relaxation time and the boundary scattering relaxation time is given by The Matthiessen rule is commonly used to combine phonon scattering mechanisms, but, to our knowledge, its accuracy has never been directly tested.Rewriting Eqs.͑10͒ and ͑11͒ in the spirit of Eqs.͑8͒ and ͑9͒, we find where In Fig. 2, F and F M are plotted against ␦ for p = 0 and p = 0.5.Both F and F M increase monotonically with p from a minimum at p = 0 to a maximum at p = 1, where F = F M = 1 for all values of ␦.The Matthiessen rule overpredicts the effect of boundary scattering ͑F M Յ F for all p and ␦͒.This overprediction is greatest when p = 0 and monotonically decreases as p increases.Plotted in the inset of Fig. 2 is the ratio ͑F − F M ͒ / F for p = 0. From the plot, we see that F M deviates from F by at most 12%.Thus, using the Matthiessen rule will result in an underprediction of the in-plane thermal conductivity by at most 12%.In our LD model, there is no cost advantage of using one boundary scattering expression over the other.Therefore, we opt to use the more accurate expressions, Eqs.͑8͒ and ͑9͒, to compute the effective relaxation times used in Eq. ͑1͒.
Boundary scattering is commonly modeled using expressions similar to Eqs. ͑8͒ and ͑9͒ ͑Refs.13, 15, 16, and 32-34͒ or Eqs.6][37] In contrast to previous work, however, our expressions for F and b do not use the isotropic approximation.We instead use our ability to compute the properties of all phonons in the Brillouin zone to differentiate between the in-plane and cross-plane motions of each phonon mode.The LD framework does not, however, indicate the value of the specularity parameter, which, in general, depends on the phonon mode and the atomic structure at the boundary.We assume completely diffuse boundaries ͑p =0͒, rather than use p as a fitting factor, as done by others. 32,38,39e base this assumption on the idea that the reconstruction of free argon and silicon surfaces 12 disrupts phonons traveling in the cross-plane direction, leading to a high probability that phonons incident on the boundary will scatter diffusely.Similarly, for thin films bounded by an amorphous material ͑such as a silicon thin film bounded by silica layers͒, we believe that the transition from a crystalline to an amorphous material presents a large disruption to the phonon propagation and will diffusely scatter the majority of incident phonons.

A. Thermal conductivity predictions
First, consider LJ argon at a temperature of 20 K.The LD-based predictions are made using phonon properties in the classical limit 19 so that they may be compared to thermal conductivity predictions made with the Green-Kubo method using classical MD simulation. 40The bulk thermal conductivities predicted by these two methods are given in Table I.The two predictions are in good agreement because the small atomic displacement assumption is valid at this temperature. 19In Fig. 3, the thin film to bulk thermal conductivity ratio is plotted against film thickness.We use the thermal conductivity ratio to minimize differences between the different prediction methods.The thermal conductivity ratios predicted by the LD and Green-Kubo MD methods are in excellent agreement for films thicker than 4.3 nm, suggesting that the assumption of diffuse boundary scattering is valid for LJ argon thin films.The LD-based predictions for films thinner than 4.3 nm diverge from the Green-Kubo MD results.We believe that this disagreement is due to differences between the thin film and bulk phonon DOS.
As evidence for this assertion, the phonon DOS for a 2.1 nm thick LJ argon film at a temperature of 0 K computed for a structure with free surfaces that is relaxed using MD is plotted in Fig. 4. 12 Also plotted is the DOS for a 2.1 nm thick film embedded in bulk ͓i.e., limited to bulk phonons that satisfy Eq. ͑7͔͒.There are notable differences between the DOS of the free and embedded films, particularly between frequencies of 0.5 and 0.8 rad/ps, where the DOS of the free film is larger than the DOS of the embedded film.As the film thickness increases, the DOS of the free and embedded films converge to the bulk DOS.For films thinner than 2.1 the DOS of the free and embedded structures become more distinct.This qualitative comparison between the DOS of the free and embedded films suggests that our LD-based model for in-plane thin film thermal conductivity, which uses the phonon population of the embedded film, cannot be applied to LJ argon films with thicknesses of 2.1 nm ͑eight atomic layers͒ or less.This assessment is confirmed by the predicted thin films thermal conductivities, where, for films smaller than 4.3 nm ͑16 atomic layers͒, the predicted thermal conductivities diverge from the Green-Kubo results.Similarly, from the DOS of free and embedded SW silicon films and the predicted thermal conductivities, we find our LD-based approach to be appropriate when the film thickness is 17.4 nm ͑128 atomic layers͒ or greater.
For the SW silicon thin films, the LD-based results at a temperature of 400 K are computed in the classical limit so that they may be compared to Green-Kubo MD predictions. 17The predicted bulk thermal conductivities are given in Table I.The LD-based bulk thermal conductivity prediction is 52% higher than the corresponding Green-Kubo MD prediction.This difference is due to the lowtemperature approximations inherent in the LD technique.As shown in Fig. 5͑a͒, our LD model underpredicts the SW silicon thin film thermal conductivity ratio found from Green-Kubo MD by 30% on average.The majority of this underprediction is likely caused by our assumption of completely diffuse boundary scattering.Assuming partially specular boundary scattering would increase the thermal conductivity.The underprediction may also be partially due to the low-temperature approximations inherent in the LD methods, leading to an overprediction of the phonon-phonon relaxation times and an underprediction of F.
We compare the available experimental measurements at a temperature of 300 K ͑Refs.13, 14, 16, and 46͒ to quantum LD predictions at 300 K.The large difference between the LD-predicted and experimental bulk thermal conductivities, seen in Table I, is due to the combined effects of the approximations in the LD techniques and the inability of the SW potential to accurately model thermal transport in silicon. 47he four sets of experimental results plotted in Fig. 5͑b͒ represent two p-type 14,46 and one n-type 13 silicon films sandwiched between silica layers and one suspended bridge of pure silicon. 16We use the same bulk thermal conductivity value to compute the thermal conductivity ratio for each film. 45The large scatter in the experimental data makes it difficult to discern an accurate trend.We can say, however, that the thermal conductivity ratio predicted with our LD model is in good agreement with the lower bound of the experimental results, consistent with our assumption of diffuse boundary scattering.

B. Insights from phonon properties
Using our LD results, we compute the average phonon MFP, ⌳, from Comparing Eq. ͑14͒ to Eqs. ͑1͒ and ͑8͒, we see that ⌳ is a weighted average of ͉v g,x ͉ p-p F. The classical effective MFP for bulk LJ argon at a temperature of 20 K is 4.03 nm.For SW silicon at a temperature of 300 K ͑400 K͒, the bulk quantum ͑classical͒ effective MFP is 243 nm ͑125 nm͒.The reductions in the effective MFPs in the thin films are the same as the reductions in the thermal conductivity ͑i.e., ⌳ film / ⌳ bulk = k film / k bulk ͒.The reductions in the thermal conductivity and MFPs for the thin films are driven by the boundary scattering ͓Eq.͑9͔͒, which, under the assumption of diffuse scattering, does a good job of accounting for the altered phonon dynamics.In Fig. 6, the thermal conductivity ratio is plotted versus L z / ⌳ bulk for the LJ argon film and the SW silicon film at a temperature of 300 K. Also plotted in Fig. 6 are the two proposed simple models 25,26 for predicting the reduction in the in-plane thermal conductivity using the bulk MFP.Neither expression can predict the differing rates of convergence to bulk of the LJ argon and SW silicon thin film thermal conductivity ratios.
The frequency-dependent contribution to the thermal conductivity for bulk SW silicon and for the 556 and 34.8 nm thin films at a temperature of 300 K are plotted in Fig. 7.The area under each curve is proportional to the total thermal conductivity.Phonon modes with frequencies below 50 rad/ps dominate thermal transport.The film boundaries have the largest impact on these phonon modes.The phonon modes with higher frequencies, which contribute little to the total thermal conductivity, are less sensitive to the boundaries.These results can be explained by noting that highfrequency phonons generally have small group velocities and phonon-phonon relaxation times.In the limit of zero groupvelocity or phonon-phonon relaxation time, F is unity ͓see Eq. ͑9͔͒ and the boundaries do not affect the phonon proper- LD (Quantum) Experiment [13]  Experiment [14]  Experiment [46]  Experiment [16]   FIG. 5. ͑Color online͒ Silicon thin film thermal conductivity normalized by the bulk thermal conductivity as found from ͑a͒ classical LD calculations and the Green-Kubo MD method at a temperature of 400 K and ͑b͒ quantum LD calculations and experiment at a temperature of 300 K. ties.For thinner films, however, the boundaries affect even the high-frequency phonons, as illustrated by the reduction in the peak at 86 rad/ps in the 34.8 nm film.The information contained in Fig. 7 can be used to devise strategies to alter the in-plane thermal transport properties of silicon thin films.9][50] In bulk, the most effective nanoparticles would be those that scatter phonons with wavelengths that correspond to frequencies around 10 rad/ps.In thin films, however, the most effective nanoparticles would scatter phonons with frequencies closer to 30 rad/ps.
The in-plane thermal conductivity of a thin film can be increased by reducing the probability that phonons will scatter diffusely at its boundary.Careful fabrication and intelligent selection of the materials bonded to the surfaces of the thin film ͑i.e., controlling the substrate material and growth process͒ will increase the thermal conductivity.High quality interfaces encourage specular phonon scattering and the transmission of phonons between the thin film and the adjacent material, particularly for long wavelength phonon modes, which tend to have low frequencies. 10

IV. CONCLUSION
The in-plane thermal conductivity of a thin film is reduced from the corresponding bulk value.The thinner the film, the more severe the reduction.Efforts to explain the mechanisms behind this reduction are complicated by the routine use of fitting parameters and the isotropic approximation when modeling phonon transport.To eliminate these complications, we developed a methodology for predicting the in-plane thermal conductivity of thin films.This method uses bulk phonon properties obtained from LD calculations along with the Boltzmann transport equation under the relaxation time approximation.The effects of the boundaries are included by considering only those phonons with wavelengths that fit within the film and adjusting the relaxation times to account for mode-dependent, diffuse boundary scattering ͓see Eqs.͑8͒ and ͑9͔͒.Using our boundary scattering expression, we showed that using the Matthiessen rule will result in an underprediction of the in-plane thin film thermal conductivity of at most 12%.
Given that bulk phonon properties are being used, we do not expect our LD model to work for very thin films, where the phonon population deviates from that of a bulk system.We find that the use of bulk phonon properties is appropriate for LJ argon films thicker than 4.3 nm and SW silicon films thicker than 17.4 nm.Our LD-based thermal conductivity predictions for argon thin films at a temperature of 20 K are in agreement with predictions made using the Green-Kubo MD method ͑see Fig. 3͒.This agreement supports our assumption that phonon-boundary scattering is diffuse in LJ argon films.
When compared to Green-Kubo MD predictions for SW silicon at a temperature of 400 K, our LD model overpredicts the thermal conductivity reduction by an average of 30% ͑see Fig. 5͒, suggesting the existence of partially specular phonon-boundary scattering.The LD-predicted reductions in the thermal conductivity are in good agreement with the lower bound of the available experimental measurements.This last result, combined with the spread in the experimental data, indicates that the quality of the thin film boundary affects the phonon-boundary scattering, creating boundaries that can be partially specular to completely diffuse.
Upon computing the bulk MFPs in the argon and silicon systems, we found that the reduction in the in-plane thin film thermal conductivity cannot be explained by this parameter alone.Instead, the properties of the individual phonon modes are required.Knowledge of the phonon properties can also be used to gain insight into the mechanisms of phonon transport through bulk crystals and nanostructures and to quantify the importance of individual phonon modes.This insight can potentially be used to devise and tune techniques to adjust the thermal conductivity of nanostructured materials.

TABLE I .
Predicted and/or measured bulk thermal conductivities for argon and silicon.We estimate the uncertainty in all the LD predictions ͑due to the finite Brillouin zone resolution͒ to be less than Ϯ5% ͑Ref.19͒.The uncertainty in the argon MD predictions is less than Ϯ5% ͑Ref.21͒.Based on typical predictions for SW silicon, the uncertainty in the silicon Green-Kubo MD result is estimated to be Ϯ20% ͑Refs.17, 41, and 42͒.
b References 43 and 44.c Reference 45.FIG.3.͑Color online͒ Thin film thermal conductivity as a fraction of the bulk thermal conductivity for LJ argon at a temperature of 20 K predicted from LD calculations and MD simulations using the Green-Kubo method.
FIG.6.͑Color online͒ Thin film thermal conductivity ratio ͑k film / k bulk ͒ versus L z / ⌬ bulk for argon and silicon ͑300 K͒ thin films and estimations derived from theory. .7. ͑Color online͒ Contribution to the phonon thermal conductivity vs frequency for SW silicon in bulk and 556 and 34.8 nm thin films at a temperature of 300 K.The area under each curve is proportional to the total thermal conductivity.