Thermal conductivity accumulation in amorphous silica and amorphous silicon

We predict the properties of the propagating and nonpropagating vibrational modes in amorphous silica (a-SiO2) and amorphous silicon (a-Si) and, from them, thermal conductivity accumulation functions. The calculations are performed using molecular dynamics simulations, lattice dynamics calculations, and theoretical models. For a-SiO2, the propagating modes contribute negligibly to thermal conductivity (6%), in agreement with the thermal conductivity accumulation measured by Regner et al. [Nat. Commun. 4, 1640 (2013)]. For a-Si, propagating modes with mean-free paths up to 1 μm contribute 40% of the total thermal conductivity. The predicted contribution to thermal conductivity from nonpropagating modes and the total thermal conductivity for a-Si are in agreement with the measurements of Regner et al. The accumulation in the measurements, however, takes place over a narrower band of mean-free paths (100 nm–1 μm) than that predicted (10 nm–1 μm).


I. INTRODUCTION
The vibrational modes in disordered solids (e.g., alloys, amorphous materials) can be classified as propagons (propagating and delocalized, i.e., phononlike), diffusons (nonpropagating and delocalized), and locons (nonpropagating and localized) [1,2].Diffusons contribute to thermal conductivity by harmonic coupling with other modes due to the disorder.Locons do not contribute significantly to thermal transport in three-dimensional systems [3].
Experimental measurements, estimates based on experiments, and modeling predictions have demonstrated that propagating modes contribute significantly to the thermal conductivity of amorphous silicon (a-Si) [4][5][6][7][8][9][10] and amorphous silicon nitride [11], but not to that of amorphous silica (a-SiO 2 ) [5,10,[12][13][14][15][16][17][18].Notably, using broadband frequency domain thermoreflectance, Regner et al. measured how the thermal conductivity of a-SiO 2 and a-Si thin films at a temperature of 300 K change with the thermal penetration depth associated with the heating laser, which identifies the depth normal to the surface at which the temperature amplitude is 1/e of its surface amplitude [10].Adopting the convention of Koh and Cahill [19], they interpret the measured thermal conductivity at a given thermal penetration depth to be representative of the phonons with mean-free paths (MFP) less than that value, allowing for the construction of the thermal conductivity accumulation function [20,21].For a-SiO 2 , the thermal conductivity of a 1000-nm-thick film did not vary for thermal penetration depths between 100 and 1000 nm.This result suggests that any propagating modes that contribute to thermal conductivity have MFPs below 100 nm.For a-Si, they find that the thermal conductivities of films with thicknesses of 500 and 2000 nm vary by 40% between thermal penetration depths of 100 and 1000 nm.This result suggests that propagating modes with MFPs in this range contribute significantly to thermal conductivity.
Interpreting the results of Regner et al. requires knowledge of the MFPs of the propagating modes and the contribution to thermal conductivity from the nonpropagating modes.
The objective of this work is to investigate the propagating and nonpropagating contributions to the thermal conductivity of a-SiO 2 and a-Si by predicting the MFPs and thermal conductivity accumulation functions for realistic models and comparing the predictions to experimental measurements.The paper is organized as follows.The theoretical formulation and modeling framework are discussed in Sec.II.The sample preparation for the a-SiO 2 and a-Si bulk models and the computational details are discussed in Sec.III.In Secs.IV A-IV C, harmonic lattice dynamics calculations are performed to predict the vibrational density of states, the plane-wave character of the vibrational modes, and the group velocity of the low-frequency propagating modes (i.e., the sound speed).The vibrational mode lifetimes are predicted using the molecular-dynamics-based normal-mode decomposition (NMD) method in Sec.IV D. Using the sound speeds and lifetimes, the vibrational mode diffusivities (i.e., the product of the square of the group velocity and the lifetime) are calculated and compared with predictions from the AF theory in Sec.IV E. Using this comparison, a cutoff frequency between propagating and nonpropagating modes is specified.The properties of the propagating and nonpropagating modes are then used to predict the total thermal conductivity in Sec.V A. The thermal conductivity accumulation functions are predicted in Sec.V B, where the results are compared with experimental measurements.

II. THEORETICAL FORMULATION OF VIBRATIONAL THERMAL CONDUCTIVITY
We calculate the total vibrational thermal conductivity k vib of an amorphous solid from where k pr is the contribution from propagating modes [38] and k AF is the contribution from nonpropagating modes predicted by the AF theory [4].Mode-level properties obtained from molecular dynamics (MD) simulations and lattice dynamics calculations will be used as inputs.Equation ( 1) has been used in previous studies of amorphous materials, leading to predictions that while k pr is a negligible fraction of k vib for a-SiO 2 (<10%) [12,13,15,18], it is non-negligible for a-Si (20-80%) [4][5][6][7][8][9].The propagating contribution is modeled as [4,6] where V is the system volume, ω cut is the maximum frequency of propagating modes, DOS(ω) is the vibrational density of states (DOS), C(ω) is the mode specific heat, and D pr (ω) is the mode diffusivity.When using mode properties obtained from calculations on finite-sized systems, it is common to write Eq. ( 2) as a summation over the available modes [4,6].We choose the integral form because the required use of finite-sized simulation cells limits the lowest-frequency modes that can be accessed.An extrapolation must be made to the zero-frequency limit that is more easily handled with the integral [4][5][6][7][8]15,18].Equation ( 2) is obtained by using the single-mode relaxation-time approximation to solve the Boltzmann transport equation for a phonon gas [38].
In the derivation of Eq. (2), the system is assumed to be isotropic (valid for an amorphous material) and have a single polarization, making the mode properties only a function of frequency.The choice of a single polarization (i.e., an averaging of the transverse and longitudinal branches) does not significantly change the results predicted in this work or in that of others [4][5][6][7][8]18].We will evaluate Eq. ( 2) under the Debye approximation, which assumes isotropic and linear dispersion such that the DOS is where v s is an appropriate sound speed.The specific heat in the classical, harmonic limit is k B , where k B is the Boltzmann constant.Taking this classical limit allows for a direct comparison between the latticedynamics-based predictions and those from the classical MD simulations.The harmonic approximation has been found to be valid for classical systems ranging from Lennard-Jones (LJ) argon [39] to crystalline Stillinger-Weber silicon [40] at temperatures below half the melting temperature.The full quantum expression for the specific heat is [38] where is the Planck constant divided by 2π .The quantum specific heat will be used for the nonpropagating modes to compare the k AF predictions to experimental measurements in Secs.V A and V B.
The diffusivity of the propagating modes is where τ (ω) is the frequency-dependent mode lifetime and (ω) is the MFP, defined as (ω) = v s τ (ω).The lifetimes will be modeled using where B is a constant coefficient that incorporates the effect of temperature.By using a constant sound speed, the lifetime and diffusivity frequency scalings will be the same.For amorphous materials, the scaling exponent n has been found experimentally and numerically to be between two and four [6][7][8][9][25][26][27][28][41][42][43].A value of two corresponds to anharmonic scattering [44], while a value of four corresponds to Rayleightype scattering from point defects [45].Combined with Eq. (3), choosing n 2 ensures that the thermal conductivity evaluated from Eq. ( 2) is finite.Choosing n > 2 causes the thermal conductivity to diverge in the zero-frequency limit, which can be fixed using additional anharmonic [4,6] or boundary scattering terms [5,7,8].
The AF diffuson contribution to thermal conductivity is [4,6] where ω i is the frequency of the ith diffuson mode, C(ω i ) is the diffuson specific heat, and D AF (ω i ) is the diffuson diffusivity.Equation ( 7) is written as a sum because there are enough high-frequency diffusons in the finite-size systems studied here to ensure a converged value.The AF diffusivities are calculated from [1] where δ is the Dirac delta function [46].The heat current operator S ij , which measures the thermal coupling between vibrational modes i and j based on their frequencies and spatial overlap of eigenvectors, can be calculated from harmonic lattice dynamics theory.For Eq. ( 8), S ij is directionally averaged because amorphous materials are isotropic.

