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

Abstract

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, or through quasi-Lie brackets augmented by dissipative terms. Quasi-Lie brackets possess the unique feature that, while conserving the energy (which the Noether theorem links to time-translation symmetry), they violate the time-translation symmetry of their algebra. This fact can be heuristically understood in terms of the dynamics of the open quantum subsystem. We then describe an example in which a quantum subsystem is embedded in a bath of classical spins, which are described by non-canonical coordinates. In this case, it has been shown that an off-diagonal open-bath geometric phase enters into the propagation of the quantum-classical dynamics. Next, we discuss how non-Hamiltonian dynamics may be employed to generate the constant-temperature evolution of phase space degrees of freedom coupled to the quantum subsystem. Constant-temperature dynamics may be generated by either a classical Langevin stochastic process or a Nosé–Hoover deterministic thermostat. These two approaches are not equivalent but have different advantages and drawbacks. In all cases, the calculation of the operator-valued quasi-probability function allows one to compute time-dependent statistical averages of observables. This may be accomplished in practice using a hybrid Molecular Dynamics/Monte Carlo algorithms, which we outline herein.

1. 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 [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. Such systems often exhibit an interplay between classical and quantum effects, allowing them to be modeled by means of hybrid quantum-classical methods.

It has been known for a long time, that the dynamics and statistical mechanics of a quantum subsystem coupled to classical-like DOF can be formulated in terms of operator-valued quasi-probability functions in phase space [27,28,29,30,31,32]. For example, the dynamics of nano-mechanical oscillators has been previously described by one of the authors in terms of operator-valued quasi-probability functions [33]. Such functions evolve through quasi-Lie brackets [34,35,36,37,38,39,40,41,42,43], which can also be augmented by dissipative terms when the energy is not conserved [44,45]. When the bath is described by canonically conjugate variables (and only in this case), a hybrid quantum-classical formalism may be derived. Starting from a fully quantum representation of the subsystem and bath DOF, one can perform a partial Wigner transform [46] (over the bath DOF) and then take its semiclassical limit [47]. The resulting equation of motion is commonly known as the quantum-classical Liouville equation (QCLE) [48,49,50,51,52,53,54,55,56,57,58,59,60]. The QCLE has been used to study a wide variety of problems [61,62,63,64,65,66,67,68,69,70,71,72,73,74,75] and a number of in-depth reviews on the basic formulation of the theory exist [76,77,78,79,80,81,82,83,84,85,86,87,88,89]. The mathematical structure underlying the QCLE is dictated by a quasi-Lie bracket [42,43,90,91]. Quasi Lie brackets are known within the community of classical molecular dynamics simulators as non-Hamiltonian brackets [92,93,94]. Mathematicians have also studied very similar structures known as almost Poisson brackets or quasi-Lie algebras [95,96,97,98,99]. It is interesting to note that the quasi-Lie (or non-Hamiltonian) structure of the QCLE [30,31,34,35,36,37,38,39,40,41,42,43] has both favorable and unfavorable aspects associated with it. Because the antisymmetry of the quasi-Lie bracket ensures energy conservation, one is able to verify the stability of numerical integration algorithms. However, because the quasi-Lie algebra is not invariant under time translation, the initially classical DOF acquire a quantum character as time flows, implying that one never has a true dynamical theory of quantum and classical DOF but only an approximated dynamics of a full quantum system [100]. This is somewhat paradoxical because energy conservation is linked to time-translation symmetry through the Noether theorem; nevertheless, quasi-Lie brackets break the time-translation symmetry of the algebra (which can be seen as a signature of the effect of the classical bath on the quantum subsystem).

This review deals with situations where the bath DOF are described in terms of non-canonical coordinates [101,102] or non-Hamiltonian coordinates [92,93,94], and situations where dissipation must be taken explicitly into account [44,45]. 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 [90,91], as an example of a bath described by non-canonical coordinates [101,102]. It has been shown that for such a bath an off-diagonal [103] open-path [104,105,106] geometric phase [107,108,109] 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 [110] (e.g., Langevin dynamics [45]), or by means of non-Hamiltonian fictitious coordinates acting as deterministic thermostats (e.g., the Nosé–Hoover thermostat [111,112]). Both Langevin and Nosé–Hoover deterministic time evolutions are examples of non-Hamiltonian dynamics. However, only Nosé–Hoover dynamics is defined solely in terms of a quasi-Lie bracket [42,43]. Instead, explicit dissipative dynamics requires that diffusive terms be added to the bracket.

The quantum-classical equations of motion herein discussed can be implemented in silico using a variety of simulation algorithms [78,113,114,115,116,117,118,119,120,121,122,123]. We will sketch out one such integration algorithm, which unfolds the quantum-classical dynamics of the operator-valued quasi-probability function in terms of piecewise-deterministic trajectories evolving on the adiabatic energy surfaces of the system under study [78,113].

The structure of this review is as follows. In Section 2, 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 3, 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 4, we show how the formalism allows us to also treat stochastic classical-like baths undergoing Langevin dynamics. Finally, in Section 5, 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 6.

2. Quasi-Lie Brackets and Hybrid Quantum-Classical Systems

Classical and quantum dynamics share the same algebraic structure [124,125], 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 [102,126]. Both Poisson brackets and commutators define Lie algebras. In terms of commutators, a Lie algebra possesses the following properties:

$$\begin{array}{ccc}\hfill [{\widehat{\chi }}_1,{\widehat{\chi }}_2]& =& -[{\widehat{\chi }}_2,{\widehat{\chi }}_2],\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \left[{\chi }_1{\widehat{\chi }}_2,{\widehat{\chi }}_3\right]& =& {\widehat{\chi }}_1[{\widehat{\chi }}_2,{\widehat{\chi }}_3]+[{\widehat{\chi }}_1,{\widehat{\chi }}_3]{\widehat{\chi }}_2,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \left[c,{\widehat{\chi }}_j\right]& =& 0,\hfill \end{array}$$

(3)

where c is a so-called c-number and ${\widehat{\chi }}_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

$$\mathcal{J}=[{\widehat{\chi }}_1,[{\widehat{\chi }}_2,{\widehat{\chi }}_3]]+[{\widehat{\chi }}_3,[{\widehat{\chi }}_1,{\widehat{\chi }}_2]]+[{\widehat{\chi }}_2,[{\widehat{\chi }}_3,{\widehat{\chi }}_1]]=0.$$

(4)

The time-translation invariance of the commutator algebra follows from the Jacobi relation, which therefore states an integrability condition. If $\widehat{H}$ is not explicitly time-dependent, the antisymmetry of the commutator (1), arising from the antisymmetry of the symplectic matrix $\mathcal{B}$, ensures that the energy is a constant of motion: $d\widehat{H}/dt=i\widehat{\mathcal{L}}\widehat{H}=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 $(\widehat{q},\widehat{p})=\widehat{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

$$\begin{array}{ccc}{\widehat{H}}_{\mathrm{W}}(X)& =\hfill & {\displaystyle \frac{P^2}{2M}}+{\displaystyle \frac{{\widehat{p}}^2}{2m}}+V_{\mathrm{W}}(\widehat{q},Q)\hfill \\ & =\hfill & {\displaystyle \frac{P^2}{2M}}+{\widehat{h}}_{\mathrm{W}}(Q),\hfill \end{array}$$

(5)

where m and M are the masses of the subsystem and bath DOF, respectively, and $V_{\mathrm{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 ${\widehat{h}}_{\mathrm{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 $\widehat{W}(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

$$\langle \widehat{\chi }\rangle (t)={\mathrm{Tr}}^{\prime }\int dX\widehat{W}(X,t){\widehat{\chi }}_{\mathrm{W}}(X),$$

(6)

where ${\mathrm{Tr}}^{\prime }$ denotes the partial trace involving a complete set of states of the quantum subsystem.

The operator-valued quasi-probability function in phase space evolves according to

$$\frac{\partial }{\partial t}\widehat{W}(X,t)=-\frac{i}{\hslash }\left[\begin{array}{cc}{\widehat{H}}_{\mathrm{W}}& \widehat{W}(X,t)\end{array}\right]\mathcal{D}\left[\begin{array}{c}{\widehat{H}}_{\mathrm{W}}\\ \widehat{W}(X,t)\end{array}\right]=-\frac{i}{\hslash }{\left[{\widehat{H}}_{\mathrm{W}},\widehat{W}(X,T)\right]}_{\mathcal{D}}=-i{\widehat{\mathcal{L}}}^{\mathcal{D}}\widehat{\chi },$$

(7)

where $\mathcal{D}$ is an antisymmetric matrix super-operator defined by

$$\mathcal{D}=\left[\begin{array}{cc}0& 1-\frac{\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }}{2i{\hslash }^{-1}}\\ -\left(1-\frac{\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }}{2i{\hslash }^{-1}}\right)& 0\end{array}\right],$$

(8)

with $\nabla =(\partial /\partial Q,\partial /\partial P)=\partial /\partial X$, and

$$\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }=\sum _{I,J=1}^{2N}{\overleftarrow{\nabla }}_I{\mathcal{B}}_{IJ}{\overrightarrow{\nabla }}_J$$

(9)

denotes the Poisson bracket operator. The last equality on the right-hand side of Equation (7) defines the quantum-classical Liouville operator $i{\widehat{\mathcal{L}}}^{\mathcal{D}}$. Equation (7) is the QCLE [48,49,50,51,52,53,54,55,56,57,58,59,60] of the system.

The QCLE in Equation (7) is founded upon a quasi-Lie bracket, which we may write explicitly as

$${[{\widehat{\chi }}_1(X),{\widehat{\chi }}_2(X)]}_{\mathcal{D}}=\left[\begin{array}{cc}{\widehat{\chi }}_1(X)& {\widehat{\chi }}_2(X)\end{array}\right]\mathcal{D}\left[\begin{array}{c}{\widehat{\chi }}_1(X)\\ {\widehat{\chi }}_2(X)\end{array}\right],$$

(10)

where $\mathcal{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):

$${\mathcal{J}}_{\mathcal{D}}={\left[{\widehat{\chi }}_1(X),{\left[{\widehat{\chi }}_2(X),{\widehat{\chi }}_3(X)\right]}_{\mathcal{D}}\right]}_{\mathcal{D}}+{\left[{\widehat{\chi }}_3(X),{\left[{\widehat{\chi }}_1(X),{\widehat{\chi }}_2(X)\right]}_{\mathcal{D}}\right]}_{\mathcal{D}}+{\left[{\widehat{\chi }}_2(X),{\left[{\widehat{\chi }}_3(X),{\widehat{\chi }}_1(X)\right]}_{\mathcal{D}}\right]}_{\mathcal{D}}\ne 0.$$

(11)

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

$$\begin{array}{c}\hfill e^{it{\mathcal{L}}^{\mathcal{D}}}\left[{\widehat{\chi }}_1(X,0),{\widehat{\chi }}_2(X)\right]\ne \left[e^{it{\mathcal{L}}^{\mathcal{D}}}{\widehat{\chi }}_1(X),e^{it{\mathcal{L}}^{\mathcal{D}}}{\widehat{\chi }}_2(X)\right].\end{array}$$

(12)

On the other hand, the quasi-Lie bracket conserves the energy $e^{it{\mathcal{L}}^{\mathcal{D}}}\widehat{H}{}_{\mathrm{W}}(X)=\widehat{H}{}_{\mathrm{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.

2.1. Derivation of the QCLE through a Partial Wigner Transform

When the bath DOF are described by canonically conjugate variables (and only in this case), the hybrid quantum can be derived by performing a partial Wigner transform of the quantum Liouville equation (QLE) over the bath DOF and taking a semiclassical limit of the resulting equations. To this end, let us consider the fully quantum counterpart to the Hamiltonian in Equation (5):

$$\widehat{H}=\frac{{\widehat{P}}^2}{2M}+\frac{{\widehat{p}}^2}{2m}+V(\widehat{q},\widehat{Q}).$$

(13)

The quantum statistical state of the system is described by the density matrix (or statistical operator) $\widehat{\rho }(t)$. The time dependence of the density matrix is dictated by the QLE:

$$\begin{array}{ccc}\hfill \frac{d}{dt}\widehat{\rho }(t)& =& -\frac{i}{\hslash }\left[\widehat{H},\widehat{\rho }(t)\right]=-\frac{i}{\hslash }\left[\begin{array}{cc}\widehat{\rho }& \widehat{H}\end{array}\right]\mathcal{B}\left[\begin{array}{c}\widehat{\rho }\\ \widehat{H}\end{array}\right],\hfill \end{array}$$

(14)

where $[...,...]$ denotes the commutator, and $\mathcal{B}$ is the symplectic matrix [102,126]:

$$\mathcal{B}=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right].$$

(15)

The average of an operator $\widehat{\chi }$ defined on the same Hilbert space of the system is calculated by

$$\langle \widehat{\chi }\rangle (t)=\mathrm{Tr}\left(\widehat{\rho }(t)\widehat{\chi }\right),$$

(16)

where $\mathrm{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 $\widehat{\rho }$ over the $\widehat{X}$’s:

$$\begin{array}{ccc}\hfill \widehat{W}(X,t)& =& \frac{1}{2\pi \hslash }\int dZ{e}^{iP·Z/\hslash }\langle Q-\frac{Z}{2}|\widehat{\rho }(t)|Q+\frac{Z}{2}\rangle .\hfill \end{array}$$

(17)

The symbol $\widehat{W}$ 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 the $\widehat{q}$’s and a function of the bath coordinates X. The partial Wigner transform of an arbitrary operator $\widehat{\chi }$ is analogously given by

$$\begin{array}{ccc}\hfill {\widehat{\chi }}_{\mathrm{W}}(X)& =& \int dZ{e}^{iP·Z/\hslash }\langle Q-\frac{Z}{2}|\widehat{\chi }|Q+\frac{Z}{2}\rangle .\hfill \end{array}$$

(18)

Taking the partial Wigner transform of Equation (16) leads to the expression for the average of $\widehat{\chi }$ 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 in ℏ, one arrives at the QCLE

$$\begin{array}{ccc}{\displaystyle \frac{\partial }{\partial t}}\widehat{W}(X,t)& =\hfill & -{\displaystyle \frac{i}{\hslash }}\left[{\widehat{H}}_{\mathrm{W}},\widehat{W}(X,t)\right]+{\displaystyle \frac{1}{2}}{\widehat{H}}_{\mathrm{W}}\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }\widehat{W}(X,t)-{\displaystyle \frac{1}{2}}\widehat{W}(X,t)\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }{\widehat{H}}_{\mathrm{W}}\hfill \\ & =\hfill & -i\mathcal{L}\widehat{W}(X,t),\hfill \end{array}$$

(19)

where the last equality defines the quantum Liouville operator $i\mathcal{L}=(i/\hslash )[{\widehat{H}}_{\mathrm{W}},·]-(1/2)({\widehat{H}}_{\mathrm{W}}\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }·)+(1/2)(·\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }{\widehat{H}}_{\mathrm{W}})$. To arrive at Equation (19), we have used the partial Wigner transform of a product of operators,

$$\begin{array}{ccc}\hfill {\left({\widehat{\chi }}_1{\widehat{\chi }}_2\right)}_{\mathrm{W}}(X)& =& {\widehat{\chi }}_{1,\mathrm{W}}(X)e^{\frac{i\hslash }{2}\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }}{\widehat{\chi }}_{2,\mathrm{W}}(X),\hfill \end{array}$$

(20)

and truncated the exponential after first order in ℏ, i.e.,

$$e^{\frac{i\hslash }{2}\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }}\approx 1+\frac{i\hslash }{2}\overleftarrow{\nabla }\mathcal{B}\overrightarrow{\nabla }.$$

