A. Algorithm for BKNH dynamics The Hamiltonian can be split as The measure-preserving splitting [19] of the Liouville operator yields 13 H BKNH 1 = V ( q ) (31a) H BKNH 2 = p 2 2 m (31b) H BKNH 3 = k$_{B}$T ζ (31c) H BKNH 4 = k$_{B}$T ξ (31d) H BKNH 5 = K$_{1}$ ( p$_{ζ}$ ) m$_{ζ}$ (31e) H BKNH 6 = K$_{2}$ ( p$_{ξ}$ ) m$_{ξ}$ (31f) H BKNH 7 = p 2 η 2 m$_{η}$ (31g) H BKNH 8 = 2 k$_{B}$T η (31h) L$_{α}$ = B BKNH ij ∂H BKNH α ∂x$_{j}$ ∂ ∂x$_{i}$ (32) L BKNH 1 = − ∂V ∂q ∂ ∂p + G$_{2}$ ∂V ∂q ∂ ∂p$_{ξ}$ (33a) L BKNH 2 = p m ∂ ∂q + G$_{1}$ p m ∂ ∂p$_{ζ}$ (33b) L BKNH 3 = − k$_{B}$T ∂G$_{1}$ ∂p ∂ ∂p$_{ζ}$ (33c) L BKNH 4 = − k$_{B}$T ∂G$_{2}$ ∂q ∂ ∂p$_{ξ}$ (33d) L BKNH 5 = 1 m$_{ζ}$ ∂G$_{1}$ ∂p ∂K$_{1}$ ∂p$_{ζ}$ ∂ ∂ζ − G$_{1}$ m$_{ζ}$ ∂K$_{1}$ ∂p$_{ζ}$ ∂ ∂p + p$_{ζ}$ m$_{ζ}$ ∂K$_{1}$ ∂p$_{ζ}$ ∂ ∂p$_{η}$ (33e) L BKNH 6 = − G$_{2}$ m$_{ξ}$ ∂K$_{2}$ ∂p$_{ξ}$ ∂ ∂q + 1 m$_{ξ}$ ∂G$_{2}$ ∂q ∂K$_{2}$ ∂p$_{ξ}$ ∂ ∂ξ + p$_{ξ}$ m$_{ξ}$ ∂K$_{2}$ ∂p$_{ξ}$ ∂ ∂p$_{η}$ (33f) L BKNH 7 = p$_{η}$ m$_{η}$ ∂ ∂η − p$_{η}$ m$_{η}$ p$_{ζ}$ ∂ ∂p$_{ζ}$ − p$_{η}$ m$_{η}$ p$_{ξ}$ ∂ ∂p$_{ξ}$ (33g) L BKNH 8 = − 2 k$_{B}$T ∂ ∂p$_{η}$ . (33h)