III. CALCULATION DETAILS A. Sample preparation
The three smallest a-SiO 2 samples are the same as those used in Ref. [47] and contain 288, 576, and 972 atoms at a density of 2350 kg/m 3 .The atomic interactions are modeled using the modified Beest-Kramer-van Santen (BKS) potential [48,49] from Ref. [47], except that the 24-6 LJ potential [50] is changed to a 12-6, which has a negligible effect on the predictions.The LJ potentials use a cutoff of 8.5 Å and the Buckingham potentials use a cutoff of 10 Å.The electrostatic interactions are handled using the Wolf direct summation method with a damping parameter of 0.223 Å−1 and a cutoff of 12 Å [51].The simulated density of 2350 kg/m 3 , which corresponds to zero pressure, is 7% larger than the experimental value of 2220 kg/m 3 [52].At the experimental density, the BKS potential generates a structure in tension at a pressure of 2-3 GPa [47].We use a density of 2350 kg/m 3 to be consistent with previous simulations.
Larger systems of 2880, 4608, and 34 562 atoms were created by tiling the smaller samples, melting at a temperature of 10 000 K, and quenching instantaneously to a temperature of 300 K at constant volume.The melt-quench procedure and subsequent MD simulations were performed using the MD package LAMMPS and a time step of 0.905 fs [53].The resulting a-SiO 2 structure is built from a network of rigidly bonded SiO 4 tetrahedral subunits that are weakly bonded via shared oxygen atoms, as shown in Fig. 1(a).The radial distribution function g(r) for the 4608-atom sample is shown in Fig. 1(a) along with an experimental measurement [54].which compares well with our sample.At least 99.5% of the atomic coordinations are at the expected values (4 for Si and 2 for O) [47].
For a-Si, we use samples with 216, 1000, 4096, and 100 000 atoms, generated from the modified Wooten-Winer-Weaire (WWW) algorithm from Ref. [55].The resulting a-Si structure is a rigid, predominantly tetrahedrally bonded network and is shown in Fig. 1(b).A larger sample was created from the 100 000-atom sample by tiling it twice in all directions to create an 800 000-atom sample with a side length of 24.81 nm.All a-Si structures have a density of 2330 kg/m 3 , equivalent to the perfect crystal with a lattice constant of 5.43 Å.The Stillinger-Weber (SW) potential is used to model the atomic interactions [56].The MD simulations are performed using LAMMPS with a time step of 0.5 fs.The radial distribution function for the 4096-atom sample is shown in Fig. 1(b).Also shown in Fig. 1(b) is an experimental measurement [57] which compares well with our current sample, but with a slight broadening in the first peak.
Amorphous materials may have many different atomic configurations with nearly equivalent potential energies, leading to potential metastability during MD simulations [6,9,37].This metastability can cause errors when predicting vibrational lifetimes using NMD (see Sec. IV D).To remove metastability, all a-SiO 2 and a-Si samples were annealed at a temperature of 1100 K for 10 ns [6,9].The removal of metastability is demonstrated by a decrease and plateau of the sample's potential energy during the annealing.

