Size effects in molecular dynamics thermal conductivity predictions

We predict the bulk thermal conductivity of Lennard-Jones argon and Stillinger-Weber silicon using the Green-Kubo GK and direct methods in classical molecular dynamics simulations. While system-sizeindependent thermal conductivities can be obtained with less than 1000 atoms for both materials using the GK method, the linear extrapolation procedure Schelling et al., Phys. Rev. B 65, 144306 2002 must be applied to direct method results for multiple system sizes. We find that applying the linear extrapolation procedure in a manner consistent with previous researchers can lead to an underprediction of the GK thermal conductivity e.g., by a factor of 2.5 for Stillinger-Weber silicon at a temperature of 500 K . To understand this discrepancy, we perform lattice dynamics calculations to predict phonon properties and from these, length-dependent thermal conductivities. From these results, we find that the linear extrapolation procedure is only accurate when the minimum system size used in the direct method simulations is comparable to the largest mean-free paths of the phonons that dominate the thermal transport. This condition has not typically been satisfied in previous works. To aid in future studies, we present a simple metric for determining if the system sizes used in direct method simulations are sufficiently large so that the linear extrapolation procedure can accurately predict the bulk thermal conductivity.


I. INTRODUCTION
0,[16][17][18][19][20][21][22] The Green-Kubo method is an equilibrium technique based on the fluctuation-dissipation theorem.Cubic simulation cells with periodic boundary conditions are typically used and size-independent thermal conductivities can be obtained for most materials using fewer than 10 000 atoms.Interpreting the results of the Green-Kubo method can be challenging for high thermal conductivity materials and systems with complex unit cells, such as superlattices. 6,10,11The direct method is a nonequilibrium steady-state technique in which a heat flux is applied to a simulation cell along the direction of interest.Using the imposed heat flux and the resulting steady-state temperature gradient, the Fourier law is applied to calculate the thermal conductivity.For many materials, it is computationally prohibitive to obtain sample size-independent thermal conductivity predictions using the direct method.For such cases, it becomes necessary to make thermal conductivity predictions for several sample lengths and then perform a postprocessing extrapolation procedure.
The extrapolation procedure commonly used with the direct method can be derived using the kinetic theory expression for thermal conductivity and the Matthiessen rule, 5 resulting in a predicted linear dependence of 1 / k on 1 / L, where L is the sample length.Thermal conductivities for multiple sample lengths are plotted as 1 / k versus 1 / L and a linear fit to the data is extrapolated to the 1 / L → 0 ͑i.e., bulk͒ limit.Using classical MD simulation and empirical interatomic potentials, a linear relationship between 1 / k and 1 / L has been observed for bulk silicon, 5 argon, 7 SiGe alloys, 10 and diamond, 18 Lennard-Jones ͑LJ͒ ͑Ref.6͒ and Si/Ge ͑Ref.
10͒ superlattices, and silicon 19,20 and SiC ͑Ref.21͒ nanowires.Studying gallium nitride ͑GaN͒ samples up to 128 nm in length, Zhou et al. 8 found that a linear relationship reasonably captures the 1 / k dependence on 1 / L at a temperature of 300 K.At a temperature of 800 K, however, they found a nonlinear dependence.For carbon nanotubes, Thomas et al. 22 obtained sample size-independent thermal conductivity predictions for tube lengths on the order of 1 m.A nonlinear relationship is observed when their 1 / k data is plotted versus 1 / L.
Studying bulk LJ argon 7 and LJ superlattices, 6 Landry, McGaughey and co-workers found agreement between Green-Kubo and direct method thermal conductivities to within the prediction uncertainties for temperatures between 20 and 80 K ͓the melting temperature of LJ argon is ϳ87 K ͑Ref.23͔͒.For Stillinger-Weber ͑SW͒ silicon, Schelling et al. 5 found agreement between the two methods at a temperature of 1000 K.In Sec.II, however, we will show that a seemingly accurate application of the direct method extrapolation procedure to SW silicon at a temperature of 500 K results in an underprediction of the Green-Kubo thermal conductivity by a factor of 2.5.This result, along with the nonlinear results of Zhou et al. 8 and Thomas et al., 22 call into question the general validity of the linear extrapolation procedure.
In Sec.III, we assess the validity of the linear extrapolation procedure by considering two systems: ͑i͒ LJ argon, where the two methods are in agreement for all temperatures and ͑ii͒ SW silicon, where the two methods are in better agreement at higher temperatures.To do so, sample lengths that generate length-independent thermal conductivities ͑tens of micron for SW silicon͒ must be accessed.To predict the thermal conductivity of such systems using the direct method is computationally prohibitive.To overcome this challenge, we use lattice dynamics calculations to predict phonon properties and thermal conductivities.This ability allows us to resolve nonlinearities in the 1 / k versus 1 / L trend that only become apparent for SW silicon samples longer than a few microns.We use these results to identify when the linear extrapolation procedure will accurately predict the bulk thermal conductivity.
In Sec.IV, we develop an analytical model for the length dependence of thermal conductivity.We then use this model to derive a metric that can be used to assess whether the sample sizes used in a direct method MD simulation are sufficiently large to accurately predict the bulk thermal conductivity using the linear extrapolation procedure.

