## Ergodicity of a Time-Reversibly Thermostated Harmonic Oscillator and the 2014 Ian Snook Prize

William Graham Hoover and Carol Griswold Hoover Ruby Valley Research Institute Highway Contract 60, Box 601 Ruby Valley, Nevada 89833 (Dated: August 4, 2014)

## Abstract

Shuichi Nos´ opened up a new world of atomistic simulation in 1984. e He formulated a Hamiltonian tailored to generate Gibbs' canonical distribution dynamically. This clever idea bridged the gap between microcanonical molecular dynamics and canonical statistical mechanics. Until then the canonical distribution was explored with Monte Carlo sampling. Nos´'s dynamical Hamiltonian e bridge requires the 'ergodic' support of a space-filling structure in order to reproduce the entire distribution. For sufficiently small systems, such as the harmonic oscillator, Nos´'s dynamical ape proach failed to agree with Gibbs' sampling and instead showed a complex structure, partitioned into a chaotic sea, islands, and chains of islands, that is familiar textbook fare from investigations of Hamiltonian chaos. In trying to enhance small-system ergodicity several more complicated 'thermostated' equations of motion were developed. All were consistent with the canonical Gaussian distribution for the oscillator coordinate and momentum. The ergodicity of the various approaches has undergone several investigations, with somewhat inconclusive ( contradictory ) results. Here we illustrate several ways to test ergodicity and challenge the reader to find even more convincing algorithms or an entirely new approach to this problem.

## PACS numbers:

Keywords: Ergodicity, Chaos, Algorithms, Dynamical Systems

## I. INTRODUCTION

In 1989 Shuichi Nos´ won IBM-Japan's science prize for his 1984 discovery e 1,2 of a dynamical Hamiltonian approach to Gibbs' canonical ensemble. He added a new 'time-scaling' variable s ( and its conjugate momentum p s ) to the usual list of # Cartesian degrees of freedom, coupling the kinetic energy to a logarithmic temperature-dependent potential for s :