(21)

It should be noted that Equation (21) is exact for Hamiltonians with quadratic bath terms and bilinear coupling between the $\widehat{x}$ and X DOF. In Ref. [47], it is shown how the linear expansion can be performed in terms of the parameter $\mu =\sqrt{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).

2.2. Integration Algorithm

A number of algorithms, which depend on the basis representation, exist for approximately solving the QCLE [50,51,54,55,56,60,78,113,114,115,116,117,118,119,120,121,122,123]. Herein, we illustrate the so-called Sequential Short-Time Propagation (SSTP) algorithm [78,113], 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

$${\widehat{h}}_{\mathrm{W}}|\alpha ;Q\rangle =E_{\alpha }(Q)|\alpha ;Q\rangle .$$

(22)

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 $\tau =t/n$. If one is able to calculate the propagation over a single $\tau$, the dynamics over the whole interval can be reconstructed by sequential iteration of the procedure. Let us then consider the quantum-classical propagator over a small step $\tau$ for the matrix elements of the operator-valued quasi-probability function $\widehat{W}(X)$ in the adiabatic basis. Such a propagator is written as

$${\left(e^{-i\tau \mathcal{L}}\right)}_{\alpha \alpha ,\beta {\beta }^{\prime }}\approx {\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}{e}^{-i{\int }_0^{\tau }ds{\omega }_{\alpha {\alpha }^{\prime }}(s)}{e}^{-i\tau L_{\alpha {\alpha }^{\prime }}}\left(1+\tau {\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\right).$$

(23)

On the right-hand side of Equation (23), we have introduced ${\omega }_{\alpha {\alpha }^{\prime }}$, the Bohr frequency defined in Equation (A3), $i{L}_{\alpha {\alpha }^{\prime }}$ is a classical-like Liouville operator, defined in Equation (A5), and ${\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}$ is the transition operator defined in Equation (A7). The SSTP dynamics of the matrix elements of $\widehat{W}(X,t)$ is given by

$$W_{\alpha \alpha }(X,\tau )=\sum _{\beta {\beta }^{\prime }}{\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}{e}^{-i{\int }_0^{\tau }ds{\omega }_{\alpha {\alpha }^{\prime }}(s)}{e}^{-i\tau L_{\alpha {\alpha }^{\prime }}}\left(1+\tau {\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\right)W_{\beta \beta }(X).$$

(24)

When $\tau$ 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 [113].

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 [43,62,64,65,66,67,68,69,70,71,72,73,74,75,78,82,85,87,88,89,90,91,113,118,119,120,121,122] the following expressions for the probabilities of making a transition (jump) and not-making a transition, respectively:

$$\begin{array}{ccc}\hfill {\mathcal{P}}_{\mathrm{J}}& =& \frac{|\tau \frac{P}{M}·d_{\alpha \beta }|}{1+|\tau \frac{P}{M}·d_{\alpha \beta }|},\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill Q_{\mathrm{NO}-\mathrm{J}}& =& \frac{1}{1+|\tau \frac{P}{M}·d_{\alpha \beta }|}.\hfill \end{array}$$

(26)

Another important technical ingredient of the algorithm is the approximation of the transition operator in Equation (A7) with its momentum-jump form:

$$\begin{array}{ccc}\hfill {\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathrm{MJ}}& =& {\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}\frac{P}{M}·d_{\alpha \beta }{e}^{(E_{\alpha }-E_{\beta })M\partial /\partial {(P·{\widehat{d}}_{\alpha \beta })}^2}+{\delta }_{\alpha \beta }\frac{P}{M}·d_{{\alpha }^{\prime }{\beta }^{\prime }}^{\ast }{e}^{(E_{\alpha }^{\prime }-E_{\beta }^{\prime })M\partial /\partial {(P·{\widehat{d}}_{{\alpha }^{\prime }{\beta }^{\prime }}^{\ast })}^2},\hfill \end{array}$$

(27)

where ${\widehat{d}}_{\alpha \beta }$ is the normalized coupling vector. Within the momentum-jump approximation [77,78], the action of the transition operator on the bath momenta can be easily obtained in closed form:

$$\begin{array}{ccc}\hfill e^{(E_{\alpha }-E_{\beta })M\partial /\partial {(P·{\widehat{d}}_{\alpha \beta })}^2}P& =& P-P\left(P·{\widehat{d}}_{\alpha \beta }\right)+{\widehat{d}}_{\alpha \beta }\sqrt{{\left(P·{\widehat{d}}_{\alpha \beta }\right)}^2+M\left(E_{\alpha }-E_{\beta }\right)}.\hfill \end{array}$$

(28)

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 [127] 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 [78,113] 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 [123]. The second problem is not just associated with the SSTP algorithm, but it is common to all Monte Carlo approaches to the calculation of quantum averages: the infamous sign-problem. The sign-problem is one of the major unsolved problems in the physics of quantum systems. Within the SSTP algorithm, it manifests itself both through the oscillating phase factors associated with the propagation on mean-energy surfaces and through the accumulation of fluctuating weights associated with the Monte Carlo sampling of the quantum transitions. In practice, upon analyzing the results obtained by means of this algorithm [43,62,64,65,66,67,68,69,70,71,72,73,74,75,78,82,85,87,88,89,90,91,113,115,116,117,118,119,120,121,122], we can conclude that the more quantum is the character of the bath the greater is the error in the calculation of the averages.

The mapping of the calculation of averages via the SSTP algorithm onto a stochastic process is reminiscent of the approach to open quantum system dynamics provided by the Stochastic Liouville Equation (SLE) [128,129,130,131]. 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 [123] does not involve any stochastic process whatsoever.

3. Classical Spin Baths

Contrary to what some books in quantum mechanics state (in the authors’s knowledge, an exception is Schulman’s book [132]), the concept of spin can be defined in an entirely classical way [132,133,134,135,136]. In practice, spinors provide a more fundamental representation of the rotation group than that given by tensors [132,133,134,135,136]. 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 [137], whose Hamiltonian is

$$H_{CHS}=\sum _{a=x,y,z}\sum _{I,J}^N{S}_a^I{\mathcal{C}}_{IJ}^a{S}_a^J,$$

(29)

where ${\mathbf{S}}^I$ are N classical vectors obeying the constraint

$${\left(S_x^I\right)}^2+{\left(S_y^I\right)}^2+{\left(S_z^I\right)}^2=1,$$

(30)

for $I=1,...,N$, and the ${\mathcal{C}}_{IJ}^a$ 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 $\mathbf{S}$, with components $S_a$, $a=x,y,z$, and Hamiltonian $H^{\mathrm{S}}(\mathbf{S})$. Let us define the spin gradient as ${\nabla }^{\mathbf{S}}=\partial /\partial \mathbf{S}$, which in terms of the spin components is written as ${\nabla }_a^S=\partial /\partial S_a$, with $a=x,y,z$. The equations of motion of the spin are then written as

$$\dot{\mathbf{S}}={\mathcal{B}}^{\mathbf{S}}{\nabla }^{\mathbf{S}}{H}_{\mathrm{S}},$$

(31)

where

$${\mathcal{B}}^{\mathbf{S}}=\left[\begin{array}{ccc}0& S_z& -S_y\\ -S_z& 0& S_x\\ S_y& -S_x& 0\end{array}\right].$$

(32)

One can also adopt the compact form ${\mathcal{B}}_{ab}^{\mathbf{S}}={\sum }_{c=x,y,z}{ϵ}_{abc}{S}_c$ and $a,b=x,y,z$ of the antisymmetric matrix ${\mathcal{B}}^{\mathbf{S}}$, where ${ϵ}_{abc}$ is the Levi–Civita pseudo-tensor. The Casimir $C_2=\mathbf{S}·\mathbf{S}$ is preserved by the equations of motion (31), independently of the form of the spin Hamiltonian $H^{\mathrm{S}}(\mathbf{S})$. In addition, the dynamics has a zero phase space compressibility ${\kappa }^{\mathrm{S}}=\nabla \mathbf{S}·\dot{\mathbf{S}}=0$. The classical phase space flow of the spin is defined through the non-canonical bracket

$$\sum _{a,b}A(\mathbf{S}){\overleftarrow{\nabla }}_a^{\mathbf{S}}{\mathcal{B}}_{ab}^{\mathrm{S}}{\overrightarrow{\nabla }}_b^{\mathbf{S}}B(\mathbf{S})=A(\mathbf{S}){\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathbf{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}B(\mathbf{S}),$$

(33)

where $A=A(\mathbf{S})$ and $B=B(\mathbf{S})$ are arbitrary functions of the spin DOF.

Consider now the hybrid quantum-classical Hamiltonian of a quantum subsystem coupled to the classical spin

$$\begin{array}{ccc}\widehat{\mathcal{H}}(\mathbf{S})& =\hfill & \widehat{H}(\left\{\widehat{\chi }\right\})+V_{\mathrm{C}}(\left\{\widehat{\chi }\right\},\mathbf{S})+H^{\mathbf{S}}(\mathbf{S})\hfill \\ & =\hfill & {\widehat{h}}_{\mathbf{S}}(\mathbf{S})+H^{\mathbf{S}}(\mathbf{S}),\hfill \end{array}$$

(34)

describing a quantum subsystem in terms of the Hamiltonian $\widehat{H}(\left\{\widehat{\chi }\right\})$, depending on the operators $\left\{\widehat{\chi }\right\}$, $V(\left\{\widehat{\chi }\right\},\mathbf{S})$ is the subsystem-spin interaction potential, and the second line of the equation defines the adiabatic Hamiltonian ${\widehat{h}}_{\mathbf{S}}$. The quantum-classical dynamics of the operator-valued quasi-probability function (defined in the spinor space of the total system), ${\widehat{W}}^{\mathbf{S}}(\mathbf{S},t)$, is dictated by the spin-bath QCLE [90,91]

$$\begin{array}{ccc}{\displaystyle \frac{\partial }{\partial t}}{\widehat{W}}^{\mathbf{S}}(\mathbf{S},t)& =\hfill & -{\displaystyle \frac{i}{\hslash }}\left[\begin{array}{cc}\widehat{\mathcal{H}}(\mathbf{S})& {\widehat{W}}^{\mathbf{S}}(\mathbf{S},t)\end{array}\right]{\mathcal{D}}^{\mathrm{S}}\left[\begin{array}{c}\widehat{\mathcal{H}}(\mathbf{S})\\ {\widehat{W}}^{\mathbf{S}}(\mathbf{S},t)\end{array}\right]\hfill \\ & =\hfill & -{\displaystyle \frac{i}{\hslash }}{\left[\widehat{\mathcal{H}}(\mathbf{S}),{\widehat{W}}^{\mathbf{S}}(\mathbf{S},t)\right]}_{{\mathcal{D}}^{\mathrm{S}}},\hfill \end{array}$$

(35)

where

$$\begin{array}{ccc}\hfill {\mathcal{D}}^{\mathrm{S}}& =& \left[\begin{array}{cc}0& 1+\frac{i\hslash }{2}\overleftarrow{\nabla }{\mathcal{B}}^{\mathrm{S}}\overrightarrow{\nabla }\\ -1-\frac{i\hslash }{2}\overleftarrow{\nabla }{\mathcal{B}}^{\mathrm{S}}\overrightarrow{\nabla }& 0\end{array}\right].\hfill \end{array}$$

(36)

We next set out to represent Equation (35) in the adiabatic basis $|\alpha ;\mathbf{S}\rangle$ defined by the eigenvalue equation

$${\widehat{h}}_{\mathbf{S}}(\mathbf{S})|\alpha ;\mathbf{S}\rangle =E_{\alpha }(\mathbf{S})|\alpha ;\mathbf{S}\rangle .$$

(37)

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 $\mathbf{S}$. In this basis, Equation (35) becomes

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial t}}{\widehat{W}}_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}=& -i{\omega }_{\alpha {\alpha }^{\prime }}{W}_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}-H^{\mathrm{S}}{\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathbf{S}}\langle \alpha |{\overrightarrow{\nabla }}^{\mathbf{S}}{\widehat{W}}^{\mathbf{S}}|{\alpha }^{\prime }\rangle \hfill \\ & +{\displaystyle \frac{1}{2}}\langle \alpha |{\widehat{h}}_{\mathbf{S}}{\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathbf{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}{\widehat{W}}^{\mathbf{S}}|{\alpha }^{\prime }\rangle -{\displaystyle \frac{1}{2}}\langle \alpha |{\widehat{W}}^{\mathbf{S}}{\overleftarrow{\nabla }}^{\mathrm{S}}{\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}{\widehat{h}}_{\mathbf{S}}|{\alpha }^{\prime }\rangle ,\hfill \end{array}$$

(38)

where ${\omega }_{\alpha {\alpha }^{\prime }}=E_{\alpha }(\mathbf{S})-E_{{\alpha }^{\prime }}(\mathbf{S})/\hslash$ is the Bohr frequency. Defining the spin coupling vector

$$\begin{array}{ccc}\hfill d_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}& =& \langle \alpha ;\mathbf{S}|{\overrightarrow{\nabla }}^{\mathbf{S}}|{\alpha }^{\prime };\mathbf{S}\rangle ,\hfill \end{array}$$

