Influence of Boundaries and Thermostatting on Nonequilibrium Molecular Dynamics Simulations of Heat Conduction in Solids

TL;DR

This paper examines how boundary conditions and thermostat choices affect the accuracy of nonequilibrium molecular dynamics simulations of heat conduction in solids, emphasizing the importance of proper temperature profile analysis and thermostat selection.

Contribution

It demonstrates that nonlinear temperature profiles should be included in thermal transport calculations and shows Langevin thermostats outperform Nosé-Hoover thermostats in nanoscale heat transport simulations.

Findings

Proper thermostat choice affects simulation accuracy.

Nonlinear temperature profiles are important in calculations.

Langevin thermostats reduce artifacts in simulations.

Abstract

Nonequilibrium molecular dynamics (NEMD) has been extensively used to study thermal transport at various length scales in many materials. In this method, two local thermostats at different temperatures are used to generate a nonequilibrium steady state with a constant heat flux. Conventionally, the thermal conductivity of a finite system is calculated as the ratio between the heat flux and the temperature gradient extracted from the linear part of the temperature profile away from the local thermostats. Here we show that, with a proper choice of the thermostat, the nonlinear part of the temperature profile should actually not be excluded in thermal transport calculations. We compare NEMD results against those from the atomistic Green's function method in the ballistic regime, and those from the homogeneous nonequilibrium molecular dynamics method in the ballistic-to-diffusive regime. These comparisons suggest that in all the transport regimes, one should directly calculate the thermal conductance from the temperature difference between the heat source and sink and, if needed, convert it to the thermal conductivity by multiplying it with the system length. Furthermore, we find that the Langevin thermostat outperforms the Nos\'{e}-Hoover (chain) thermostat in NEMD simulations because of its stochastic and local nature. We show that this is particularly important for studying asymmetric carbon-based nanostructures, for which the Nos\'{e}-Hoover thermostat can produce artifacts leading to unphysical thermal rectification. Our findings are important to obtain correct results from molecular dynamics simulations of nanoscale heat transport as the accuracy of the interatomic potentials is rapidly improving.

Equations

$$\kappa (L)=\frac{Q}{S|\nabla T|}.$$

$$G(L)=\frac{Q}{S\Delta T},$$

$$G(L)\equiv \frac{\kappa (L)}{L}.$$

$$\kappa (L)=\frac{Q}{S(\Delta T/L)}.$$

$$G={\int }_0^{\infty }\frac{d\omega }{2\pi }G(\omega ),$$

$$G(\omega )=G_{\mathrm{c}}(\omega )=\frac{k_{\mathrm{B}}}{S}\mathcal{T}(\omega )$$

$$G(\omega )=G_{\mathrm{q}}(\omega )=\frac{k_{\mathrm{B}}}{S}\frac{x^2{e}^x}{(e^x-1{)}^2}\mathcal{T}(\omega )$$

$$\mathcal{T}(\omega )=\mathrm{Tr}[\mathcal{G}(\omega ){\mathrm{\Gamma }}_{\mathrm{L}}(\omega ){\mathcal{G}}^{†}(\omega ){\mathrm{\Gamma }}_{\mathrm{R}}(\omega )],$$

$$\mathcal{G}(\omega )=\frac{1}{{\omega }^2-D-{\Sigma }_{\mathrm{L}}(\omega )-{\Sigma }_{\mathrm{R}}(\omega )}$$

$${\Gamma }_{L/R}(\omega )=i\left[{\Sigma }_{L/R}(\omega )-{\Sigma }_{L/R}^{†}(\omega )\right].$$

$$D_{\mu \nu }=\frac{1}{\sqrt{m_{\mu }{m}_{\nu }}}\frac{{\partial }^{2}U}{\partial u_{\mu }\partial u_{\nu }},$$

$$D_{ij}^{ab}=\frac{1}{\sqrt{m_i{m}_j}}\frac{{\partial }^{2}U}{\partial u_i^a\partial u_j^b},$$

$$D_{ij}^{ab}=-\frac{1}{\sqrt{m_i{m}_j}}\frac{\partial F_i^a}{\partial u_j^b}.$$

$$D_{ij}^{ab}\approx \frac{1}{\sqrt{m_i{m}_j}}\frac{F_i^a(r_j^b-\delta )-F_i^a(r_j^b+\delta )}{2\delta }.$$

$$\frac{d{\mathbf{r}}_i}{dt}=\frac{{\mathbf{p}}_i}{m_i};\frac{d{\mathbf{p}}_i}{dt}={\mathbf{F}}_i-\frac{{\pi }_0}{Q_0}{\mathbf{p}}_i,$$

$$\frac{d{\mathbf{r}}_i}{dt}=\frac{{\mathbf{p}}_i}{m_i};\frac{d{\mathbf{p}}_i}{dt}={\mathbf{F}}_i-\frac{{\mathbf{p}}_i}{\tau }+{\mathbf{f}}_i,$$

$${\mathbf{F}}_i^{\mathrm{ext}}=E_i{\mathbf{F}}_{\mathrm{e}}+\sum _{j\ne i}\left(\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}\otimes {\mathbf{r}}_{ij}\right)\cdot {\mathbf{F}}_{\mathrm{e}},$$