B. Simulation details
Before data collection, all MD simulations are first equilibrated in an NVT (constant number of atoms, volume, and temperature) ensemble for 10 6 time steps at a temperature of 300 K. Data are then collected from simulations in the NVE (constant number of atoms, volume, and total energy) ensemble for 2 21 time steps, where the atomic trajectories are sampled every 2 8 time steps.Ten MD simulations with different initial conditions are run and the predictions are ensemble averaged.
The Green-Kubo (GK) method [58] is used to predict a thermal conductivity k GK without using Eq. ( 1) using the first-avalanche method to specify the converged value of the integral of the heat current autocorrelation function (Sec.V A) [59].For system sizes of 4608 (a-SiO 2 , supercell side length of 4.026 nm) and 4096 (a-Si, supercell side length of 4.344 nm) atoms, the trajectories from the MD simulations are also used to predict the vibrational mode lifetimes using the NMD method (Sec.IV D).
For an amorphous supercell, the only allowed wave vector is the gamma point (i.e., κ κ κ = 0), where κ κ κ is the wave vector and there are 3N a polarization branches labeled by ν, where N a is the number of atoms.Specification of the vibrational modes at the gamma point requires the eigenvalue solution of a dynamical matrix of size (3N a ) 2 that scales as [(3N a ) 2 ] 3 , limiting the system sizes that can be considered to 4608 (a-SiO 2 ) and 4096 (a-Si) atoms.The eigenvalue solution is also required to predict the vibrational DOS (Sec.IV A) and structure factors (Sec.IV B), and to perform the NMD calculations (Sec.IV D) and the AF calculations (Sec.IV E).The frequencies and eigenvectors were computed using harmonic lattice dynamics calculations with GULP [60].The calculation of the AF thermal diffusivities [Eq.( 8)] is performed using GULP and a Lorentzian broadening of 14δω avg for a-SiO 2 and 5δω avg for a-Si, where δω avg is the average mode frequency spacing [δω avg = 1.8 × 10 10 rads/s (a-SiO 2 ) and 1.0 × 10 10 rads/s (a-Si)].Varying the broadening by 10% around these values does not change k AF within its uncertainty (Sec.V A).

A. Density of states
The vibrational DOS is computed from where a unit step function of width 100δω avg is used to broaden δ(ω i − ω).The results for a-SiO 2 and a-Si are plotted in Fig. 2.
The DOS for a-Si is similar to that of crystalline silicon [61], with peaks at mid and high frequencies.The DOS for a-SiO 2 is constant over most of the frequency range, with a gap that The DOS for a-Si has two peaks similar to the DOS of the crystalline phase [62].The DOS for a-SiO 2 is flat over most of the spectrum, with a high-frequency gap that separates the modes involving Si-O interactions [47].
separates the high-frequency Si-O interactions.There is a clear ω −2 scaling for both a-Si and a-SiO 2 at the lowest frequencies.The onset of this scaling occurs at a higher frequency for a-Si (∼1.5 ×10 13 rads/s) than a-SiO 2 (∼4.5 ×10 12 rads/s).This low-frequency scaling is predicted by the Debye model [Eq.( 3)] and suggests that these modes may be propagating (i.e., phononlike).

