can be calculated exactly. Details of this way of visualizing the phase space sampling have already been given in [14, 15]. Figure 2, displaying the comparison between the theoretical and the calculated value of the radial probability in phase space, clearly shows that the BK dynamics is not able to produce canonical sampling. A look at the inset of Fig. 2, showing the phase space distribution for the harmonic mode, also immediately shows that the dynamics is not ergodic. The same analysis has been carried out for BKNH and BKNHC phase space flows, and these are displayed in Fig 3 and Fig 4, respectively. Within numerical errors, both BKNH and BKNHC thermostats are able to produce the correct canonical distribution for the stiff harmonic modes. Introduction of a single, global Nos´ e-type variable in the BKNH thermostat effectively introduces additional coupling between the two demon variables. The effectiveness of the BKNH thermostat is consistent with our findings (results not discussed here) that introduc- tion of explicit coupling between demons in BK thermostat dynamics also leads to efficient thermostatting of the harmonic oscillator. 24