Quasi-Lie Brackets and the Breaking of Time-Translation Symmetry for Quantum Systems Embedded in Classical Baths

Many open quantum systems encountered in both natural and synthetic situations are embedded in classical-like baths. Often, the bath degrees of freedom may be represented in terms of canonically conjugate coordinates, but in some cases they may require a non-canonical or non-Hamiltonian representation. Herein, we review an approach to the dynamics and statistical mechanics of quantum subsystems embedded in either non-canonical or non-Hamiltonian classical-like baths which is based on operator-valued quasi-probability functions. These functions typically evolve through the action of quasi-Lie brackets and their associated Quantum-Classical Liouville Equations.


I. INTRODUCTION
A growing community of physicists is interested in both monitoring and controlling the time evolution of small numbers of quantum degrees of freedom (DOF) that are embedded in noisy and uncontrollable environments [1][2][3]. A specific case of such a system is encountered when the environment is classical-like in nature. This situation is one of fundamental importance because, ultimately, we and our experimental tools behave classically, at least from a coarse-grained perspective. In recent years, we have also witnessed a rising interest in nano-mechanical, opto-mechanical and other types of hybrid quantum systems . Such systems often exhibit an interplay between classical and quantum effects, allowing them to be modeled by means of hybrid quantum-classical methods.
This review deals with situations where the bath DOF are described in terms of non-canonical coordinates [102,103] or non-Hamiltonian coordinates [93][94][95], and situations where dissipation must be taken explicitly into account [45,46]. In all these cases, we will see that the operator-valued probability functions will develop new functional dependences and novel definitions of the quasi-Lie brackets will have to be introduced. In particular, we will first describe the case of a classical spin bath [91,92], as an example of a bath described by non-canonical coordinates [102,103]. It has been shown that for such a bath an off-diagonal [104] open-path [105][106][107] geometric phase [108][109][110] enters into the propagation of the quantum-classical dynamics. We will then describe the case of a non-Hamiltonian bath, which arises when the bath coordinates coupled to the quantum subsystem are also coupled to a large bath (which does not directly interact with the quantum subsystem and whose detailed dynamics is not of interest). In such cases, the secondary bath acts as a thermal reservoir and can be described either by means of stochastic processes [111] (e.g., Langevin dynamics [46]), or by means of non-Hamiltonian fictitious coordinates acting as deterministic thermostats (e.g., the Nosé-Hoover thermostat [112,113]). Both Langevin and Nosé- Nosé-Hoover dynamics is defined solely in terms of a quasi-Lie bracket [43,44]. Instead, explicit dissipative dynamics requires that diffusive terms be added to the bracket.
The structure of this review is as follows. In Section II, we illustrate the algebraic approach used to formulate the dynamics of a quantum subsystem embedded in a classical-like environment with canonically conjugate coordinates. In Section III, we show how this formalism can be generalized to the case of a bath described by non-canonical variables, namely a collection of classical spins.
Here, we will also show how an off-diagonal open-path geometric phase enters into the time evolution of the operator-valued quasi-probability function of the system. In Section IV, we show how the formalism allows us to also treat stochastic classical-like baths undergoing Langevin dynamics.
Finally, in Section V, we shed light on the quasi-Lie algebra established by the quantum-classical brackets and show how their antisymmetric structure is exploited to achieve thermal control of the bath DOF by means of deterministic thermostats such as the Nosé-Hoover and Nosé-Hoover chain thermostats. Our conclusions and perspectives are given in Section VI.