(39)

one finds the two identities

$$\begin{array}{ccc}\hfill \langle \alpha ;\mathbf{S}|\left({\overrightarrow{\nabla }}^{\mathbf{S}}{\widehat{W}}^{\mathbf{S}}(\mathbf{S})\right)|{\alpha }^{\prime };\mathbf{S}\rangle & =& {\overrightarrow{\nabla }}^{\mathbf{S}}{W}_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S})+\sum _{\beta }{d}_{\alpha \beta }^{\mathbf{S}}{W}_{\beta {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S})-\sum _{{\beta }^{\prime }}{W}_{\alpha {\beta }^{\prime }}^{\mathbf{S}}(\mathbf{S})d_{{\beta }^{\prime }{\alpha }^{\prime }}^{\mathbf{S}}\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill \langle \alpha ;\mathbf{S}|\left({\overrightarrow{\nabla }}^{\mathbf{S}}{\widehat{h}}_{\mathbf{S}}(\mathbf{S})\right)|{\alpha }^{\prime };\mathbf{S}\rangle & =& {\overrightarrow{\nabla }}^{\mathbf{S}}{h}_{\mathbf{S}}^{\alpha {\alpha }^{\prime }}-\Delta E_{\alpha {\alpha }^{\prime }}{d}_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}\hfill \end{array}$$