B. Structure factor
Calculating the structure factors of the disordered modes of the supercell at the gamma point is a method to test for their propagating (i.e., plane-wave) character at a particular wave vector and polarization.This approach has been previously used to predict effective dispersion curves of disordered and amorphous materials experimentally [18,26,28,63] and numerically [41,42,[64][65][66].The structure factor at a wave vector κ κ κ is defined as [ where the summation is over the gamma modes, E T refers to the transverse polarization and is defined as and E L refers to the longitudinal polarization and is defined as In Eqs. ( 11) and ( 12), the b summations are over the atoms in the disordered supercell, r r r 0 ( l = 0 b ) refers to the equilibrium atomic position of atom b, l labels the unit cells (l = 0 for the supercell), α labels the Cartesian coordinates, and κ κ κ is a unit vector.The vibrational mode shape is contained in the 3N a components of its eigenvector e( κ κ κ = 0 0 0 b ν α ).The transverse and longitudinal structure factors are plotted in Figs.3(a) and 3(b) for a-SiO 2 and a-Si for wave vectors along the [100] direction of the supercells.Because amorphous structures are isotropic, the structure factors are direction independent.Mode frequencies ω 0 (κ κ κ) and linewidths (κ κ κ) can be predicted by fitting each structure factor peak to a Lorentzian function of the form where C 0 (κ κ κ) is a constant related to the DOS [65].A dispersion relation is identified by plotting the ω 0 (κ κ κ) values in the middle panels of Figs.The wave vectors are normalized by κ max = 2π/a, where a is 4.8 Å (a-SiO 2 ) and 5.43 Å (a-Si), based on the lattice constants of the crystalline phases [48,56].
κ max = 2π/a and a is the lattice constant of crystalline silicon (5.43 Å) [56].For a-SiO 2 , the coefficients of determination are greater than 0.8 for |κ κ κ|/κ max 0.2 and less than 0.7 for larger wave vectors, where the structure factors peaks are less than an order of magnitude larger than the background.To evaluate κ max for a-SiO 2 , we use a lattice constant of 4.8 Å, which corresponds to the a direction of quartz [48].
For a-Si, the extracted dispersion is nearly linear at small wave vectors with a slight decrease in slope at the largest values [4,6].For a-SiO 2 , the dispersion is concave down for the smallest wave vectors considered, transitioning to a strong concave-up dispersion at intermediate wave vectors.For the intermediate wave vectors, the longitudinal dispersion for a-SiO 2 is well described by the so-called "dispersion law for diffusons," where ω ∝ κ 2 [65].This large concave-up dispersion has been observed in experimental measurements and numerical models of amorphous materials including a-SiO 2 [18,41,42,63,64].We note that at frequencies lower than 10 12 rads/s, experimental measurements of a-SiO 2 recover a linear dispersion [18,26,28,63,67].This frequency range is not accessible with the a-SiO 2 models studied in this work.
The atomic structures of a-SiO 2 and a-Si play an important role in determining the differences in the low-frequency mode properties.The weakly bonded network of tetrahedra in a-SiO 2 [48,49] results in a Debye scaling of the DOS that occurs at a lower frequency than in a-Si (Fig. 2), which is a network of strongly bonded tetrahedra [2,55,56].The lower-frequency onset of the Debye scaling of the DOS for a-SiO 2 leads to the strong nonlinear dispersion seen in Fig. 3(a).The behavior of the DOS and structure factors demonstrate a clear difference in the properties of the low-frequency modes for our models of a-SiO 2 and a-Si, which is further investigated in the following sections.

C. Sound speed
For a disordered solid, except for the transverse and longitudinal sound speeds, there is not an accepted method to predict the group velocity of individual vibrational modes.While the structure factor gives the frequency spectrum needed to construct a propagating state with pure wave vector κ κ κ, the individual mode spectra E T ( κ κ κ ν ) and E L ( κ κ κ ν ) predict the plane-wave character of each mode [2].It is not generally possible to assign a unique wave vector to individual modes, even at low frequency [2], which makes predicting their group velocities challenging.While attempts have been made to predict individual mode group velocities [9,61,68,69], there is no theoretical basis for the proposed methods.
We now use the DOS and structure factors predicted in Secs.IV A and IV B to predict the group velocities of the lowfrequency modes for a-SiO 2 and a-Si (i.e., the sound speeds).By fitting the DOS from Fig. 2 to Eq. (3), a sound speed is obtained and is reported in Table I.Because the DOS is a mixture of transverse and longitudinal modes, only a single sound speed can be predicted.Both longitudinal and transverse sound speeds can be predicted from the structure factor peaks by forward differencing the dispersion relation as where κ min is 0.1κ max for a-SiO 2 and 0.125κ max for a-Si.The results are provided in Table I.The transverse and longitudinal sound speeds can also be predicted from the material's bulk (G) and shear (K) moduli from [60] and Using the bulk and shear moduli defined in terms of the elastic constants according to the Voigt convention [60], the corresponding sound speeds are reported in Table I.
The longitudinal and transverse sound speeds for a-SiO 2 predicted using the moduli are 10%-20% lower than predictions made by Horbach et al. using a linear fit to the peaks of the current correlation function for a model with 8016 atoms using the BKS potential [3568 m/s (transverse) and 5937 m/s (longitudinal)] [41].The smaller values predicted by the structure factors and DOS result from the concavedown dispersion seen at low wave vector (i.e., we are not able to reach the linear portion of the dispersion curve).Experimental measurements of the sound speeds of a-SiO 2 using Brillouin light and inelastic x-ray scattering range between 3800 to 4000 m/s (transverse) and 6000 to 6400 m/s (longitudinal) [63,67,70].Differences between our predictions and experimental measurements may be related to limitations of the BKS potential.
The effect of the concave-down dispersion is less pronounced for a-Si than for a-SiO 2 , where the sound speeds predicted by all three methods are within 5% of each other.Our sound speed predictions for a-Si using all three methods are within 10% of predictions made using the elastic moduli [71,72] and structure factor [42] from models created by the original WWW algorithm [73].The 4096-atom model created by the modified WWW algorithm [55] predicted a longitudinal sound speed of 7670 m/s from the structure factor [43], within 5% of our prediction.In an attempt to explain the anomalously high longitudinal sound speed (8300 m/s) and thermal conductivity measurements in Ref. [7], three 1000-atom a-Si models relaxed using a tight-binding electron structure method predicted an average of 4740 m/s (transverse) and 7830 m/s (longitudinal) [7].By annealing our structures to remove metastability, the sound speeds predicted by the elastic moduli are increased, but not by the amount reported in Ref. [7].Experimental transverse sound speed measurements using Rayleigh wave scattering are 3420 and 4290 m/s for sputtered and ion-bombarded a-Si thin films [35], which is within 15% of the predictions from our models.It is clear that the experimentally measured sound speeds for a-Si show a wide range, which depends on the deposition method and impurity concentration [7,8,35].
The sound speed v s,DOS will be used for both a-SiO 2 and a-Si for the rest of this work, allowing for the use of a single polarization for the propagating contribution [Eq.( 2)].By comparing the sound speeds in Table I, it is clear that the low-frequency DOS of our models for a-Si and a-SiO 2 are dominated by transverse modes, which is expected due to their degeneracy and lower frequencies compared to the longitudinal modes.The transverse sound speed predicted for our model of a-SiO 2 is 85% of that predicted by the other methods (Table I) and that measured by experiment [63,67,70], which is likely related to the larger density of our samples (see Sec. III A).While using a smaller transverse sound speed leads to an underprediction of the mode diffusivities [Eq.( 5)], it leads to an overprediction of the DOS [Eq.( 3)].Holding all other input parameters in Eq. ( 1) constant, a smaller sound speed leads to a larger k pr because the DOS scales as 1/v 3 s .We can thus regard our k pr prediction as an upper bound.

D. Lifetimes
We now predict the lifetimes of all vibrational modes in our models of a-SiO 2 and a-Si using the MD simulation-based NMD method, which explicitly includes the disorder in the supercell [9,66,74].In NMD, the atomic trajectories from an MD simulation are first mapped onto the vibrational mode coordinate time derivatives Here, m b is the mass of the b th atom in the supercell, uα is the α component of the atomic velocity, and t is time.Because the supercells of a-SiO 2 and a-Si are disordered, the NMD method can only be performed at the gamma point (κ κ κ = 0 0 0).The spectral energy of each vibrational mode (ν,ω) is calculated from We choose the frequency-domain representation of the normal mode energy because we find it to be less sensitive to any remaining metastability of the amorphous structure than the time-domain representation.The vibrational mode frequency and lifetime are predicted by fitting each mode's spectral energy to a Lorentzian function where the constant C 0 (ν) is related to the average energy of each mode.This expression is valid when the linewidth (ν) ω 0 (ν) [40].The mode lifetime is [74] The NMD-predicted lifetimes are plotted in Figs.4(a) and 4(b) for a-SiO 2 and a-Si.Also plotted are the time scales extracted from the structure factor linewidths 1/[2 (κ)] (Sec.IV B).For a-SiO 2 , the NMD lifetimes are larger than the Ioffe-Regel (IR) limit τ = 2π/ω [64], and are bounded by this limit at low frequencies.Similarly for a-Si, the IR limit is a lower limit for the lifetimes predicted by NMD.While lifetimes predicted near the IR limit do not satisfy the constraint (ν) ω 0 (ν), only a limited number of these  (20)] and the time scales extracted from the structure factors [Eq.( 13)] (a) a-SiO 2 and (b) a-Si.For a-Si, a clear ω −2 scaling is observed at low frequencies, while the lifetimes plateau at higher frequencies and over a wider range of frequencies than for a-SiO There is no clear evidence for an ω −2 scaling in a-SiO 2 , which would correspond to propagating modes.midfrequencies, the NMD lifetimes are approximately constant and there is peak near 2 ×10 14 rads/s, which corresponds to the peak in the DOS (see Fig. 2).The time scales predicted from the structure factor fall below the NMD-predicted lifetimes and the IR limit by up to one order of magnitude.These low values result the structure factors for a-SiO 2 are evaluated for wave vectors where the resulting wave packets are formed by nonpropagating modes [4,6].
For a-Si, the NMD lifetimes show a clear ω −2 scaling at low frequency.The lifetimes plateau at higher frequencies, over a wider range of frequencies than for a-SiO 2 , with two peaks corresponding to peaks in the DOS (see 2).A similar plateau of lifetimes at high frequencies has reported for lattices [66,75] and in another study of a-Si [9].The transition from the low-frequency scaling to the plateau region occurs near 10 13 rads/s, which to where the DOS first peaks Fig. 2. Similar behavior has been observed for models of disordered lattices [66].The time scales predicted by the structure factors are in good agreement with those predicted by NMD at low frequencies.Similar agreement has been reported in other models of amorphous [6,[76][77][78].The agreement between the NMD-predicted lifetimes and the structure factor time scales for a-Si at low frequencies indicates  5) and ( 20) with the DOS sound speed from Table I] and the AF theory [Eq.( 8)].Also shown are extrapolations based on an ω −2 scaling with Eqs. ( 5) and ( 6) for a-SiO 2 and a-Si, and an additional ω −4 scaling for a-Si.For both systems, the diffusivities are larger than the high-scatter limit [Eq.( 21)] except at high frequencies, where the modes are localized.that these modes are plane-wave-like and that the wave packets formed by these modes are propagating [4,6].
The NMD-predicted lifetimes for a-Si range from 0.5 to 10 ps and are similar in magnitude to those predicted for previous WWW-generated models of a-Si [77][78][79][80].We note that one previous study of a-Si modeled using the Tersoff potential predicted vibrational lifetimes on the order of 100 ps [9], an order of magnitude larger than the values reported here and in previous studies [77][78][79][80].There are several issues to consider when comparing our results to those of He et al. [9].The a-Si models Ref. [9] were prepared using a melt-quench technique that may lead to differences.When we applied the Tersoff potential (as used by He et al.) to the WWW a-Si models, we predict similar lifetimes to those from the SW potential.Furthermore, in Ref. [9] the NMD analysis was performed in the time domain, where the effects of metastability can be more strongly pronounced.Finally, we note that the a-Si bulk thermal conductivity predicted by He et al. using the Green-Kubo method is 40% larger than our prediction (3 W/m-K versus 2.1 W/m-K).

E. Diffusivities
Using the sound speeds predicted from the DOS (Table I), the NMD-predicted lifetimes for a-SiO 2 and a-Si are used to predict the mode diffusivities with Eq. ( 5).The results are plotted in Figs.5(a) and 5(b).The mode diffusivities are predicted from the NMD lifetimes for the low-frequency modes where the DOS scales as ω 2 (Fig. 2).The AF theory is used to predict the mode diffusivities for all frequencies and the results are also plotted in Figs. 5

(a) and 5(b).
For a-SiO 2 , the mode diffusivities predicted by NMD and AF agree well at low frequency.The AF diffusivities at the highest frequencies show a sharp decrease, which is an indication that these modes are localized [4].The low-and mid-frequency diffusivities are above the high-scatter limit which assumes that all vibrational modes travel with the sound speed and scatter over a distance of the lattice constant [14].
In evaluating Eq. ( 21), we use the lattice constant of the crystalline phase (see Sec. IV B).The low-frequency NMD diffusivities do not show a definitive scaling.Based on the results in Ref. [18], we choose a propagating/nonpropagating cutoff frequency of 4.55 × 10 12 rads/s, which is at the onset of the Debye scaling of the DOS (Fig. 2).The constant B in Eq. ( 6) for n = 2 is then fit to the AF-predicted diffusivities for frequencies below the cutoff by dividing the diffusivities by v s,DOS .The fit value is B = 5.65 × 10 13 rads 2 s −1 .
For a-Si, the mode diffusivities predicted by NMD show a clear ω −2 scaling.The NMD-predicted diffusivities are larger and show less scatter than those predicted by the AF theory, which is due to the finite-size and the broadening is required to evaluate (8) [4].By using a larger broadening (100δω avg ), the scatter in the AF-predicted diffusivities at low frequency can be smoothed, but at the cost of decreasing the diffusivities at intermediate and high frequencies, which affects the predicted diffuson contribution to thermal conductivity (see Sec. V A).It is possible that a frequency-dependent broadening may be necessary for a-Si and the AF theory, but determining this dependence is not necessary for interpreting our results.For a-Si, the AF diffusivities are larger than the high-scatter limit [Eq.( 21)], except for the highest-frequency modes, which are localized [4].
For a-Si, we choose ω cut and B so that Eq. ( 5) is equal to the average AF-predicted diffusivity at the cutoff frequency.The resulting values are ω cut = 1.16 × 10 13 rads/s (which is at the onset of the Debye scaling of the DOS, Fig. 2) and B = 2.76 × 10 14 rads 2 s −1 .This choice allows Eq. ( 5) to pass reasonably well through both the AF-and NMD-predicted diffusivities.
While experiments on a-SiO 2 show that there is a crossover region for the low-frequency lifetime scaling from ω −2 to ω −4 [25] and then back to ω −2 [25][26][27][28], our present model is not large enough to investigate the mode properties in this crossover region.Because experiments are limited for a-Si thin films [24], we also consider an ω −4 scaling for Eq.(6).Because this scaling is not clear from the data in Fig. 5(b), we use a cutoff frequency of 1.52 ×10 13 rads/s (which is at the onset of the Debye scaling of the DOS, Fig. 2) based on Refs.[4,5] and choose B = 2.07 × 10 40 rads 4 s −3 so that Eq. ( 5) is equal to the average AF-predicted diffusivity at the cutoff frequency.
Both a-SiO 2 and a-Si have a region at higher frequencies where the AF-predicted mode diffusivities are relatively constant.This behavior has been reported for model disordered systems such as disordered lattices [65,66,75] and jammed systems [81,82].While diffusons are nonpropagating modes whose MFPs are not well defined [4], a diffuson MFP can be calculated from where τ (ω i ) is the NMD-predicted lifetime for that mode.Using this definition, AF (ω i ) for both a-SiO 2 and a-Si is found to vary between the crystal lattice constant (∼0.5 nm) and the supercell size (∼5 nm) for modes with frequency above the cutoff.Similar MFPs have been estimated for diffusons in a-Si in previous studies [4,6].For modes with frequency below the cutoff, the NMD-predicted MFPs range up to 16 nm (a-SiO 2 ) and 43 nm (a-Si).This result is in contrast to the MFPs estimated in Ref. [9] for a-Si, which ranged up to 500 nm.We believe that the origin of the large MFPs in Ref. [9] is a combination of the predicted lifetimes (see Sec. IV D) and the method used to estimate the mode group velocities.

A. Bulk
To predict the bulk thermal conductivity for our models of a-SiO 2 and a-Si, we use both Eq. ( 1) and the GK method.The GK method is computationally inexpensive compared to the NMD and AF methods so that larger system sizes can be accessed.The GK-predicted thermal conductivities for a-SiO 2 and a-Si are plotted in Fig. 6 versus the inverse of the length of the simulation cell.For a-SiO 2 , there is no system-size dependence.The bulk thermal conductivity is estimated to be 2.1 ± 0.2 W/m-K by averaging over all the samples.This prediction is in agreement with the GK predictions in Ref. [47] within the uncertainties, but larger than the MD-based 6. (Color online) Thermal conductivities of a-SiO 2 and a-Si predicted using the GK method and Eq. ( 1).For a-SiO 2 , the GKpredicted thermal conductivity is size independent, indicating that there is not an important contribution from propagating modes.For a-Si, there is a clear size dependence, indicating the importance of propagating modes.direct-method predictions in Ref. [83].Shenogin et al. predicted the total thermal conductivity of a-SiO 2 using nonequilibrium MD simulations of the same small structures used in this work [84].They find 2.0 W/m-K for their largest system, which was based on a 972-atom model tiled six times in one direction.Our GK-predicted value is larger than experimental measurements, which range between 1.3 and 1.5 W/m-K [10,14,16,17], which may be due to the classical nature of the MD simulation and/or the suitability of the BKS interatomic potential for modeling thermal transport in a-SiO 2 [47,83].Quantum statistical effects are considered later in this section.
For a-Si, there is a clear system-size dependence of thermal conductivity.Because the low-frequency DOS has the form of Eq. ( 3) and the diffusivities scale as ω −2 , the thermal conductivity will scale as the inverse of the system size.The bulk value can be found by extrapolating to an infinite system size [40,85,86].The extrapolation is performed using the three largest system sizes [87], leading to a bulk value of 2.0 ± 0.2 W/m-K, where the uncertainty is estimated from the ensemble averaging for each system size.Our extrapolated bulk value is in reasonable agreement with experimental values for a wide range of thin-film thicknesses (see Fig. 7 in Sec.V B).We note that a-Si can be only prepared experimentally as a thin film, where voids and other inhomogeneities are unavoidable [4,7,8,35,36] and can influence the vibrational structure at low frequencies.
To predict thermal conductivity from Eq. ( 1), we use the parameters B and ω cut specified in Sec.IV E assuming an ω −2 scaling below ω cut and the AF-predicted diffusivities.For a-SiO 2 , the propagating, nonpropagating, and total thermal conductivities are 0.10 ± 0.05, 1.9 ± 0.1, and 2.0 ± 0.1 W/m-K (see Table II).The uncertainties are estimated by varying ω cut and the AF broadening by 10%.The total value agrees with the GK value within the uncertainties.For the propagating contribution, using an expression similar to Eq. ( 2), Baldi et al. [18] estimated 0.1 W/m-K and Love and Anderson [15] estimated 0.03 W/m-K.
By using the ω −2 diffusivity scaling for a-Si, the propagating, nonpropagating, and total thermal conductivities are 0.6 ± 0.1, 1.2 ± 0.1, and 1.8 ± 0.2 W/m-K.This value for total thermal conductivity is in agreement with the GK-predicted bulk value within the uncertainties.Earlier studies using similar models of a-Si found that k pr is less than half of k vib [4,6], in agreement with our results.A recent study of a-Si modeled using the Tersoff potential found k pr ≈ k AF [9].Regner 1000 nm [10] Expt.[16] Expt.[17] (a) a-SiO 2  24)] for a-SiO 2 compared with experimental broadband frequency domain reflectance measurements by Regner et al. [10] and thin-film measurements from Refs.[16,17].The predicted thermal conductivity accumulation demonstrates that the propagating contribution is negligible in our model, which is in accord with the experimental measurements.(b) Predicted thermal conductivity accumulation function for a-Si compared with experimental measurements by Regner et al. and thin films fabricated by sputtering (Experiment A) [5,31,32] and chemical vapor deposition (Experiment B) [7,8,30,33].The predicted thermal conductivity accumulation demonstrates that the propagating contribution is significant for a-Si.We note that thermal conductivities as high as 6 W/m-K (not plotted) have been measured for a-Si thin films deposited using hot-wire chemical vapor deposition [8].
Estimates based on experimental measurements have found k pr to be as low as 20% [5,6] and as high as 80% of k vib [7,8].
If an ω −4 lifetime scaling is assumed for a-Si, the thermal conductivity diverges at low frequency.We bound the thermal conductivity by assuming the sample to be a thin film of thickness t f and modify the lifetimes using the Matthiessen rule [38] 1 Using the largest film thickness from the experimental literature (80 μm) [7] gives a propagating contribution to thermal conductivity of 3.0 ± 0.4 W/m-K, which is significantly larger than the GK-predicted value.Using the ω −2 scaling and this film thickness gives a propagating contribution of 0.6 W/m-K (i.e., there is no change from the bulk value).While predictions for k pr for a-Si vary based on the assumed scaling of the low-frequency vibrational lifetimes, all evidence supports that k pr is a significant fraction of the total thermal conductivity [4][5][6][7][8][9][10].
To this point we approximated the specific heat of the propagating and nonpropagating modes by the classical, harmoniclimit value of k B .At a temperature of 300 K, the quantum heat capacity [Eq.( 4)] at the largest cutoff frequency for either a-SiO 2 or a-Si is 0.98k B , justifying the use of the classical specific heat in the propagating term in Eq. ( 2).For the AF contribution, however, the effect of the quantum specific heat is important.At the highest frequency in each of a-SiO 2 and a-Si, the specific heat is 0.073k B and 0.47k B .Using Eq. ( 4) in Eq. ( 7) gives AF thermal conductivities of 1.4 ± 0.1 and 1.0 ± 0.1 W/m-K for a-SiO 2 and a-Si (Table II).This correction brings the estimate of k vib for a-SiO 2 into good agreement with experimental measurements [10,14,16,17].For a-Si, the modified k AF is 20% lower than the classical-limit value.