A. Green-Kubo method
In the Green-Kubo method, the thermal conductivity, k GK , is predicted using the equilibrium fluctuations of the heat current vector, S, via the fluctuation-dissipation theorem.For a cubically isotropic material, the diagonal components of the thermal conductivity tensor can be averaged to give 2 where k B is the Boltzmann constant, V and T are the system volume and absolute temperature, t is the time, and ͗S͑t͒ • S͑0͒͘ is the heat current autocorrelation function ͑HCACF͒.The details of our Green-Kubo methodology can be found in Refs.9 and 11.Two challenges are encountered when applying the Green-Kubo method to predict thermal conductivity.The first challenge is to accurately specify the converged value of the HCACF integral, which is proportional to the thermal conductivity through Eq. ͑1͒.The HCACFs ͑normalized by their initial values͒ and their integrals for SW silicon at temperatures of 500 and 1000 K are shown in Fig. 1.For both cases, the simulation cells contain 6 3 conventional diamond unit cells ͑N = 1728 atoms͒ and the data correspond to the average of ten independent simulations of 1600 ps ͑3 ϫ 10 6 time steps͒.While the HCACF integral has converged after 400 ps for both temperatures, it begins to drift after 800 ps due to noise in the HCACF. 5For SW silicon, we specify the converged value of the HCACF integral by averaging its value between times of 400 and 800 ps ͑the shaded region of the inset in Fig. 1͒.Specifying the converged value of the HCACF integral for LJ argon, where the HCACF has less noise, 7 is less challenging.
To ensure that sufficient data are available to accurately predict the converged value of the HCACF integral, the total simulation time should be many times greater than the largest phonon relaxation times that dominate the thermal conductivity.For this reason, it can be challenging to apply the Green-Kubo method for high thermal conductivity materials, which typically have large relaxation times.For some complex unit-cell materials ͑e.g., superlattices͒, the amount of noise in the HCACF makes it impossible to specify a converged HCACF integral, even after increasing the total simulation time and averaging over independent simulations. 6,10,11][14] The second challenge is to address the effect of the finite simulation cell size.The thermal conductivity will depend on the size of the simulation cell if there are not enough phonon modes to accurately reproduce the phonon-phonon scattering in the associated bulk material ͓i.e., the Brillouin-zone ͑BZ͒ resolution is too coarse͔.This effect can be removed by increasing the simulation cell size until the thermal conductivity reaches a size-independent value.The effect of the simulation cell size on the thermal conductivity of bulk SW silicon at a temperature of 500 K is provided in Table I.The prediction uncertainties are the 95% confidence intervals based on the results of ten independent simulations of 1600 ps.The uncertainty is larger for SW silicon ͑ϳ20%͒ than for LJ argon ͑ϳ5%͒. 9There is no discernible size dependence for simulation cells containing 4 3 to 8 3 conventional unit cells, a finding that is in agreement with the results of Schell- The shaded region in the inset indicates the time range over which the HCACF integral is averaged to predict the thermal conductivity, k GK , which is represented by a dashed line for each temperature.Note that the time scale of the inset is four times that of the body.ing et al. 5 The thermal conductivity predicted using 8 3 conventional unit cells will be used when comparing to the direct method predictions for SW silicon at a temperature of 500 K.To reduce computational load, 6 3 conventional unit cells are used to generate the other Green-Kubo predictions.The Green-Kubo predicted thermal conductivities for SW silicon and LJ argon ͑at a temperature of 40 K͒ are provided in Table II.The SW silicon prediction at a temperature of 1000 K is in agreement with the predictions of Schelling et al. 5 and Goicochea et al. 25

