A. Algorithm for BK Dynamics In order to derive a measure preserving algorithms, the first step, following Eq. (5), is to introduce a splitting of H $^{BK}$: A measure-preserving splitting of the Liouville operator then follows from Eq. (6): Upon choosing a symmetric Trotter factorization of the BL Liouville operator based on the decomposition a measure-preserving algorithm can be produced in full generality. In practice, a choice of K$_{1}$ , K$_{2}$ , G$_{1}$ , G$_{2}$ must be made in order obtain explicit formulas. L BK = 8 ∑ α =1 L BK α (19) H BK 1 = V ( q ) (17a) H BK 2 = p 2 2 m (17b) H BK 3 = k$_{B}$T ζ (17c) H BK 4 = k$_{B}$T ξ (17d) H BK 5 = K$_{1}$ ( p$_{ζ}$ ) m$_{ζ}$ (17e) H BK 6 = K$_{2}$ ( p$_{ξ}$ ) m$_{ξ}$ . (17f) L BK 1 = − ∂V ∂q ∂ ∂p + G$_{2}$ ∂V ∂q ∂ ∂p$_{ξ}$ (18a) L BK 2 = p m ∂ ∂q + G$_{1}$ p m ∂ ∂p$_{ζ}$ (18b) L BK 3 = − k$_{B}$T ∂G$_{1}$ ∂p ∂ ∂p$_{ζ}$ (18c) L BK 4 = − k$_{B}$T ∂G$_{2}$ ∂q ∂ ∂p$_{ξ}$ (18d) L BK 5 = 1 m$_{ζ}$ ∂G$_{1}$ ∂p ∂K$_{1}$ ∂p$_{ζ}$ ∂ ∂ζ − G$_{1}$ m$_{ζ}$ ∂K$_{1}$ ∂p$_{ζ}$ ∂ ∂p (18e) L BK 6 = − G$_{2}$ m$_{ξ}$ ∂K$_{2}$ ∂p$_{ξ}$ ∂ ∂q + 1 m$_{ξ}$ ∂G$_{2}$ ∂q ∂K$_{2}$ ∂p$_{ξ}$ ∂ ∂ξ . (18f) 8