$$\kappa (t)=\frac{⟨J(t){⟩}_{\mathrm{ne}}}{TV{F}_{\mathrm{e}}},$$

$$\mathbf{J}=\sum _i{\mathbf{v}}_i{E}_i+\sum _i\sum _{j\ne i}{\mathbf{r}}_{ij}\left(\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}\cdot {\mathbf{v}}_i\right),$$

$$Q=-\sum _{i\in A}\sum _{j\in B}⟨\left(\frac{\partial U_i}{\partial {\mathbf{r}}_{ij}}\cdot {\mathbf{v}}_j-\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}\cdot {\mathbf{v}}_i\right)⟩.$$

$$K(t)=-\sum _{i\in A}\sum _{j\in B}⟨\left(\frac{\partial U_i}{\partial {\mathbf{r}}_{ij}}(0)\cdot {\mathbf{v}}_j(t)-\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}(0)\cdot {\mathbf{v}}_i(t)\right)⟩$$

$$\tilde{K}(\omega )={\int }_{-\infty }^{\infty }dt{e}^{i\omega t}K(t).$$

$$K(t)={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }{e}^{-i\omega t}\tilde{K}(\omega ).$$

$$Q=K(t=0)={\int }_0^{\infty }\frac{d\omega }{2\pi }\left[2\tilde{K}(\omega )\right],$$

$$G(\omega )=\frac{2\tilde{K}(\omega )}{S\Delta T}\text{with}G={\int }_0^{\infty }\frac{d\omega }{2\pi }G(\omega ).$$

$$\kappa (\omega )=\frac{2\tilde{K}(\omega )}{ST{F}_{\mathrm{e}}}\text{with}\kappa ={\int }_0^{\infty }\frac{d\omega }{2\pi }\kappa (\omega ).$$

$$k^{\prime }(L)=\frac{Q}{S(\Delta T^{\prime }/L)}.$$

$$G^{\prime }(L)=\frac{Q}{S\Delta T^{\prime }}.$$

$$\kappa (L)={\int }_0^{\infty }\frac{d\omega }{2\pi }\frac{\kappa (\omega )}{1+\lambda (\omega )/L}.$$

$$\eta =\frac{G_{\mathrm{H}}-G_{\mathrm{L}}}{G_{\mathrm{L}}}\times 100\%.$$

Full text

I Introduction

Molecular dynamics (MD) is the most versatile and complete classical method to study heat transport at the nanoscale, which is vital for many technological applications Cahill et al. (2014); Volz et al. (2016); Donadio (2018) such as thermoelectric energy conversion and thermal management of electronic devices. As the interatomic interactions used in MD simulations have become inreasingly accurate by using quantum mechanical density functional based equilibrium MD Marcolongo et al. (2015); Carbogno et al. (2017); Kang and Wang (2017), nonequilibrium MD Stackhouse et al. (2010) and approach-to-equilibrium MD Puligheddu et al. (2017); Martin et al. (2018), it is crucial to develop a deeper understanding of the MD methods used for heat transport studies. In this work, we focus on one of the most popular MD methods for heat transport: the nonequilibrium MD (NEMD) method.

In NEMD simulations Stackhouse and Stixrude (2010); Shiomi (2014), one usually calculates the thermal conductivity $\kappa(L)$ of a finite system with length $L$ form the (presumably constant) temperature gradient $\nabla T$ and the heat flux $Q/S$ determined from a steady state according to Fourier’s law:

[TABLE]

Here, $Q$ is the thermal power across a cross-section area $S$ perpendicular to the transport direction and $|\nabla T|$ is the magnitude of the temperature gradient. However, at the nanoscale transport is not diffusive and the conventional concept of conductivity as a materials property becomes invalid Datta (1995). For example, thermal transport at the nanoscale, especially in materials with high thermal conductivity such as graphene Balandin et al. (2008); Ghosh et al. (2008); Xu et al. (2014), is almost ballistic with a length-dependent $\kappa(L)$ increasing with increasing $L$. In this situation, a more appropriate quantity for describing heat transport is the thermal conductance per unit area $G(L)$, which is defined as

[TABLE]

where $\Delta T>0$ is the temperature difference between the heat source and the sink. This quantity is constant in the ballistic regime and only weakly dependent on the system length in the nanoscale transport regime. In a system with a uniform cross section, the length-dependent conductivity and the conductance are related by the following equation:

[TABLE]

A question then arises as to whether the conductivity and conductance calculated using Eqs. (1) and (2) are consistent with each other. Clearly, if the temperature gradient is replaced by $\Delta T/L$, $\kappa$ and $G$ as calculated from Eqs. (1) and (2) imply Eq. (3). However, in most previous works using NEMD simulations, the temperature gradient was not calculated as $\Delta T/L$, but was instead determined as the slope of the so-called linear region of the temperature profile, ignoring the nonlinear parts of the temperature profile near the thermal baths. This practice assumes that transport is diffusive, i.e. in accordance with Fourier’s law that predicts a linear temperature profile at steady state conditions. It has been usually argued that nonlinearities are artifacts related to e.g., the local thermostats in the NEMD simulation setup and should be excluded. However, recent studies suggest that such nonlinearities have a physical origin, related to transport in finite size systems Allen (2014); Cepellotti and Marzari (2017). If the nonlinear parts need to be excluded, then $|\nabla T|\neq\Delta T/L$, leading to an inconsistency among Eqs. (1)-(3).