Overview
In the direct method, a known heat flux, q, is applied across a sample of finite length and the resulting temperature gradient, ‫ץ‬T / ‫ץ‬z, is calculated.The thermal conductivity, k DM , is then determined using the Fourier law, The simulation cell configuration used in our direct method simulations is shown in Fig. 2. The simulation cell consists of a sample region of length L bordered by hot and cold reservoirs in the z direction.Periodic boundary conditions are applied in the x and y directions.Fixed boundaries bound the reservoirs to prevent the sublimation of reservoir atoms.Details related to our direct method methodology are available elsewhere. 7,11

Simulation cell size effects
The size dependence of the thermal conductivity predicted in a direct method simulation can arise in two ways.The first way is the same phenomenon that occurs in the Green-Kubo method.For samples that are too small, the BZ resolution is too coarse and too few phonon modes are available to accurately reproduce the phonon-phonon scattering present in the bulk material.Because direct method simulation cells are much larger than those used for Green-Kubo simulations, this effect is negligible.
The second mechanism by which size dependence arises is through phonon scattering at the sample/reservoir boundaries.This effect is significant when the sample length is smaller than the bulk phonon mean-free paths, ⌳ ϱ ͑ , ͒ ͑here, the subscript ϱ refers to bulk and each phonon mode is identified by its wave vector, , and dispersion branch, ͒.If the sample length is smaller than ⌳ ϱ ͑ , ͒, a phonon of mode ͑ , ͒ emitted from the hot reservoir can travel across the sample without scattering ͑i.e., ballistic transport͒.As the sample length is increased, phonons will have a greater likelihood of scattering before they reach the other sample/ reservoir boundary, resulting in a transition to diffusive ͑i.e., bulklike͒ transport.Phonons traveling ballistically contribute less to the thermal conductivity than in bulk as their meanfree paths are reduced to the sample length.Thermal conductivity is therefore a length-dependent property when ballistic effects are present.
To describe this length dependence, we now present a model based the Boltzmann transport equation ͑BTE͒ and the Matthiessen rule.The derivation begins with an expression for the thermal conductivity in the z direction that is obtained by solving the BTE under the relaxation-time approximation and using the Fourier law 26,27 The summation is over all phonon modes, the volumetric specific heat of each mode, c ph , is k B / V for the classical systems considered here, and v g,z ͑ , ͒ is the z component of the group-velocity vector, v g ͑ , ͒.The phonon transport is described using a set of mode-and system-size-dependent relaxation times, ͑ , , L͒, defined as the average time be- tween successive scattering events, i.e., ͑ , , L͒ = ⌳͑ , , L͒ / ͉v g ͑ , ͉͒.
We apply the Matthiesen rule, which assumes that different scattering mechanisms are independent, to model the length dependence of ͑ , , L͒ such that where ϱ ͑ , ͒ and b ͑ , , L͒ are the intrinsic scattering and boundary scattering relaxation times.We assume that the crystal contains no defects, no free electrons, and no internal interfaces so that ϱ ͑ , ͒ is equal to the relaxation time associated with phonon-phonon scattering.The boundary scattering relaxation time is taken to be the average time between boundary scattering events in the absence of intrinsic scattering, i.e., b ͑,,L͒ = L/2 ͉v g,z ͑,͉͒ .͑5͒ Substituting Eq. ͑4͒ into Eq.͑3͒ and applying Eq. ͑5͒ leads to an expression that describes the length dependence of thermal conductivity,

͑6͒
As L approaches infinity, the bracketed term in Eq. ͑6͒ approaches unity and k z ͑L͒ approaches the bulk value, k ϱ .