II. QUASI-LIE BRACKETS AND HYBRID QUANTUM-CLASSICAL SYSTEMS
Classical and quantum dynamics share the same algebraic structure [125,126], which is realized by means of Poisson brackets in the classical case and commutators in the quantum theory. Poisson brackets have a symplectic structure that is easily represented in matrix form [103,127]. Both Poisson brackets and commutators define Lie algebras. In terms of commutators, a Lie algebra possesses the following properties: [c, where c is a so-called c-number andχ j , j = 1, 2, 3 are quantum operators. In order to have a Lie algebra, together with Equations (1)-(3), the Jacobi relation must also hold The time-translation invariance of the commutator algebra follows from the Jacobi relation, which therefore states an integrability condition. IfĤ is not explicitly time-dependent, the antisymmetry of the commutator (1), arising from the antisymmetry of the symplectic matrix B, ensures that the energy is a constant of motion: dĤ/dt = iLĤ = 0. Energy conservation under time-translation is a fundamental property shared by the algebra of Poisson brackets and the algebra of commutators that is in agreement with Noether theorem. Now, let us consider a hybrid quantum-classical system, in which the quantum subsystem, described by a few canonically conjugate operators (q,p) =x is embedded in a classical bath with many DOF, described by many canonically conjugate phase space coordinates, X = (Q, P ). We will assume that the Hamiltonian of this hybrid system has the form where m and M are the masses of the subsystem and bath DOF, respectively, and V W is the potential energy describing the interactions among the subsystem DOF, among the bath DOF, and between these two sets of DOF. The last equality on the right-hand side of Equation (5) defines the adiabatic Hamiltonianĥ W (Q) of the system. It has been known for many years that the statistical mechanics of such hybrid quantum-classical systems may be formulated in terms of an operator-valued quasi-probability functionŴ (X, t) [27][28][29][30][31][32]. Specifically, the statistical average of hybrid quantum-classical operators, representing a dynamical property of the system, may be calculated according to where Tr ′ denotes the partial trace involving a complete set of states of the quantum subsystem.
The QCLE in Equation (7) is founded upon a quasi-Lie bracket, which we may write explicitly where D is the antisymmetric matrix operator defined in Equation (8). However, in contrast to the Lie brackets of quantum and classical mechanics, the quasi-Lie bracket defined in Equation (10) violates the Jacobi relation (4): The failure of the Jacobi implies that the algebra of quasi-Lie brackets is not invariant under time-translation. For example, it can be generally proven that e itL D [χ 1 (X, 0),χ 2 (X)] = e itL Dχ 1 (X), e itL Dχ 2 (X) .
On the other hand, the quasi-Lie bracket conserves the energy e itL DĤ W (X) =Ĥ W (X). Hence, the dynamics generated by the QCLE displays energy conservation and lack of time-translation invariance of the bracket algebra. The situation is surprising because one does not expect a broken time-translation invariance symmetry in an isolated system. However, while a total hybrid quantum-classical system is closed from the point of view of energy conservation, the quasi-Lie bracket describes the irreversible transfer of quantum information from the subsystem to the classical DOF, which acquire a quantum character as the time flows. In this sense, one can heuristically argue that the lack of time-translation invariance or the algebra is a mere consequence of the open dynamics of the quantum subsystem.
The quantum statistical state of the system is described by the density matrix (or statistical operator)ρ(t). The time dependence of the density matrix is dictated by the QLE: where [..., ...] denotes the commutator, and B is the symplectic matrix [103,127]: The average of an operatorχ defined on the same Hilbert space of the system is calculated by where Tr denotes the trace operation. Now, in order to derive a classical-like description of the bath, one introduces the partial Wigner transform of the density matrixρ over theX's: The symbolŴ denotes an operator-valued Wigner function (also known as the partially-Wigner transformed density matrix), which is both an operator in the Hilbert space of theq's and a function of the bath coordinates X. The partial Wigner transform of an arbitrary operatorχ is analogously given byχ Taking the partial Wigner transform of Equation (16) leads to the expression for the average ofχ given in Equation (6). The partial Wigner transform of the Hamiltonian in Equation (13) is given in Equation (5).
Upon taking the partial Wigner transform of the QLE, Equation (14), and truncating the resulting equation after first order inh, one arrives at the QCLE where the last equality defines the quantum Liouville . To arrive at Equation (19), we have used the partial Wigner transform of a product of operators, and truncated the exponential after first order inh, i.e., It should be noted that Equation (21) is exact for Hamiltonians with quadratic bath terms and bilinear coupling between thex and X DOF. In Ref. [48], it is shown how the linear expansion can be performed in terms of the parameter µ = m/M , which is small in cases where the bath DOF are much more massive than those of the subsystem. Equation (19) is exactly equivalent to Equation (7).

B. Integration Algorithm
A number of algorithms, which depend on the basis representation, exist for approximately solving the QCLE [51, 52, 55-57, 61, 79, 114-124]. Herein, we illustrate the so-called Sequential Short-Time Propagation (SSTP) algorithm [79,114], which offers a good compromise between accuracy and simplicity of implementation. The SSTP algorithm is based on the representation of the QCLE in the adiabatic basis, which is defined by the eigenvalue equation The representation of the QCLE in the adiabatic basis is sketched in Appendix A. In the adiabatic basis, the QCLE is given by Equation (A1) and the quantum-classical Liouville super-operator matrix elements are given in Equation (A4).
To derive the SSTP algorithm, we divide the time interval t into n equal small steps τ = t/n.
If one is able to calculate the propagation over a single τ , the dynamics over the whole interval can be reconstructed by sequential iteration of the procedure. Let us then consider the quantumclassical propagator over a small step τ for the matrix elements of the operator-valued quasiprobability functionŴ (X) in the adiabatic basis. Such a propagator is written as On the right-hand side of Equation (23), we have introduced ω αα ′ , the Bohr frequency defined in Equation (A3), iL αα ′ is a classical-like Liouville operator, defined in Equation (A5), and T αα ′ ,ββ ′ is the transition operator defined in Equation (A7). The SSTP dynamics of the matrix elements ofŴ (X, t) is given by When τ is infinitesimal, the right-hand side of Equations (23) and (24), become essentially equal to the left-hand side, as can be seen from the Dyson identity [114].
The transition operator is purely off-diagonal. Its action generates quantum transitions in the subsystems and changes the bath momenta accordingly. Upon setting the transition operator to zero, we obtain an adiabatic expression for the propagator. If the non-adiabatic effects are not too strong, they may be treated in a perturbative fashion by sampling the action of the transition operator in a stochastic fashion. Typically, researchers have used [44, 63, 65-76, 79, 83, 86, 88-92, 114, 119-123] the following expressions for the probabilities of making a transition (jump) and not-making a transition, respectively: Another important technical ingredient of the algorithm is the approximation of the transition operator in Equation (A7) with its momentum-jump form: whered αβ is the normalized coupling vector. Within the momentum-jump approximation [78,79], the action of the transition operator on the bath momenta can be easily obtained in closed form: Considering Equations (6) and (24), together with its SSTP implementation just described, one can see that the solution of the QCLE can be obtained from an ensemble of classical-like trajectories, where each trajectory (whose initial conditions arise from a Monte Carlo sampling [128] of the X's), involves deterministic evolution segments on a given adiabatic energy surfaces interspersed with stochastic quantum transitions, caused by the momentum-jump operator in Equation (27).
The SSTP algorithm [79,114] maps the calculation of averages through the QCLE (19) onto a stochastic process. It is a hybrid Molecular Dynamics/Monte Carlo procedure suffering from two main problems. The first is given by the momentum-jump approximation, which is not valid in general. One can avoid this approximation by devising different integration schemes, but usually at the expense of other approximations [124]. The second problem is not just associated with the SSTP algorithm, but it is common to all Monte Carlo approaches to the calculation of  [129][130][131][132]. However, in contrast to the SLE, the QCLE is a deterministic equation that explicitly takes into account all the DOF of the system without approximating the memory of the total hybrid quantum-classical system. The stochastic process only enters through the specific hybrid Molecular Dynamics/Monte Carlo implementation provided by the SSTP algorithm. Indeed, a recently proposed scheme of integration [124] does not involve any stochastic process whatsoever.

III. CLASSICAL SPIN BATHS
Contrary to what some books in quantum mechanics state (in the authors's knowledge, an exception is Schulman's book [133]), the concept of spin can be defined in an entirely classical way [133][134][135][136][137]. In practice, spinors provide a more fundamental representation of the rotation group than that given by tensors [133][134][135][136][137]. Hence, one can think of a collection, e.g., a bath, of DOF comprising classical spinors (or, for brevity, spins): a classical spin-bath. An example of a classical spin baths is given by the Classical Heisenberg Model [138], whose Hamiltonian is where S I are N classical vectors obeying the constraint for I = 1, ..., N , and the C a IJ are coupling constants. However, since the generalization to baths with many spins is straightforward, in the following, we will illustrate the theory using a bath comprising a single classical spin. Consider a classical spin vector S, with components S a , a = x, y, z, and Hamiltonian H S (S). Let us define the spin gradient as ∇ S = ∂/∂S, which in terms of the spin components is written as ∇ S a = ∂/∂S a , with a = x, y, z. The equations of motion of the spin are then written asṠ where One can also adopt the compact form B S ab = c=x,y,z ǫ abc S c and a, b = x, y, z of the antisymmetric matrix B S , where ǫ abc is the Levi-Civita pseudo-tensor. The Casimir C 2 = S · S is preserved by the equations of motion (31), independently of the form of the spin Hamiltonian H S (S). In addition, the dynamics has a zero phase space compressibility κ S = ∇S ·Ṡ = 0. The classical phase space flow of the spin is defined through the non-canonical bracket where A = A(S) and B = B(S) are arbitrary functions of the spin DOF.
Consider now the hybrid quantum-classical Hamiltonian of a quantum subsystem coupled to the classical spinĤ where We next set out to represent Equation (35) in the adiabatic basis |α; S defined by the eigenvalue equationĥ It should be noted that, in contrast to the case of canonically conjugate phase space coordinates which depends only on the positions Q and not on the conjugate momenta P , this adiabatic basis depends on all the non-canonical spin coordinates S. In this basis, Equation (35) becomes where ω αα ′ = E α (S) − E α ′ (S)/h is the Bohr frequency. Defining the spin coupling vector one finds the two identities (40) and (41), the spin-bath QCLE may be rewritten where we have defined the classical-like spin-Liouville operator with the average adiabatic Hamiltonian The transition operator for the spin bath is given by The limit d S αα → 0 of the spin transition operator in Equation (45) provides the form of the standard transition operator for canonical conjugate coordinates, given in Equation (A7). Finally, because of the spin nature of the bath, one finds a higher order transition operator (which does not appear in the case of canonical conjugate bath coordinates): The adiabatic limit of the spin-bath QCLE in (42) can be taken by setting to zero the offdiagonal elements of d αα ′ , which appear in the operators in Equations (45) and (46). This is physically reasonable whenever the coupling between the different adiabatic energy surfaces is negligible. One obtains The geometric phase has been introduced exploiting the purely imaginary character of d S αα . Similarly, the higher order transition operator becomes Putting everything together, the adiabatic approximation of the spin-bath QCLE may be written as In Equation (50), the phase ω αα ′ has a dynamical nature while the phase φ S αα is of a geometric origin and it can be considered an instance of the famous Berry phase [108][109][110]. Interestingly, Equation (35) predicts that the geometric phase φ S αα can be non-zero also for open paths of the classical spins of the bath (open-path Berry phases were discussed in Ref. [105]). Moreover, the phase factor φ S αα −φ S α ′ α ′ is purely off-diagonal (off-diagonal Berry phases for environments described by canonically conjugate variables were discussed in Refs. [104,106,107]). It is worth mentioning that the geometric phase φ S αα is predicted also for non-adiabatic dynamics. When the total Hamiltonian is time-independent, as the one in Equation (34), the adiabatic evolution of the matrix elements of the spin-bath operator-valued quasi-probability function, given by Equation (50), can be rewritten as Using the Dyson identity, one can obtain the following form forŴ S (S, t) in terms of the adiabatic propagator: Equation (52) provides a convenient starting point for devising numerical integration schemes based on the SSTP propagation scheme [114].
In Ref. [92], the following model Hamiltonian was considered: where Ω, c 1 , and c 2 are real parameters, b is the z component of the magnetic field B = (0, 0, b), while σ = (σ x , σ y , σ z ) is a vector having the Pauli matrices σ x , σ y , and σ z as components. The SSTP algorithm was applied to Equation (52) and the action of the classical like Liouville operator S was evaluated using time reversible integration algorithms based on the symmetric break-up of the Liouville propagator [139][140][141].

IV. STOCHASTIC CLASSICAL BATHS
Consider a quantum-classical system comprising a quantum subsystem and a classical environment whose classical phase space coordinates are partitioned into two sets: one set X = (Q, P ) interacts directly with the quantum subsystem while the second set X ′ = (Q ′ , P ′ ) interacts only with the coordinates X (and therefore is not directly coupled to the quantum subsystem). We assume that the detailed dynamics of the coordinates X ′ is not interesting: their function is just that of working as a thermal bath, leading to dissipative dynamics [45].
An equation of motion for the hybrid quantum-classical system composed of the quantum subsystem and the classical DOF X only has been derived using projection operator methods [45].
It takes the form, where ∇ P = ∂/∂P , ζ is the friction constant, k B is the Boltzmann constant, and T is the temperature of the bath. The Hamiltonian in Equation (55) is defined in Equation (5). However, in the present case, we must interpret V W (q, Q) as the potential of mean force arising from the average over the primed bath variables Q ′ . The Liouville operator iL D , defined on the right-hand side of Equation (55), determines the dissipative dynamics of the system. This Fokker-Planck-like operator and the potential of mean force make the dissipative quantum-classical Liouville operator in Equation (55) different from that describing an isolated quantum-classical system [48]. In particular, the term ζ − → ∇ P (P/M ) + k B T − → ∇ P directly breaks the time-translation symmetry leading to diffusive motion and energy dissipation.
The dissipative Liouville operator can be written in the adiabatic basis as where we have defined the Kramers operator as The quantum-classical average of any operator or dynamical variableχ(X) can be written as where iL DB β ′ βα ′ α is the backward operator, defined as The backward Kramers iL KB αα ′ operator is written as According to the classical theory of random processes [111], the time evolution under the backward Kramers operator iL KB αα ′ ββ ′ can be unfolded it via an average over realizations of stochastic Langevin trajectories. In such a picture, the classical trajectory segments obey the Langevin equations of motion,Q where R(t) is a Gaussian white noise process with the properties, To Equations (61) and (62), one can associate a time-dependent Langevin-Liouville operator and a time-ordered propagator In order to generate the stochastic Langevin trajectories, we can use a total time-dependent Langevin-Liouville super-operator and the associated propagator Within such a Langevin picture, the quantum-classical average of any operatorχ(X) can be calculated as where the over-line denotes an average over an ensemble of stochastic Langevin trajectories.
Since they are independent from each other, the order in which the average over phase space and the average over the stochastic Langevin process are performed can be permuted. Hence, one can write Equation (70) allows one to calculate averages in a quantum-classical dissipative system as phase space weighted averages over many Langevin trajectories.
In Ref. [46], a quantum subsystem with two energy levels interacting with a dissipative classical quartic oscillator was considered. The Hamiltonian of the hybrid quantum-classical system readŝ where V q (Q) = a 4 R 4 − b 2 R 2 , Ω, a, b, and γ 0 are real parameters, M is the mass of the quartic oscillator, andσ x andσ z are Pauli matrices.
The calculation of quantum-classical averages using the dynamics defined by the time-dependent Langevin-Liouville propagator U L ss ′ (t) in Equation (68) is no more complicated than that for deterministic quantum-classical dynamics. The momentum-jump approximation [78,79] and a simple generalization of the SSTP algorithm [79,114] to the time dependent propagator were used in Ref. [46]. The explicitly time-dependent propagator U L ss ′ (t) must be defined as a time ordered product. A simple way to achieve that is to employ the decomposition scheme devised by Suzuki [142].
Details of the numerical procedures are found in Ref. [46] V. NON-HAMILTONIAN DYNAMICS IN THERMAL BATHS By exploiting the antisymmetric structure of the quantum-classical commutator, arising from the matrix operator D given in Equation (8), one can impose the thermodynamic constraints of constant temperature on the classical-like DOF [43,44]. Following Refs. [93][94][95] As in the classical case, the Nosé variables are where Q η and P η are the Nosé coordinate and momentum. The Nosé quantum-classical Hamiltonian is obtained by adding the Nosé kinetic energy P 2 η /2M η and potential energy N k B T Q η toĤ W in Equation (5) where M η is the Nosé inertial parameter, k B is the Boltzmann constant, T is the constant temperature, and N is the number of Q coordinates. Using the matrix B N in Equation (B4), the classical phase space quasi-Hamiltonian bracket of two variables A 1 and A 2 can be defined as The explicit form of the matrix operator, which defines the quantum-classical bracket and the law of motion through Equation (19), is then given by The Nosé-Hoover QCLE for the operator-valued quasi-probability functionŴ N (X N , t) is given by The presence of the term −κ N (X N )W N (X N , t) in the left-hand side of Equation (76) derives from the passage from the Heisenberg to the Schrödinger picture, as it is explained in Appendix C.
Upon considering the term in the right-hand side of (76), one obtainŝ Finally, using the above result, the Nosé-Hoover QCLE reads d dtŴ In the adiabatic states defined in Equation (22), Equation (78) reads where We have used the definition of the Bohr frequency ω αα ′ in Equation (A3) and of the transition The existence of the stationary operator-valued Nosé quasi-probability functionŴ N,e (X N ) is discussed in Appendix C.

A. Nosé-Hoover Chain Thermal Baths
The Nosé-Hoover thermostat suffers from lack of ergodic dynamics when the bath has high frequencies of motion. The Nosé-Hoover chain [143] is a more general non-Hamiltonian thermostat that solves the ergodicity problems suffered by the standard Nosé-Hoover thermostat in the case of stiff variables. The Nosé-Hoover chain thermostat can also be formulated in a quantum-classical framework with minimal changes with respect to what is shown in Section V. To this end, considering for simplicity a chain of just two thermostat coordinates, one can define the classical phase space point as where M η 1 and M η 2 are the inertial parameters of the thermostat variables. As shown in Ref. [93,94], one can define an antisymmetric matrix The matrix B NHC can be used to define the quasi-Hamiltonian bracket according to Equation (9). The Nosé-Hoover chain classical equations of motion in phase space [93] are then given Quantum-classical dynamics is then introduced using the matrix super-operator The quantum-classical equations of motion can then be written as The equations of motion can be represented using the adiabatic basis obtaining the Liouville super-operator where with F Qη 2 = (P 2 The proof of the existence of stationary density matrix in the case of Nosé-Hoover chains follows the same logic of the simpler Nosé-Hoover case. In the adiabatic basis, the density matrix stationary up to order bar has the same form as that given in Equations (C19) and (C21). One has just to replace Equation (C19) for the order zero term with with an obvious definition of Z NHC .

VI. CONCLUSIONS AND PERSPECTIVES
In this review, we discussed how to mathematically describe the dynamics and statistical mechanics of quantum subsystems embedded in classical baths. The formalism is founded on an operator-valued quasi-probability function evolving through a QCLE defined in terms of a quasi-Lie bracket. It is worth emphasizing that the QCLE is a fully deterministic equation that takes into account explicitly all the DOF of the system, i.e., it describes the quantum and classical DOF of the total hybrid system. Hence, the QCLE generates a unitary dynamics, conserving both the system's probability and energy. However, the time-translation invariance of the quasi-Lie bracket algebra is broken. This situation is surprising: one does not expect a broken time-translation invariance symmetry in an isolated system when all its degrees of freedom are taken into account.
This can be seen as a signature of the effect of the classical bath on the quantum subsystem, and of the back-reaction of the subsystem onto the bath. In other words, the total hybrid system is closed from the point of view of energy and probability conservation but, because of the above mentioned back-reaction, it is also open: the quasi-Lie bracket describes the irreversible transfer of quantum information onto the classical DOF. We also reviewed how the hybrid quantum-classical theory can be derived from a partial Wigner transform and a semiclassical limit of the QLE only in the case when the bath is described by canonically conjugate coordinates. After this, we discussed how to treat quantum subsystems embedded in both non-canonical and non-Hamiltonian bath. In all cases, the mathematical object representing the state of the system is an operatorvalued quasi-probability function that depend on the coordinates of the bath and whose equation These methods scales favorably in terms of bath DOF but, to date, have been limited to relatively short time intervals and Markovian systems. When the dynamics is non-Markovian, the memory function, i.e., the autocorrelation function of the random force [3,111], cannot be approximated by a delta function. The memory function of the bath can be expected to become more and more different from a delta function as the quantum character of the bath becomes more pronounced (for example, at low temperature) and as the subsystem-bath coupling grows in strength.
The QCLE discussed herein constitutes an approach to open quantum system dynamics (in the case of hybrid quantum-classical systems) that is both distinct and complementary to that given by master equations [3,111]. The QCLE-based approach to quantum dynamics in classical baths has proven to be successful in modeling a variety of quantum processes in the condensed phase. Nevertheless, the currently algorithms also present significant challenges, necessitating the need for further improvements and developments. In light of the above, we hope that this review will attract the attention of a broad community of researchers and spur further work along this direction. In addition to further algorithm developments, we are interested in broadening the scope of applications studied by this approach. For example, based on preliminary results, we believe that this approach can be successfully applied to studying the interplay between quantum and classical fluctuations in hybrid nanoscale devices.
A. * Funding: A.S. and R.G. acknowledge support by research funds in memory of Francesca Palumbo, difc 3100050001d08+, University of Palermo.

Appendix A: Representation in the Adiabatic Basis
In the adiabatic basis, Equation (19) reads where are the matrix elements of the density matrix. Upon defining the Bohr frequency as the Liouville super-operator may be written as We have also introduced a classical-like Liouville operator where is the Hellmann-Feynman force.
In Equation (A4), the transition operator T αα ′ ,ββ ′ is defined as In turn, the transition operator is defined in terms of the shift vector and of the coupling vector Appendix B: The Nosé-Hoover Thermostat The Nosé-Hoover thermostat was originally formulated in Refs. [112,113]. Herein, we follow Refs. [93][94][95]. The Hamiltonian of the subsystem with phase space coordinates (R, P ) is: where V (R) is the potential energy. One can introduce an extended system comprised by the coordinates of the original subsystem augmented with the additional variables Q η and conjugate momentum P η . The dimension of such an extended phase space is obviously 2N + 2, which is computationally tractable whenever N is computationally tractable. As a consequence, the phase space point of the extended system is while the energy reads: where M η is a fictitious mass associated with the additional degree of freedom, k B is Boltzmann constant, and T the bath constant temperature. In order to define time evolution, we abandon the Hamiltonian structure of the theory. To this end, using the general formalism of Refs. [93][94][95], we introduce the antisymmetric matrix: so that Nosé's equations of motion can be written aṡ where the first equality on the right-hand side of Equations (B5) introduces the Nosé bracket, while the extended phase space gradient is denoted as ∇ N J = ∂/∂X N J . We remark here that the Nosé bracket does not satisfy the Jacobi relation [93][94][95], and thus defines a quasi-Hamiltonian algebra. The Liouville equation for the Nosé distribution function is where the compressibility of the phase space reads: As implied by Equation (B7), Nosé's phase space flow has a non-zero compressibility (however, this does not always occur for a quasi-Hamiltonian dynamics). In terms of the Nosé bracket, the equilibrium Liouville equation for Nosé distribution function reads: By direct substitution, one can verify that the solution of Equation (B8) is: where w is defined by the equation dw/dt = κ N . Equations (B5) can be written explicitly in the form:Ṙ In order to write explicitly the Nosé distribution function, it is useful to introduce the following extended phase space function: Using the equations of motion, one finds which is related to the compressibility by At this point, we have all the ingredients that are needed to prove that extended phase space averages of functions of the subsystem coordinates (R, P ) can be written as canonical averages.
We start by considering The integral is calculated by using the identity where the sum runs over the zeros Q η 0 of f (Q η ). Upon identifying f (Q η ) = E − H N , one gets with the above results, the integral over Q η becomes a trivial Gaussian integral over P η : Finally, one obtains: Hence, averages in the canonical ensemble can be calculated by letting the trajectories evolve according to Nosé's dynamics.
The quasi-Hamiltonian Nosé dynamics is a well-established tool of molecular dynamics simulations. In practice, it is adopted whenever one wants to calculate dynamical properties at constant temperature and/or study phase transitions. Discussions and pointers to the relevant literature on the subject can be found in Ref. [128].

Appendix C: Stationary Operator-Valued Nosé Quasi-Probability Function
The quantum average of any operatorŴ N (X N ), in a dynamics where the temperature of the X degrees of freedom is controlled by the Nosé-Hoover thermostat can be calculated as The action of exp iL N t can be transferred fromχ(X N ) toŴ N (X N ) by using the cyclic invariance of the trace and integrating by parts the terms coming from the classical brackets. One can write The action of iL N on an arbitrary operatorχ(X N ) is defined by when integrating by parts the right-hand side, one obtains a term proportional to the compress- As a result, the quantum Liouville operator, partially depending on phase space variables, is non-Hermitian The average value can then be written as The operator-valued Nosé quasi-probability function evolves under the equation: The stationary operator-valued Nosé quasi-probability functionŴ N,e is defined by To find the explicit expression, one can follow Ref. [42]: the density matrix is expanded in powers ofhŴ and an explicit solution in the adiabatic basis is searched for. On such a basis, the Nosé-Liouville operator is expressed by Equation (80) and the Nosé Hamiltonian is given by One obtains an infinite set of equations corresponding to the various power ofh ) * , T αα,ββ ′ = T * αα,β ′ β and the fact that T αα,ββ = 0 when a real basis is chosen, it is useful to re-write Equation (C11) in the form One has [93] (−iL N αα − κ N ) † = iL N αα . The right-hand side of this equation is expressed by means of the generalized bracket in Equation (74): H α N and any general function f (H α N ) are constants of motion under the action of iL N αα . The phase space compressibility κ N associated with the generalized bracket in the case of Nosé dynamics is where N is the number of classical momenta P in the Hamiltonian.
To ensure that a solution to Equation (C12) exists, one must invoke the theorem of Fredholm alternative, requiring that the right-hand side of Equation (C12) is orthogonal to the null space of (iL N αα ) † = −iL N αα − κ N [144]. The null-space of this operator is defined by the equation (iL N αα + κ N )G α (X) = 0, with G α (X) = f (H N α ) exp(−w N α ). Hence, the condition to be satisfied is The fact that 2 exp(−w α )Re T αα,ββ ′ W (C16) Equations (C15) and (C16) allow one to calculate W N,e αα ′ to all orders inh once W N,e,(0) αα ′ is given.
This order zero term is obtained by the solution of (iL N αα + κ N )W N,e,(0) αα = 0. All higher order terms are obtained by the action of H N αα ′ , the imaginary unit i and T αα ′ ββ ′ (involving factors of d αα ′ , P and derivatives with respect to P . Hence, one can conclude that functional dependence of W (0)αα Ne on the Nosé variables Q η and P η is preserved in higher order terms W where Z N is One can now prove that, when calculating averages of quantum-classical operators depending only on physical phase space variables, G α (R, P ), the canonical form of the stationary density is obtained. It can be noted that it will suffice to prove this result for thē h 0 term since, as discussed before, the differences with the standard case are contained therein.
Indeed, when calculating considering the integral of the delta function over Nosé variables, one has where the property δ(f (s)) = [df /ds] −1 s=s 0 δ(s − s 0 ) has been used (s 0 is the zero of f (s)).