In this work we take graphene, a material with an exceptionally high lattice thermal conductivity and long phonon mean free paths Balandin et al. (2008); Ghosh et al. (2008), as an example to explore the interpretation of NEMD results both in the ballistic and in the ballistic-to-diffusive transport regimes. In the ballistic limit, we compare NEMD (with spectral decomposition Sääskilahti et al. (2014, 2015); Fan et al. (2017a)) against the standard atomistic Green’s function (AGF) method Mingo and Yang (2003); Wang et al. (2008); Sadasivam et al. (2014). We show that in order to reach an agreement with the AGF method, one needs to use Eq. (2) to calculate the conductance and regard $\Delta T$ as the temperature difference between the hot and cold thermostats, not excluding any local nonlinear regions of the temperature profile. If one first calculates a conductivity using Eq. (1) from the linear region of the temperature profile and then converts it to the conductance using Eq. (3), one obtains a ballistic conductance much larger than the correct one.

On the other hand, although the NEMD method has been shown to be fully equivalent to the equilibrium MD (EMD) and the homogeneous nonequilibrium MD (HNEMD) methods in the diffusive limit Fan et al. (2017a); Dong et al. (2018a); Fan et al. (2019); Xu et al. (2018); Dong et al. (2018b); Xu et al. (2019), the influence of the simulation details in the NEMD method on the results in the ballistic-to-diffusive crossover regime have not been addressed. Here, we clarify the interpretation of the NEMD results by comparing them against those from the HNEMD method (also with spectral decomposition) Fan et al. (2019). Again, we show that one should not calculate the thermal conductivity using Eq. (1) from a purely linear region of the temperature profile. Instead, the correct way is to calculate the conductance from Eq. (2) with $\Delta T$ being the temperature difference between the heat source and sink and then convert it to the conductivity using Eq. (3). That is, one should calculate the conductivity as

[TABLE]

Finally, we address thermal rectification in asymmetric graphene nanostructures, showing that the choice of thermostatting method and related parameters is of crucial importance, not only to calculate $G$ in the ballistic regime, but also to estimate correctly the thermal rectification efficiency of these systems.

This paper is organized as follows. In Sec. II, we review the various numerical phonon thermal transport methods used in this work, including the AGF method based on harmonic force constants (Sec. II.1), the NEMD method (Sec. II.2), the HNEMD method (Sec. II.3), and the spectral decomposition method (Sec. II.4). After briefly presenting the AGF results in Sec. III.1, we compare the NEMD results against the AGF results in Sec. III.2. Then, we examine the temperature profiles in the NEMD simulations in Sec. III.3. In Sec. III.4, we present the NEMD results for the ballistic-to-diffusive transport regimes and compare them with the HNEMD results. Connection between the NEMD method and the Boltzmann transport equation method is discussed in Sec. III.5. Finally in Sec. III.6 we show the application of NEMD to thermal rectification in asymmetric graphene devices. In Sec. IV, we give a summary and the main conclusions of this work.

II Methods

II.1 Atomistic Green’s function method

Following the early work by Mingo and Yang Mingo and Yang (2003), the AGF method has become a standard tool to study ballistic phonon transport in the harmonic approximation. In this method, one can calculate the phonon transmission $\mathcal{T}(\omega)$ between two leads as a function of the phonon frequency $\omega$. The ballistic conductance can be calculated through the Landauer expression as Imry and Landauer (1999); Young and Maris (1989):

[TABLE]

where the spectral conductance $G(\omega)$ is

[TABLE]

using classical statistics, and

[TABLE]

using quantum statistics. Here, $x=\hbar\omega/k_{\rm B}T$, where $\hbar$ is Planck’s constant, $k_{\rm B}$ is Boltzmann’s constant, and $T$ is the system temperature.

As in the case of electron transport Datta (1995), there are many equivalent representations of the phonon transmission Wang et al. (2008). Here, we adopt Caroli’s formula Caroli et al. (1971):

[TABLE]

where

[TABLE]

is the retarded Green’s function for the system, $\mathcal{G}^{\dagger}(\omega)$ is the advanced Green’s function, $D$ is the dynamical matrix of the system, and $\Sigma_{\rm L/R}(\omega)$ is the self energy matrix of the left (right) lead. The coupling matrices $\Gamma_{\rm L}(\omega)$ and $\Gamma_{\rm R}(\omega)$ are the imaginary part of the self energy matrices:

[TABLE]

The dynamical matrix is the mass-weighted Hessian matrix of the empirical potential as used in our MD simulations:

[TABLE]

where $U$ is the total potential energy of the system and $u_{\mu}$ is a vibrational degree of freedom in the system with a mass of $m_{\mu}$. The self energy matrices can be obtained from the dynamical matrices for the semi-infinite leads by using an iterative method Sancho et al. (1985). In numerical calculations, each atom contributes three degrees of freedom. The elements of the dynamical matrix can also be expressed as

[TABLE]

where $i$ and $j$ run over the atom indices, while $a$ and $b$ run over the three directions ($a,b=x,y,z$). From the relationship between force and potential, $F_{i}^{a}=-\partial U/\partial u_{i}^{a}$, we can write the above equation as