(41)

where $\Delta E_{\alpha {\alpha }^{\prime }}=E_{\alpha }-E_{{\alpha }^{\prime }}$. Using Equations (40) and (41), the spin-bath QCLE may be rewritten as

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}{W}_{\alpha {\alpha }^{\prime }}^{\mathrm{S}}(\mathbf{S},t)& =& -\sum _{\beta {\beta }^{\prime }}\left(i{\omega }_{\alpha {\alpha }^{\prime }}{\delta }_{\alpha \beta }{\delta }_{\alpha {\alpha }^{\prime }}+i{L}_{\alpha {\alpha }^{\prime }}{\delta }_{\alpha \beta }{\delta }_{\alpha {\alpha }^{\prime }}+{\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathrm{S}}+{\mathcal{S}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\right)W_{\beta {\beta }^{\prime }}^{\mathrm{S}}(\mathrm{S},t),\hfill \end{array}$$

(42)

where we have defined the classical-like spin-Liouville operator

$$\begin{array}{ccc}\hfill i{L}_{\alpha {\alpha }^{\prime }}& =\hfill & H_{\mathrm{S}}{\overleftarrow{\nabla }}^{\mathrm{S}}{\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}+{\displaystyle \frac{1}{2}}{E}_{{\alpha }^{\prime }}{\overleftarrow{\nabla }}^{\mathrm{S}}{\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}+{\displaystyle \frac{1}{2}}{E}_{\alpha }{\overleftarrow{\nabla }}^{\mathrm{S}}{\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}\hfill \\ & =\hfill & \left({\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}{H}_{\alpha {\alpha }^{\prime }}^{\mathrm{S}}\right)·{\overrightarrow{\nabla }}^{\mathrm{S}},\hfill \end{array}$$

(43)

with the average adiabatic Hamiltonian

$$H_{\alpha {\alpha }^{\prime }}^{\mathrm{S}}=H_{\mathrm{S}}+\frac{1}{2}\left(E_{\alpha }+E_{{\alpha }^{\prime }}\right).$$

(44)

The transition operator for the spin bath is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathbf{S}}=& d_{\alpha \beta }^{\mathbf{S}}·\left({\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}{H}_{\mathrm{S}}\right){\delta }_{{\beta }^{\prime }{\alpha }^{\prime }}+{\displaystyle \frac{1}{2}}\Delta E_{\alpha \beta }{d}_{\alpha \beta }^{\mathbf{S}}·\left({\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}\right){\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}\hfill \\ & +d_{{\alpha }^{\prime }{\beta }^{\prime }}^{\mathbf{S}\ast }·\left({\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}{H}_{\mathrm{S}}\right){\delta }_{\alpha \beta }+{\displaystyle \frac{1}{2}}\Delta E_{{\alpha }^{\prime }{\beta }^{\prime }}{d}_{{\alpha }^{\prime }{\beta }^{\prime }}^{\mathbf{S}\ast }·\left({\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathrm{S}}\right){\delta }_{\alpha \beta }.\hfill \end{array}\end{array}$$

