Vibrational mean free paths and thermal conductivity of amorphous silicon from non-equilibrium molecular dynamics simulations

The frequency-dependent mean free paths (MFPs) of vibrational heat carriers in amorphous silicon are predicted from the length dependence of the spectrally decomposed heat current (SDHC) obtained from non-equilibrium molecular dynamics simulations. The results suggest a (frequency)$^{-2}$ scaling of the room-temperature MFPs below 5 THz. The MFPs exhibit a local maximum at a frequency of 8 THz and fall below 1 nm at frequencies greater than 10 THz, indicating localized vibrations. The MFPs extracted from sub-10 nm system-size simulations are used to predict the length-dependence of thermal conductivity up to system sizes of 100 nm and good agreement is found with separate molecular dynamics simulations. Weighting the SDHC by the frequency-dependent quantum occupation function provides a simple and convenient method to account for quantum statistics and provides reasonable agreement with the experimentally-measured trend and magnitude.


I. INTRODUCTION
Compared to heat transfer by phonons in crystalline materials, heat transfer in amorphous materials is complicated by the existence of three regimes of vibrational modes. 1 Low-frequency propagons are delocalized and have a well-defined wave vector and group velocity, 2 similar to phonons in crystals, while high-frequency locons are localized and contribute negligibly to thermal conduction. 3 Diffusons have intermediate frequencies and are delocalized, but do not have well-defined wave vectors or group velocities. The contribution of diffusons to thermal conduction can be notable, however, as they occupy the majority of the vibrational spectrum. 2 From kinetic theory, 4 the contribution of an individual phonon or propagon mode to thermal conductivity is proportional to its mean free path (MFP). Because diffusons do not have a well-defined group velocity, it is not clear if they have a MFP or how it can be defined. Their contribution to thermal conductivity can be predicted using their diffusivity, which is well-defined and can be calculated from Allen-Feldman theory. 2,5 Nevertheless, it would be insightful to identify a frequencydependent length scale for propagons and diffusons describing the decay of the heat flux at each vibrational frequency. Such a definition would lift the (fundamentally) arbitrary distinction between propagons and diffusons and enable a unified description of heat transfer at all vibrational frequencies.
In this paper, we apply the spectrally-decomposed MFP method 6 to probe the non-equilibrium MFPs of vibrational heat carriers in amorphous silicon (a-Si). This method is based on calculating the spectrallydecomposed heat current (SDHC) 7 in systems of differ-ent lengths using non-equilibrium molecular dynamics (NEMD) simulations. The MFPs are determined from the variation of the SDHC as a function of system length at each vibrational frequency. 8 We previously used the spectrally-decomposed MFP method to calculate the non-equilibrium MFPs in lowdimensional systems such as carbon nanotubes 6 and anharmonic chains. 8 We demonstrated that the nonequilibrium MFPs transparently describe the ballistic-todiffusive transition in the length-dependence of thermal conductivity. Compared to previous calculations for a-Si, the spectrally-decomposed MFP method has several advantages. Unlike in modal life-time calculations, 9 we do not need to estimate the group velocities of individual modes to calculate their MFPs. We also do not need to distinguish between propagons and diffusons 10 nor resort to the harmonic approximation. 2 In contrast to recent calculations studying the spectral conductivity of a-Si in fixed-size systems, 11,12 we focus on the MFPs of heat carriers and the system-size dependence of thermal conductivity.
The rest of the paper is organized as follows. The calculation methods are presented in Sec. II and the numerical results are discussed in Sec. III. We also introduce a simple method for the quantum correction of thermal conductivity from classical MD simulations, based on weighting the SDHC by the quantum occupation function. Because this quantum correction method operates at the frequency level, it is more reasonable than quantum-correction methods that operate at the system level (see Ref. 13 and references therein) and allows us to compare our predictions to experimental measurements.

