New insights into thermal conductivity by non-equilibrium molecular dynamics

Non-equilibrium molecular dynamics (NEMD) is often used to simulate thermal conductivity ($\kappa$). A steady state heat current and corresponding temperature gradient are created computationally over a simulation cell of thousands of atoms. We advocate a variation that gives directly $\kappa(q)$, the Fourier transform of the non-local $\kappa(x-x^\prime)$ that relates $J(x)$ to $\nabla T(x^\prime)$. The algorithm is tested on the Lennard-Jones liquid and crystal, and is efficient for extraction of the macroscopic $\kappa={\rm Lim}_{q\rightarrow0}\kappa(q)$. Peierls-Boltzmann theory gives (in relaxation-time approximation) a closed-form expression for $\kappa(q)$ that can be used to study the $q$-dependence in the small $q$ limit, and how it depends on simulation cell dimensions in NEMD. The frequency-dependent relaxation rate $1/\tau_Q \propto \omega_Q^2$ was chosen for detailed comparison with simulation. For an isotropic cell ($N_x=N_y=N_z$), the behavior is $\kappa(q)=\kappa-A*q^{1/2}$. For the more typical anisotropic cell with one length ($N_z$) large compared to the others, there is an additional term $\propto q^{-1/2}/N_xN_y$. This divergent contribution disappears in the bulk limit. Strategies for extrapolation of simulations are suggested.