Atomistic Simulations of Spectral Phonon Properties of Solids

Achieving sustainable energy and effective thermal management has been one of the greatest challenges facing the society. One way to achieve sustainable energy is by using thermoelectrics, which requires materials of low lattice thermal conductivities to achieve a high figure of merit (ZT) or efficiency. On the other hand, effective thermal management, which requires materials of high thermal conductivities to aid the heat dissipation, has become increasingly important to ensure the functionality and reliability of electronics, as the key components shrink down to micro-/nano-scales. Regardless of the applications, however, a fundamental understanding of the thermal transport mechanisms at the micro-/nano-scales is a prerequisite for the relevant analysis or design. We seek to develop new methodologies to predict thermal transport properties, and use these methods to gain such an understanding through atomistic simulations. We have predicted the thermal conductivities of freestanding silicene and silicene supported on amorphous SiO2 substrates in the temperature range from 300 to 900 K using nonequilibrium molecular dynamics (NEMD) simulations. The thermal conductivities of freestanding and supported silicene both decrease with increasing temperature, but that of the supported silicene has a weaker temperature dependence. The presence of a SiO2 substrate results in a large reduction to the thermal conductivity of silicene, up to 78.3% at 300 K, which is mainly due to the reduction in the phonon relaxation times. We also find that 97% of thermal conductivities of freestanding and supported silicene are contributed by phonons with mean free paths (MFPs) less than 100 and 20 nm, respectively. We have examined the geometry and discretization of the first Brillouin zone of face-centered cubic (FCC) crystals, and proposed a way of reducing the total number k points by using the symmetry k points. A systematic procedure is presented to determine the coordinates and degeneracy of the symmetry k points. By using the symmetry k points, the number of k points that need to be explicitly resolved could be reduced by as much as 97.92%. We have also proposed a formula for calculating the lattice thermal conductivity of FCC crystals by using the symmetry k points, which is validated by calculating the thermal conductivity of solid argon in the temperature range from 10 to 80 K. We have conducted equilibrium molecular dynamics (EMD) and NEMD simulations as well as phonon normal mode analysis to study the domain size effects of the thermal conductivities of silicon at 1000 K, graphene at 300 K, and silicene at 300 K. The EMD-predicted thermal conductivity of silicon shows a normal domain size effect, which could be attributed to the dominating effect of more phonon modes contributing to the thermal conductivity over the effect of more phonon scattering processes, as the domain size increases. The EMD-predicted thermal conductivities of graphene and silicene show an anomalous domain size effect, which could be due to the dominating effect of more phonon scattering processes (particularly those associated with the low-frequency flexural phonons) over the effect of more phonon modes contributing to the thermal conductivity, as the domain size increases. The NEMD-predicted thermal conductivities of the three material systems all show a normal domain size effect, which agrees well with the phonon MFP accumulation profile. We have studied the uncertainty quantification of EMD-predicted thermal conductivities by using solid argon, silicon, and germanium as model materials systems. The uncertainty increases with the upper limit of the correlation time, tcorre, UL, and decreases with the total simulation time, ttotal, whereas the velocity initialization seed, simulation domain size, temperature, and type of material have minimal effects. The uncertainty can be quantified by using a “universal” square-root relation, as σkx/kx,ave = 2(ttotal/tcorre, UL) –0.5. With this relation, it is possible to predict the uncertainty of the thermal conductivities from EMD simulations based on the chosen simulation parameters, even before the simulations are done. From statistical analysis, we have derived a formula that correlates the relative error bound (Q), confidence level (P), t corre, UL, ttotal, and number of independent simulations (N). We have provided recommendations on choosing tcorre, UL, ttotal, and N to achieve a desired relative error bound and confidence level. We have evaluated the performance of the Rescaling, Berendsen, Langevin, and Nosé-Hoover thermostats in controlling temperatures, reproducing the MB speed distributions, and predicting thermal conductivities, by conducting NEMD simulations of solid argon. The Langevin thermostat provides the best temperature and exchanged energy profiles, but the four thermostats provide similar results at a small thermostat size (e.g., a few lattice constants). When the thermostat time constant is small ( e.g.,, 4 fs), all the four thermostats are able to maintain temperatures to the desired values. When the thermostat time constant is large ( e.g.,, 40 ps), the Langevin and Nosé-Hoover thermostats provide slightly better temperature profiles, but none of the four thermostats is able to maintain temperatures to the desired values. The four thermostats show similar performance in reproducing the Maxwell-Boltzmann speed distribution and predicting the thermal conductivity of solid argon at 60 K. For large thermostated regions, the Langevin thermostat is recommended because it provides the best temperature and exchanged energy profiles. Otherwise, any of the four thermostats could be used in NEMD simulations. The results from this study offer new understandings of the physics of thermal transport at the micro-/nano-scales and of the numerical methods used in the relevant simulations. They are potentially useful for similar studies in this area and ultimately for the design of better thermoelectrics and thermal management solutions