$$\sigma \. \\ \mathcal { H } _ { N } ( q, p, s, p _ { s } ) = \sum ^ { \# } ( p ^ { 2 } / 2 m s ^ { 2 } ) + \Phi ( q ) + ( p _ { s } ^ { 2 } / 2 M ) + ( \# + 1 ) k T \ln ( s ) \.$$

Multiplying the resulting equations of motion by s ( 'scaling the time' ) gives a set of (2# + 2) dynamical motion equations fully consistent with Gibbs' canonical ensemble ,

$$f ( q, p ) \in e ^ { - \mathcal { H } / k T }$$

provided that we replace ( p/s ) by p in the scaled equations of motion.

Hoover 3 soon pointed out the lack of ergodicity for this approach, applied to a harmonic oscillator, and at the same time simplified the derivation of a set of (2# + 1) equations of motion only slightly different to Nos´'s : e

$$\{ \, \dot { q } = ( p / m ) \ ; \, \dot { p } = F - \zeta p \, \} \ ; \, \dot { \zeta } = \sum ^ { \# } [ \ ( p ^ { 2 } / m k T ) - 1 \ ] / \tau ^ { 2 } \.$$

The friction coefficient ζ in these 'Nos´ e-Hoover' equations is proportional to Nos´'s e p s . These equations are much better behaved numerically and also make the extra s variable redundant. Carl Dettmann showed 4 that these same equations follow from a Hamiltonian like Nos´'s provided that the Dettmann Hamiltonian is arbitrarily set equal to zero and the e number of degrees of freedom is # rather than # + 1 ,

$$\mathcal { H } _ { D } \equiv s \mathcal { H } _ { N } \equiv 0 \.$$

Hoover and Harald Posch and Franz Vesely 5,6 investigated the details of these dynamics for the harmonic oscillator and discovered an infinite variety of periodic, toroidal, and chaotic solutions, with the sum of all of these disparate parts equal to the canonical distribution, the product of three Gaussian functions, in q , in p , and in ζ for the oscillator.

Bauer, Bulgac, and Kusnezov 7,8 generalized the thermostating approach to include two or more control variables, even managing to reproduce Brownian motion by using three of them. Soon after, Martyna, Klein, and Tuckerman 9 suggested the use of a chain of

thermostat variables to enhance ergodicity. Their equations of motion ( for the shortest chain and with all the free parameters equal to unity ) are :

$$\{ \ \dot { q } = p \ ; \ \dot { p } = - q - \zeta p \ ; \ \dot { \zeta } = p ^ { 2 } - 1 - \xi \zeta \ ; \ \dot { \xi } = \zeta ^ { 2 } - 1 \ \} \. \ [ \ M K T \ ]$$

Presently, Hoover and Holian suggested a different dual-control approach 10 , fixing both the second and fourth velocity moments :

$$\{ \ \dot { q } = p \ ; \ \dot { p } = - q - \zeta p - \xi p ^ { 3 } \ ; \ \dot { \zeta } = p ^ { 2 } - 1 \ ; \ \dot { \xi } = p ^ { 4 } - 3 p ^ { 2 } \ \} \. \ [ \ H H \ ]$$

Equations controlling the sixth or higher moments tend to be too stiff for practical use.

Very recently Patra and Bhattacharya 11 took a different direction, stimulated by work on the 'configurational temperature', as opposed to the usual kinetic one. For the oscillator the configurational temperature 12 is proportional to the potential energy, 〈 q 2 〉 . Patra and Bhattacharya had the clever idea of imposing both temperatures simultaneously on the oscillator. Their feedback equations use one friction coefficient, ζ , to control kinetic temperature, and the other, ξ , to control configurational temperature , 〈 q 2 〉 :

$$\{ \ \dot { q } = p - \xi q \ ; \ \dot { p } = - q - \zeta p \ ; \ \dot { \zeta } = p ^ { 2 } - 1 \ ; \ \dot { \xi } = q ^ { 2 } - 1 \ \} \. \ [ \ P B \ ]$$

We wish to emphasize that all three approaches followed here, Martyna-Klein-Tuckerman, Hoover-Holian, and Patra-Bhattacharya, if ergodic, provide exactly the same fourdimensional Gaussian distribution ,

$$( 2 \pi ) ^ { 2 } f ( q, p, \zeta, \xi ) \equiv e ^ { - ( q ^ { 2 } / 2 ) } e ^ { - ( p ^ { 2 } / 2 ) } e ^ { - ( \zeta ^ { 2 } / 2 ) } e ^ { - ( \xi ^ { 2 } / 2 ) } \.$$

To prove that this is so it is only necessary to evaluate Liouville's continuity equation for the flow in phase space, showing that the four-dimensional Gaussian distribution is stationary :

$$( \partial f / \partial t ) \equiv - f [ \, ( \partial \dot { q } / \partial q ) + ( \partial \dot { p } / \partial p ) + ( \partial \dot { \zeta } / \partial \zeta ) + ( \partial \dot { \xi } / \partial \xi ) \ ]$$

$$& ( \partial f / \partial t ) \equiv - f [ \ ( \partial \dot { q } / \partial q ) + ( \partial \dot { p } / \partial p ) + ( \partial \dot { \zeta } / \partial \zeta ) + ( \partial \dot { \xi } / \partial \xi ) \ ] \\ & \quad - \dot { q } ( \partial f / \partial q ) - \dot { p } ( \partial f / \partial p ) - \dot { \zeta } ( \partial f / \partial \zeta ) - \dot { \xi } ( \partial f / \partial \xi ) \equiv 0 \.$$

Let us turn to the question of establishing ergodicity for the three dynamical models, MKT, HH, and PB . Our own approach, because we know it best, is computational. It is possible that a more convincing 'rigor mortis' approach could be developed from a mathematical standpoint. At the moment any 'proof' of ergodicity is an 'open problem'. In seven brief sections we describe some of the methods that have been applied to this problem 13 , ending with a challenge for the reader.

(q,p) Projections for the MKT, HH, and PB Oscillator Models

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0256v1-figure-3fd0f46e-6404-4784-b530-2fc31dce6ac0.jpg)

FIG. 1: Chaotic sea for the Martyna-Klein-Tuckerman, Hoover-Holian, and Patra-Bhattacharya oscillators. The initial condition for these trajectories was ( q, p, ζ, ξ ) = (0 5 0 0) in all three cases. , , , 100,000 equally spaced ( q, p ) points are plotted for each model. The Patra-Bhattacharya oscillator has clear deviations from hyperspherical symmetry. The fourth-order Runge-Kutta timestep is 0.001 for all of the computations in our work. Here ten million timesteps were used. Note the holes in the MKT and HH projections and the lack of circular symmetry in the PB projection.

## II. TIME AVERAGES OF THE MOMENTS

Because kinetic temperature is the second velocity moment it is usual to confirm that the first few velocity moments agree with the canonical distribution. A short fourth-order Runge-Kutta calculation, adding ( p , p 2 4 , p 6 ) to the list of righthandsides being integrated, shows that the second, fourth, and sixth moments of the MKT and HH equations match the values expected from a Gaussian, (1, 3, 15) to three or four significant figures while the PB model gives instead (1.000, 3.8, 23.5) . These results are conclusive evidence that only the MKT and HH equations are candidates for ergodicity. Nevertheless we will apply each of our numerical tests to all three sets of motion equations, in order to see the difference between the difficult and the easy.

## III. PROJECTIONS

The Projections of the four-dimensional dynamics onto any plane, such as the ( q, p ) , ( q, ζ ) , ( q, ξ ) , ( p, ζ ) , ( p, ξ ) , ( ζ, ξ ) planes, should look Gaussian without any holes . For simplicity we use the initial condition p = 5 with the other variables initially zero, promoting chaos for all three sets of equations. The MKT and HH equations both appear to have a small hole near the origin of the ( qp ) plane. The PB equations produce a very different pattern, showing a strong positive correlation between the two variables, another indication that these equations though sometimes chaotic, are not ergodic. See Figure 1 for a comparison of the three models in the ( q, p ) planar projection. These results suggest two further checks, a look at the density near the ( q, p ) origin and a search for periodic orbits, which would be necessary in order that any holes form in the Gaussian distributions.

## IV. DENSITY

The probability density for the four-dimensional Gaussian at a radius of 0 1 has decreased . from its maximum value at the (0 0 0 0) origin by a factor of , , ,

$$[ \ \rho ( r = 0. 1 ) / \rho ( r = 0. 0 ) \ ] = e ^ { - r ^ { 2 } / 2 } = e ^ { - 0. 0 0 5 } \simeq 1. 0 0 0 0 0 - 0. 0 0 5 0 0 + 0. 0 0 0 0 1 = 0. 9 9 5 0 1$$

so that the densities determined from trajectories should hardly vary within that fourdimensional sphere provided that the flow is ergodic . Measuring density requires considerably longer runs because the time interval between visits near the origin increases at least as rapidly as ( r -4 ) . We measured the probability densities for the MKT and HH algorithms inside ( q, p ) circles of area 0 01 . π × (1 , 1 / 2 1 , / 4 1 , / 8 1 , / 16) using a billion timesteps, and found no significant difference in density over that range. Evidently the MKT and HH algorithms behave as though they were ergodic. The difference between those densities and the PB chaotic density was approximately a factor of two, indicating that approximately half the PB measure is chaotic and the other half quasiperiodic, and so omitted from this chaotic simulation.

It is worth mentioning that the MKT equations have two fixed points, though in the end neither one is actually attractive. In the special case that q and p both vanish so that ˙ q and ˙ p are likewise zero, the remaining equations of motion are :

$$\dot { \zeta } = - 1 - \xi \zeta \ ; \ \dot { \xi } = \zeta ^ { 2 } - 1 \.$$

These equations describe a flow from the unstable two-dimensional fixed point ( -1 +1) , to the stable ( only in two dimensions ) fixed point (+1 , -1) . These flow equations are isomorphic to those of a one-dimensional particle in a constant field, with ζ playing the rˆle o of momentum and ξ acting as the friction coefficient. Patra and Bhattacharya's assertion that the MKT equations are not ergodic relies on analyses in the vicinity of these two fixed points. That type of work is complicated by the fact that the four-dimensional flow becomes more and more intermittent as one nears either of these two points. Nevertheless we believe that our analysis in this Section and particularly in Section VII below casts doubt on their claim ( hence our offering of the 2014 Ian Snook Prize on that subject, as is described in Section IX ) .

## V. PERIODIC ORBITS

Choosing an initial condition somewhat closer to the origin

$$\{ \ q, p, \zeta, \xi \ \} = \{ \ 0. 6 8, 0. 6 8, 0, 0 \ \}$$

gives ( q, p, ζ ) projection plots which 'look' spherical for the MKT and HH equations. The PB equations, on the other hand, provide a clearly-defined torus, shown in Figure 2 . The PB dynamics has been projected into three-dimensional ( q, p, ζ ) space. The PB projection shows that the phase-space distribution has at least one large cavity in the chaotic sea and that the cavity contains a family of nested tori.

## VI. ENSEMBLE TESTS

Because a lack of ergodicity implies segregated regions in the ( q, p, ζ, ξ ) space it is a tempting idea to study ensembles of initial conditions, expecting to find two or more distinct longtime-average values rather than a Gaussian Central-Limit-Theorem distribution around a single ergodic average value. In addition to moments, and their correlations, the largest Lyapunov exponent should be a particularly 'good' property to follow. Although some complicated attractors ( such as Rayleigh-B´nard flows ) have more than a single chaotic e strange-attractor region such a situation is implausible for the simple oscillator. Therefore

FIG. 2: Patra-Bhattacharya torus with initial conditions ( q, p, ζ, ξ ) = (0 68 0 68 0 0). . , . , , Choosing (0.5,0.5,0.0,0.0) generates a very thin torus close to the periodic orbit at the center of the larger torus shown in the figure. 50,000 points are plotted.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0256v1-figure-6303a60d-40a0-4aef-9d96-9201aed76595.jpg)

it seems likely that a random or a grid-based selection of initial points would produce at least a bimodal distribution of values for a nonergodic set of equations. Some ensemble tests of this kind were carried out in References 6 and 13.

## VII. MANY LYAPUNOV EXPONENTS

Taking up the ensemble idea we choose 1000 different initial conditions in the fourdimensional hypercube of sidelength four centered on the origin, following each of them for 10 6 timesteps. Though there is no reason for the HH and MKT exponents to agree there turns out to be good agreement between them, as Figure 3 shows. On the other hand the PB equations' Lyapunov exponents separate into two distinct values, 0.0000 and 0.14, characteristic of their separate portions of the phase space. We view the density and Lyapunov tests as the most convincing evidence that both the MKT and the HH algorithms are ergodic.

## VIII. HISTORY OF THE EXTREMAL LYAPUNOV EXPONENTS

To nail this conclusion down we took the two extreme cases ( the maximum and minimum Lyapunov exponent out of 1000 simulations ) and ran them 1000 times longer. Figure 3 shows the behavior of the ensemble of initial conditions for the Martyna-Klein-Tuckerman and Hoover-Holian models. Figure 4 shows a comparison of the maximum and minimum

## Exponents for Oscillator Models

FIG. 3: Time-averaged Lyapunov exponents for the Martyna-Klein-Tuckerman, Hoover-Holian, and Patra-Bhattacharya oscillators. 1000 trajectories are traced from points located randomly within a 4 × 4 × 4 × 4 hypercube centered at the coordinate origin. Although there is no reason that the MKT and HH results should agree, or nearly so, the two models have similar Lyapunov exponents. The Patra-Bhattacharya motion equations lead to two clumps of exponents at long times, one at 0.0000 and the other at 0.14.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0256v1-figure-c5e09640-0d4b-46a6-a53b-a001ced4a7f2.jpg)

Lyapunov-exponent calculations, carried out for a time 1000 times longer. The MKT and HH data show no significant difference between them, strongly suggesting that the entire set of random trajectories started in a 4 × 4 × 4 × 4 hypercube are part of the same chaotic sea, filling all of space in a Gaussian manner for the Martyna-Klein-Tuckerman and HooverHolian oscillators.

## IX. FURTHER TESTS

Patra and Bhattacharya also determined phase-space density numerically, in one- or twodimensional cross sections 15 . They also computed mean deviations of long-time-averaged distributions from the canonical one as averages. Their final conclusion is that the MartynaKlein-Tuckerman algorithm is not ergodic. Our own evidence for that same problem points in

## Exponents for Oscillator Models

FIG. 4: Time-averaged Lyapunov exponents for two runs 1000 times longer than those in Figure 3. The initial conditions were chosen to match those of the minimum and maximum values shown at the righthand margin of Figure 3. The convergence of the two outliers toward a common value is a strong indicator of ergodicity.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0256v1-figure-ec5197d3-ea05-473d-abfb-143599f2539b.jpg)