(45)

The limit $d_{\alpha \alpha }^{\mathbf{S}}\to 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):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{S}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}=& {\displaystyle \frac{1}{2}}\Delta E_{\alpha \sigma }{d}_{\alpha \sigma }^{\mathbf{S}}{\mathcal{B}}^{\mathrm{S}}{d}_{\sigma \beta }^{\mathbf{S}}{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+{\displaystyle \frac{1}{2}}\Delta E_{\alpha \beta }{d}_{\alpha \beta }^{\mathbf{S}}{\mathcal{B}}^{\mathrm{S}}{d}_{{\alpha }^{\prime }{\beta }^{\prime }}^{\mathbf{S}\ast }\hfill \\ & +{\displaystyle \frac{1}{2}}\Delta E_{{\alpha }^{\prime }{\sigma }^{\prime }}{d}_{{\alpha }^{\prime }{\sigma }^{\prime }}^{\mathbf{S}\ast }{\mathcal{B}}^{\mathrm{S}}{d}_{{\sigma }^{\prime }{\beta }^{\prime }}^{\mathbf{S}\ast }{\delta }_{\alpha \beta }+{\displaystyle \frac{1}{2}}\Delta E_{{\alpha }^{\prime }{\beta }^{\prime }}{d}_{{\alpha }^{\prime }{\beta }^{\prime }}^{\mathbf{S}\ast }{\mathcal{B}}^{\mathrm{S}}{d}_{\alpha \beta }^{\mathbf{S}}\hfill \\ & -{\displaystyle \frac{1}{2}}\left(E_{\alpha }+E_{{\alpha }^{\prime }}\right){\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathrm{S}}·d_{\alpha \beta }^{\mathbf{S}}{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}-{\displaystyle \frac{1}{2}}\left(E_{\alpha }+E_{{\alpha }^{\prime }}\right){\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathrm{S}}·d_{{\alpha }^{\prime }{\beta }^{\prime }}^{\mathrm{S}\ast }{\delta }_{\alpha \beta }.\hfill \end{array}\end{array}$$

(46)

The adiabatic limit of the spin-bath QCLE in (42) can be taken by setting to zero the off-diagonal elements of $d_{\alpha {\alpha }^{\prime }}$, 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

$$\begin{array}{c}\hfill \begin{array}{ccc}{\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathbf{S},\mathrm{ad}}\hfill & =\hfill & \left(d_{\alpha \alpha }^{\mathbf{S}}+d_{{\alpha }^{\prime }{\alpha }^{\prime }}^{\mathbf{S}\ast }\right){\mathcal{B}}^{\mathbf{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}{H}_{\mathrm{S}}{\delta }_{\alpha \beta }{\delta }_{{\beta }^{\prime }{\alpha }^{\prime }}\hfill \\ & =\hfill & -i\left({\varphi }_{\alpha \alpha }^{\mathbf{S}}-{\varphi }_{{\alpha }^{\prime }{\alpha }^{\prime }}^{\mathbf{S}}\right){\mathcal{B}}^{\mathbf{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}{\delta }_{\alpha \beta }{\delta }_{{\beta }^{\prime }{\alpha }^{\prime }}.\hfill \end{array}\end{array}$$

(47)

The geometric phase

$${\varphi }_{\alpha \alpha }^{\mathrm{S}}=-i{d}_{\alpha \alpha }^{\mathbf{S}}$$

(48)

has been introduced exploiting the purely imaginary character of $d_{\alpha \alpha }^{\mathbf{S}}$. Similarly, the higher order transition operator becomes

$$\begin{array}{ccc}\hfill {\mathcal{S}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathrm{ad}}& =& -\frac{i}{2}\sum _{I,J}\left({\varphi }_{\alpha \alpha }^{\mathbf{S}}-{\varphi }_{{\alpha }^{\prime }{\alpha }^{\prime }}^{\mathbf{S}}\right){\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}\left(E_{\alpha }+E_{{\alpha }^{\prime }}\right){\delta }_{\alpha \alpha }{\delta }_{{\alpha }^{\prime }{\alpha }^{\prime }}\hfill \end{array}$$

(49)

Putting everything together, the adiabatic approximation of the spin-bath QCLE may be written as

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}{W}_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S},t)& =& \left[-i{\omega }_{\alpha {\alpha }^{\prime }}-i\left({\varphi }_{\alpha \alpha }^{\mathbf{S}}-{\varphi }_{{\alpha }^{\prime }{\alpha }^{\prime }}^{\mathbf{S}}\right)\mathcal{B}{\overrightarrow{\nabla }}^{\mathbf{S}}{H}_{\mathbf{S}}^{\alpha {\alpha }^{\prime }}-H_{\mathrm{S}}^{\alpha {\alpha }^{\prime }}{\overleftarrow{\nabla }}^{\mathbf{S}}\mathcal{B}{\overrightarrow{\nabla }}^{\mathbf{S}}\right]W_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S},t).\hfill \end{array}$$

(50)

In Equation (50), the phase ${\omega }_{\alpha {\alpha }^{\prime }}$ has a dynamical nature while the phase ${\varphi }_{\alpha \alpha }^{\mathbf{S}}$ is of a geometric origin and it can be considered an instance of the famous Berry phase [107,108,109]. Interestingly, Equation (35) predicts that the geometric phase ${\varphi }_{\alpha \alpha }^{\mathbf{S}}$ can be non-zero also for open paths of the classical spins of the bath (open-path Berry phases were discussed in Ref. [104]). Moreover, the phase factor ${\varphi }_{\alpha \alpha }^{\mathbf{S}}-{\varphi }_{{\alpha }^{\prime }{\alpha }^{\prime }}^{\mathbf{S}}$ is purely off-diagonal (off-diagonal Berry phases for environments described by canonically conjugate variables were discussed in Refs. [103,105,106]). It is worth mentioning that the geometric phase ${\varphi }_{\alpha \alpha }^{\mathbf{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

$$\frac{\partial }{\partial t}{W}_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S},t)=\left[-i{\omega }_{\alpha {\alpha }^{\prime }}-\left(\langle \alpha ,S|\frac{d}{dt}|\alpha ,S\rangle -\langle {\alpha }^{\prime },S|\frac{d}{dt}|{\alpha }^{\prime },S\rangle \right)-H_{\mathrm{S}}^{\alpha {\alpha }^{\prime }}{\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathbf{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}\right]W_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S},t).$$

(51)

Using the Dyson identity, one can obtain the following form for ${\widehat{W}}^{\mathbf{S}}(\mathbf{S},t)$ in terms of the adiabatic propagator:

$$\begin{array}{cc}\hfill W_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S},t)=& exp\left[-i{\int }_{t_0}^{t}d{t}^{\prime }{\omega }_{\alpha {\alpha }^{\prime }}(t^{\prime })\right]exp\left[-{\int }_{t_0}^{t}d{t}^{\prime }\left(\langle \alpha ,S|\frac{d}{d{t}^{\prime }}|\alpha ,S\rangle -\langle {\alpha }^{\prime },S|\frac{d}{d{t}^{\prime }}|{\alpha }^{\prime },S\rangle \right)\right]\hfill \\ & \times exp\left[-(t-t_0)H_{\alpha {\alpha }^{\prime }}^{\mathrm{S}}{\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}\right]W_{\alpha {\alpha }^{\prime }}^{\mathbf{S}}(\mathbf{S},t_0).\hfill \end{array}$$

(52)

Equation (52) provides a convenient starting point for devising numerical integration schemes based on the SSTP propagation scheme [113].

In Ref. [91], the following model Hamiltonian was considered:

$$\begin{array}{ccc}\hfill \widehat{H}(\mathbf{S})& =& -\Omega {\widehat{\sigma }}_x-c_{1}b{\widehat{\sigma }}_z-\mu \mathbf{S}·\mathbf{\sigma }-c_{2}b{S}_z+\frac{S_z^2}{2}\hfill \end{array}$$

(53)

$$\begin{array}{ccc}\hfill \widehat{H}(\mathbf{S})& =& {\widehat{h}}_{\mathbf{S}}(\mathbf{S})-c_{2}b{S}_z+\frac{S_z^2}{2},\hfill \end{array}$$

