Long-time instability in the Runge-Kutta algorithm for a Nosé-Hoover heat bath of a harmonic chain and its stabilization

In this paper, we investigate the Runge-Kutta algorithm for the Nosé-Hoover heat bath of a harmonic chain. The Runge-Kutta algorithm is found to be unstable in long-time calculations, with the system temperature growing exponentially. The growth rate increases if time step size is chosen larger. By analyzing the Fourier spectra in both space (wave number) and time (frequency), we discover that the growth is caused by spurious energy accumulation, particularly at the largest wave number. Such accumulation may be explained by von Neumann analysis for an infinite chain, with the nonlinear heat bath being ignored. Furthermore, we propose to add a filter to remove excessive energy, which effectively stabilizes the algorithm.

Figure 1

Runge-Kutta algorithm with step size $h=0.01$: (a) system temperature, (b) local temperature, (c) time average of wave number spectrum, (d) instantaneous wave number spectrum at $\xi =\pi$, (e) frequency spectrum, and (f) instantaneous frequency spectrum at $\omega =2$.

Figure 2

Runge-Kutta algorithm with step size $h=0.05,0.02,0.01,0.008,0.006$ from top to bottom. Stars represent numerical growth rate, and solid lines represent theoretical growth rate.

Figure 3

Runge-Kutta algorithm with step size $h=0.005$: (a) System temperature, (b) Local temperature.

Figure 4

The stabilized Runge-Kutta algorithm with a filtering time $\mathrm{\Delta }{t}_f=4\times {10}^5$: (a) system temperature, (b) local temperature, (c) time average of wave number spectrum, and (d) instantaneous wave number spectrum at $\xi =\pi$.

Figure 5

The stabilized Runge-Kutta algorithm with different filtering time: (a) local temperature with $\mathrm{\Delta }{t}_f=1\times {10}^5$, (b) time average of wave number spectrum with $\mathrm{\Delta }{t}_f=1\times {10}^5$, (c) local temperature with $\mathrm{\Delta }{t}_f=8\times {10}^5$, and (d) time average of wave number spectrum with $\mathrm{\Delta }{t}_f=8\times {10}^5$.