B. Accumulation function
In their broadband frequency domain thermoreflectance measurements, Regner et al. [10], adopting the convention of Koh and Cahill [19], interpret the measured thermal conductivity at a given thermal penetration depth to be representative of the thermal conductivity accumulation function at a MFP equal to the thermal penetration depth.Their results are plotted in Fig. 7(a) for a 1000-nm-thick film of a-SiO 2 and in Fig. 7(b) for 500-nm-and 2000-nm-thick films of a-Si.The vertical coordinate of any point on the accumulation function represents the thermal conductivity that comes from phonons with MFPs less than the horizontal coordinate at that point.Also plotted in Figs.7(a) and 7(b) are experimental measurements of thin-film thermal conductivities.For a-Si, the experimental measurements are broadly grouped by sample preparation technique: (A) chemical vapor deposition [7,8,33] and (B) sputtering [5,31,32].
Based on the results in Sec.IV E, we build thermal conductivity accumulation functions for a-SiO 2 and a-Si from where cut is the MFP at the cutoff frequency, * is the maximum MFP considered in the thermal conductivity accumulation, k( ) is the thermal conductivity as a function of MFP [21], and the propagating mode MFPs are calculated using lifetimes from Eq. ( 23).The nonpropagating contribution k AF is evaluated using the quantum specific heat (see Sec. V A).The results are plotted for a-SiO 2 in Fig. 7(a) using an infinite film thickness and for a-Si in Fig. 7(b) using a film thickness of 80 μm [88].
The predicted thermal conductivity accumulation function for a-SiO 2 saturates at a MFP of 10 nm, which is on the order of the finite size of our model.This result is in good quantitative agreement with the thermal penetration depth-independent thermal conductivity measurements using broadband frequency domain thermoreflectance [10] and experimental measurements that show minimal film-thickness dependence [16,17].
For a-Si, the low-MFP plateau of thermal conductivity in the measurements of Regner et al. is consistent with our predicted k AF .The propagating contribution to the accumulation is predicted using ω −2 and ω −4 lifetime scalings, which have both been inferred from thin-film experiments [4][5][6][7][8]23,34].Predictions for both the ω −2 and ω −4 scalings pass reasonably through the thin-film thermal conductivity measurements, particularly for thicknesses in the 50-2000 nm range.Overall, the film-thickness-dependent measurements show a large variation which results from the deposition method and impurity concentration [7,8,35,36].The measurements of Regner et al. show sharper accumulations than either the ω −2 or ω −4 scalings, particularly for the 2000-nm film.For the ω −2 scaling, which best matches our model [see Fig. 4(b)], the thermal conductivity accumulation saturates at 1 μm, in good agreement with where the measurements of Regner et al. saturate for their 500-nm film.The 2000-nm film accumulation shows no sign of saturation.