(54)

where $\Omega$, $c_1$, and $c_2$ are real parameters, b is the z component of the magnetic field $\mathbf{B}=(0,0,b)$, while $\mathit{\sigma }=({\sigma }_x,{\sigma }_y,{\sigma }_z)$ is a vector having the Pauli matrices ${\sigma }_x$, ${\sigma }_y$, and ${\sigma }_z$ as components. The SSTP algorithm was applied to Equation (52) and the action of the classical like Liouville operator $H_{\alpha {\alpha }^{\prime }}^{\mathrm{S}}{\overleftarrow{\nabla }}^{\mathbf{S}}{\mathcal{B}}^{\mathrm{S}}{\overrightarrow{\nabla }}^{\mathbf{S}}$ was evaluated using time reversible integration algorithms based on the symmetric break-up of the Liouville propagator [138,139,140].

4. 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^{\prime }=(Q^{\prime },P^{\prime })$ 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^{\prime }$ is not interesting: their function is just that of working as a thermal bath, leading to dissipative dynamics [44].

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 [44]. It takes the form,

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}\widehat{W}(X,t)& =& -\frac{i}{\hslash }\left[\begin{array}{cc}{\widehat{H}}_{\mathrm{W}}& \widehat{W}(X,t)\end{array}\right]\mathcal{D}\left[\begin{array}{c}{\widehat{H}}_{\mathrm{W}}\\ \widehat{W}(X,t)\end{array}\right]\hfill \\ & +& \zeta {\overrightarrow{\nabla }}_P\left(\frac{P}{M}+k_{B}T{\overrightarrow{\nabla }}_P\right)\widehat{W}(X,t),=-i{\widehat{\mathcal{L}}}^{\mathrm{D}}\widehat{W}(X,t),\hfill \end{array}$$

(55)

where ${\nabla }_P=\partial /\partial P$, $\zeta$ 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_{\mathrm{W}}(\widehat{q},Q)$ as the potential of mean force arising from the average over the primed bath variables $Q^{\prime }$. The Liouville operator $i{\widehat{\mathcal{L}}}^{\mathrm{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 [47]. In particular, the term $\zeta {\overrightarrow{\nabla }}_P\left[(P/M)+k_{B}T{\overrightarrow{\nabla }}_P\right]$ 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

$$i{\widehat{\mathcal{L}}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}^{\mathrm{D}}=\left(i{\omega }_{\alpha {\alpha }^{\prime }}(R)+i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{K}}\right){\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+{\mathcal{T}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }},$$

(56)

where we have defined the Kramers operator as

$$\begin{array}{ccc}\hfill i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{K}}& =& \left[\frac{P}{M}{\overrightarrow{\nabla }}_Q+\frac{1}{2}\left(F_W^{\alpha }+F_W^{{\alpha }^{\prime }}\right){\overrightarrow{\nabla }}_P-\zeta {\overrightarrow{\nabla }}_P\left(\frac{P}{M}+k_{B}T{\overrightarrow{\nabla }}_P\right)\right].\hfill \end{array}$$

(57)

The quantum-classical average of any operator or dynamical variable $\widehat{\chi }(X)$ can be written as

$$\begin{array}{c}\hfill \begin{array}{ccc}\langle \widehat{\chi }\rangle (t)& =\hfill & \sum _{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}\int dX{\chi }_{{\alpha }^{\prime }\alpha }(X)exp[-i{\mathcal{L}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}^{D}t]W^{\beta {\beta }^{\prime }}(X)\hfill \\ & =\hfill & \sum _{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}\int dX{W}^{\beta {\beta }^{\prime }}(X)exp\left[i{\mathcal{L}}_{{\beta }^{\prime }\beta {\alpha }^{\prime }\alpha }^{\mathrm{DB}}t\right]{\chi }_{{\alpha }^{\prime }\alpha }(R,P),\hfill \end{array}\end{array}$$

(58)

where $i{\mathcal{L}}_{{\beta }^{\prime }\beta {\alpha }^{\prime }\alpha }^{\mathrm{DB}}$ is the backward operator, defined as

$$\begin{array}{ccc}\hfill i{\widehat{\mathcal{L}}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}^{DB}& =& \left(i{\omega }_{\alpha {\alpha }^{\prime }}(R)+i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{KB}}\right){\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+{\mathcal{T}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}\hfill \end{array}$$

(59)

The backward Kramers $i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{KB}}$ operator is written as

$$\begin{array}{ccc}\hfill i{L}_{\alpha \alpha }^{\mathrm{KB}}& =& \left[\frac{P}{M}{\overrightarrow{\nabla }}_Q+\frac{1}{2}\left(F_{\mathrm{W}}^{\alpha }+F_{\mathrm{W}}^{{\alpha }^{\prime }}\right){\overrightarrow{\nabla }}_P-\zeta \left(\frac{P}{M}-k_{B}T{\overrightarrow{\nabla }}_P\right){\overrightarrow{\nabla }}_P\right]{\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}.\hfill \end{array}$$

(60)

According to the classical theory of random processes [110], the time evolution under the backward Kramers operator $i{\mathcal{L}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}^{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,

$$\begin{array}{ccc}\hfill \dot{Q}& =& \frac{P}{M},\hfill \end{array}$$

(61)

$$\begin{array}{ccc}\hfill \dot{P}& =& -\frac{\zeta }{M}P+\frac{1}{2}\left(F_{\mathrm{W}}{W}^{\alpha }+F_{\mathrm{W}}^{{\alpha }^{\prime }}\right)+\mathcal{R}(t),\hfill \end{array}$$

(62)

where $\mathcal{R}(t)$ is a Gaussian white noise process with the properties,

$$\begin{array}{ccc}\hfill \langle \mathcal{R}(t)\rangle & =& 0,\hfill \end{array}$$

(63)

$$\begin{array}{ccc}\hfill \langle \mathcal{R}(t)\mathcal{R}(t^{\prime })\rangle & =& 2k_{B}T\zeta \delta (t-t^{\prime }).\hfill \end{array}$$

(64)

To Equations (61) and (62), one can associate a time-dependent Langevin–Liouville operator

$$i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{L}}(t)=\frac{P}{M}{\overrightarrow{\nabla }}_Q+\left(-\frac{\zeta }{M}P+\frac{1}{2}(F_{\mathrm{W}}^{\alpha }+F_{\mathrm{W}}^{\alpha })+\mathcal{R}(t)\right){\overrightarrow{\nabla }}_P,$$

(65)

and a time-ordered propagator

$$U_{\alpha {\alpha }^{\prime }}^{\mathrm{L}}(t,0)=\mathcal{T}exp\left[{\int }_0^{t}d{t}^{\prime }i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{L}}(t^{\prime })\right].$$

(66)

In order to generate the stochastic Langevin trajectories, we can use a total time-dependent Langevin–Liouville super-operator

$$i{\widehat{\mathcal{L}}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}^{\mathrm{L}}(t)=\left(i{\omega }_{\alpha {\alpha }^{\prime }}(Q)+i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{L}}(t)\right){\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+{\mathcal{T}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}$$

(67)

and the associated propagator

$${\mathcal{U}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}^{\mathrm{L}}(t,0)=\mathcal{T}exp\left[{\int }_0^{t}d{t}^{\prime }i{\mathcal{L}}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}^{\mathrm{L}}(t^{\prime })\right].$$

(68)

Within such a Langevin picture, the quantum-classical average of any operator $\widehat{\chi }(X)$ can be calculated as

$$\begin{array}{ccc}\hfill \langle \widehat{\chi }\rangle (t)& =& \sum _{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}\int dX{W}^{\beta {\beta }^{\prime }}(Q)\overline{{\mathcal{U}}_{\beta {\beta }^{\prime }\alpha {\alpha }^{\prime }}^{\mathrm{L}}(t){\chi }_{{\alpha }^{\prime }\alpha }(Q)}\hfill \end{array}$$

(69)

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

$$\begin{array}{ccc}\hfill \langle \widehat{\chi }(X,t)\rangle & =& \sum _{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}\overline{\int dRdP{W}^{\beta {\beta }^{\prime }}(X){\mathcal{U}}_{\beta {\beta }^{\prime }\alpha {\alpha }^{\prime }}^{\mathrm{L}}(t){\chi }_{{\alpha }^{\prime }\alpha }^{\prime }(XP)}.\hfill \\ \hfill \end{array}$$