[TABLE]

This can be evaluated using a finite displacement $\delta$:

[TABLE]

Here $F_{i}^{a}(r_{j}^{b}\pm\delta)$ is the total force on atom $i$ in the $a$ direction caused by displacing the position of atom $j$ in the $b$ direction by an amount of $\pm\delta$ (keeping all the other atoms at their relaxed positions). Here we use $\delta=0.005$ Å.

II.2 NEMD method

There are many variants of the NEMD method, both in terms of boundary conditions in the transport direction and methods of generating the temperature difference. Periodic boundary conditions in the transport direction have been a popular choice, perhaps due to the ease of computer implementation. However, this choice cannot be directly compared to the AGF calculation setup, and is usually different from experimental situations, although it can give equivalent results as obtained using fixed boundary conditions as long as the sample length $L$ is carefully defined Xu et al. (2018). We thus consider fixed boundary conditions in the transport direction. Regarding the methods of generating temperature difference, there are many algorithms, including the constant heat current method Ikeshoji and Hafskjold (1994); Jund and Jullien (1999), the momentum swapping method Müller-Plathe (1997); Kuang and Gezelter (2010), and methods based on thermostats. In both the constant heat current method and the momentum swapping method, the temperature difference cannot be precisely controlled. We therefore use thermostatting methods to generate the temperature difference.

In our NEMD simulations, we fixed a few layers of atoms at the two ends of the sample in the transport direction to achieve the fixed boundary conditions. Next to the two fixed layers, atoms within a length of $L_{\rm th}$ were coupled to a hot and a cold thermal bath, respectively. The distance between the two thermal baths defines the system length $L$, see Fig. 1 for an illustration. The whole simulation cell was first equilibrated at the target temperature using the Berendsen thermostat Berendsen et al. (1984) for $1$ ns, after which we switched off the global thermostat and switched on the local thermostats to realize the hot and cold thermal baths. The local thermostats were applied for $11$ ns, while steady state can be well achieved within the first $1$ ns in all our simulated systems. We used the data within the last $10$ ns to determine the temperature profile and the nonequilibrium heat current. Each NEMD simulation has been repeated three times and the statistical errors are very small and are thus omitted in the relevant figures.

To realize the local thermostats, we consider both the Nosé-Hoover chain Nosé (1984); Hoover (1985); Martyna et al. (1992) and the Langevin methods Lepri et al. (2003); Dhar (2008). In the Nosé-Hoover chain method , the equations of motion for the particles (with position $\boldsymbol{r}_{i}$, momentum $\boldsymbol{p}_{i}$, force $\boldsymbol{F}_{i}$, and mass $m_{i}$) in the thermostatted region (with $N_{\rm f}$ degrees of freedom) are

[TABLE]

where $Q_{0}=N_{\rm f}k_{\rm B}T\tau^{2}$ is the “mass” of the thermosat variable directly coupled to the system, with $\tau$ being a time parameter, and $\pi_{0}$ is the corresponding “momentum”. In the Langevin method, the equations of motion for the particles in the thermostatted region are

[TABLE]

where $\tau$ is a time parameter and $\boldsymbol{f}_{i}$ is a random force with a variation determined by the fluctuation-dissipation relation to recover the canonical ensemble distribution.

There are many implementations of these thermostats and here we use a velocity-Verlet integrator Martyna et al. (1996); Bussi and Parrinello (2007). In both thermostatting methods, there is a time parameter $\tau$ that dictates the coupling between the thermostat and the system. In the Nosé-Hoover chain method, the thermostat mass is proportional to $\tau^{2}$ and a larger $\tau$ would, in principle, decouple adiabatically the degrees of freedom of the thermostat from those of the system. However, when $\tau$ is too large the large inertia of the thermostat degrees of freedom produces unphysical fluctuations in the kinetic energy of the system Bussi et al. (2007). In the Langevin method, $\tau$ is the inverse of the coefficient of friction and a larger $\tau$ also gives a weaker temperature control. The major difference between them is that in the Nosé-Hoover chain method, the temperature in the thermal bath region is adjusted globally and deterministically by re-scaling the velocities of the atoms with a common factor, while in Langevin dynamics, the velocities of the atoms are adjusted locally and stochastically.

II.3 HNEMD method

Although the focus of this work is the NEMD method, we also use some results from the HNEMD method for comparison. This method is physically equivalent to the Green-Kubo method but is computationally much faster Fan et al. (2019). In this method, an external force of the form Fan et al. (2019)

[TABLE]

is added to each atom $i$, driving the system out of equilibrium. Here, $E_{i}$ and $U_{i}$ are the total and potential energies of atom $i$, respectively, $\boldsymbol{r}_{ij}\equiv\boldsymbol{r}_{j}-\boldsymbol{r}_{i}$, and $\boldsymbol{r}_{i}$ is the position of particle $i$. The parameter $\boldsymbol{F}_{\rm e}$ is of the dimension of inverse length and should be small enough to keep the system within the linear response regime. A global thermostat should be applied to control the temperature of the system. The driving force will induce a nonequilibrium heat current $\langle\boldsymbol{J}\rangle_{\rm ne}$ linearly proportional to $\boldsymbol{F}_{\rm e}$. For a given transport direction, this linear relation provides a way to compute the thermal conductivity in this direction:

[TABLE]

where $T$ is the system temperature, $V$ is the system volume, $J=|\boldsymbol{J}|$, and $F_{\rm e}=|\boldsymbol{F}_{\rm e}|$. For a many-body potential, the heat current $\boldsymbol{J}$ is given by Fan et al. (2015)

[TABLE]

where $\boldsymbol{v}_{i}$ is the velocity of atom $i$. Because there is no boundary scattering in this method, the calculated thermal conductivity can be considered as that for an infinitely long system, as long as a sufficiently large periodic simulation cell is used. For graphene, a rectangular cell of dimension $25\times 25$ nm2 is large enough to eliminate the finite-size effects Fan et al. (2019).

II.4 Spectral decomposition method

In both the NEMD and the HNEMD methods, a nonzero nonequilibrim heat current exists and can be spectrally decomposed. Considering an imaginary interface separating two groups of atoms $A$ and $B$ as schematically shown in Fig. 1, this nonequilibrim heat current can be expressed as Fan et al. (2017a)

[TABLE]

For a spatially homogeneous system, $Q/S=J/V$, where $J$ is the magnitude of the heat current given in Eq. (19).

In the spectral decomposition method developed by Sääskilahti et al. Sääskilahti et al. (2014, 2015), one first defines the force-velocity correlation function, which can be expressed as

[TABLE]

for a general many-body potential Fan et al. (2017a). This correlation function reduces to the nonequilibrium heat current in Eq. (20) at zero correlation time. The spectrally decomposed heat current $\tilde{K}(\omega)$ can be obtained from a Fourier transform of the force-velocity correlation function

[TABLE]

The inverse Fourier transform is

[TABLE]

Because

[TABLE]

we obtain the following spectral decomposition of the thermal conductance in NEMD simulations:

[TABLE]

Similarly, the thermal conductivity from HNEMD simulations can be spectrally decomposed as:

[TABLE]

II.5 Details on the numerical calculations

The NEMD, HNEMD, and spectral decomposition methods have been implemented in the highly efficient Graphics Processing Units Molecular Dynamics (GPUMD) package Fan et al. (2013, 2017b); gpu (2017), which is the code we used for most of the MD simulations. The simulations of thermal rectification in asymmetric graphene devices were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package Plimpton (1995). We used the Tersoff potential Tersoff (1989) optimized for graphene systems Lindsay and Broido (2010). For multilayer graphene the van der Waals interactions among different layers were modeled with the Lennard-Jones potential with $\varepsilon=3.296$ eV and $\sigma=3.55$ Å Robertson et al. (2015). A time step of $0.5$ fs, which ensures good energy conservation, was used in all the MD simulations. The AGF calculations were performed by using a Matlab code, which is also publicly available AGF (2019). More details on simulations of thermal rectification are presented in Sec. III.6.

III Results and Discussion

III.1 Ballistic thermal conductance from AGF calculations

We first calculate the ballistic conductance of the graphene sheet (with the same width as in the NEMD simulations) using the AGF method. The ballistic spectral phonon conductance, obtained by using either classical or quantum statistics, is shown in Fig 2. Because anharmonicity is totally absent in this method, the classical conductance (which is essentially the transmission) exhibits many quantized plateaus. By integrating the spectral conductance with respect to the frequency, we can get the total thermal conductance, which is about 12.2 and 4.2 GWm*-2K-1* at room temperature in the cases of classical and quantum statistics, respectively. The quantum thermal conductance obtained here agrees with that reported by Serov et al. Serov et al. (2013) where the same empirical potential was used.

III.2 Comparison between NEMD and AGF in the ballistic regime

To obtain the ballistic conductance using the NEMD method, we set the system length as $L=10$ nm and the simulation temperature as $T=300$ K. The temperatures of the hot and cold thermal baths were set to $T\pm\Delta T/2$ with $\Delta T=10$ K. For such a short sample, anharmonicity caused by phonon-phonon scattering can be largely ignored (except for the optical phonons with the highest frequencies) and the phonon transport can be essentially regarded as ballistic. Here, we calculate the temperature difference from the actual average temperatures in the thermal baths. In the next subsection, we will examine the temperature profile and its implications in the interpretation of the results.

Figure 3 compares the spectral conductance obtained from the NEMD simulations (with spectral decomposition) using both the Nosé-Hoover chain and the Langevin thermostatting methods, against that from the AGF method. For the results from Fig. 3(a) to 3(f), the parameter $\tau$ in both thermostatting methods is increased from $0.05$ ps to $2$ ps. Because the MD simulations are classical, we need to compare with the classical ballistic thermal conductance from the AGF method.

The spectral conductance obtained by using the Nosé-Hoover chain thermostatting method is significantly larger than the reference AGF value, while that obtained by using the Langevin thermostatting method agrees with the AGF reference value well for all the $\tau$ values considered, particularly for intermediate $\tau$ values (0.1 ps to 1 ps). Chen et al. Chen et al. (2010) also recommended using intermediate $\tau$ values in NEMD simulations of heat conduction. This is reasonable, as too small a $\tau$ results in too strong a perturbation on the system, and too large a $\tau$ results in too weak a control of the temperatures in the thermal baths.

