## VI. NUMERICAL RESULTS 

In its simplicity, the dynamics of a harmonic mode in one dimension is a paradigmatic example for checking the chaotic (ergodic) properties of constant-temperature phase space flows and the correct sampling of the canonical distribution. It is well known that it is necessary to generalize basic Nos´ e-Hoover dynamics [1, 8, 9] to thermostats such as the Nos´ e-Hoover chain [32, 37] in order to produce correct constant-temperature averages for systems such as the harmonic oscillator. 

Some time ago, BK dynamics was devised to provide a deterministic thermostat for systems such as classical spins [23, 24]. To ensure efficient thermostatting, BK found it necessary to introduce several ‘demons’ per thermostatted degree of freedom, where each demon was taken to have a different and in principle complicated coupling to the system degree of freedom [23, 24]. In the present work, we keep the form of the system-thermostat coupling as simple as possible, in order to facilitate the formulation of explicit, reversible and measure-preserving integrators [19]. It is then of interest to investigate the ability of our BK-type thermostats to produce the correct canonical sampling in the case of the harmonic oscillator. Interest in harmonic modes is also justified by the possibility of devising models of condensed matter systems in terms of coupled spins and harmonic modes, as already done in quantum dynamics with so-called spin-boson models [38]. We therefore investigate the performance of our integration schemes on the simple one-dimensional harmonic oscillator. 

For the particular calculations reported here, the oscillator angular frequency, all masses and k$_{B}$T were taken to be unity. The time step in all cases was τ = 0 . 0025, and all runs were calculated for 10 6 time steps, starting from the same initial conditions: harmonic oscillator coordinate q = 0 . 3, all other phase space variables zero at t = 0. 

The measure-preserving algorithms derived here result in stable numerical integration for all the three cases treated: BK, BKNH, and BKNHC chain dynamics. Figure 1 shows the three extended Hamiltonians (normalized by their respective initial time value) versus time. All three Hamiltonians are numerically conserved by the corresponding measure-preserving algorithm to very high accuracy (which is maintained in all the three cases). 

However, the basic BK phase space flow is not capable of producing the correct canonical sampling for a harmonic mode. This can be easily checked since the canonical distribution function of the harmonic oscillator is isotropic in phase space and its radial dependence 