(70)

Equation (70) allows one to calculate averages in a quantum-classical dissipative system as phase space weighted averages over many Langevin trajectories.

In Ref. [45], a quantum subsystem with two energy levels interacting with a dissipative classical quartic oscillator was considered. The Hamiltonian of the hybrid quantum-classical system reads

$$\begin{array}{ccc}\hfill {\widehat{H}}_W(X)& =& \frac{P^2}{2M}+V_q(Q)-\hslash \Omega {\widehat{\sigma }}_x-\hslash {\gamma }_{0}Q{\widehat{\sigma }}_z,\hfill \end{array}$$

(71)

where $V_q(Q)=\frac{a}{4}{R}^4-\frac{b}{2}{R}^2$, $\Omega$, a, b, and ${\gamma }_0$ are real parameters, M is the mass of the quartic oscillator, and ${\widehat{\sigma }}_x$ and ${\widehat{\sigma }}_z$ are Pauli matrices.

The calculation of quantum-classical averages using the dynamics defined by the time-dependent Langevin–Liouville propagator ${\mathcal{U}}_{s{s}^{\prime }}^L(t)$ in Equation (68) is no more complicated than that for deterministic quantum-classical dynamics. The momentum-jump approximation [77,78] and a simple generalization of the SSTP algorithm [78,113] to the time dependent propagator were used in Ref. [45]. The explicitly time-dependent propagator ${\mathcal{U}}_{s{s}^{\prime }}^L(t)$ must be defined as a time ordered product. A simple way to achieve that is to employ the decomposition scheme devised by Suzuki [141]. Details of the numerical procedures are found in Ref. [45]

5. Non-Hamiltonian Dynamics in Thermal Baths

By exploiting the antisymmetric structure of the quantum-classical commutator, arising from the matrix operator $\mathcal{D}$ given in Equation (8), one can impose the thermodynamic constraints of constant temperature on the classical-like DOF [42,43]. Following Refs. [92,93,94], constant-temperature dynamics for the classical bath coordinates, as defined through the non-Hamiltonian Nosé–Hoover equations of motion, can be introduced by modifying the matrix $\mathcal{B}$ and augmenting in a minimal way the dimension of the phase space bath. The classical Nosé–Hoover thermostat is briefly discussed in Appendix B.

As in the classical case, the Nosé variables are

$$X^{\mathrm{N}}\equiv (Q,Q_{\eta },P,P_{\eta }),$$

(72)

where $Q_{\eta }$ and $P_{\eta }$ are the Nosé coordinate and momentum. The Nosé quantum-classical Hamiltonian is obtained by adding the Nosé kinetic energy $P_{\eta }^2/2M_{\eta }$ and potential energy $N{k}_{B}T{Q}_{\eta }$ to ${\widehat{H}}_{\mathrm{W}}$ in Equation (5)

$$H^{\mathrm{N}}=\frac{P^2}{2M}+\frac{P_{\eta }^2}{2M_{\eta }}+N{k}_{\mathrm{B}}T{Q}_{\eta }+{\widehat{h}}_{\mathrm{W}}(Q),$$

(73)

where $M_{\eta }$ is the Nosé inertial parameter, $k_{\mathrm{B}}$ is the Boltzmann constant, T is the constant temperature, and N is the number of Q coordinates. Using the matrix ${\mathcal{B}}^{\mathrm{N}}$ in Equation (A13), the classical phase space quasi-Hamiltonian bracket of two variables $A_1$ and $A_2$ can be defined as

$$A_1{\overleftarrow{\nabla }}^{\mathrm{N}}{\mathcal{B}}^{\mathrm{N}}{\overrightarrow{\nabla }}^{\mathrm{N}}{A}_2=\sum _{I,J=1}^{2(N+1)}{A}_1{\overleftarrow{\nabla }}_I^{\mathrm{N}}{\mathcal{B}}_{IJ}^{\mathrm{N}}{\overrightarrow{\nabla }}_J^{\mathrm{N}}{A}_2.$$

(74)

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

$${\mathcal{D}}^{\mathrm{N}}=\left[\begin{array}{cc}0& 1-\frac{{\overleftarrow{\nabla }}^{\mathrm{N}}{\mathcal{B}}^{\mathrm{N}}{\overrightarrow{\nabla }}^{\mathrm{N}}}{2i{\hslash }^{-1}}\\ -\left(1-\frac{{\overleftarrow{\nabla }}^{\mathrm{N}}{\mathcal{B}}^{\mathrm{N}}{\overrightarrow{\nabla }}^{\mathrm{N}}}{2i{\hslash }^{-1}}\right)& 0\end{array}\right].$$

(75)

The Nosé–Hoover QCLE for the operator-valued quasi-probability function ${\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)$ is given by

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \frac{d}{dt}{\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)& =\hfill & -i{\mathcal{L}}^{\mathrm{N}}{W}^{\mathrm{N}}(X^{\mathrm{N}},t)-{\kappa }^{\mathrm{N}}(X^{\mathrm{N}})W^{\mathrm{N}}(X^{\mathrm{N}},t)\hfill \\ & =\hfill & -\frac{i}{\hslash }\left[\begin{array}{cc}{\widehat{H}}^{\mathrm{N}}& {\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)\end{array}\right]·{\mathcal{D}}^{\mathrm{N}}·\left[\begin{array}{c}{\widehat{H}}^{\mathrm{N}}\\ {\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)\end{array}\right]-{\kappa }^{\mathrm{N}}(X^{\mathrm{N}})W^{\mathrm{N}}(X^{\mathrm{N}},t).\hfill \end{array}\end{array}$$

(76)

The presence of the term $-{\kappa }^{\mathrm{N}}(X^{\mathrm{N}})W^{\mathrm{N}}(X^{\mathrm{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 obtains

$$\begin{array}{ccc}\hfill {\widehat{H}}^{\mathrm{N}}{\overleftarrow{\nabla }}^{\mathrm{N}}{\mathcal{B}}^{\mathrm{N}}{\overrightarrow{\nabla }}^{\mathrm{N}}\widehat{\chi }(X^{\mathrm{N}},t)-\widehat{\chi }(X^{\mathrm{N}},t){\overleftarrow{\nabla }}^{\mathrm{N}}{\mathcal{B}}^{\mathrm{N}}{\overrightarrow{\nabla }}^{\mathrm{N}}{\widehat{H}}_{\mathrm{N}}& =& \frac{\partial \widehat{V}}{\partial Q}\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial P}+\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial P}\frac{\partial \widehat{V}}{\partial Q}\hfill \\ & -& 2F_{Q_{\eta }}\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial P_{\eta }}-2\frac{P}{M}\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial Q}\hfill \\ & -& 2\frac{P_{\eta }}{M_{\eta }}\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial Q_{\eta }}+2\frac{P_{\eta }}{M_{\eta }}P\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial P},\hfill \end{array}$$

(77)

where $F_{Q_{\eta }}=\frac{P^2}{M}-N{k}_{B}T$. Finally, using the above result, the Nosé–Hoover QCLE reads

$$\begin{array}{c}\hfill \begin{array}{ccc}\frac{d}{dt}{\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)\hfill & =\hfill & -\frac{i}{\hslash }\left(H^{\mathrm{N}}{\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)-\widehat{\chi }(X^{\mathrm{N}},t)H^{\mathrm{N}}\right)+\frac{1}{2}\left(\frac{\partial {\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)}{\partial P}\frac{\partial \widehat{V}}{\partial Q}+\frac{\partial \widehat{V}}{\partial Q}\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial P}\right)\hfill \\ & & -\frac{P}{M}\frac{\partial {\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)}{\partial Q}-\frac{P_{\eta }}{M_{\eta }}\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial Q_{\eta }}+\frac{P_{\eta }}{M_{\eta }}P\frac{\partial \widehat{\chi }(X^{\mathrm{N}},t)}{\partial P}-F_{Q_{\eta }}\frac{\partial {\widehat{W}}^{\mathrm{N}}(X^{\mathrm{N}},t)}{\partial P_{\eta }}.\hfill \end{array}\end{array}$$

(78)

In the adiabatic states defined in Equation (22), Equation (78) reads