the opposite direction, toward ergodicity. We thank Puneet Patra, Baidurya Bhattacharya, and Clint Sprott for an exchange of hundreds of stimulating emails on this general subject, and specially appreciated a personal visit from Professor Bhattacharya this summer.

## X. CONCLUSIONS - IAN SNOOK PRIZE FOR 2014

The disagreement between our own investigations past and present ( which agree with those of Martyna, Klein, and Tuckerman ) and those of Patra and Bhattacharya are both thought-provoking and stimulating. Last year we offered a prize to honor the memory of our Australian Colleague Ian Snook. We asked arXiv readers for a time-reversed version of a simple pseudorandom number generator with a cycle length of 2 22 . Within 24 hours Professor Federico Ricci-Tersenghi ( University of Rome ) won the prize with the elegant solution described in his arXiv.1305.1805 contribution 16 , 'The Solution to the Challenge in 'Time-Reversible Random Number Generators' by Wm. G. Hoover and Carol G. Hoover'. His generator can be used, for instance, to reverse Brownian dynamics trajectories, a demonstration problem in our forthcoming book on the management and control of nonequilibrium systems.

FIG. 5: Shuichi Nos´ ( 1951-2005 ) and Ian Snook ( 1945-2013 ) e

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0256v1-figure-36260867-8b92-4f1d-a474-19e5e5508518.jpg)