III.3 Temperature profile

To better understand the results above, we plot the detailed temperature profiles in the thermal baths as well as in the middle of the system with different coupling constants $\tau$ in Fig. 4. Obviously, the Langevin thermostat gives much better temperature control than the Nosé-Hoover chain thermostat. When $\tau\leq 0.5$ ps, the temperatures in the thermal baths under the action of the Langevin thermostat are close to the target values. In contrast, the Nosé-Hoover chain thermostat cannot maintain a constant temperature in the thermal baths for any value of $\tau$. The temperatures close to the fixed boundaries overshoot the target source temperature or undershoot the target sink temperature. The reason is that Nosé-Hoover is a “global” thermostat, i.e. it rescales the velocities of all the particles by the same amount at every MD step and it can guarantee only that the average kinetic temperature of the thermostatted region is the target one. In contrast Langevin is “local” and it ensures that the whole thermostatted region is at the same temperature, as long as a strong coupling (small $\tau$) is enforced. The use of Nosé-Hoover or any other global thermostat in NEMD results in an effective temperature difference that is larger than that calculated from the average temperatures in the thermal baths and an overestimated thermal conductance as shown in Fig. 3.

Let us now focus on the temperature profile in the middle of the system (the $L=10$ nm part). To this end, we represent the thermal baths as single points with the average temperatures within them; see Fig. 5. We can see that there are abrupt temperature jumps between the thermal baths and the system in the middle. Previously, it has been frequently argued that one should apply Fourier’s law to the linear region of the temperature profile only. This motivates to fit the middle part of the temperature profile using a linear function and extract a temperature gradient $|\nabla T^{\prime}|=\Delta T^{\prime}/L$, where $\Delta T^{\prime}$ (c.f. Fig. 5) is the temperature difference as determined by the interception between the fitted line for the middle part of the temperature profile and the vertical lines at the source and sink. The (effective) thermal conductivity for the system with length $L$ is then calculated as

[TABLE]

According to the relationship between the conductivity and conductance given in Eq. (3), this amounts to using the temperature difference $\Delta T^{\prime}$ to calculate the thermal conductance:

[TABLE]

Figure 6 shows the total thermal conductance (integrated over the frequency) values obtained by using Eqs. (2) and (28) and different thermostatting methods. When Eq. (2) is used, the thermal conductance obtained with both thermostatting methods are relatively close to the AGF value, while the underestimation using the Langevin thermostat is due to the phonon-phonon scattering for the highest-frequency phonons, and the overestimation using the NHC thermostat is due to the nonuniform temperature in the source and sink regions. In contrast, using Eq. (28), the obtained thermal conductance is several times higher than the AGF value, for both thermostatting methods. We thus reach the most important conclusion in this study: Eqs. (27) and (28) are incorrect and should not be used. The correct way is to calculate the thermal conductance using Eqs. (2), taking $\Delta T$ as the temperature difference between the averaged temperatures in the source and sink. Correspondingly, the (effective) thermal conductivity should be calculated using Eq. (4), as we will discuss in Sec. III.4.

III.4 Comparison between NEMD and HNEMD in the ballistic-to-diffusive regime

We now move from the ballistic regime to the ballistic-to-diffusive regime. To this end, we consider systems with different lengths: $L=25$, $50$, $100$, $200$, $500$, $1000$, and $2000$ nm. Other simulation parameters are fixed: $W=10$ nm, $\tau=0.2$ ps, $L_{\rm th}=40$ nm.

Figure 7 shows the temperature profiles for two different system lengths obtained by using the Langevin or the Nosé-Hoover chain thermostat. The difference between the two definitions of temperature difference, $\Delta T$ and $\Delta T^{\prime}$, decreases with increasing $L$, which is still significant when $L=50$ nm, but becomes much smaller when $L=2000$ nm. Due to this difference, the conductance as calculated from Eq. (2) and that from Eq. (28) are different, as can be seen from Fig. 8(a). From the comparison between the NEMD and AGF results in Sec. III.3, we know that Eq. (28) is wrong in the ballistic regime. To show that the corresponding expression for thermal conductivity Eq. (27) is also wrong away from the ballistic regime, we convert the conductance in Fig. 8(a) into the conductivity as shown in Fig. 8(b) using Eq. (3).

There is a very efficient way to calculate the length-dependent thermal conductivity in the ballistic-to-diffusive regime based on the HNEMD-based spectral decomposition method developed recently Fan et al. (2019). In this method, one can calculate the spectral thermal conductivity $\kappa(\omega)$ using Eq. (26). Then one can obtain the spectral phonon mean free path $\lambda(\omega)=\kappa(\omega)/G(\omega)$, from which one can calculate the length-dependent thermal conductivity as Fan et al. (2019):

[TABLE]