$$\frac{d}{dt}{\widehat{W}}_{\alpha {\alpha }^{\prime }}^{\mathrm{N}}(X^{\mathrm{N}},t)=-\sum _{\beta {\beta }^{\prime }}i{\mathcal{L}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathrm{N}}{\widehat{W}}_{\beta {\beta }^{\prime }}^{\mathrm{N}}(X^{\mathrm{N}},t),$$

(79)

where

$$i{\mathcal{L}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathrm{N}}=i{\omega }_{\alpha {\alpha }^{\prime }}{\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+{\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{N}}+{\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}.$$

(80)

We have used the definition of the Bohr frequency ${\omega }_{\alpha {\alpha }^{\prime }}$ in Equation (A3) and of the transition operator ${\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}$ in Equation (A7) in Appendix A. We have introduced a classical-like Nosé–Liouville operator

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill i{\widehat{L}}_{\alpha {\alpha }^{\prime }}^{\mathrm{N}}=& {\displaystyle \frac{P}{M}}\frac{\partial }{\partial Q}+{\displaystyle \frac{1}{2}}\left(F^{\alpha }+F^{{\alpha }^{\prime }}\right)\hfill \\ & -P{\displaystyle \frac{P_{\eta }}{M_{\eta }}}\frac{\partial }{\partial P}+{\displaystyle \frac{P_{\eta }}{M_{\eta }}}\frac{\partial }{\partial Q_{\eta }}+F_{Q_{\eta }}{\displaystyle \frac{\partial }{\partial P_{\eta }}}{\displaystyle \frac{\partial }{\partial P}}.\hfill \end{array}\end{array}$$

(81)

The existence of the stationary operator-valued Nosé quasi-probability function ${\widehat{W}}^{\mathrm{N},\mathrm{e}}(X^{\mathrm{N}})$ is discussed in Appendix C.

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 [142] 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 5. To this end, considering for simplicity a chain of just two thermostat coordinates, one can define the classical phase space point as

$$X^{\mathrm{NHC}}=(R,Q_{{\eta }_1},Q_{{\eta }_2},P,P_{{\eta }_1},P_{{\eta }_2}),$$

(82)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{H}}^{\mathrm{NHC}}=& {\displaystyle \frac{{\widehat{p}}^2}{2m}}+{\displaystyle \frac{P^2}{2M}}+{\displaystyle \frac{P_{{\eta }_1}^2}{2M_{{\eta }_1}}}+{\displaystyle \frac{P_{{\eta }_2}^2}{2M_{{\eta }_2}}}\hfill \\ & +\widehat{V}(\widehat{q},R)+N{k}_{B}T{Q}_{{\eta }_1}+k_{B}T{Q}_{{\eta }_2},\hfill \end{array}\end{array}$$

(83)

where $M_{{\eta }_1}$ and $M_{{\eta }_2}$ are the inertial parameters of the thermostat variables. As shown in Ref. [92,93], one can define an antisymmetric matrix

$${\mathcal{B}}^{\mathrm{NHC}}=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ -1& 0& 0& 0& -P& 0\\ 0& -1& 0& P& 0& -P_{{\eta }_1}\\ 0& 0& -1& 0& P_{{\eta }_1}& 0\end{array}\right].$$

(84)

The matrix ${\mathcal{B}}^{\mathrm{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 [92] are then given by

$$\dot{X}=-X^{\mathrm{NHC}}{\overleftarrow{\nabla }}^{\mathrm{NHC}}{\mathcal{B}}^{\mathrm{NHC}}{\overrightarrow{\nabla }}^{\mathrm{NHC}}{\widehat{H}}^{\mathrm{NHC}}.$$

(85)

Quantum-classical dynamics is then introduced using the matrix super-operator

$${\mathcal{D}}^{\mathrm{NHC}}=\left[\begin{array}{cc}0& 1-\frac{{\overleftarrow{\nabla }}^{\mathrm{NHC}}{\mathcal{B}}^{\mathrm{NHC}}{\overrightarrow{\nabla }}^{\mathrm{NHC}}}{2i{\hslash }^{-1}}\\ -\left(1-\frac{{\overleftarrow{\nabla }}^{\mathrm{NHC}}{\mathcal{B}}^{\mathrm{NHC}}{\overrightarrow{\nabla }}^{\mathrm{NHC}}}{2i{\hslash }^{-1}}\right)& 0\end{array}\right].$$

(86)

The quantum-classical equations of motion can then be written as

$$\begin{array}{ccc}\hfill \frac{d\widehat{\chi }}{dt}& =& \frac{i}{\hslash }\left[\begin{array}{cc}{\widehat{H}}^{\mathrm{NHC}}& \widehat{\chi }\end{array}\right]·{\mathcal{D}}^{\mathrm{NHC}}·\left[\begin{array}{c}{\widehat{H}}^{\mathrm{NHC}}\\ \widehat{\chi }\end{array}\right].\hfill \end{array}$$

(87)

The equations of motion can be represented using the adiabatic basis obtaining the Liouville super-operator

$$\begin{array}{ccc}\hfill i{\mathcal{L}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}^{\mathrm{NHC}}& =& (i{\omega }_{\alpha {\alpha }^{\prime }}+i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{NHC}}){\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}-{\mathcal{T}}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }},\hfill \end{array}$$

(88)

where

$$\begin{array}{ccc}\hfill i{L}_{\alpha {\alpha }^{\prime }}^{\mathrm{NHC}}& =& \frac{P}{M}\frac{\partial }{\partial R}+\frac{1}{2}(F^{\alpha }+F^{{\alpha }^{\prime }})\frac{\partial }{\partial P}+\sum _{k=1}^2(\frac{P_{{\eta }_k}}{M_{{\eta }_k}}\frac{\partial }{\partial Q_{{\eta }_k}}+F_{Q_{{\eta }_k}}\frac{\partial }{\partial P_{{\eta }_k}})-\frac{P_{{\eta }_2}}{M_{{\eta }_2}}{P}_{{\eta }_1}\frac{\partial }{\partial P_{{\eta }_1}},\hfill \end{array}$$

(89)

with $F_{Q_{{\eta }_2}}=(P_{{\eta }_1}^2/M_{{\eta }_1})-k_{B}T$. 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 (A50) and (A52). One has just to replace Equation (A50) for the order zero term with

$$\begin{array}{ccc}\hfill W{\alpha \alpha }^{\mathrm{NHC},\mathrm{e},(0)}& =& \frac{1}{Z^{\mathrm{NHC}}}{e}^{-\beta \left[\frac{P^2}{2M}+E_{\alpha }(R)+{\sum }_{k=1}^2\left(\frac{P_{{\eta }_k}^2}{2M_{{\eta }_k}}\right)+N{k}_{B}T{Q}_{{\eta }_1}+k_{B}T{Q}_{{\eta }_2}\right]}\hfill \end{array}$$

(90)

with an obvious definition of $Z^{\mathrm{NHC}}$.

6. 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 operator-valued quasi-probability function that depend on the coordinates of the bath and whose equation of motion depends on the specific case under study. It is explained how classical spin baths are described in terms of non-canonical coordinates and how this fact leads to the appearance of an off-diagonal open-path geometric phase in the dynamics of the operator-valued quasi-probability function of the system. We then discussed how the effect of thermal baths can be implemented by means of a stochastic, quantum-classical Langevin dynamics and by means of a deterministic, non-Hamiltonian Nosé–Hoover thermostatted dynamics. The formulation of the dynamics in both the spin and Nosé–Hoover case was achieved by generalizing the quasi-Lie bracket of the canonical case.

The formalisms were presented in such a way to shed light on practical implementation via computer simulation algorithms. The particular class of algorithms upon which we focused is based on the unfolding of the evolution of the operator-valued quasi-probability function in terms of piecewise-deterministic trajectories evolving on the adiabatic energy surfaces of the system. 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,110], 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,110]. Within the QCLE approach, the degrees of freedom of the bath are not integrated out of the dynamics but are explicitly taken into account at every time step. Hence, there is no memory function to be approximated and bath properties can be calculated with the same ease with which subsystem properties are computed. The limitations of the QCLE approach are mostly numerical in character and arise in the SSTP algorithm, herein discussed, from the momentum-jump approximation and the accumulation of fluctuating statistical weights associated with the Monte Carlo sampling of the quantum transitions of the subsystem.

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.