Interpretation of apparent thermal conductivity in finite systems from equilibrium molecular dynamics simulations

We propose a way to properly interpret the apparent thermal conductivity obtained for finite systems using equilibrium molecular dynamics simulations (EMD) with fixed or open boundary conditions in the transport direction. In such systems the heat current autocorrelation function develops negative values after a correlation time which is proportional to the length of the simulation cell in the transport direction. Accordingly, the running thermal conductivity develops a maximum value at the same correlation time and eventually decays to zero. By comparing EMD with nonequilibrium molecular dynamics (NEMD) simulations, we conclude that the maximum thermal conductivity from EMD in a system with domain length $2L$ is equal to the thermal conductivity from NEMD in a system with domain length $L$. This facilitates the use of nonperiodic-boundary EMD for thermal transport in finite samples in close correspondence to NEMD.

Figure 1

Schematic illustration of the setups for the different MD simulations: (a)–(c) EMD with open boundary conditions in the transport direction, EMD with fixed boundary conditions in the transport direction ([100] direction), and NEMD with fixed boundary conditions in the transport direction, respectively, for 3D silicon. (d)–(f) similarly for quasi-1D $(10,10)$ CNT, with the transport being along the tube. (g)–(i) similarly for 2D graphene, with the transport being along the zigzag direction.

Figure 2

(a) Steady-state temperature profile and (b) energies of the thermostats coupled to the heat source and sink regions as a function of the time. The heat transfer rate $dE/dt$ is calculated as the average of the absolute slopes in the two curves. The energy in the thermostat coupled to the heat source (sink) region is decreasing (increasing) because it releases (absorbs) energy to maintain the higher (lower) temperature in the thermostatted region.

Figure 3

(a) Normalized HCACF and (b) running thermal conductivity $\kappa \left(\tau \right)$ as a function of the correlation time for different sample lengths $L$ from EMD simulations with open boundary conditions in the transport direction. For the running thermal conductivity, the red thick line represents the mean values from ten independent runs (the gray thin lines) and the black dashed lines indicate the standard error. For the HCACF, only the average values are shown for clarity. The blue dashed vertical lines correspond to the time ${\tau }_{\max }$ at which $\kappa \left(\tau \right)$ reaches a maximum, and equivalently after which the HCACF develops negative values. The systems here are 3D silicon crystals at 300 K and zero pressure, with the length $L$ for each sample written in the corresponding panel. Note that the positive and negative areas enclosed by the HCACF and the time axis, which are equal, appear to be unequal due to the logarithmic scale of the time axis.

Figure 4

Similar to Fig. 3, but for EMD simulations with fixed boundary conditions, instead of open boundary conditions, in the transport direction.

Figure 5

(a) The maximum thermal conductivity ${\kappa }_{\max }$ from EMD simulations with open and fixed boundary conditions in the transport direction and the apparent thermal conductivity from NEMD simulations as a function of the domain length $L$ for 3D crystal silicon at 300 K and zero pressure. (b) Similar to (a) but using $L/2$ as the horizontal axis for the EMD data.

Figure 6

Normalized HCACF in silicon crystal (300 K and zero pressure) in different conditions. The blue dashed line represents the HCACF obtained from EMD simulations with periodic boundary conditions in the transport direction. The red dot-dashed line represents the HCACF obtained from EMD simulations with open boundary conditions in the transport direction. The black solid line represents the difference between the above two, which is induced by phonon-boundary scattering.

Figure 7

The correlation time ${\tau }_{\max }$ at which the running thermal conductivity $\kappa \left(\tau \right)$ attains the maximum value ${\kappa }_{\max }$ as shown in Fig. 3 and Fig. 4 against the EMD simulation domain length $L$.

Figure 8

Validation of Eq. (4) for (a) $(10,10)$ CNT with different lengths, (b) graphene sheet with different lengths, and (c) $(10,10)$ CNT with fixed lengths but at different temperatures $T$.