VI. SUMMARY
We investigated the contributions of propagating (k pr ) and nonpropagating (k AF ) modes to the total vibrational thermal conductivity (k vib ) of a-SiO 2 and a-Si using the NMD method (Sec.IV D), AF theory (Sec.IV E), and the GK method (Sec.V A).The atomic structures of a-SiO 2 and a-Si play an important role in determining the mode-level properties needed to predict the propagating and nonpropagating contributions.The propagating regime ends at a lower frequency for a-SiO 2 , which is evident from the DOS (Fig. 2) and the effective dispersion extracted from the structure factors [Fig.3(a)].This smaller maximum frequency of propagating modes is due, in part, to the weak bonding that exists between the SiO 4 tetrahedra in a-SiO 2 , while a-Si is formed by a network of strongly bonded tetrahedra.The structural differences are also apparent in the low-frequency scalings of the mode lifetimes (Fig. 4) which show a clear ω −2 dependence (i.e., phononlike) for a-Si, but not for a-SiO 2 .The combined effect of all the mode-level properties results in a significant difference in the propagating and nonpropagating contributions to thermal conductivity for a-SiO 2 and a-Si (Table II).
For our model of a-SiO 2 , the contribution from propagating modes is negligible (∼6%).Our predictions align with experimental measurements of the film-thickness independence of thermal conductivity [16,17] and thermal penetration-depth independence in the measurements of Regner et al. [10].While the finite size of our model makes it difficult to identify a clear scaling of the low-frequency lifetime scaling, experiments show that both ω −2 and ω −4 scalings exist in a-SiO 2 [25,26,28].Still, in previous studies the propagating contribution to thermal conductivity has been found to be negligible [15][16][17].
For our model of a-Si, the thermal conductivity has a significant (∼35%) contribution from propagating modes that are best described by a lifetime scaling of ω −2 .Our predicted nonpropagating thermal conductivity contribution is in good agreement with the plateau at low-MFP for both films studied by Regner et al.For both films, the thermal conductivities accumulate much faster than our predictions.The large range of thermal conductivity measurements on a-Si thin films suggest that a comprehensive experimental study using thermoreflectance techniques on varying film thicknesses and preparation techniques is necessary.It may be particularly helpful to perform the experiments at temperatures less than 10 K, where the propagating contribution dominates for both a-SiO 2 and a-Si and the low-frequency lifetime scaling, which is still under debate, can be better resolved.