Motivating the linear extrapolation procedure
Motivated by Eq. ͑6͒, we can describe the length dependence of thermal conductivity using 1 / k as a function of 1 / L, i.e.,

͑7͒
Here, X is an unknown function of 1 / L that converges to 1 / k ϱ as 1 / L → 0. The bulk thermal conductivity can be estimated using X and its derivatives with respect to 1 / L ͑XЈ, XЉ, etc.͒ at a finite sample length and extrapolating to X͑0͒ using a Taylor-series expansion, i.e., To obtain the derivatives of X, one could use the direct method to predict the thermal conductivities for a range of sample sizes and take the derivatives of a polynomial fit to 1 / k ͑i.e., X͒ versus 1 / L. It is difficult, however, to accurately predict the second-order and higher-order derivatives with this approach due to uncertainty in the direct method predictions.Additionally, to predict an Nth order derivative requires at least N + 1 sample sizes to be considered, which can be computationally expensive.To avoid these challenges, 1 / k ϱ can be approximated by truncating Eq. ͑8͒ after the first-order term.Under this first-order approximation, the procedure for estimating the bulk thermal conductivity reduces to plotting 1 / k versus 1 / L for a range of sample lengths and extrapolating a linear fit to the data to 1 / L =0.This procedure is equivalent to that proposed by Schelling et al., 5 who derived it by assuming that all the phonon modes in Eq. ͑6͒ have the same, averaged properties.We will demonstrate the conditions under which the linear extrapolation procedure ͑i.e., the first-order approximation͒ is valid in Sec.III B.

Direct method thermal conductivity predictions
Our direct method thermal conductivity predictions are shown in Figs.3͑a͒ and 3͑b͒.The linearly extrapolated bulk values, k DM e , are provided in Table II.The discrete data points in Fig. 3͑a͒ represent predictions for SW silicon samples with lengths ranging from 40 to 80 nm.In Fig. 3͑b͒ e e e FIG. 3. ͑Color online͒ Length-dependent thermal conductivities for ͑a͒ SW silicon at temperatures of 500 and 1000 K and ͑b͒ LJ argon at a temperature of 40 K predicted from the direct method and MD simulations.The dashed lines are linear fits to the discrete data.
For SW silicon, the bulk thermal conductivities extrapolated from the direct method data underpredict the values predicted using the Green-Kubo method at both temperatures ͑see Table II͒.The underprediction is more severe at a temperature of 500 K ͑k DM e / k GK = 0.40͒ than at a temperature of 1000 K ͑k DM e / k GK = 0.66͒.Schelling et al. 5 report 1000 K direct method results using larger sample lengths ͑up to 209 nm͒ that are in better agreement with their Green-Kubo result ͑k DM e / k GK = 0.95͒.For LJ argon, the extrapolated bulk thermal conductivity from the direct method data is within the uncertainty of the Green-Kubo prediction.As reported by Turney et al., 7 the predictions of these two methods are within the prediction uncertainties for temperatures ranging from 20 to 80 K. To understand why the Green-Kubo and direct method results are in better agreement for LJ argon than for SW silicon and why the effect is temperature dependent for SW silicon, we next present a carrier-level analysis that uses Eq. ͑6͒ and phonon properties predicted from lattice dynamics calculations.

A. Phonon properties
The advantage of using Eq.͑6͒ to predict thermal conductivity is that length scales that are computationally prohibitive with direct method MD simulations can be studied.To predict the phonon properties required to evaluate Eq. ͑6͒, we use harmonic and anharmonic lattice dynamics calculations.Harmonic lattice dynamics calculations use the second-order derivatives of the interatomic potential to calculate the phonon frequencies ͓͑ , ͔͒.In an anharmonic lattice dynamics calculation, the harmonic solution is perturbed by the third-order and fourth-order derivatives of the interatomic potential, providing the anharmonic shifts in the phonon frequencies and the phonon-phonon relaxation times.The phonon group velocities are found from ‫ץ‬ / ‫.ץ‬The phonon properties predicted from lattice dynamics calculations can depend on the resolution of the BZ.To eliminate this size effect, the BZ resolution is increased until sizeindependent phonon properties and bulk thermal conductivity are attained.To meet these criteria, we use the conventional unit cell and a grid of 12ϫ 12ϫ 12 wave vectors for both LJ argon and SW silicon.A description of lattice dynamics calculations can be found elsewhere. 7,28,29he bulk thermal conductivities calculated using the lattice dynamics-predicted phonon properties are presented in Table II as k LD .Because the lattice dynamics techniques do not include the full anharmonicity of the atomic interactions ͑i.e., they neglect the fifth-order and higher-order derivatives of the interatomic potential͒, k LD is larger than k GK .We further discuss this effect in Sec.III B.