The length-dependent thermal conductivity calculated using this method as well as the corresponding thermal conductance [using Eq. (3)] are shown as the solid lines in Fig. 8. One can see that the HNEMD results agree well with the NEMD values obtained by using Eqs. (2) and (4). The Langevin thermostat gives better agreement with the HNEMD results, but the difference between the results from the two thermostatting methods are quite small. The NEMD results obtained by using Eqs. (28) and (27), on the other hand, deviate from the HNEMD results significantly. The relative errors caused by using Eqs. (28) and (27) decrease with increasing system length $L$. If one focuses on the diffusive regime with relatively long systems, using Eq. (27) will not result in large errors. This is why previous works Fan et al. (2017a); Dong et al. (2018a) using Eq. (27) can get agreement between NEMD and other methods in the diffusive regime. However, when $L$ is relatively short (compared to the average phonon mean free path), using Eq. (27) in NEMD simulations will result in large errors.

III.5 Relation to the Boltzmann transport equation

The temperature drop between the thermal baths and the middle system is not an artifact of NEMD simulations, as it has also been observed in Boltzmann transport equation (BTE) calculations Allen (2014). To give a simple demonstration of such an effect, we analytically solved the gray 1D BTE for graphene, assuming an average phonon mean free path of 800 nm Ghosh et al. (2008); Feng et al. (2015). To be consistent with the NEMD (with the Langevin thermostat) and AGF simulations, the left and right boundaries are assumed to be at constant temperatures of 305 K and 295 K, respectively.

The temperature profiles for different domain lengths $L$ are shown in Fig. 9. Note that in BTE the constant temperature (or thermalizing) boundary condition is implemented in the way that all the outgoing phonons leave the boundary unaffected, while all the incoming phonons have an intensity corresponding to equilibrium distribution at the given temperature Zahiri et al. (2018), as schematically shown in the inset of Fig. 9. Such an implementation is equivalent to having two boundaries in contact with infinitely large external thermal baths, and therefore is consistent with the Langevin thermostat as implemented in our NEMD simulations and also the AGF calculations. The temperature discontinuity is straightforward in the BTE picture. For example, in the ballistic limit, since phonons are not thermalized in the middle region, the phonons traveling from left to right have the same energy as the left boundary (305 K), while the phonons traveling from right to left have the same energy as the right boundary (295 K). Under such a non-equilibrium condition, the “effective” temperature in the simulation domain is the average temperature of the two boundaries. Therefore, there is a temperature discontinuity at the boundary. As the system length increases, the phonons within the domain experience stronger phonon-phonon scatterings, so that the discontinuity gradually decreases and eventually diminishes. This is similar to what has been demonstrated in Fig. 7. Note that for BTE simulation, it is well known that the conductance in this case should be calculated using the temperature difference between the boundaries Chen (2005) instead of that in the middle region. Therefore, it further confirms that the conductance and conductivity should be calculated using Eq. (2) and Eq. (4), instead of Eq. (28) and Eq. (27).

Note that here we have only considered a simplified gray BTE model to demonstrate that the temperature profiles in BTE calculations are qualitatively similar to those in NEMD simulations. Gu et al. Gu et al. (2019) have recently performed BTE calculations with a full iterative scheme, considering both three- and four-phonon scattering processes as well as boundary scatterings in graphene sheets. In such calculations, one can directly obtain length dependent thermal conductivity. They found that good agreement between NEMD and BTE regarding the length-dependence of thermal conductivity can only be obtained if Eq. (4) is used to compute the thermal conductivity in NEMD simulations.

III.6 Thermal rectification in asymmetric graphene-based systems

The discovery of thermal rectification in low-dimensional non-linear lattice models Terraneo et al. (2002); Li et al. (2004) fostered significant efforts to identify efficient thermal diodes Chen et al. (2015); Benenti et al. (2016); Chen et al. (2018), which would lay the cornerstone of nanophononic circuitry Li et al. (2012); Volz et al. (2016). In this context NEMD has been extensively used to probe thermal rectification in asymmetric graphene-based nanodevices, including branched nanoribbons, triangular patches and multilayer junctions  Hu et al. (2009); Yang et al. (2009); Zhong et al. (2011); Wang et al. (2017); López-Suárez et al. (2018). Using the definition of thermal conductance in Eq. (2), the rectification factor is defined in terms of the difference between high and low thermal conductances $G_{\rm H}$ and $G_{\rm L}$ obtained by swapping the temperature bias of the thermal diode:

[TABLE]

$G_{\rm H/L}$ may be replaced by the heat currents $J_{\rm H/L}$, provided that the temperature difference between the thermal reservoirs remains the same when the temperature bias is inverted. The above mentioned simulation works predict extremely high $\eta$, up to 350$\%$, for carbon-based devices, but experimental measurements show much smaller thermal rectification, if any at all Chang et al. (2006); Wang et al. (2017); Rojo et al. (2019).

Here we perform NEMD simulations of thermal rectification in both large trapezoid (LT) and small trapezoid (ST) monolayer graphene patches, where the large trapezoid compares to the smallest trapezoid studied in Ref. Wang et al., 2017, which was reported to have large rectification, and in multilayer graphene junctions: bilayer to monolayer (BTM), trilayer to monolayer (TTM) and quadlayer to monolayer (QTM), where the multilayers are of the same geometry as those in Zhong et al. (2011). All the layers were thermalized in order to compare with the previous numerical studies  Zhong et al. (2011). In all our thermal rectification simulations, we fixed two layers of atoms at the two ends of the sample in the transport direction to achieve the fixed boundary conditions. Next to the two fixed layers, atoms within a length of $L_{\rm th}$ ($1.7$ nm for LT, $0.5$ nm for ST, and $0.8$ nm for other systems) were coupled to a hot and a cold thermal bath, respectively. We obtain the rectification factor of these nanodevices by computing the heat current in two separate NEMD simulations for each system, in which the hot and cold reservoirs are swapped, so to probe the high-conductance and the low-conductance conditions. We compare the rectification factor obtained using either Nosé-Hoover, Nosé-Hoover chain, or Langevin dynamics to set the temperatures at 350 and 250 K for the hot and cold thermal reservoirs. After equilibration at 300 K, NEMD simulations are run for 4 ns (8 million time steps), and the thermal conductance is computed from the last 2 ns of these runs.