FIG. 1 .
FIG. 1. (Color online) (a) Radial distribution function g(r) for the 4608-atom a-SiO 2 structure created from a melt-quench technique.The radial distribution function compares well with the experimental measurement from Ref. [54].Inset: Small sample of the a-SiO 2 structure showing the Si-O tetrahedral bond network.Bond lengths range between 1.6 and 1.8 Å.(b) Radial distribution function for the 4096-atom a-Si structure created by the modified WWW algorithm.The radial distribution function compares well with the experimental measurement from Ref. [57].Inset: Small sample of the a-Si structure.Bond lengths range between 2.3 and 2.7 Å.

2 FIG. 2 .
FIG. 2. (Color online)Vibrational DOS of a-SiO 2 and a-Si plotted on a log-log scale.Both models show an ω 2 scaling at low frequency.The dashed lines indicate an extrapolation of the DOS based on this scaling.The DOS for a-Si has two peaks similar to the DOS of the crystalline phase[62].The DOS for a-SiO 2 is flat over most of the spectrum, with a high-frequency gap that separates the modes involving Si-O interactions[47].

2 FIG. 4 .
FIG.4.online) Vibrational mode lifetimes predicted by NMD [Eq.(20)] and the time scales extracted from the structure factors [Eq.(13)] (a) a-SiO 2 and (b) a-Si.For a-Si, a clear ω −2 scaling is observed at low frequencies, while the lifetimes plateau at higher frequencies and over a wider range of frequencies than for a-SiO 2 .

