VII. CONCLUSIONS We have formulated Bulgac-Kusnezov [23, 24], Nos´ e-Hoover controlled Bulgac-Kusnezov, and Bulgac-Kusnezov-Nos´ e-Hoover chain thermostats in phase space by means of non- Hamiltonian brackets [14, 15]. We have derived time-reversible measure-preserving algo- rithms [19] for these three cases and showed that additional control by a single Nos´ e-Hoover thermostat or independent Nos´ e-Hoover thermostats is necessary to produce the correct canonical distribution for a stiff harmonic mode. Measure-preserving dynamics of the kind discussed here is associated with equilibrium simulations (where, for example, there is a single temperature parameter T ). Stationary phase space distributions associated with non-equilibrium situations are much more com- plicated than the smooth equilbrium densities analyzed in the present paper [11, 39, 40]. Nonequilibrium simulations of heat flow could be carried out by extending the present ap- proach to multimode systems (e.g., a chain of oscillators) coupled to BK-type demons with associated NH thermostats corresponding to two different temperatures [41–43]. The techniques presented here for derivation and implementation of thermostats have been shown to be efficient and versatile. We anticipate that analogous approaches can be usefully applied to systems of classical spins coupled to both harmonic and anharmonic modes. 25