We first test the effect of using different types of thermostats on the rectification factor of the LT graphene patch with 21.6 nm and 2 nm bases and height of 17 nm. This system is made of 8049 atoms and former NEMD simulations using the Nosé-Hoover thermostat predicted a rectification factor $\sim 95\%$ Wang et al. (2017). Our simulations with the standard Nosé-Hoover thermostat with $\tau=0.1$ ps are in accord with these former results, predicting $\eta=111(\pm 2)\%$. Adding more degrees of freedom to the thermostat, using Nosé-Hoover chains does not significantly change the estimate for $\eta$, which remains $\geq 100\%$ (Fig. 10). However, when Langevin dynamics is employed to fix the temperature of the thermal baths and produce a stationary flux, $\eta$ is considerably reduced which matches the results of a previous study Chen et al. (2010). Furthermore, $\eta$ depends on the relaxation time used in Eq. (16). The smaller the $\tau$, i.e. the stronger the coupling, the lower the rectification factor and its uncertainty. For $\tau\leq 0.1$ ps $\eta$ can be considered statistically zero. As we have shown in the previous sections the NEMD Langevin dynamics provides more accurate results than Nosé-Hoover chain, and thus we can argue that the LT graphene patch considered here should not exhibit any significant thermal rectification.

To unravel the origin of the large discrepancy between the two simulation methods, we calculated the power spectrum of the thermal baths in the NEMD simulations of trapezoid graphene at direct and reverse bias conditions (Fig. 11). When the hot bath is on the large base of the trapezoid the hot and cold power spectra overlap, except for the high frequency peak, which does not contribute significantly to heat transport. This condition corresponds to high conductance and the two thermostatting methods give similar results. With reverse bias, however, the hot bath exhibits an extremely intense peak at 12 THz in the case of the Nosé-Hoover chain thermostat, indicating that vibrational modes at this frequency are overpopulated. The energy accumulated in these specific modes cannot transfer through the device, leading to low conductance. Such overpopulation of a specific mode is an artifact of the Nosé-Hoover chain thermostat, occurring in non-ergodic system with a small number of degrees of freedom Patra and Bhattacharya (2014). This effect leads to a disbalance between the power pumped by the thermostat into the hot bath and that removed from the cold bath, which should be equal at stationary conditions. With Nosé-Hoover chain these conditions are not achieved over the typical run time of several ns. When Langevin dynamics is used, the coupling parameter $\tau$ determines how rapidly stationary conditions with constant flux are attained, with shorter $\tau$ providing faster thermalization.

These results suggest that Nosé-Hoover is unsuitable to study thermal rectification, unless one carefully addresses the very long thermalization time. Hence we use Langevin dynamics with $\tau=0.01$ ps to calculate the rectification factor of the other junctions considered. Results in Fig. 12 show that for either trapezoidal or multilayer graphene junctions $\eta<15\%$, in contrast with the high rectification efficiencies predicted in NEMD simulations carried out using Nosé-Hoover. Whereas multilayer graphene junctions may attain measurable rectification, all these systems are unsuitable for practical applications as thermal diodes.

IV Summary and conclusions

In summary, we have carefully and systematically revisited some critical issues in the frequently used NEMD methods for thermal transport calculations. By comparing with the AGF method in the ballistic regime and the HNEMD method in the ballistic-to-diffusive regime, we found that the nonlinear part of the temperature profile in NEMD simulations should not be excluded in the calculations of the thermal conductivity and conductance. We also found that the Langevin thermostatting method controls the local temperatures better and is more reliable than the Nosé-Hoover (chain) thermostatting method for NEMD simulations. This is particularly important for studying asymmetric nanostructures, for which the Nosé-Hoover thermostat can produce artifacts leading to unphysical thermal rectification. Based on our results, we recommend an intermediate value ($\tau=0.1-1$ ps) for the time parameter in the Langevin thermostat for thermal conductivity calculations and a small value ($\tau\leq 0.1$ ps) for thermal rectification studies, where one usually considers large temperature differences.

Acknowledgements.

ZL and NW were supported by the National Natural Science Foundation of China (Grant No 11502217). SX acknowledges the support from the National Natural Science Foundation of China (Grant No. 11804242) and the Jiangsu Provincial Natural Science Foundation (Grant No. BK20160308). YH and HB acknowledge the support from the National Natural Science Foundation of China (Grant No. 51676121). ZF and TA-N acknowledge support from the Academy of Finland Centre of Excellence program QTF (Project 312298) and the computational resources provided by Aalto Science-IT project and Finland’s IT Center for Science (CSC).