B. Assessing the linear extrapolation procedure
We predict a length-dependent thermal conductivity, k LD ͑L͒, using Eq.͑6͒ and the phonon properties predicted from lattice dynamics calculations.The inverse of the normalized thermal conductivity, k LD / k LD ͑L͒, is plotted as a solid line in Fig. 4͑a͒ ͑SW silicon at T = 500 K͒ and Fig. 4͑b͒ ͑LJ argon at T =40 K͒ as a function of the inverse of the sample length, 1 / L. The discrete data points in Figs.4͑a͒ and  4͑b͒ correspond to the sample lengths used in the direct method calculations shown in Figs.3͑a͒ and 3͑b͒   than for SW silicon.To explain this result, we must assess when neglecting the second-order and higher-order terms in Eq. ͑8͒ is valid.Using Eq. ͑6͒, it can be shown that 1 / k and 1 / L are linearly dependent if each phonon property is well approximated by an average value.Since we consider classical systems, the phonon specific heat is equal to k B / V for all modes.To check if this condition is also true for the group velocities and relaxation times, we plot the bulk thermal conductivity contribution dependence on the phonon mean-free path ͓⌳͑ , ͒ = ͉v g ͑ , ͉͒ ϱ ͑ , ͔͒ in Fig. 5.For LJ argon, thermal conductivity is dominated by phonons with meanfree paths that span 0.5-3 nm.This small spread suggests that using an average mean-free path to describe the phonon transport is reasonable.The near linear dependence of k LD / k LD ͑L͒ on 1 / L in Fig. 4͑b͒ and the accuracy of the bulk thermal conductivity predicted using the linear extrapolation procedure are therefore not surprising.
For SW silicon, phonons with mean-free paths that span several orders of magnitude ͑10 2 -10 4 nm͒ contribute significantly to the bulk thermal conductivity, thus generating the strong deviation from linearity found when k LD / k LD ͑L͒ is plotted versus 1 / L in Fig. 4͑a͒.Using a first-order Taylorseries expansion ͑i.e., the linear extrapolation procedure͒ to model such a nonlinear trend is only accurate around the point where the Taylor series is evaluated.Accurately predicting a bulk thermal conductivity thus requires consideration of sample lengths that correspond to the largest bulk mean-free paths that dominate the thermal transport.The largest bulk mean-free paths for SW silicon are two orders of magnitude larger than the sample lengths we considered.This difference makes the linear extrapolation procedure underpredict k LD by a factor of 2.1 at a temperature of 500 K.When sample lengths comparable to the largest bulk meanfree paths are considered, the error associated with the firstorder Taylor-series approximation decreases.This trend is seen in Fig. 6, where thermal conductivities predicted using sample lengths between 4 and 8 m are plotted and the linear extrapolation procedure ͑dashed line͒ predicts a bulk thermal conductivity that is within 10% of k LD .
For SW silicon at a temperature of 1000 K ͑not shown͒, a nonlinear trend is also found between k LD / k LD ͑L͒ and 1 / L. As temperature is increased, the phonon-phonon-scattering rates increase due to increased anharmonicity, resulting in a reduction in the mean-free paths.Thus, the difference between the mean-free paths and sample lengths also decreases, resulting in more accurate predictions of the linear extrapolation procedure at higher temperatures ͑k LD e / k LD = 0.48 at T = 500 K while k LD e / k LD = 0.57 at T = 1000 K͒.This increase in accuracy is not as significant as that seen in the MD results ͑k DM e / k GK = 0.40 at T = 500 K while k DM e / k GK = 0.66 at T = 1000 K͒.As temperature increases, the fifth-order and higher-order anharmonic terms neglected in the lattice dynamics calculations become more important.Since MD simulations include these higher-order terms, the difference between the mean-free paths and sample lengths decreases at a faster rate than in a lattice dynamics calculation as temperature increases.