This year we invite the readers of Computational Methods in Science and Technology to consider thoughtfully this interesting problem in classical statistical mechanics and dynamical systems theory. Specifically the Ian Snook Prize for 2014 will be awarded to whomever submits the most convincing treatment of this problem to us prior to 1 January 2015. The question to be discussed is 'To what extent are trajectory-based solutions of the MartynaKlein-Tuckerman harmonic-oscillator motion equations ergodic?'

It is our intention to reward the most convincing entry received ( or submitted to the arXiv or to Computational Methods in Science and Technology ) prior to 1 January 2015. The 2014 Ian Snook prize of five hundred United States dollars will be presented to the winner in January 2015. We would be very grateful for your contributions. We dedicate this work to the memories of our two colleagues shown in Figure 5 .

## XI. ACKNOWLEDGMENTS

We thank Puneet Patra and Baidurya Bhattacharya for bringing this problem to our attention and additionally thank Clint Sprott for his insights into a wide variety of fascinating problems in dynamical systems theory and for his thoughts on the current prize offer.

- 1 S. Nos´, 'A Unified Formulation of the Constant Temperature Molecular Dynamics Methods', e Journal of Chemical Physics 81 , 511-519 (1984).
- 2 S. Nos´ e, 'Constant Temperature Molecular Dynamics Methods', Progress in Theoretical Physics Supplement 103 , 1-46 (1991).
- 3 Wm. G. Hoover, 'Canonical Dynamics: Equilibrium Phase-Space Distributions', Physical Review A 31 , 1695-1697 (1985).
- 4 C. P. Dettmann and G. P. Morriss, 'Hamiltonian Formulation of the Gaussian Isokinetic Thermostat', Physical Review E 54 , 2495-2500 (1996).
- 5 H. A. Posch, W. G. Hoover, and F. J. Vesely, 'Canonical Dynamics of the Nos´ e Oscillator: Stability, Order, and Chaos', Physical Review A 33 , 4253-4265 (1986).
- 6 H. A. Posch and Wm. G. Hoover, 'Time-Reversible Dissipative Attractors in Three and Four Phase-Space Dimensions', Physical Review E 55 , 6803-6810 (1997).
- 7 D. Kusnezov, A. Bulgac, and W. Bauer, 'Canonical Ensembles from Chaos', Annals of Physics 204 , 155-185 (1990).
- 8 D. Kusnezov and A. Bulgac, 'Canonical Ensembles from Chaos: Constrained Dynamical Systems', Annals of Physics 214 , 180-218 (1992).
- 9 G. J. Martyna, M. L. Klein, and M. Tuckerman, 'Nos´-Hoover Chains-the Canonical Ensemble e via Continuous Dynamics', Journal of Chemical Physics 97 , 2635-2643 (1992).
- 10 Wm. G. Hoover and B. L. Holian, 'Kinetic Moments Method for the Canonical Ensemble Distribution', Physics Letters A 211 , 253-257 (1996).
- 11 P. K. Patra and B. Bhattacharya, 'A Deterministic Thermostat for Controlling Temperature using All Degrees of Freedom', Journal of Chemical Physics 140 , 064106 (2014).
- 12 K. P. Travis and C. Braga, 'Configurational Temperature and Pressure Molecular Dynamics: Review of Current Methodology and Applications to the Shear Flow of a Simple Fluid',