2 .
FIG.4.online) Vibrational mode lifetimes predicted by NMD [Eq.(20)] and the time scales extracted from the structure factors [Eq.(13)] (a) a-SiO 2 and (b) a-Si.For a-Si, a clear ω −2 scaling is observed at low frequencies, while the lifetimes plateau at higher frequencies and over a wider range of frequencies than for a-SiO 2 .lifetimes are used to determine the coefficient of the lowfrequency scaling [see Figs.5(a) and 5(b)].There is no clear evidence for an ω −2 scaling in a-SiO 2 , which would correspond to propagating modes.midfrequencies, the NMD lifetimes are approximately constant and there is peak near 2 ×10 14 rads/s, which corresponds to the peak in the DOS (see Fig.2).The time scales predicted from the structure factor fall below the NMD-predicted lifetimes and the IR limit by up to one order of magnitude.These low values result the structure factors for a-SiO 2 are evaluated for wave vectors where the resulting wave packets are formed by nonpropagating modes[4,6].For a-Si, the NMD lifetimes show a clear ω −2 scaling at low frequency.The lifetimes plateau at higher frequencies, over a wider range of frequencies than for a-SiO 2 , with two peaks corresponding to peaks in the DOS (see 2).A similar plateau of lifetimes at high frequencies has reported for lattices[66,75] and in another study of a-Si[9].The transition from the low-frequency scaling to the plateau region occurs near 10 13 rads/s, which to where the DOS first peaks Fig.2.Similar behavior has been observed for models of disordered lattices[66].The time scales predicted by the structure factors are in good agreement with those predicted by NMD at low frequencies.Similar agreement has been reported in other models of amorphous[6,[76][77][78].The agreement between the NMD-predicted lifetimes and the structure factor time scales for a-Si at low frequencies indicates FIG. 5. (Color online) Vibrational mode diffusivities predicted from NMD [using Eqs.(5) and(20) with the DOS sound speed from TableI] and the AF theory [Eq.(8)].Also shown are extrapolations based on an ω −2 scaling with Eqs.(5) and (6) for a-SiO 2 and a-Si, and an additional ω −4 scaling for a-Si.For both systems, the diffusivities are larger than the high-scatter limit [Eq.(21)] except at high frequencies, where the modes are localized.

FIG. 7 .
FIG. 7. (Color online) (a) Predicted thermal conductivity accumulation function [Eq.(24)] for a-SiO 2 compared with experimental broadband frequency domain reflectance measurements by Regner et al.[10] and thin-film measurements from Refs.[16,17].The predicted thermal conductivity accumulation demonstrates that the propagating contribution is negligible in our model, which is in accord with the experimental measurements.(b) Predicted thermal conductivity accumulation function for a-Si compared with experimental measurements by Regner et al. and thin films fabricated by sputtering (Experiment A)[5,31,32] and chemical vapor deposition (Experiment B)[7,8,30,33].The predicted thermal conductivity accumulation demonstrates that the propagating contribution is significant for a-Si.We note that thermal conductivities as high as 6 W/m-K (not plotted) have been measured for a-Si thin films deposited using hot-wire chemical vapor deposition[8].

TABLE II .
Thermal conductivities for bulk a-SiO 2 and a-Si