C. Discussion
Previous applications of the direct method assumed that 1 / k and 1 / L are linearly dependent.If this assumption is applied to the discrete data points in Fig. 4͑a͒, which correspond to lengths that span less than one order of magnitude ͑40-80 nm͒, one would conclude that the linear fit reasonably describes the 1 / k dependence on 1 / L and that a linear extrapolation procedure can be used to predict the bulk thermal conductivity.By analyzing results for sample lengths that span multiple orders of magnitude ͑10-10 7 nm, see Fig. 6͒, however, it becomes apparent that the observed linear relationship is a construct of using too small a range of sample lengths ͑i.e., any nonlinear function will appear linear if a small enough region is analyzed͒.0,[16][17][18][19][20][21][22] For example, using samples lengths ranging from 29 to 128 nm, Zhou et al. 8 found that a linear relationship reasonably captures the 1 / k dependence on 1 / L for SW GaN at T = 300 K.But at a temperature of 800 K, where the phonon mean-free paths are reduced, they observed a nonlinear trend similar to that shown in Fig. 4͑a͒.To explain this observation, they suggest that their direct method data at T = 800 K correspond to a nonlinear-transport regime ͑i.e., a breakdown of the Fourier law͒.Based on the results presented in Sec.III B for SW silicon, however, we believe that the nonlinear trend observed at T = 800 K is the true 1 / k dependence on 1 / L, and that the linear trend found at 300 K is a result of using too small a range of sample lengths.Zhou et al. state that their observed nonlinear behavior calls into question the linear extrapolation procedure, a statement that is supported by this work.
From the results presented here, it is clear that caution should be used when attempting to predict bulk thermal conductivities using the linear extrapolation procedure.In Sec.III B, we showed that knowledge of the mode-dependent phonon mean-free paths is required to predict what sample lengths need to be considered for the linear extrapolation procedure to be appropriate.In practice, however, the linear extrapolation procedure will be applied to results from direct method MD simulations, where the phonon mean-free paths are not predicted.A more practical metric for determining when the linear extrapolation procedure is appropriate is thus needed.

IV. LENGTH-DEPENDENT THERMAL CONDUCTIVITY MODEL
We now develop an analytical model for the length dependence of thermal conductivity.We will use this model to develop a metric to check if the sample sizes considered in direct method MD simulations are sufficiently large to accurately predict the bulk thermal conductivity using the linear extrapolation procedure.
We begin by converting the summation over all phonon modes in Eq. ͑3͒ to an integral over the first BZ, resulting in where the relaxation time is length dependent.Assuming that optical phonons do not contribute to thermal conductivity and transforming the integral to spherical coordinates gives ϫ 2 sin ddd, ͑10͒ where BZ ͑ , ͒ is the magnitude of the wave vector at the first BZ boundary along the and directions and the summation is restricted to acoustic modes ͑ac͒.Note that multiplying the magnitude of the velocity vector by cos gives its component in the ͓001͔ direction.We now make the Debye approximation for the phonon dispersion, which assumes a single acoustic branch ͑i.e., no distinction between longitudinal and transverse polariza-tions͒, such that = v ac .Under this approximation for a classical system, Eq. ͑10͒ can be simplified to k͑L͒ = 3 where the integral over wave vector has been changed to an integral over frequency and c ph has been replaced with k B / V.The Debye frequency, D , is defined as 26,27 where ⍀ is the volume of the primitive cell and is equal to a 3 / 4 for face-centered cubic and diamond crystal lattices, where a is the lattice constant.The Matthiessen rule ͑see Sec.II B 2͒ is then used to combine the phonon-phonon and phonon-boundary scatterings so that The phonon-phonon relaxation times are modeled using a relationship derived by Callaway for low frequencies, ϱ = A / 2 .This form is in agreement with the lattice dynamicspredicted phonon-phonon relaxation times for SW silicon at frequencies below 3 THz. 30,31The constant A is calculated in the bulk limit and is Substituting Eqs.͑12͒-͑14͒ into Eq.͑11͒ and integrating over , , and yields