II. SIMULATION SETUP AND METHODS
All simulations are carried out using the LAMMPS package 14 with a time step of 2.5 fs. The Si-Si interactions are modeled by the Stillinger-Weber potential 15 with the parameters of Ref. 16. The NEMD simulation geometry is shown in Fig. 1. The atomic coordinates for a-Si are generated by following the melt-quench procedure of Ref. 17 and the final density is 2,291 kg/m 3 . After the equilibration of the quenched system, atoms located within a distance L bath = 5 nm from the left and right edges of the structure are coupled to Langevin heat baths at temperatures T H = T + ∆T /2 and T C = T − ∆T /2 with bath relaxation times of 1 ps. To prevent sublimation, atoms at the far left and right edges are fixed to their equilibrium positions. Periodic boundary conditions are applied at the boundaries transverse to the current flow. The width of the system cross-section is 7 nm. System lengths L (i.e., the region between the baths) between 1 and 10 nm at intervals of 1 nm are considered.
The SDHC is calculated through the plane of decomposition located halfway between the hot and cold baths (dashed line in Fig. 1). 6,8 The SDHC q i→j (ω) between particles i and j located on opposite sides of this plane is given by the pair-wise SDHC equation 7 where t simu is the simulation time, ω is the angular frequency, and the interatomic force constant K αβ ij is defined as The velocitiesv α i (ω) andv β j (ω) are the discrete Fourier transforms of the atomic velocities v α i (t) =u α i (t) and v β j (t) =u β j (t) (the exact definitions are in Ref. 6), where u α i and u β j are the displacements of atoms i and j from their equilibrium positions in directions α,β ∈ {x, y, z}. In Eq. (2), V is the interatomic potential energy function. The spectral flux through the plane of decomposition is obtained from Eq. (1) by summing over all pairs of atoms (one of the left side, denoted byL, and one on the right side, denoted byR) within the potential cut-off distance of each other and dividing by the interface area A: While Eq. (1) is the first-order approximation to the inter-particle SDHC, 7 we have confirmed that the contribution of higher-order terms is negligible for a-Si by comparing the integral of Eq. (3), which we denote as Q, to the total flux determined from the work done by the heat baths. The two results agree within 4%. We attribute this good agreement to the stiffness of the interatomic bonds in a-Si, which ensures that the first-order term in the current (proportional to the harmonic force constants K αβ ij ) dominates the higher-order terms that are related to anharmonic force constants. This restriction to the first-order current term at the plane of decomposition does not, however, mean that anharmonic scattering in the bulk is neglected, because all anharmonic effects are included in the NEMD simulations. 6,7 Frequency-dependent MFPs Λ(ω) are calculated by determining q(ω, L) for different system lengths L and fitting the length-dependent q(ω, L) to the equation 6 where q 0 (ω) is the spectral flux when the baths are in contact (L → 0 + ). Both q 0 (ω) and Λ(ω) are determined from the fitting procedure. While q 0 (ω) depends on the details of the heat baths, the MFPs extracted from the length-dependence are not expected to depend on the bath details. 8 The frequency-dependent MFPs determined from Eq. (4) are mode-averaged and projected along the direction of heat transfer. 6 The MFPs are independent of the system length and therefore correspond to the bulk values. 6 We note that the MFPs determined from Eq. (4) correspond to the decay length of the heat flux. This definition does not necessarily coincide with the conventional definition of the MFP as the decay length of a wave packet. 4 This distinction is important for diffusons, which do not have a well-defined wave-vector so that the traditional definition cannot be applied.
Once the spectral MFPs are determined, the thermal conductivity κ for length L can be determined from 6 The length-dependence of the thermal conductivity can be intuitively understood by writing Eq. (5) in the equivalent form where the "effective" MFP Λ eff (ω) has been introduced, which is similar to the well-known Matthiessen rule. 18 The effective MFP accounts for boundary scattering through the additional L/2 term and is limited to below this value. Finally, Eq. (5) allows for a simple quantum correction to the thermal conductivity prediction, because the contributions of different frequencies can be weighted by the vibrational mode energy and occupation, as in the Landauer-Büttiker formalism [19][20][21] . We define the quantum corrected thermal conductivity as where f BE (ω, T ) = [exp( ω/k B T ) − 1] −1 is the Bose-Einstein distribution function, k B is the Boltzmann constant, and is the Planck constant divided by 2π. By defining the dimensionless, length-dependent bath-tobath transmission function as Eq. (7) can be written in the familiar Landauer-Büttiker form as 19 The proposed quantum correction accounts for the quantum specific heat of the modes at each frequency, but does not account for quantum effects in the dynamics. The method is thus similar to the one recently introduced by Lv and Henry, 11 who weight the modal contributions to the equilibrium Green-Kubo thermal conductivity by the quantum population function. All our NEMD simulations are performed at a mean temperature of T = 300 K with temperature bias ∆T = 100 K. Choosing a relatively large temperature bias allows for very good signal-to-noise ratio in the spectral heat flux, suppressing the statistical noise. We checked that halving the bias to ∆T = 50 K does not change the spectral MFPs. In addition, we checked that the heat flux is not sensitive to the exact arrangement of atoms arising from the melt-and-quench procedure, which we attribute to the large cross-section of the system giving rise to spatial averaging in the spectral currents. Therefore, we performed a single melt-quench for each system length.

