Quantifying Modal Thermal Conductivity in Amorphous Silicon
11 Quantifying Modal Thermal Conductivity in Amorphous Silicon
Yanguang Zhou Department of Mechanical and Aerospace Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China
Abstract
While there are several methods, e.g., anharmonic lattice dynamics and normal mode decomposition, to compute the modal lattice vibrational information in perfect crystals, the modal information of vibrations, e.g., vibrational relaxation time, group velocity and mean free path, in amorphous solids are still challenge to be captured. By systematically analyzing the normal mode decomposition and structure factor methods, we conclude that the vibrational dispersion can be calculated by applying effective wave vectors in the structure factor method, while the vibrational relaxation time calculated by the normal mode decomposition method is questionable since the group velocity cannot be defined on the Gamma point. We also show that the anharmonicity caused by the system temperature has little effect on the relaxation times of the propagating modes in amorphous materials, and therefore, the corresponding modal and total thermal conductivity is temperature independent when all the vibrations are assumed to be excited. The non-propagating modes, i.e., diffusons, conduct heat via thermal coupling between different vibrational modes, and can be calculated by harmonic lattice dynamics using Allen-Feldman theory. As a result, the thermal * Author to whom all correspondence should be addressed. E-Mail: [email protected] (Y.Z.) conductivity contributed from diffusons is also temperature independent when all the vibrational modes are activated which is the situation in molecular dynamics simulations. The total thermal conductivity concerning both propagons (50%) and diffusons (50%) agree quite well with the results computed using Green-Kubo equilibrium molecular dynamics. By correcting the excitation state of the vibrations in amorphous solids, the thermal conductivity calculated by the structure factor method and Allen-Feldman theory can fully capture the experimentally measured temperature-dependent thermal conductivity. Our results here show that the phonon-gas theory is still valid in amorphous materials by assuming the system with the effective possible wave vectors, and the inaccuracy of scattering process in amorphous systems due to Boltzmann distribution can be ignored. Our study provides a fundamental understanding of the temperature-dependent thermal excitations at play in amorphous materials
I. INTRODUCTION
Amorphous material, e.g., amorphous silicon, has many promising physical, chemical and mechanical properties such as extremely low thermal conductivity [1–10], exceptional charge capacity [11–13] and high strength [14,15], comparing to its crystalline counterpart due to the local variation of the topology of atomic connectivity. Understanding the contributions to thermal conductivity from various heat carriers is important because the means by which their contributions can be manipulated can be explored once the dominant heat carriers in the amorphous systems are characterized. While it is well-known that phonons, i.e., an energy quantum of propagating lattice waves, are the main heat carriers in perfect crystals [16], the heat carriers in amorphous materials are much more complicit and difficult to be identified. The heat carriers in amorphous solid systems, e.g., amorphous silicon and amorphous silica [1], can be classified as propagating and delocalized phonon-like modes (phonons), non-propagating and delocalized modes (diffusons), and non-propagating and localized modes (locons) depending on the degree of delocalization of atomistic vibrations and their corresponding vibrational mean free paths (MFPs) [17–19]. To quantitatively predict the contribution to thermal conductivity from various heat carries, one requires the diffusivity of the non-propagating vibrational modes [17], i.e., diffusons, and the group velocities and relaxation times of the propagating vibrational modes [2], i.e., propagons. Locons do not contribute to thermal transport significantly in the three-dimensional systems [20]. The contribution to the thermal transport from diffusons can be calculated using Allen-Feldman (AF) theory [17], which are coupled through the harmonic vibrations. As a result, the thermal conductivity contributed from the diffusons should be temperature independent when all the vibrations are excited. Molecular-dynamics-based normal mode decomposition (NMD) method [1,21] and structure factor (SF) approach [2,22–24] are the two most popular computational methods to calculate the vibrational relaxation times. Unlike the group velocity that cannot be well defined in the NMD method, the effective group velocity in the SF approach can be computed from the SF spectrum. Therefore, the modal thermal conductivity of the propagating vibrational modes in the SF method can be obtained via the phonon-gas theory v g C v [2,4,10] in which v C is the volumetric heat capacity, g v is the group velocity and is the vibrational relaxation time. Meanwhile, using NMD analysis, the vibrational relaxation time of these propagating modes shows a strong reduction when the system temperature is increasing [25], and therefore the corresponding thermal conductivity should be decreasing with temperature. This is in contradiction with the experimental observations at relatively high temperatures, i.e., above 300 K, when all the heat carriers are thought to be excited [26]. Moreover, modal thermal conductivity calculated by Green-Kubo modal analysis (GKMA) [3] and the total thermal conductivity computed using Green-Kubo equilibrium molecular dynamics (GKEMD) show weak temperature dependences in amorphous silicon and amorphous silica [25]. Most recently, Moon et. al. [2,22] numerically and experimentally shows the vibrational relaxation times of the propagating modes in amorphous silicon are temperature independent based on the SF analysis. Despite these efforts in characterizing the modal thermal conductivity in amorphous materials, a systematic and clear understanding of the thermal transport mechanisms of various heat carriers in amorphous remains challenging and poorly documented. In this paper, we investigate the modal propagating and non-propagating modes’ contributions to the thermal conductivity, and the temperature-dependent thermal conductivity of amorphous silicon using normal mode analysis which is a combination of normal-mode decomposition, dynamical structure factor analysis and Allen-Feldman theory. We emphasize that the relaxation time in amorphous silicon calculated by the normal mode decomposition method using only Gamma wave vector is questionable since the group velocity cannot be well defined. The paper is organized as follows. The computational details are given in Sec. II. The modal analysis methods on normal mode decomposition, structure factor and Allen-Feldman theory are discussed in Sec. III. In Sec. IV, we present the modal relaxation time, effective mean free path and the corresponding accumulative thermal conductivity, and the temperature-dependent thermal conductivity calculated using GKEMD and the approach in combination of structure factor and Allen-Feldman theory. We discuss the possible reason for the failure of the NMD method in amorphous systems in Sec. V. Conclusions are given in Sec. IV. II. COMPUTATIONAL DETAILS
A silicon crystal including 512 atoms with the interactions depicted by Tersoff potential [27] was firstly equilibrated at a high-temperature liquid state, i.e., 4000 K, using NVT (constant particle number, volume and temperature) ensemble for 250 ps and then cool down to 300 K in the NPT (constant particle number, pressure and temperature) ensemble within 2.5 ns. The structure is finally equilibrated at 300 K for another 250 ps in the NVT ensemble to reduce the metastabilities. The generated amorphous silicon was annealed from 30 K to 800 K in the NVT ensemble to simulate the amorphous systems at various temperatures. The heat current Q in the NVE (constant particle number, volume and energy) ensemble can be calculated for 1 ns to 10 ns via [28]: , ; , , ; N N Ni i ij i ij ijk ij iki i j i j i j k i j k
Q e v F v r F v r r , (1) in which i e refers to the ionic energy, i v is the atomic velocity, ij r the distance between two ions i and j and ij F and ijk F represent the two-body and three-body forces, respectively. Then, the total thermal conductivity GK can be obtained using the computed heat current based on the Green-Kubo theory [29] GK b
Q t Q dtk VT , (2) in which, b k is the Boltzmann constant, V the system volume, T represents the system temperature and t denotes the autocorrelation time. The angular bracket indicates ensemble averaging. For each case, 30 independent runs are performed to obtain a stable averaged value of GK . The correlation time considered in our simulation is from 400 ps to 4.5 ns, which is long enough to obtain the converged thermal conductivity ( Figure 1a ). Here, we emphasize that the long autocorrelation time as shown in
Figure 1a at low temperatures, e. g., 30 K, is to make sure the ergodicity of the systems. We also test the size effects on our results. It is shown the thermal conductivity will increase around 0.2 W/mK when we increase the number of the atoms in the system from 512 to 4096 (
Figure 1b ). Considering the computational efficiency, we thereafter select the models with 512 atoms in all our later calculations. The time step in all the simulations is 0.5 fs, and periodic boundary conditions were applied in all three directions. The thermal conductivity calculated using GKEMD in our simulations is from 1.51 W/mK to 1.77 W/mK at 300 K, which is in accordance with the values from Moon el. al. (1.47 W/mK) [2] and Lv and Herny (1.7 W/mK) [25] using the same potential. For the NMD calculations (Sec. IIIA), the atomistic velocity which are used as the inputs, are dumped every 0.02 ps in a 400 ps NVE simulation, following the 250 ps NVT runnings. For structure factor calculations [Eqs. (10) and (11) in Sec. IIIB], all the simulations start with 250 ps NVT relaxation to allow the systems to reach an equilibrium state, and then the atomistic velocity which are the inputs for Eqs. (10) and (11), are dumped every 0.02 ps in the next 20 ps NVE simulation.
III. NORMAL MODE VIBRATION ANALYSIS METHODS
A. Normal mode decomposition analysis
Based on the phonon gas theory, the thermal conductivity of the crystals can be calculated under the relaxation time approximation via [30]
2( , ) ( , ) ( , ) ( , ) ph V gk C q v q q (3)where ( , ) V C q , ( , ) g v q and ( , ) q are the volumetric heat capacity, the group velocity and the relaxation time of phonon mode with wavevector q and branch . In classical molecular dynamics simulation, the volumetric specific heat can be obtained by , ( , ) / ph V b C q k V . For the amorphous solids, Larkin and MacGauhey [1] suggested (0, 0, 0) q and therefore the thermal conductivity contributed by the propagating vibrational modes can be calculated via ( ) ( ) ( ) propagon V g C v , (4)in which ( ) g v may be obtained by extrapolating the dispersions near the Gamma point as suggested by He et. al. [7], and the relaxation time ( ) can be calculated by fitting the spectral energy density [1] tt X t i t dtt , (5)using the Lorentzian function with center at ( ) and line width equal to ( ) , one can obtain ( ) /( ) ( ) [ ( ) ] ( ) C , (6) where ( ) C is a mode-dependent constant. The lifetime of the vibrational mode ( ) is then given by ( ) 1 / 2 ( ) . ( ; ) X t in Eq. (4) is time-dependent normal velocity and can be calculated via ( ; ) ( , ) ( , ) ji X t m e j v j t , (7) in which, j m is the atomic mass, ( , ) v j t is the velocity of atom j at a time t , and ( , ) e j is the eigenvector of the basic atom j and branch , which can be obtained by diagonalizing the dynamic matrix. Here, we emphasize that it is difficult to calculate the modal level thermal conductivity by the NMD analysis since the group velocity in it cannot be well defined when only the Gamma point is taken into consideration. Furthermore, the calculated relaxation time shows a strong temperature dependence which implies the thermal conductivity should also depend on the temperature, is contrary to the GKEMD results (see Sec. IIIA for details). B. Structure factor analysis
Another optional approach to calculate the relaxation time of vibrations in amorphous solids is fitting the structure factor S which can be computed via [23,24] (| |, ) (| |, ) bT or L T or L k TS q E qm , (8) where refers to the angular frequency of the lattice vibrations. q corresponds to the effective wave vector, and is computed using
22 2 yx z nn nq i j kLx Ly Lz , (9) in which, xor yor z n is the integer and should be smaller than / xor yor z min L a where xor yor z L and min a are the length of system along x or y or z direction and the minimal interatomic space, respectively. T E and L E stand for the atomic velocities related to transverse and longitudinal polarizations, respectively, which can be computed via
21 1 ˆ(| |, ) ( ) exp( (0)) ( )
N NL j jj
E q q v iq r , (10)
21 1 ˆ(| |, ) ( ) exp( (0)) ( )
N NT j jj
E q q v iq r , (11) where ˆ / | | q q q , and (0) j r refers to the equilibrium position of atom j . ( ) j v is the atomic velocity of the th j atom for th eigenmode, which can be calculated via
1( ) ( , ) ( ; ) j j v e j X tm , (12)We also calculate the T or L E only considering the harmonic vibrations through
21 1 ˆ(| |, ) ( ) exp( (0)) ( )
N NHarmonicL j jj
E q q e iq r , (13)
21 1 ˆ(| |, ) ( ) exp( (0)) ( )
N NHarmonicL j jj
E q q e iq r , (14)The calculation of the Dirac delta function in Eq. (10), (11), (13) and (14), is performed using a Gauss broadening of in , where in is the frequency interval of 0.04 THz. Finally, the relaxation time of the propagating vibrations in amorphous materials by fitting the structure factor using the Lorentzian function [Eq. (6)], and the corresponding thermal conductivity can be obtained using Eq. (4). It is worth noticing that the structure factor is a directionally averaged because the effective wavevectors are computed using Eq. (9) and the amorphous materials are isotropic. The structure factor has been successfully used to depict the dispersion relation in many disorder solids and liquids [2,10,22–24,31]. C. Allen-Feldman theory
The thermal conductivity contributed by the non-propagating modes, i.e., diffusons, can be computed using the AF theory [17]. In AF theory, the thermal conductivity can be expressed in the form of ( ) ( ) diffuson V AF C D , (15)in which, ( ) AF D is the diffuson diffusivity of the frequency with th diffuson mode. The AF diffusivities are calculated using [17] ( ) | | ( ) AF VD S , (16)where S is the heat current operator, which measures the thermal coupling between vibrational mode and based on their frequencies and spatial overlap of eigenvectors, and can be computed through ;( , ) ( )2 g S vV , (17) ;( , ) ( ; ) (0) (0) ( ; )2 g jj jjjj iv e j D r e j , (18)in which, (0) (0) (0) jj j j r r r and (0) jj D is the dynamic matrix and can be calculated using the second-order force constants (0) jj (0) (0) / jj jj j j D m m , (19)The Dirac delta function in Eq. (16) is approximated using Lorentzian broadening with the width of ave suggested by Larkin and MacGaughey [1], where ave is the average frequency spacing. In AF theory, the calculated thermal conductivity is also directionally averaged [Eq. (16)] since amorphous materials are isotropic. IV. RESULTS A.Modal relaxation time of the lattice vibrations
The NMD-predicted vibrational spectrums at three typical temperatures, i.e., 50 K, 300 K and 800 K, are plotted in
Figure 2a . It is clearly shown that the broadening of the spectrum becomes wider with the increase of the temperature, and therefore the vibrational relaxation time calculated using the NMD approach will decrease largely with the enhancement of temperature as shown in
Figure 2b . We also find the vibrational frequency will shift to a lower value with the increase of temperature due to the phonon softening (
Figure 2a ), which has been previously proved Feng et. al. [32] in crystal Si. A plateau of vibrational relaxation times at high frequencies as reported in our study is also observed in other disorder lattices [33] and other studies of a-Si [1,7]. The NMD-predicted vibrational relaxation time is ranging from several picoseconds to around 80 picoseconds at 300 K, which is following the calculated results of Lv and Herny [25] and Larkin and MacGaughey [1]. However, He et. al . [7] reported the vibrational relaxation time of amorphous Si at 300 K can be as high as several hundred of picoseconds using the same potential used in our paper. The difference of the vibrational relaxation time may come from the structural difference. We also note the total thermal conductivity predicted by He et. al. [7] using GKEMD is much higher than ours (3 W/mK versus 1.77 W/mK). Moreover, the vibrational relaxation time in our study drops significantly as temperature increases from 300 K to 800 K, which is also proved by Lv and Herny [25]. Based on the phonon gas theory [Eq. (4)], the thermal conductivity contributed by the propagons should decrease largely since the group velocity and heat capacity are temperature independent when all the vibrations are excited. The total thermal conductivity of the amorphous Si should also decrease with temperature since diffusons conduct energy through harmonic coupling, which contradicts with the phenomena that the total thermal conductivity is almost temperature independent calculated by GKEMD ( Figure 1 ). Therefore, the phonon gas theory is invalid in amorphous systems if one treats the branch at Gamma points as the principal gas particle. The possible reason is that the vibrational relaxation time defined in NMD is questionable because of the ill-definition of group velocity at Gamma points (see Sec. VI for detailed discussions). On the other hand, the vibrational relaxation time can be calculated by fitting the structure factor [Eq. (8)]. Figure 3a shows the vibrational spectrum at -1 | | 3 nm q calculated by the structure factor approach using the atomistic velocity at 800 K and the harmonic eigenvector. The corresponding results show there is a little temperature effect on the vibrational spectrum calculated by the SF method. As a consequence, the vibrational relaxation time via fitting the SF using the Lorentzian function (see Sec. IIIB for details) in amorphous Si has a weak temperature dependence as shown in Figure 3b . The computed relaxation time is from 0.01 ps to 1.3 ps, which agrees well with the results calculated by Larkin and McGaughey [1] using the same method. Our results here also show that the anharmonicity resulting from the temperature effect has essentially little effect on the vibrational relaxation times, and the broadening is dominated by the discontinuous of the structure of the amorphous systems (also see the analysis of the Sec. IVB for details). This phenomenon is also observed and proved by Moon et. al. [2] in their SF calculations. Thus, the corresponding thermal conductivity calculated using the phonon gas theory [Eq. (4)] has a weak temperature independence, which is consistent with the outputs of GKEMD. It means the phonon gas theory should be valid for the propagating modes in the amorphous systems once the vibrational group velocities or dispersions are well defined. B. Mean free paths of the propagating modes
In this section, we are studying the mean free path of the propagating mode (| |, ) q in amorphous silicon using the SF method ; (| |, ) (| |, ) (| |, ) Lor T L or T g Lor T q q v q , (20)where, Lor T is the longitudinal or transverse vibrational relaxation time, and ; g L or T v is the longitudinal or transverse group velocity which can be calculated via fitting the dispersions in the vibrational spectrums ( Figure 4 ). Figures 4a and show the normalized harmonic longitudinal and transverse SF spectrums, respectively. We fit the vibrational dispersions using the data below 12 THz which show clear propagating behaviors, and the vibrational modes with frequency larger than 12 THz are fitted by the extrapolations of the fitting lines (black dot lines in Figure 4a and ). The fitting longitudinal and transverse dispersions are q and q , respectively. We also plot the normalized SF spectrums at 800 K ( Figure 4b and ). The results clearly show that the effect of the anharmonicity due to the temperature on the vibrations’ scattering can be ignored in amorphous materials, which is in accordance with the conclusion of Moon et. al. [2]. With knowing the vibrational dispersions, the group velocity of the vibrations in amorphous systems can then be calculated using / | | g v d d q . The corresponding MFPs, without considering anharmonicity and at 300 K and 800 K, are plotted in
Figure 5a . Our results show the MFP in amorphous silicon is ranging from around 0.01 nm to several nanometers. The maximum calculated MFP in our simulations is around 5 nm and is lower than Moon’s reported value of 10 nm. The reason may because of the different fitting models, i.e., Lorenztian function in our calculations and harmonic oscillator in their case, and the difference of the initial structures. In their study, they use the constant group velocity to calculate the MFP, which may overestimate the mean free path at the high-frequency region as shown in our calculations (
Figure 5a ). We also notice that the MFPs computed at different temperatures are generally identical, which once again verifies that the anharmonicity has little influence on the vibrational scattering in amorphous materials as we discussed above. We then plot accumulative thermal conductivity function versus mean free path in
Figure 5b . The results show the vibrations with MFPs smaller than 3 nm contribute around 80% of the total thermal conductivity. Meanwhile, increasing the temperature of systems shows little effect on the accumulative thermal conductivity function. To further confirm that the vibrations in amorphous materials are mainly scattered by the discontinuities of the structure rather than the anharmonicity caused by the temperature, we record the decay process of a specific vibration mode, i.e., , in both amorphous and crystalline Si using trigger wave method (calculation details can be found in Ref. [2]) ( Figure 6 ). To perform these calculations, we create an amorphous model by repeating the 512-atom cell 100 times in one direction, and a corresponding crystalline model with the same size. Periodic boundary conditions are applied all three directions, and the system temperature is set as 0.1 K. Figure 6a shows the vibration will vibrate as a harmonic oscillator in the crystalline Si. The vibration in amorphous Si will be scattered quickly because of the discontinuities of the structure (
Figure 6b ), e.g., the amplitude of the vibration is only half of the original value at a position of 5 nm. If we identify the location at which the wave amplitude has decreased to e of its original value as the MFP of the corresponding vibration as suggested by Moon et. al. [2], the MFP from the trigger wave method is larger than that calculated using the SF method, i.e., 14 nm versus 5 nm. The reason may because we use a sinusoidal wave with amplitude 0.001 nm along the longitudinal direction to describe all the atoms in the first slab, which is not following the real situation, the amplitude of atoms is determined by its eigenvector. Nevertheless, our results in this section show anharmonicity caused by the temperature has little influence on the vibrational scattering process in amorphous materials. C. Temperature-dependent thermal conductivity
To calculate the total thermal conductivity, we also need to calculate the thermal conductivity of the non-propagating modes, i.e., diffusons. Here, we use the AF theory (see Sec. IIC for details) to calculate the modal thermal conductivity of diffusons (
Figure 7a ). The modal thermal conductivities are temperature independent since we assume all the vibrations are excited here, i.e., ( ) / diffuson b c k V . The modal thermal conductivities at low frequency show a scaling, which agrees with Larkin and MacGauhey’s calculations [1]. We then calculate the temperature-dependent thermal conductivity of the diffusons or propagons by correcting the excitation of the vibrations using / 2( ) sinh( / 2 ) bdiffuson or propagon V b b k TC k k T , (21)in which, is the reduced Planck constant, T is the concerning temperature and sinh is the hyperbolic sine function. Figure 7b shows the temperature-dependent thermal conductivities of propagons, diffusons and their resulting summations. By assuming the vibrations’ excited states following the distribution of Bose-Einstein function [Eq. (21)], which should be the situation for the bosons, the thermal conductivity at low temperatures should much lower those values at the high-temperature region because only a part of the vibrations are excited at low temperatures. With the increase of the temperature of the system, more vibrations will be excited, and the resulting thermal conductivity will saturate to a constant value as shown in
Figure 7b . Because all the vibrations are excited in MD simulations, the thermal conductivity calculated by GKEMD is temperature independent (black dots in
Figure 7b ), and the thermal conductivity is found to be identical to the summation of the saturated thermal conductivities contributed from propgagons and diffusons in which all vibrations are excited. We also plot the measured temperature-dependent thermal conductivities of amorphous Si adopted from Ref. [26] in
Figure 7b . The experimental measurement shows the same trend with our corrected thermal conductivities, which indicates only the correction of excited states, i.e., the heat capacity V C , is needed to obtain the correct thermal conductivity in amorphous materials using MD simulations. What is also worth noting that, our calculated thermal conductivities are generally lower than the experimental measurements, which should be due to the choice of the potential and the size of the system we used to do our calculations. V. DISCUSSIONS
Before closing, we discuss the possible reasons for why the NMD method is invalid to calculate the relaxation time of vibrations in amorphous materials. We start from the expression of the harmonic heat current , ,, ,
1( ) ( ) ( ) ( )21 1 1( ; ) ( ; ) ( ; ) ( ; )2 x x harmonic jj j j jjj jxjj jjj j j j Q t v t u t r tr e j X t e j X tm m , (22)in which, we derivate the finial harmonic heat current expression using Eq. (12). When , Eq. (22) is the heat current carried by the diffusons as proposed by Allen and Feldman [17]. The heat current resulting from the propagating modes are the harmonic heat current contributed by vibrational modes in terms of , we then have * *, , x xharmonic jj ijj j j j Q t e j e j r X t X tm m , (23)Based on the definition of dynamic matrix, we know ( ) ( ; ) (0) ( ; ) jjj j e j D e j , (24)and ( )( ) ( )2 ( ) ( ; ) ( ; )( )( ; ) ( ) ( ; )1( ; ) (0) ( ; ) jjx x xj jqq qijjj xj j qxjj jjj j j j D qe j e jq q qq ri e j D q e jqi e j r e jm m , (25)Combing Eq. (23) and (25), we have * 0 0 ( ) ( )( ; ) ( ) ( ; ) ( ; ) ( ) xharmonic x xq q Q t i X t X t Eq q , (26)where * ( ) ( ) ( ; ) ( ; ) E i X t X t is the energy of vibrational mode . Combing Eq. (2) and (26), we can know that
1( ; ) ( ; ) ( ; 0)1 ( ) ( ; ) ( ; 0)( ) ( ) ( ; ) ( ; 0)( ; 0) ( ; 0)( ) ( ; ) ( ; 0( ) xx x xharmonic harmonicb xb qb xb qV x q t Q t Q dtk VT E t E dtk VT qk T E t E dtk VT q E EE t EC q )( ; 0) ( ; 0) dtE E , (27)From Eq. (27), it clearly shows that the vibrational relaxation time in amorphous materials calculated using NMD, i.e., ( ; ) ( ; 0) / ( ; 0) ( ; 0) E t E E E dt , will be the true vibrational relaxation time only if ;0 ( ) / x xgq q v . However, in amorphous materials, the definition of vibrational group velocity at the Gamma points is invalid since the vibrational dispersion is not existing, which is also proved by Lv and Herny [25]. Therefore, the vibrational relaxation time in amorphous materials computed using NMD is questionable, and cannot capture the behaviors of the corresponding vibrations, i.e., the temperature independence of vibrational thermal conductivity. VI. CONCLUSIONS
In conclusion, we investigated the relaxation times of the propagating vibrational modes om amorphous silicon using both the normal mode decomposition method and the structure factor approach. We found the relaxation time calculated using the normal-mode decomposition is strongly dependent on anharmonicity caused by the temperature of the systems. The resulting total thermal conductivity then should be also strongly related to the temperature when all the vibrations are excited since the thermal conductivity contributed by the non-propagating modes are harmonic, which is contradictory to the investigations using equilibrium molecular dynamics simulations. We owe this to the ill-definition of group velocity at Gamma point in the normal mode decomposition method. On the opposite, the vibrational relaxation time calculated using both dynamical and static structure factor is almost the same, which indicates that it is suitable to calculate the vibrational dispersions using the effective wave vectors. The trigger wave method also shows that the propagating vibrational modes in amorphous materials are mainly scattered by the discontinuous of the structure rather than anharmonicity caused by temperature as in the perfect crystals. We then studied the modal thermal conductivity contributed from both the propagating vibrational modes and the non-propagating vibrational modes using structure factors and Allen-Feldman theory assumed all the vibrations are excited, i.e., the case in molecular dynamics simulations. The propagating modes with a mean free path smaller than 3 nm contribute to around 80% of the total thermal conductivity resulting from the propagons, and the thermal conductivity contributed from the propagons around is around 0.7 W/mK. The Allen-Feldman shows the thermal conductivity stemming from the diffusons is around 0.7 W/mK when all the vibrations are excited. Therefore, the total thermal conductivity contributed by propagons and diffusons is around 1.4 W/mK, which agrees quite well with the results computed using GKEMD. By correcting the excitations of the vibrations, i.e., calculating the volumetric heat capacity using the Bose-Einstein distribution function, the experimentally observed temperature dependence of the thermal conductivity in amorphous materials can be reproduced. Our results also imply that the inaccuracy of vibrational scattering rates caused by the distribution functions for the amorphous materials in molecular dynamics simulations, i.e., Boltzmann distributions in molecular dynamics simulations versus Bose-Einstein distributions for the bosons, can be ignored. This work provides a systematic analysis of the long-standing problem of thermal transport in disorder systems. Acknowledgments
Y.Z. thanks startup fund from The Hong Kong University of Science and Technology (HKUST). References [1] J. M. Larkin and A. J. H. McGaughey, Phys. Rev. B , 144303 (2014). [2] J. Moon, B. Latour, and A. J. Minnich, Phys. Rev. B , 024201 (2018). [3] W. Lv and A. Henry, New J. Phys. , 013028 (2016). [4] D. G. Cahill, S. K. Watson, and R. O. Pohl, Phys. Rev. B , 6131 (1992). [5] Y. Pan, J. Zhou, and G. Chen, Phys. Rev. B , 144203 (2020). [6] K. T. Regner, D. P. Sellan, Z. Su, C. H. Amon, A. J. H. McGaughey, and J. A. Malen, Nat. Commun. , 1640 (2013). [7] Y. He, D. Donadio, and G. Galli, Appl. Phys. Lett. , 144101 (2011). [8] Y. Zhou and M. Hu, Nano Lett. , 6178 (2016). [9] Y. Zhou, X. Zhang, and M. Hu, Phys. Rev. B , 195205 (2015). [10] Y. Zhou and S. Volz, arXiv:2006.02130 (2020). [11] S. Bourderau, T. Brousse, and D. M. Schleich, J. Power Sources , 233 (1999). [12] H. Jung, M. Park, Y. G. Yoon, G. B. Kim, and S. K. Joo, J. Power Sources , 346 (2003). [13] J. Sakabe, N. Ohta, T. Ohnishi, K. Mitsuishi, and K. Takada, Commun. Chem. , 24 (2018). [14] F. Fan, S. Huang, H. Yang, M. Raju, D. Datta, V. B. Shenoy, A. C. T. Van Duin, S. Zhang, and T. Zhu, Model. Simul. Mater. Sci. Eng. , 074002 (2013). [15] S. Patinet, D. Vandembroucq, and M. L. Falk, Phys. Rev. Lett. , 045501 (2016). [16] C. Kittel, Introduction to Solid State Physics (Wiley, New York, 1976). [17] P. B. Allen and J. L. Feldman, Phys. Rev. B , 12581 (1993). [18] H. R. Seyf and A. Henry, J. Appl. Phys. , 025101 (2016). [19] Y. Zhou, S. Xiong, X. Zhang, S. Volz, and M. Hu, Nat. Commun. , 4712 (2018). [20] D. M. Leitner, Phys. Rev. B , 094201 (2001). [21] J. M. Larkin and A. J. H. McGaughey, J. Appl. Phys. , 023507 (2013). [22] J. Moon, R. P. Hermann, M. E. Manley, A. Alatas, A. H. Said, and A. J. Minnich, Phys. Rev. Mater. , 065601 (2019). [23] H. Shintani and H. Tanaka, Nat. Mater. , 870 (2008). [24] Y. M. Beltukov, C. Fusco, D. A. Parshin, and A. Tanguy, Phys. Rev. E , 023006 (2016). [25] W. Lv and A. Henry, Sci. Rep. , 37675 (2016). [26] D. G. Cahill, M. Katiyar, and J. R. Abelson, Phys. Rev. B , 6077 (1994). [27] J. Tersoff, Phys. Rev. B , 5566 (1989). [28] R. J. Hardy, Phys. Rev. , 168 (1963). [29] R. Kubo, J. Phys. Soc. Japan , 570 (1957). [30] J. Ziman, Electrons and phonons: the theory of transport phenomena in solids (Oxford University Press, Oxford, 2001). [31] T. Iwashita, D. M. Nicholson, and T. Egami, Phys. Rev. Lett. , 205504 (2013). [32] T. Feng, X. Yang, and X. Ruan, J. Appl. Phys. , 145101 (2018). [33] P. Sheng and M. Zhou, Science , 539 (1991). Figure 1 . (a) the thermal conductivity versus the autocorrelation time for amorphous silicon with 512 atoms at 30 K and 50 K. (b) Thermal conductivity versus the autocorrelation time for amorphous silicon with 512 to 4096 atoms at 100 K, 300 K and 600 K. Figure 2 . (a) The normalized spectral energy density (SED) calculated using normal mode decomposition method at 50 K, 300 K and 800 K for frequency at 1.77 THz, the data are fitted using Lorentzian function. (b) The frequency-dependent vibrational relaxation time calculated using normal mode decomposition method at 50 K, 300 K and 800 K. Figure 3 . (a) The normalized dynamical and harmonic structure factor calculated using Eqs. (10), (11), (13) and (14). The data are fitted by the Lorentzian functions. (b) The frequency-dependent vibrational relaxation time calculated using structure factor calculation at static state, 300 K and 800 K. The results show the anharmonicity on the vibrational relaxation time caused by the temperature can be ignored.8
Figure 4 . The longitudinal harmonic (a) and dynamical (b) structure factor spectrum, and the transverse harmonic (c) and dynamical (d) structure factor spectrum. The dynamical structure factor are calculated at 800 K. The black dot lines are the fitting vibrational dispersions. Figure 5 . (a) The mean free path calculated using harmonic and dynamical structure factor method. (b) The corresponding accumulative thermal conductivity functions concerning the mean free path. Figure 6 . (a) The wave decay process of a specific vibration mode, i.e., , in (a) crystalline silicon and (b) amorphous silicon using the trigger wave method. 1 Figure 7 . (a) The modal thermal conductivity of diffusons. (b) The temperature dependent thermal conductivity of propagons, diffusons and their summation. The black dots are calculated using Green-Kubo equilibrium molecular dynamics simulations, and the experimental data (blue open dots) are adopted from Ref. [26].. (a) The modal thermal conductivity of diffusons. (b) The temperature dependent thermal conductivity of propagons, diffusons and their summation. The black dots are calculated using Green-Kubo equilibrium molecular dynamics simulations, and the experimental data (blue open dots) are adopted from Ref. [26].