Equilibrium Molecular Dynamics Study of Phonon Transport Properties at High Temperatures

Abstract

Thermal conductivity is a material property of great importance in many areas of physics including coat barriers, renewable energy, nano-electronic devices, and geological processes inside the Earth. At the microscopic level, two type of heat carriers, electronic and vibrational, regulate this macroscopic heat transport property. The lattice thermal conductivity, from lattice vibrations or phonons, dominates in semiconductors or insulators at moderate temperatures while the electronic contribution is non-negligible at ambient temperatures. While experiments are only able to measure the total conductivity, extrapolation equations are given to estimate the lattice component. Traditionally, perturbation theory predicts a temperature (T) dependence of the form 1/T at high temperatures within the Peierl’s heat flux approximation and three-phonon scattering mechanisms yet the experimental extrapolations do not have the same predictions. Even for simple crystals, such as Si, the temperature dependence of the lattice conductivity remain speculative near the melting point (1680K for Si) Alternatively, Molecular Dynamic (MD) simulations can be used to get a non-perturbative modelling of anharmonic lattice dynamics, but calculations of the thermal conductivity with equilibrium and non-equilibrium MD is susceptible to accuracy issues, such as finite-size effect, reliability of the inter-atomic potential, and requirements of large number of ensembles. In this work, we focus on the high temperature regime of diamond-lattice Silicon through Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) with the highly accurate Gaussian Approximation Potential (GAP) machine-learning inter-atomic potential. First, we calculate the thermal conductivity with renormalized phonon frequency and lifetime given the q-projected velocity Auto-Correlation Function (ACF). The damp harmonic oscillator theory is validated at high temperature and the q-space symmetry is used to reduce the necessary large number of ensembles by averaging over equivalent q-points, with a mean equivalence number of 24. A fitting algorithm, utilizing LMFIT, is developed to fix the sensitivity of the phonon lifetime parameter by implementing an ”accidental degeneracy” constraint if non-degenerate mode frequencies are too close relative to the MD resolution. When compared to similar perturbation calculations, we find that the simulation parameters of 100 ps simulation time may induce sampling issues of the acoustic mode with long lifetimes that leads to severe underestimation of the thermal conductivity. Then, we calculate the thermal conductivity through the Green-Kubo formalism with different orders of the Heat Flux Auto-Correlation Function (HFACF) based on the Hardy heat flux operator. Applying an orthogonality rule to the second-order HFACF mitigates convergence issues when integrating the correlation function, enabling us to confirm the underestimation of phonon lifetimes. A unified transport theory, based on the second-order HFACF, is discussed that combines the inter-mode, intra-mode HFACF and the cross-correlation term, where the intra-term corresponds to the Peierl’s heat flux discussed above. Our MD simulations independently confirmed that the conductivity, from the inter-mode HFACF, is comparable to Wigner transport equation results at room temperature and saturates to 0.32 W/(m-K) at 1500K. In contrast, the conductivity, from the cross-correlation term, lowers the overall prediction by about -1 W/(m-K) at all temperatures and this term is neglected by the Wigner transport equation, as well as all the recently proposed unified theory. Finally, we perform at q-projected force fitting to get the temperature-dependent higher-order Force Constant (FC) matrices, namely the third- and fourth- order anharmonic and use them for the higher-order HFACF. The effective fourth-order FC matrix with a selection rule ∆(−⃗q + ⃗q2) provided no improvements to q-projected force fitting when used in conjunction with the bare or renormalized FC2 and a different approach should be used to handle the fourth- order effects. We demonstrate that a similar orthogonality rule to the one applied to the second-order HFACF enhances convergence by integration of the third-order HFACF. Only accounting for the third-order HFACF lowers the overall conductivity prediction and an improved method of the fourth-order FC is needed to clarify the importance of the fourth-order anharmonicity at high temperatures.