III. RESULTS
The spectral heat flux q(ω) for selected system lengths L as a function of frequency f = ω/(2π) is plotted in Fig. 2. The spectral distribution of the heat flux for L = 20 nm was recently analyzed in detail by Zhou and Hu, 12 so we focus here on its length-dependence. As expected, increasing the system length reduces the heat current throughout the whole frequency range because of increased phonon-phonon scattering. The reduction is strongest at high frequencies, where the MFPs are shorter compared to low frequencies. At frequencies less than 2 THz, the spectral current is nearly independent of system length. Such nearly ballistic conduction suggest that the low frequency MFPs are notably longer than the system sizes considered here.
Equation (4) suggests that the inverse of the spectral flux will be linearly proportional to the system length, with the slope given by [2Λ(ω)] −1 . To determine the MFPs, we calculated the spectral flux for system sizes L ∈ {1, 2, . . . , 10} nm, plotted q(ω) −1 versus L and fitted a linear function at each frequency using least squares fitting. 6,8 A linear function accurately reproduces the length-dependence of q(ω) −1 (not shown), as previously also observed for other systems. 6,8 Because of the high computational cost associated with calculating the spectral heat fluxes for large systems, we limited our study to systems at most 10 nm long. This restriction precludes extracting MFPs longer than 10 nm accurately, limiting the current analysis to frequencies greater than 2 THz. The MFPs Λ(ω) extracted from the linear fitting procedure are shown in a log-log plot in Fig. 3. At frequencies below 5 THz, the MFPs obey a power-law scaling Λ(ω) ∝ ω −2 . This scaling agrees with modal life-time calculations on a-Si. 10 At frequencies greater than 5 THz, the power-law scaling breaks down and the MFPs increase with increasing frequency, giving rise to a local maximum around 8 THz. A similar maximum for a-Si has been observed in effective MFPs 9 and in lifetimes. 10 At higher frequencies, the MFPs decrease again and fall below 1 nm, which is on the order of the silicon-silicon bond length. At such high frequencies, the uncertainty is large because of the sensitivity of the spectral flux to the system size.
Larkin and McGaughey reported a propagon-diffuson transition frequency of 1.8 THz, 10 such that the frequency range considered in Fig. 3 mostly corresponds to diffuson-like modes. In the analysis below, we assume that the scaling Λ(ω) ∝ ω −2 (solid line in Fig. 3) remains valid at frequencies below 2 THz. With such scaling, the MFPs exceed 100 nm below 530 GHz and 1 µm below 170 GHz. While the ω −2 scaling may break down in real situations because of defects, boundary scattering, or even the onset of a Rayleigh-like ω −4 scaling at very low frequencies, 2 we assume it to hold for simplicity.
We now investigate the length-dependence of the thermal conductivity using Eq. (5). The calculated thermal conductivity (continuous line) is compared with that determined directly from NEMD simulations (data points) for lengths up to 100 nm in Fig. 4. In the evaluation of the integral in Eq. (5), the MFPs have been assumed to scale as Λ(ω) ∝ ω −2 at frequencies below 2 THz. These calculations were carried out without the quantum correction as the NEMD simulations are classical. Equation (5) combined with the MFP data of Fig. 3 reproduces the length-dependence of thermal conductivity up to lengths L = 100 nm to within 2%. This close agreement (i) supports the assumption of Λ(ω) ∝ ω −2 scaling at low frequencies, (ii) shows that the MFP data in Fig. 3, which were determined from simulations of systems shorter than 10 nm, can be reliably used to estimate the relative contributions of different vibrational frequencies to thermal transport in much larger systems, and (iii) provides support for the accuracy of Eq. (5) in describing the length-dependence.
Because the Debye temperature of a-Si is 530 K 22 , which is well above room temperature, we need to apply the quantum correction to compare the predicted temperature-dependence of thermal conductivity to experimental data. To do so, we use Eq. (7) and evaluate the integral as a function of temperature using the MFPs from Fig. 3, again assuming the scaling Λ(ω) ∝ ω −2 at frequencies below 2 THz. A full quantum-corrected analysis would require determining the MFPs at each temperature. For simplicity and based on the recent results of Lv and Henry 11 , we assume that the MFPs calculated at a temperature of 300 K remain valid at other temperatures.
The quantum-corrected thermal conductivity [Eq. (7)] is plotted as a function of temperature for system lengths of 50, 250, and 1000 nm in Fig. 5. As noted above, assuming finite L in Eq. (7) limits the MFPs to L/2. Experimental data from Cahill et al. for a 520 nm thick film of hydrogenated a-Si with one atomic percent hydrogen content are also plotted. 23 Because available thermal conductivity measurements for a-Si contain significant scatter, 10   MFPs are assumed to scale as Λ(ω) ∝ ω −2 below frequencies of 2 THz and to be independent of temperature. Also plotted is the thermal conductivity of a 520 nm thick hydrogenated a-Si thin film measured by Cahill et al. 23 .
NEMD simulations. The increase of thermal conductivity with increasing temperature is well described by the quantum-corrected thermal conductivity. At temperatures higher than 300 K, the experimentally measured thermal conductivity increases slightly slower as a function of temperature than our prediction, but this disagreement may be related to our approximation that the MFPs are independent of temperature. At such high temperatures, anharmonic scattering will reduce MFPs and therefore decrease the thermal conductivity. Without the quantum-correction, the predicted thermal conductivity would depend only very weakly on temperature (due to the weak temperature-dependence of the MFPs), precluding reasonable agreement with the trend of the experimental data. We caution that this quantum correction has only been examined for a-Si here and that its application to other systems warrants further investiga-tion.

