Examining the Validity of the Phonon Gas Model in Amorphous Materials

The idea of treating phonon transport as equivalent to transport through a gas of particles is termed the phonon gas model (PGM), and it has been used almost ubiquitously to try and understand heat conduction in all solids. However, most of the modes in disordered materials do not propagate and thus may contribute to heat conduction in a fundamentally different way than is described by the PGM. From a practical perspective, the problem with trying to apply the PGM to amorphous materials is the fact that one cannot rigorously define the phonon velocities for non-propagating modes, since there is no periodicity. Here, we tested the validity of the PGM for amorphous materials by assuming the PGM is applicable, and then, using a combination of lattice dynamics, molecular dynamics (MD) and experimental thermal conductivity data, we back-calculated the phonon velocities for the vibrational modes. The results of this approach show that if the PGM was valid, a large number of the mid and high frequency modes would have to have either imaginary or extremely high velocities to reproduce the experimental thermal conductivity data. Furthermore, the results of MD based relaxation time calculations suggest that in amorphous materials there is little, if any, connection between relaxation times and thermal conductivity. This then strongly suggests that the PGM is inapplicable to amorphous solids.


Simulation details of generating a-SiO 2 structure
The a-SiO 2 initial structure was generated by the melting-quenching method using Tersoff potential with parameters published by Munetoh et al. 2 . The detailed procedures of melting-quenching method have been described by Ong et al. 5 After quenching, the structure was annealed at 1100 K for 10 ns to avoid the meta-stability reported by Larkin et al. 4 . After generating the a-SiO 2 structure, we calculated the density of states (DOS), and compared with experiments 6 . The agreement is overall reasonable. The supercell is approximately cubic has 4608 atoms, with a length of ~ 40.44 A. The density for the relaxed structure is 2317 Kg/m 3 , which is 4% larger than the experimental value 2220 kg/m 37 . The time step used in the simulation is 0.1 fs.

Lattice dynamics simulation and GKMA
After we obtained the initial atomistic structures, we applied lattice dynamics (LD) at gamma point ( k = 0 ) for the supercell with periodic boundary conditions to obtain the normal mode eigen values and eigen vectors which allow one to visualize the normal mode shapes. The LD calculations were performed using the General Utility Lattice Program (GULP) 8 . Before LD calculation, we relaxed structure at 0K with zero pressure. Then one obtained the eigenvectors and harmonic frequencies.
Once all the eigen-vectors have been calculated, they are read into the molecular dynamics (MD) simulation to calculate the mode level thermal conductivity contributions from GKMA. All MD simulations were performed using the Large Atomic/Molecular Massively Parallel Simulator (LAMMPS). After equilibrated for 100 ps at 300K using NVT (constant number of atoms, volume, temperature), the heat flux, and mode heat flux were captured for another 4 ns (4 x 10 7 time steps) using NVE (constant number of atoms, constant volume, constant energy). After supplying the eigenvectors once at the beginning of the MD simulation, one is able to obtain the heat flux and kinetic energy of each mode. The integral of the heat flux autocorrelation function is cut off at 30 ps 9 since the largest relaxation time in the supercell is less than 10 ps. After the modal heat flux is computed, the modal thermal conductivity is determined by calculating the correlation between modal heat flux and the total heat flux.

GKMA formulations
Here, we recapitulate the key points of the GKMA formalism. Firstly, the normal mode eigenvectors are computed from a LD calculation. Secondly, one needs to project the atomic velocities from the MD simulation of the same supercell onto the normal mode eigenvectors. Then one can obtain the time history of each normal mode's amplitude. Thirdly, each atom's instantaneous velocity can then be decomposed into individual mode contributions based on the respective instantaneous normal mode amplitudes, whereby summing the modal contributions returns each atom's velocity. Fourthly, one substitutes the modal components of each atom's velocity into the heat flux operator 10 to obtain each mode's instantaneous contribution to the heat flux. The total heat flux can then be obtained from the sum of all individual mode contributions to the heat flux, via where n is mode index, N is total number of atoms in the super cell, V is volume of the super cell, E i is the kinetic and potential energy of atom i , Φ j denotes potential energy of atom j , !
x i (n,t) is the contribution mode n makes to the velocity of atom i and r ij is distance between atom i and j . After having access to the individual mode heat fluxes, we can substitute the summation over modes in Eq. (1) directly into the Green-Kubo expression for TC, which calculates the TC as proportional to the heat flux autocorrelation function. One then obtains the TC as a direct summation over individual mode contributions, where κ (n) stands for TC contribution of mode n , k B is Boltzmann constant, T is the temperature and V is volume. Using Eq. (2) one can calculate the TC of individual modes in any material where the atoms vibrate around stable equilibrium sites. Classical MD has been considered to be inaccurate at low temperatures (i.e., below a material's Debye temperature) because it does not reproduce the proper mode amplitudes that correspond to the quantum occupations. As a result, classical MD results in a constant heat capacity with respect to temperature, since every mode is equally excited at all temperatures. However, once each individual mode's TC is obtained, one can easily apply a quantum specific heat correction, which allows one to extend the MD based predictions to essentially any temperature.

Quantum specific heat correction
The quantum expression of volumetric specific heat, based on Bose-Einstein statistics is given by 14 and the classical volumetric specific heat is given by c m (ω ) = k B V . Thus, the quantum heat capacity correction factor is the ratio, Calculation of phonon group velocities using PGM The relaxation time τ (n) is calculated using NMA method with Tersoff potential on a-SiO 2 at six different temperatures (30K, 100K, 100K, 200K, 400K, 800K). Then one could easily calculate relaxation time τ (n,T ) ( 30K ≤ T ≤ 800K ) by interpolating between different temperatures. The specific heat from Bose-Einstein statistics 14 is also temperature dependent, as c(n, is independent of temperature. The measured thermal conductivity of a-SiO 2 is available from 30K to 750K [2]. The first step is to get temperature dependent thermal conductivity as a continuous function κ (T ) by fitting from the experiments result. Then start from 30K, we could calculate the square of phonon group velocity ν g (n) 2 from ν g (n) 2 = κ (T ) / c(T ,n)τ (T ,n) ⎡ ⎣ ⎤ ⎦ . At 30K, there are only a few modes are excited, in other words, most of modes c(n,30) is negligible. We choose c(n,T ) < 0.001 as a criteria to set mode n is excited or not at temperature T. After we calculated the phonon velocities of a few modes that are excited at 30K, we increased temperature to be slightly higher (30.072K). Since there are total 13824 modes in the system, we divided temperature into 10000 bins so that each temperature increase only excites a small number of modes. When the temperature changes, the relaxation time and specific heat are changing for all modes.
The contribution of thermal conductivity k(n,T ) n ′ n ∑ from the modes that phonon velocity has been calculated is changing with the relaxation time and specific heat. The difference between experimental value k(T ) and k(n,T )