͑15͒
where L is a nondimensional length defined as

͑16͒
In the limit L→ϱ, the bracketed term in Eq. ͑15͒ approaches unity and k͑L͒ approaches k ϱ .The assumptions underlying Eq.͑15͒ are reasonable for single-element crystals with facecentered cubic or diamond lattices ͑e.g., argon and silicon͒.Optical modes contribute less than 5% to the bulk thermal conductivity of SW silicon 30 and do not exist in crystals with a monoatomic primitive cell, such as argon.Using Eq. ͑15͒ as a model for the length dependence of the thermal conductivity and the Taylor-series analysis introduced in Sec.II B 3, we now examine the linear extrapolation procedure.The bulk thermal conductivity estimated from the linear extrapolation procedure applied to a sample of nondimensional length L, k ϱ e ͑L͒, is obtained by truncating the second-order and higher-order terms of the Taylor series in Eq. ͑8͒ such that The thermal conductivity predicted using Eq.͑15͒ ͓k͑L͔͒ and the bulk thermal conductivity estimated using the linear extrapolation procedure ͓k ϱ e ͑L͔͒ are plotted in Fig. 7.Both thermal conductivities converge to the bulk value as the sample length increases.k͑L͒ is within 10% of the bulk value for a nondimensional length of 166 while a nondimensional length of 36 is required to obtain the same level of convergence when the linear extrapolation procedure is applied.
A nondimensional length of 36 corresponds to a dimensional-and material-specific length of where we define L min to be the minimum sample length required to ensure that the bulk thermal conductivity predicted using the linear extrapolation procedure is within 10% of the true value.To evaluate L min , we take v ac to the be the average of the ͓001͔ acoustic phonon group velocities in the → 0 limit, where L and T indicate the longitudinal and transverse dispersion branches.Setting k ϱ equal to k LD and using the group velocities predicted from lattice dynamics calculations for SW silicon, L min is 3.3 m͑1.5 m͒ at a temperature of 500 K ͑1000 K͒.As shown in Fig. 6, this length is consistent with the sample lengths required to predict a bulk thermal conductivity that is within 10% of k LD using the linear extrapolation procedure and lattice dynamics-predicted phonon properties in Eq. ͑6͒.Based on this result, we believe that Eq. ͑18͒ can be used to estimate the minimum sample length required when using the linear extrapolation procedure.Since neither v ac nor k ϱ is available prior to performing direct method MD simulations, Eq. ͑18͒ cannot be applied as is.We therefore substitute v ac,L ͓001͔ = ͑C 11 / ͒ 1/2 and v ac,T

͓001͔
= ͑C 44 / ͒ 1/2 into Eq.͑19͒, where , C 11 , and C 44 are the bulk material density and elastic constants.Substituting this definition for v ac into Eq.͑18͒ and solving for k ϱ results in We define k ϱ max to be the maximum thermal conductivity that can be accurately predicted using the linear extrapolation procedure with a minimum direct method sample length, L. In other words, if the bulk thermal conductivity estimated using the linear extrapolation procedure and direct method MD data ͑k DM e ͒ exceeds k ϱ max , then the prediction is not accurate.The minimum sample length should be increased until k DM e Շ k ϱ max .Examining the minimum sample lengths used in our direct method MD simulations ͑26 nm for LJ argon and 40 nm for SW silicon͒ and using the material properties provided in Table III and k GK .This result is consistent with the findings presented in Sec.III B and those of Schelling et al., 5 who report direct method results that are in better agreement with their Green-Kubo predictions using samples lengths up to 209 nm ͑for SW silicon at 1000 K͒. ͑Color online͒ Length dependence of the thermal conductivity using Eq.͑15͒ and the bulk thermal conductivity estimated using the linear extrapolation procedure ͓Eq.͑17͔͒.