IV. CONCLUSION
We investigated vibrational heat transfer in a-Si by determining the SDHC and MFPs from NEMD simulations. The calculated MFPs directly reflect the decay of the heat flux at each vibrational frequency and do not rely on the existence of a well-defined modal wave vector, thereby avoiding the separate treatment of diffusons and propagons. As shown in Fig. 3, the MFPs exhibit ω −2 scaling at frequencies above 2 THz and below 5 THz. At frequencies higher than 10 THz, the MFPs fall below 1 nm, corresponding to strongly localized vibrations. The length-independent MFPs can be used to accurately predict the thermal conductivity in systems as long as 100 nm (Fig. 4). Weighting the SDHC by the frequencydependent quantum occupation function provides a simple method for a quantum-correction of thermal conductivity and is able to reproduce the experimentally measured temperature-dependence of thermal conductivity, as shown in Fig. 5.
In the future, it would be useful to calculate the SDHC for systems longer than those considered here, enabling direct extraction of MFPs at frequencies below 2 THz. Such an analysis could inform the ongoing discussion 10 of the low-frequency scaling of MFPs in a-Si. It will also be important to compare the non-equilibrium MFPs to those calculated from equilibrium molecular dynamics simulations.

V. ACKNOWLEDGEMENTS
This work was initiated during the research visit of K.S. to Carnegie Mellon University. The computational resources were provided by the Finnish IT Center for Science and Aalto Science-IT project. The work was partially funded by the Aalto Energy Efficiency Research Programme (AEF) and the Academy of Finland.