Computation of thermal conductivity based on Path Integral Monte Carlo methods

Abstract

The calculation of thermal conductivity in insulating solids at temperatures below the Debye temperature is problematic, due to the breakdown of classical and semi-classical approaches. In this work, we present a fully quantum methodology to compute thermal conductivity based on Path Integral Monte Carlo (PIMC) simulations combined with the Green-Kubo linear response theory. The method is applied to crystalline argon modeled by a Lennard-Jones potential, a paradigmatic system where quantum effects strongly affect both thermodynamic and transport properties. From PIMC simulations, we obtain the temperature-dependent phonon frequencies, lifetimes, and specific heat. From the imaginary time correlations of the energy current, we extract the thermal transport coefficients based on a physically motivated prior. We show that the experimentally observed increase of the thermal conductivity at low temperatures cannot be explained within a standard Peierls-Boltzmann framework or quasi-harmonic approximation using phonon lifetimes alone. Instead, a distinct transport lifetime emerges from the analysis of heat-current correlations. Our results demonstrate that quantum Monte Carlo methods provide a robust, non-perturbative framework to investigate heat transport in insulating solids, beyond the limits of classical molecular dynamics and quasi-harmonic approximations.

Figures (5)

• Figure 1: Phonon dispersion relationship, $\omega_{\mathbf{k}a}$ for solid argon, along the indicated directions of the wave-vector, $\mathbf{k}$. The results for longitudinal ($a=$ L) and transverse ($a=$ T) phonons, $\omega_{\mathbf{k}a}^{ph}$, (red circles) obtained from the normal modes correlation function for a simulated system with $N=$ 864 atoms are compared with experimental data Fujii1974. The grey dotted lines indicate the harmonic frequencies, $\omega_{\mathbf{k}a}^0$, pertaining to the ideal harmonic crystal. The vertical (red) bars quantify the broadening of the phonon lines due to their finite lifetimes, $\tau_{\mathbf{k}a}^{ph}$. The experimental line-width Fujii1974 is also indicated as a vertical (green) bar, for selected data points.
• Figure 2: Specific heat of solid argon from PIMC calculations compared to the experimental results. The black squares are obtained from the numerical derivative of the total energy, whereas $C_V^h$ (triangles) is the specific heat determined in the harmonic approximation based on the phonon frequencies, $\omega_{\mathbf{k}a}^{ph}$, obtained from the normal mode correlation functions. The experimental data from Peterson1966 ($\times$) are also shown for comparison, together with the corresponding values predicted by the Debye model with $T_D=$ 84.4 K (dashed line).
• Figure 3: Main panel: Spectral function of the harmonic current correlations at $T=20$K. Inset: The corresponding imaginary time correlations $C_{xx}^h(\tau)$, obtained from PIMC (circles) together with the reconstructed one, $C_{xx}^h(\tau |\Lambda)$, based on the prior $\Lambda(\omega,\Gamma^{tr})$, determining the inverse transport lifetime $\Gamma^{tr}$ ( solid line). For comparison, we also show the best reconstruction obtained by imposing the inverse phonon lifetime, $C_{xx}^h(\tau|\Lambda(\omega,\Gamma^{ph}))$ (dashed line), incompatible with the calculated imaginary time correlation.
• Figure 4: Thermal conductivity of solid argon obtained from the spectral reconstruction of $C_{xx}(\tau)$, compared with the result of the Peierls-Boltzmann equation, $\kappa_{BP}$ of Eq. (\ref{['KappaPB']}), based on the PIMC phonon frequencies and lifetimes. The red circles are results from the NEMD calculations reported in JLB2014. The experimental data from christensen1975 and krupski1969multiphonon are indicated in black as $\times$ and $+$, respectively.
• Figure 5: Spectral function of total current current correlations at $T=20$K, based on the corresponding reconstruction $\Lambda(\omega,\Gamma^{tr})$ of the harmonic current correlations and the additional spectral weight $\widetilde{\Lambda}(\omega)$. Notice the difference in scale between the two components of the spectral density. The inset shows the reconstructed imaginary time correlation function of the total current obtained by optimizing only the additional weight $\widetilde{\Lambda}(\omega)$ without modifying $\Lambda(\omega,\Gamma^{tr})$.