V. RECOMMENDATIONS
As summarized in Sec.III C, care must be taken when applying the linear extrapolation procedure to predict bulk thermal conductivity using the results of direct method MD simulations.If systems sizes smaller than the largest bulk mean-free paths that dominate the thermal conductivity are considered, a linear relationship between 1 / k and 1 / L may be incorrectly inferred and the thermal conductivity can be severely underestimated.Based on the results presented in Secs.II and III and our previous experiences, 6,7,9,10,15 we suggest the following procedure for using MD simulation to predict bulk thermal conductivity.
͑1͒ Use the Green-Kubo method, which is discussed in Sec.II A. The Green-Kubo method is advantageous in that: ͑i͒ it uses equilibrium simulations, ͑ii͒ the full thermal conductivity tensor is predicted, and ͑iii͒ system-size effects can typically be eliminated with less than 10 000 atoms.The disadvantage of the Green-Kubo method is that the converged value of the HCACF can be sometimes difficult, or impossible, to specify.][14] ͑2͒ If the converged value of the HCACF cannot be specified, use the direct method, which is described in Sec.II B. First, estimate k ϱ max from Eq. ͑20͒.Then, perform direct method simulations using systems of increasing length until the thermal conductivity predicted from the linear extrapolation procedure, k DM e Ӎ k ϱ max .͑3͒ If the system lengths accessible in the MD simulations cannot generate k DM e Ӎ k ϱ max , then MD simulation cannot be used to predict the bulk thermal conductivity.In such a case, we recommend using phonon properties obtained from lattice dynamics calculations and Eq.͑6͒ in the L → ϱ limit.The technique we used to predict phonon properties, described in Sec.III, is based on harmonic and third-order anharmonic lattice dynamics calculations.Because higherorder terms in the potential energy are neglected, this approach is only expected to be valid up to about half the Debye temperature. 7An alternative method that includes the full anharmonicity of the atomic interactions is to calculate the phonon properties using the spectral energy density.We have described this approach elsewhere. 35

FIG. 1 .
FIG.1.͑Color online͒ Heat current autocorrelation function ͑body͒ and its integral ͑inset͒ for SW silicon at temperatures of 500 and 1000 K.The HCACFs are normalized by their initial values.The shaded region in the inset indicates the time range over which the HCACF integral is averaged to predict the thermal conductivity, k GK , which is represented by a dashed line for each temperature.Note that the time scale of the inset is four times that of the body.
, the results for LJ argon samples with lengths ranging from 26 to

FIG. 4 .
FIG.4.͑Color online͒ Inverse of the normalized lengthdependent thermal conductivities for ͑a͒ SW silicon ͑T = 500 K͒ and ͑b͒ LJ argon ͑T =40 K͒.The squares correspond to the sample lengths used in the direct method simulations ͓see Figs.3͑a͒ and 3͑b͔͒ but are calculated using phonon properties obtained from lattice dynamics calculations.Note the difference in the scales of the vertical axes.

FIG. 5 .FIG. 6 .
FIG.5.͑Color online͒ Bulk thermal conductivity contribution dependence on mean-free path for SW silicon ͑T = 500 K͒ and LJ argon ͑T =40 K͒.The mean-free paths for each mode are sorted using a histogram with a bin width of 2 nm for SW silicon and 0.1 nm for LJ argon.
FIG.7.͑Color online͒ Length dependence of the thermal conductivity using Eq.͑15͒ and the bulk thermal conductivity estimated using the linear extrapolation procedure ͓Eq.͑17͔͒.

TABLE II .
Bulk thermal conductivities, in W/m K, for SW silicon and LJ argon found using the Green-Kubo method ͑k GK , Sec.II A͒ and the direct method ͑k DM e , Sec.II B͒ in MD simulations and from lattice dynamics calculations ͑k LD and k LD e , Sec.III͒.The direct method uncertainty is estimated to be Ϯ20% for SW silicon and Ϯ10% for LJ Argon based on the prediction repeatability.The superscript e indicates that the value was predicted using the linear extrapolation procedure.
GK inTable II is therefore not surprising.For SW silicon, using a minimum sample length of 40 nm results in k DM e Ͼ k ϱ max for both temperatures tested ͑see Table II͒, indicating that larger sample sizes are required to achieve agreement between k DM e

TABLE III .
Lattice constants, densities, and elastic constants for LJ argon and SW silicon at a temperature of 0 K.