Molecular Physics 104 , 3735-3749 (2006).

- 13 Wm. G. Hoover, Carol Hoover, and Dennis Isbister, 'Chaos, Ergodic Convergence, and Fractal Instability for a Thermostated Canonical Harmonic Oscillator', Physical Review E 63 , 35413546 (2000).
- 14 Wm. G. Hoover and C. G. Hoover, 'Time-Reversible Random Number Generators : Solution of Our Challenge by Federico Ricci-Tersenghi': arXiv.1305.0961.
- 15 Puneet Kumar Patra and Baidurya Bhattacharya, 'Non-Ergodicity of Nos´-Hoover Chain There mostat in Computationally Achievable Time': arXiv.1407.2353.
- 16 F. Ricci-Tersenghi, 'The Solution to the Challenge in 'Time-Reversible Random Number Generators' by Wm. G. Hoover and Carol G. Hoover': arXiv.1305.1805.
William Graham Hoover Carol Griswold Hoover
Ergodicity of a Time-Reversibly Thermostated Harmonic Oscillator and the 2014 Ian Snook Prize
Shuichi Nos\'e opened up a new world of atomistic simulation in 1984. He
formulated a Hamiltonian tailored to generate Gibbs' canonical distribution
dynamically. This clever idea bridged the gap between microcanonical molecular
dynamics and canonical statistical mechanics. Until then the canonical
distribution was explored with Monte Carlo sampling. Nos\'e's dynamical
Hamiltonian bridge requires the "ergodic" support of a space-filling structure
in order to reproduce the entire distribution. For sufficiently small systems,
such as the harmonic oscillator, Nos\'e's dynamical approach failed to agree
with Gibbs' sampling and instead showed a complex structure, partitioned into a
chaotic sea, islands, and chains of islands, that is familiar textbook fare
from investigations of Hamiltonian chaos. In trying to enhance small-system
ergodicity several more complicated "thermostated" equations of motion were
developed. All were consistent with the canonical Gaussian distribution for the
oscillator coordinate and momentum. The ergodicity of the various approaches
has undergone several investigations, with somewhat inconclusive (
contradictory ) results. Here we illustrate several ways to test ergodicity and
challenge the reader to find even more convincing algorithms or an entirely new
approach to this problem.