Hybrid quantum-classical dynamics of pure-dephasing systems

We consider the interaction dynamics of a classical oscillator and a quantum two-level system for different pure-dephasing Hamiltonians of the type $\widehat{H}(q,p)=H_C(q,p)\boldsymbol{1}+H_I(q,p)\widehat\sigma_z$. This type of systems represents a severe challenge for popular hybrid quantum-classical descriptions. For example, in the case of the common Ehrenfest model, the classical density evolution is shown to decouple entirely from the pure-dephasing quantum dynamics. We focus on a recently proposed hybrid wave equation that is based on Koopman's wavefunction description of classical mechanics. This model retains quantum-classical correlations whenever a coupling potential is present. Here, several benchmark problems are considered and the results are compared with those arising from fully quantum dynamics. A good agreement is found for a series of study cases involving harmonic oscillators with linear and quadratic coupling, as well as time-varying coupling parameters. In all these cases the classical evolution coincides exactly with the oscillator dynamics resulting from the fully quantum description. In the special case of time-independent coupling involving a classical oscillator with varying frequency, the quantum Bloch rotation exhibits peculiar features that escape from the hybrid description. In addition, nonlinear corrections to the harmonic Hamiltonian lead to an overall growth of decoherence at long times, which is absent in the fully quantum treatment.


Introduction
The search for models describing the interaction dynamics of quantum and classical systems goes back to the early days in the history of quantum mechanics and is typically related to the measurement problem. In that context, a quantum system is coupled to a classical system comprising a pointer and the environment, that is the rest of the measuring apparatus. While the environment is usually thought of as a heat bath inducing irreversible thermodynamical processes, the system-pointer interaction is an intrinsically unitary process known as pre-measurement [42,50]. Besides its relevance in such foundational issues, the quantum-classical coupling problem has also emerged over the years in a series of different contexts. In theoretical physics, for example, the possibility that the theory of gravity may not be quantized leads to consider coupling the classical Einstein tensor to quantum matter [12]. In this context, the theory of semiclassical gravity is based on a mean-field model that fails to retain quantum-classical correlations [10]. Furthermore, several activities on quantum-classical coupling are currently motivated by the necessity of mitigating the computational costs of fully quantum many-body simulations. For example, in chemical physics molecular dynamics algorithms often approximate nuclei as classical while retaining the full quantum effects of electron dynamics [15,35,55]. In solid state physics, recent proposals [32] suggest treating the orbital degrees of freedom as classical while leaving spins as fully quantum observables. Similar quantum-classical descriptions have also appeared in spintronics, where quantum spins are coupled to classical ferromagnets thereby reaching unexplored parameter regimes in quantum control [51].

Models of hybrid quantum-classical dynamics
Despite several attempts [5,10,17,25,33,44,57], most hybrid quantum-classical models suffer from various consistency issues [2,10,26,43,47,52]. For example, some models allow for the violation of the uncertainty principle, while some others fail to return uncoupled quantum and classical motion in the absence of an interaction potential. Following [10], we can propose five fundamental properties that quantum-classical dynamical models should possess [23]: 1. the classical system is identified by a positive probability density on phase space at all times; 2. the quantum system is identified by a positive-semidefinite density operatorρ at all times; 3. the model is covariant under both quantum unitary transformations and classical canonical transformations; 4. in the absence of an interaction potential, the model reduces to uncoupled quantum and classical dynamics; 5. in the presence of an interaction potential, the quantum purity Trρ 2 is not a constant of motion (decoherence property).
Notice that the last property rules out the standard mean-field model, which entirely neglects quantum-classical correlations thereby conserving purity. To our knowledge, the only available model satisfying all five properties is the Ehrenfest model : Here, D is the classical Liouville distribution, while the quantum density matrix isρ = Dψψ † dqdp. Also, we have denoted A = ψ, Aψ , while the hybrid Hamiltonian H = H(q, p) is a phase-space function with values in the space of Hermitian operators on the quantum Hilbert space. Despite its popularity in chemical physics [6], the model (1.1) is well known to be unreliable in reproducing the correct decoherence levels arising from the fully quantum description. In particular, the Ehrenfest model suffers from the overcoherence problem, that is it typically leads to higher purity levels than those arising from the fully quantum treatment [3]. In some cases, a certain undercoherence was also reported [30]. As we will see, the Ehrenfest model suffers also from other serious issues that become manifest in the class of Hamiltonians considered in this paper.
Several other hybrid quantum-classical models are available, although they generally fail to fulfill all the five properties listed above. In some cases this lack of consistency leads to unrealistic timescales of the quantum evolution [9]. A candidate model beyond Ehrenfest which satisfies all the five properties was recently presented in [21,22,23]. In that case, the quantum-classical interaction triggers challenging nonlinear terms and suitable closure schemes are currently under development. This nonlinear model was formulated upon revisiting a linear quantum-classical wave equation previously proposed in [9]. This equation governs the unitary evolution of a square-integrable hybrid wavefunction Υ(q, p, x), where x is the quantum position coordinate. In its simplest form, the quantum-classical wave equation reads Satisfying the properties 2)-5) above, this equation was also shown to retain property 1) when the Hamiltonian depends on the quantum degrees of freedom through a set of mutually commuting observables [24] (pure-dephasing Hamiltonian). However, the hybrid wave equation is not known to satisfy property 1) in the general case. In addition, while this equation is covariant under gauge transformations Υ → e iϕ(q,p) Υ, it is not gauge-invariant. This kind of gauge invariance was indicated as a desirable property in [10], based on the fact that classical phases cannot lead to measurable effects [50]. Indeed, this gauge invariance plays a key role in the nonlinear model from [21,22,23].
While these aspects deserve appropriate care, the hybrid wave equation from [9] seems to represent a step forward in the current state of the art of quantum-classical modeling. For example, unlike widely popular models in chemical physics [5,25,35], the hybrid wave equation preserves the Heisenberg uncertainty principle (as a consequence of property 2). Also, it reduces to uncoupled quantum and classical mechanics in the absence of an interaction potential (property 4), unlike alternative approaches based on master equations [17] or quantum hydrodynamics [20,28]. As we will see, it also overcomes the Ehrenfest model in terms of its predictions on the dynamics of the classical system. Last, the hybrid wave equation was recently shown to be Galilean-covariant [7], thereby ensuring another desirable property. Given these aspects of merit, we are motivated to present a benchmark study comparing the results obtained from the quantum-classical wave equation with those arising from fully quantum dynamics. In particular we will consider the simple case of pure-dephasing Hamiltonians, for which the classical density was shown to be always positive-definite [24]. Pure-dephasing problems provide a suitable benchmark platform for our proof-of-principle study of the mathematical model under consideration, while the development of a cost-effective computational method is left for future work. For example, the presence of characteristic curves in its Madelung form makes the quantum-classical wave equation amenable to trajectory-based methods, as opposed to the more expensive (yet more accurate) finite-volume schemes adopted here. This paper will study the case of harmonic and anharmonic oscillators coupled linearly and nonlinearly to a quantum two-level system, and focus especially on the dynamics of the quantum subsystem. In addition, we will present some results on the quantum-classical spinmomentum correlations and we will also consider the case of time-dependent parameters. While predicting slightly higher levels of decoherence, the hybrid results are in good agreement with the quantum dynamics when the hybrid Hamiltonian H(q, p) is a quadratic function of the classical coordinates. One exception is given by the particular case when the oscillator frequency varies with time. In that case, the fully quantum dynamics displays rather peculiar features that cannot be reproduced by the hybrid wave model. Furthermore, in the case of non-harmonic hybrid dynamics, the decoherence becomes particularly pronounced over time thereby leading to distinctive differences from the fully quantum evolution.

Quantum-classical wave equation
The formulation of the quantum-classical wave equation (QCWE) exploits an earlier suggestion from George Sudarshan [14,50]. The idea is to resort to Koopman's wavefunction description of classical dynamics [34,36,56] in such a way that the hybrid quantum-classical motion can be formulated on the tensor product Hilbert space comprising products of classical and quantum wavefunctions. The general theory of the resulting model for hybrid quantum-classical wavefunctions has been presented in several instances [9,22,23,24,54]. At this stage we simply present its hybrid wave equation in the general form: As mentioned previously, Υ(q, p, x) is a wavefunction depending on the classical phase-space coordinates z ≡ (q, p) as well as the quantum coordinate x. The right-hand side of (1.2) deserves some comments. First, the vector potential A = (A q , A p ) is called symplectic potential and its defining relation is where the gradient ∇ = (∂ q , ∂ p ) is computed in the phase space. The symplectic potential is defined up to a pure gradient gauge term. In the standard gauge, one has A = (p, 0) so that in the purely classical case H = H1 the parenthesis on the right-hand side of (1.2) reduces to the phase-space Lagrangian L = p∂ p H − H. Indeed, this term on the right-hand side of (1.2) has the role of incorporating the phase evolution in the classical sector [16,24]. Another relevant gauge is given by In this gauge, the QCWE reads which is the form used throughout this paper. This harmonic gauge is particularly convenient for purely quadratic Hamiltonians since it makes the right-hand side of (1.2) vanish identically.

1.4)
While this expression is generally sign-indefinite, the corresponding distribution was shown to remain positive for pure-dephasing Hamiltonians, that is Hamiltonians depending only on a set of mutually commuting quantum observables [24]. In addition, upon using the exact factorization Υ(q, p, x) = D(q, p)e iS(q,p)/ ψ(x; q, p) [1,22], we observe that the classical phase S(q, p) enters explicitly the expression (1.4) of the classical density, thereby affecting classical averages. This feature may appear counterintuitive since classical phases should not generally lead to observable effects. Indeed, we do not measure Hamilton-Jacobi functions. One can argue that in the present formalism the single terms in (1.4) do not possess any physical meaning per se, and only the full expression identifies a relevant quantity. Alternatively, one can modify the current model by enforcing a gauge principle to treat classical phases as a gauge freedom; this is the approach recently proposed in [22,23].
In the present context, both the quantum density matrix and the classical distribution may be realized as suitable projections of a hybrid von-Neumann-type operator D(q, p) so that ρ c = Tr D andρ =´ D dqdp. Explicitly, one has This quantity is needed to compute quantum-classical expectation values. For example, the total quantum-classical energy is given by Tr( D H) = H and analogously A = Tr( D A) for an arbitrary hybrid observable A(q, p).
A typical initial condition for the QCWE is given by a factorized hybrid wavefunction of the type Υ 0 (q, p, x) = χ(q, p)Ψ(x), where χ and Ψ are quantum and classical wavefunctions, respecitvely. The Koopman wavefunction χ is chosen in such a way that at the initial time ρ c is a Gaussian distribution centered at (q 0 , 0). Then, the expression (1.4) leads to solving a differential equation whose solution leads to [9] Υ 0 (q, p, x) = Ψ(x) e −i pq 0 2 1 2π (1.5) Here, β represents an inverse classical temperature. In the fully quantum description, this initial condition corresponds to the tensor product of a quantum state Ψ(x) with a Gaussian wavepacket centered at (q 0 , 0). While the function (1.5) is suitably normalized, it possesses a long tail in the phase space which arises from a decay of the type 1/H 0 . The slow decay of the Koopman wavefunction represents a challenging aspect for the numerical implementation, which especially intervenes in the computation of expectation values.

Pure-dephasing systems
This paper presents a case of study focusing on the interaction of a classical oscillator with a quantum two-level system. In the fully quantum domain, these models were widely studied in the open systems literature [11], and are often known in optics as variations of the quantum Rabi [31,58] and Jaynes-Cummings models [37]. In the corresponding quantum-classical setting, the hybrid wavefunction is a two-component spinor Υ = (Υ + , Υ − ). In particular, we focus on puredephasing Hamiltonians of the general type [27,48,49] Here, H C is a purely classical Hamiltonian function, while H I is a quantum-classical coupling term and B 0 represents an external magnetic field driving the spin variable. Also, σ z = diag(1, −1) is the third Pauli matrix and 1 is the identity matrix.

Hybrid wave equation for pure-dephasing Hamiltonians
For generic nonlinear Hamiltonians in this class, the QCWE reads Our numerical solution takes advantage of the Madelung representation Υ ± = √ D ± exp (iS ± / ). Then, the densities D ± and phases S ± obey the following equations in the phase space These are the transport equations that will be solved numerically to obtain the quantumclassical hybrid results. The chosen numerical algorithm is a grid-based split-operator technique that advances the functions D ± and S ± alternatively in position space q and in momentum space p, using a finite volume method [19]. Grid-based computational methods possess good accuracy and stability properties, usually better than methods based on trajectories, although they require a higher computational cost. More details are given in the Appendix A. Our initial condition for Υ 0± = √ D 0± exp(iS 0± / ) reads which corresponds to (1.5) with Ψ(x) replaced by the spinor Ψ = (1, 1)/ √ 2. In this Madelung formulation, the classical density is written as: dz, while the quantum density matrix is: where ∆S = S + − S − is the phase difference.

Inadequacy of the Ehrenfest model and other features
Despite their intrinsic simplicity, pure-dephasing systems represent a relevant benchmark test case for two main reasons. First, pure-dephasing systems cannot be adequately modeled by the Ehrenfest equations (1.1), since in that case the classical motion decouples entirely from the quantum dynamics. This is easy to see by noticing that the initial condition ψ| σ z ψ = 0 is preserved in time, so that the classical evolution in (1.1) simply becomes In particular, this indicates the complete absence of quantum backreaction on the classical motion, so that, when H C = (p 2 + q 2 )/2 the classical system undergoes trivial rotations in phase space, in obvious contrast with the behavior arising from the fully quantum treatment.
Another relevant property of pure-dephasing systems is that, whenever H C and H I are quadratic functions, the QCWE can be solved exactly [9,24] thereby providing important information to test the numerical results. In this particular case the evolution of the classical density (1.4) is shown to be identical to the Wigner oscillator dynamics predicted by the corresponding fully quantum problem. Indeed, since the two components Υ ± of the hybrid wavefunction decouple entirely, one shows that the two densities Notice that the densities ρ ± are transported along the characteristics of the time-independent vector field (∂ p H ± , −∂ q H ± ) and thus they remain strictly positive at all times. Then, the sign of the classical density ρ c = ρ + + ρ − is also strictly preserved, thereby indicating the absence of wavefunction nodes that may instead appear in the general case of quantum dynamics.
For the corresponding quantum problem, the Wigner function of the oscillator wavefunction satisfying is given by W = W + + W − . Here, W ± is the Wigner transform of ψ ± and obeys ∂ t W ± = {{H ± , W ± }}, where {{·, ·}} is the Moyal bracket. Then, if H ± is a quadratic function the Moyal bracket reduces to the Poisson bracket, thereby recovering exactly the same classical flow that arises from the QCWE. Notice how this result is in direct contrast with that produced by the Ehrenfest model, for which the oscillator dynamics decouples completely and becomes harmonic.
In the following sections, we will compare the results of the QCWE with those arising from the fully quantum treatment of the Hamiltonian (2.1). For more details on the quantum dynamics, see Appendix B. In particular, we will consider different cases of H C and H I corresponding to harmonic oscillator motion with linear and quadratic coupling (Sec. 3) and two types of anharmonic potentials and coupling (Sec. 4). Regimes with both constant and time-dependent parameters will be investigated.

Harmonic motion
This section presents the results obtained from the numerical implementation of the QCWE for pure-dephasing Hamiltonians of the type (2.1), in the special case when and H I is either linear or quadratic. This hybrid Hamiltonian represents the pure-dephasing interaction of a harmonic oscillator that is coupled linearly or quadratically with a two-level system. In particular, we will compare the results obtained with the QCWE to those obtained from the fully quantum evolution. In the last part of this section, we will also study the case in which either the oscillator frequency ω or the parameters in H I depend explicitly on time.
Notice that, as shown before, in the case of Hamiltonians that are at most quadratic, the classical density evolution obtained from the QCWE coincides exactly with the Wigner distribution dynamics that is obtained from the fully quantum treatment. Thus, the present study will only consider the dynamics of the Bloch vector describing the quantum two-level subsystem.

Linear coupling
Here, we consider in (2.1) a linear interaction Hamiltonian of the type and the parameter values: The initial condition is given by (1.5) with initial shift q 0 = 0 and inverse temperature β = 2. Units in which ω = = m = 1 are used throughout this paper. A fully quantum system corresponding to (2.1), (3.1), and (3.2) was discussed in [45] in the context of electron-phonon coupling. In the case B 0 = 0, the corresponding quantum Hamiltonian identifies the simplest version of the Jahn-Teller problem in chemical physics, known as E ⊗ β Jahn-Teller model [41]. Figure 1 shows the Bloch vector dynamics along with purity evolution, where purity is given by P (t) = Trρ 2 . Since in the present pure-dephasing case σ z is conserved, the motion occurs in the plane σ z = 0. While the former gives an overview of the quantum subsystem dynamics (spin dynamics), the latter quantifies decoherence. Notice the good agreement between the hybrid and fully quantum results, although we see the slightly lower purity level attained by the hybrid system. The latter feature means that the decoherence induced on the two-level subsystem is higher in the case of an interaction with a classical oscillator, compared to a fully quantum oscillator.
Another possibly relevant quantity is the spin current p σ x . Here, we are particularly interested in the quantum-classical spin-momentum correlations p σ x − p σ x whose time evolution is represented in Figure 2. The spin-momentum correlations achieve lower absolute values in the hybrid case. In other words, the correlations generated in time by the fully quantum dynamics are higher than the quantum-classical correlations achieved during the hybrid evolution. In addition, we notice that when the quantum state turns back pure, the spin-momentum correlations vanish for both quantum and hybrid dynamics. Overall, we found that the Bloch dynamics resulting from the QCWE is in good agreement with the predictions of the fully quantum theory, while the correlations between the oscillator and the two-level system reach higher values in the fully quantum case. We recall that the harmonic oscillator evolution is exactly the same in the two cases.

Quadratic coupling
This section focuses on pure-dephasing systems of the type (2.1), where H C is the harmonic oscillator Hamiltonian (3.1) and the interaction term is given by with parameters ω = 1, η = 0.5, In the absence of momentum coupling, the corresponding quantum case was considered in [45]. In the absence of both momentum coupling and magnetic field, an analogue quantumclassical system was approached by the QCWE in [9]. In the case considered here, the quantum correspondent of the coupling Hamiltonian (3.3) has made its appearance in the context of the two-photon Rabi dynamics [4]. The pure-dephasing limit was considered in [18,40,53], although in those cases the magnetic field B 0 appearing in (2.1) is absent. While the inverse temperature of the initial condition is fixed (β = 2), here we will study both cases q 0 = 0 and q 0 = 1 to observe the effects of an initial shift in position. It is perhaps interesting to mention that, if q 0 = 0, switching off the external magnetic field B 0 leads to an absence of Bloch rotation in both the quantum and the hybrid cases. This can be seen explicitly for the QCWE equation upon using the polar form Υ ± = √ D ± e iS ± / in the corresponding wave equations. Then, the QCWE leads to the expectation value expressions σ x , σ y = 2ˆ D + D − cos(∆S/ ), − sin(∆S/ ) dqdp , Figure 2: Spin-momentum correlations (linear coupling). The pale blue curves indicate the evolution of the spin current correlation p σ x − p σ x , for the quantum evolution (solid curve) and the hybrid evolution (dashed curve). Note that in this case p = 0, because no initial shift was present (q 0 = 0).
where ∆S = S + − S − . The latter quantity vanishes at all times due to the initial condition (1.5) and the particular form ∂ t Υ − {H C 1 + H I σ z , Υ} = 0 of the QCWE when B 0 = q 0 = 0. Thus, σ y = 0 remains true at all times and the Bloch vector simply oscillates in amplitude along the x−direction. The same peculiar behavior is also observed in the simulations of the fully quantum dynamics, although its analytical justification is less obvious.
The remainder of this section deals with the case of a nonzero magnetic field, B 0 = 0.3. Figure 3 illustrates the dynamics of the two-level subsystem. Independently of the presence of an initial shift, we observe essentially the same features that had already appeared in the case of linear coupling, although now the difference in purity reaches slightly higher values, meaning that the quadratic coupling leads to stronger decoherence effects in the hybrid case.
In the presence of an initial shift q 0 = 1, the spin current may again be adopted to illustrate the correlation dynamics between the oscillator and the two-level system (see Figure 4). In contrast with the case of linear coupling (Figure 2), the hybrid correlations are now higher than those observed in the full quantum dynamics. We recall, however, that the linear coupling case in Figure 2 did not involve any initial shift, so the two results are not directly comparable. Again, we notice that, when purity reaches its maximum, spin-momentum correlations vanish and this happens for both quantum and hybrid dynamics.
In summary, the case of quadratic coupling does not lead to any substantial differences with respect to the conclusions found when the coupling is linear. The main difference resides in a slightly bigger difference in purity, whose quantum predictions are however still qualitatively reproduced.

Time-dependent parameters
In this section, we keep the form of the Hamiltonians (3.1) and (3.3) and consider the separate cases where either the oscillator frequency ω or the coupling constant η are dependent on time. As we will see, while the case of a variable coupling yields again a good qualitative agreement, the case of variable frequency leads to a very peculiar quantum dynamics that escapes from the QCWE description, which instead reproduces the quantum oscillator dynamics exactly. Notice that a variable oscillator frequency is a less realistic situation than the case of a variable coupling strength.
We start by considering the following temporal expression for ω(t): with ∆ω = 0.4, t 1 = 20, and t 2 = 60. In this case, we ignore the effects of an initial displacement q 0 and we consider a much lower magnetic field B 0 = 0.02 in order to emphasize the effects of the time-dependence. The inverse temperature is still β = 2 and the coupling strength is η = 0.5. This case exhibits a very different qualitative behavior of the Bloch rotation between the hybrid and fully quantum cases, while purity is still satisfactorily reproduced by the QCWE.
In particular, the Bloch rotation in the hybrid case is entirely determined by the magnetic field, while the fully quantum case exhibits an atypical Bloch dynamics that is triggered by a nonzero relative phase of the spinor components. As shown in Figure 5, at the time when the frequency develops a derivative (t 1 = 20), the positive (counterclockwise) Bloch rotation observed in the fully quantum case reverses its sign and becomes very fast until the frequency returns to its original value. After that point in time (t 2 = 60), the Bloch rotation slows down and turns back positive. On the other hand, this remarkable effect cannot be captured by the hybrid QCWE dynamics, for which the rotation remains slow and positive at all times. Instead, the classical density D = Υ 2 resulting from the QCWE reproduces the quantum dynamics of the time-dependent oscillator exactly, as discussed in Section 2.2 In the second simulation, we take a constant oscillator frequency (ω = 1) and a timedependent coupling coefficient η(t) given by with η 0 = 0.5 and ∆η = −0.1. The results are illustrated in Figure 6. Contrarily to the previous situation (constant coupling and varying frequency), here the QCWE reproduces rather well both the Bloch rotation and the purity evolution. Once again, as in earlier cases, the purity predicted by the hybrid dynamics reaches slightly lower values compared to the fully quantum case.

Anharmonic motion
For a Hamiltonian that is at most quadratic, the oscillator dynamics is identical in the quantum and hybrid cases, as was shown in Sec. 2. In such cases, we observed a good agreement also for the quantum two-level subsystem, although distinctive differences emerged in the case of a variable oscillator frequency.
Here, we extend our study to the more challenging scenario of anharmonic motion, for which the quantum and classical dynamics differ at long times. We recall that, as discussed in Section 2.2, the Ehrenfest model (1.1) is inadequate also for anharmonic pure-dephasing dynamics since it leads to a completely decoupled classical evolution. This problem is absent in the QCWE and here we present results for two cases, where either the classical Hamiltonian H C (q, p) or the coupling term H I (q, p) involve some nonlinear terms. Two types of nonlinear correction terms are considered: a quartic term and a rational term of the type αq 4 /(1 + αq 2 ). The latter belongs to a wider class of potentials known as Mitra potentials [13,39]. Among others, Mitra potentials were used in the calculation of bond state spectra; see [46] for this and other occurrences of Mitra potentials.

Quartic potential correction
In this section, we modify the hybrid dynamics in Section 3.2 by adding a quartic potential correction, that is while the coupling terms remain quadratic: H I (p, q) = 1 2 η(q 2 − p 2 ) + B 0 . We recall that the quantum dynamics arising from a quartic potential term differs substantially from its classical counterpart [38]. Here, we use the following parameters: oscillator frequency ω = 1, quadratic coupling η = 0.3, external field B 0 = 0.1, and correction parameter = 0.01. The initial condition is still the one of Eq. (1.5), with inverse temperature β = 2.
The results for the Bloch and purity dynamics of the quantum subsystem are displayed in Fig. 7. The purity oscillates while decreasing to a lower value in the hybrid case with respect to the fully quantum case. This result should be compared to the one obtained for a purely quadratic Hamiltonian, see Fig. 3: in that case, the hybrid purity was following more closely the fully quantum one, both returning to values very close to unity. In other words, in this case of quartic potential correction we observe that the decoherence in the hybrid dynamics tends to become more and more pronounced over time, while its overall average levels do not seem to vary in the fully quantum case. Also, we notice a faster Bloch rotation in the hybrid case with respect to the fully quantum result.
Due to the slow decay of the hybrid wavefunction (1.5) and the rapid growth of the quartic potential ∼ q 4 , the implementation of the hybrid dynamics corresponding to the classical Hamiltonian (4.1) is particularly challenging and the long-time simulations do not allow a satisfactory energy conservation. In order to test the long-time hybrid dynamics, we have considered a rational correction term, which is tailored to grow asymptotically as ∼ q 2 . This is discussed in the next section.

Rational correction terms
In this case, we modify the dynamics from Section 3.2 by adding a rational correction term [13,39,46] αq 4 /(1 + 2αq 2 ) to either the classical Hamiltonian H C or the coupling H I . Hence, Figure 7: Dynamics of the quantum subsystem (quartic potential). Bloch and purity dynamics (left and right panel, respectively) corresponding to the pure-dephasing dynamics for a harmonic oscillator with a quartic potential correction.
we consider hybrid Hamiltonians (2.1) with either For this set of simulations, we use the following parameters: oscillator frequency ω = 1, quadratic coupling η = 0.5, external field B 0 = 0.1, and inverse temperature β = 2. In addition, we set α = 0.1. Compared to the quartic correction of the preceding section, which dominates at large values of q, the rational potential is strongly anharmonic at small values of q while it becomes again quadratic for q → ∞.
In the first case, we consider a correction of the classical potential ( = 0.1 and c = 0), while in the second case we insert a correction in the coupling ( = 0 and c = 0.05). The results are plotted respectively in Fig. 8 and Fig. 9, for short times (t ∈ [0, 20], top panels) and long times (t ∈ [0, 100], bottom panels). For both cases, the short-time hybrid dynamics reproduces relatively well its fully quantum counterpart. Nevertheless, we note a significant difference: for a correction in the classical potential, the Bloch rotation appears faster in the hybrid dynamics than in the quantum evolution, whereas the opposite occurs when the correction applies to the coupling term. This is similar to what was found in the case of a quartic correction. At longer times, the quantum dynamics found by the hybrid model deviates from the fully quantum predictions. The average value of purity appears to decay again, although this decay is slower compared to the case of a corrected classical potential.
Finally, we also present a comparison of the uncertainty evolution in the quantum and hybrid cases; see Fig. 10. For short times, the hybrid predictions are in excellent agreement with the fully quantum results. For long times, the behaviour is very different depending on whether the rational correction term is inserted in the potential or in the coupling. In the first case, we observe uncertainty oscillations with a much smaller amplitude in the case of the Figure 8: Dynamics of the quantum subsystem (rational potential correction). Bloch and purity dynamics (left and right panels, respectively) corresponding to the pure-dephasing dynamics for a harmonic oscillator with a rational potential correction. The top panels illustrate the short-time dynamics, while the bottom panels represent the long-time motion.
hybrid dynamics, while in the second case this difference in amplitude is less pronounced so that the hybrid model still yields a certain agreement with the fully quantum dynamics.

Conclusions
A fully quantum description of complex systems is often a fierce challenge, both for analytical developments and numerical calculations, especially for objects containing as many as hundreds or thousands of particles. In those cases, it may be necessary to separate the system under study into various subsystems, some of which are treated classically and some quantum mechanically. For instance, we may consider a quantum object immersed in a solution of water molecules that are described classically.
However, quantum and classical motion are governed by very different mathematical descriptions, i.e., point trajectories in the phase space for classical systems and wave equations Figure 9: Dynamics of the quantum subsystem (rational coupling correction). Bloch and purity dynamics (left and right panel, respectively) corresponding to the pure-dephasing dynamics for a harmonic oscillator with a rational coupling correction. The top panels illustrate the short-time dynamics, while the bottom panels represent the long-time motion.
in a Hilbert space for the quantum dynamics. Hence, coupling a quantum system to a classical one is by no means an easy task. Several approaches to couple quantum and classical systems have been proposed in the past, but each of them presents some drawbacks. Pure-dephasing Hamiltonians are an especially challenging test case for hybrid quantum-classical models and popular approaches are known to struggle. An example is given by the Ehrenfest model, in which the classical system decouples completely from the quantum degrees of freedom.
Here, we have used pure-dephasing systems as a benchmark test ground for a recently developed method, jointly proposed by one of us [9]. The method is based on the theory of Koopman wavefunctions, which is a representation of classical mechanics based on the same Hilbert-space formalism as in quantum mechanics (self-adjoint operators, unitary evolution, etc. . . ). This feature makes it a compelling candidate to represent hybrid systems that are partly quantum and partly classical. The resulting quantum-classical wave equation (QCWE) governs the dynamics of a hybrid wavefunction comprising both quantum and classical properties.
Our cases of study focussed on the hybrid dynamics of a classical oscillator coupled to a Figure 10: Uncertainty evolution (rational potential/coupling correction). Evolution of the uncertainty q 2 p 2 for the case of rational potential correction (left) and rational coupling correction (right). Once again, solid lines correspond to the fully quantum result, while dashed lines identify the hybrid dynamics.
quantum two-level system. In the fully quantum case involving a harmonic oscillator, these test cases emerge as limit cases of the Jaynes-Cummings model [8] (large-detuning approximation) or, more generally, Rabi models. See also [45] for the experimental relevance of this type of dynamics in the context of electron-phonon coupling at low temperatures. The hybrid Koopman equations were implemented into a numerical code that follows the evolution of the hybrid wave function (or rather, its amplitude and phase) in the classical phase space. The long tail of the initial condition (1.5) requires special care in the numerical implementation.
For several Hamiltonians involving a harmonic oscillator with linear or quadratic coupling, an overall good agreement was found between the QCWE model and the fully quantum description. For example, the classical evolution coincides exactly with the oscillator dynamics predicted by the fully quantum treatment. Also, the Bloch sphere motion and the purity evolution (quantifying quantum decoherence) appear to be well reproduced by the hybrid model. This good agreement persists in the case of a time-dependent coupling, while the case of a time-varying classical frequency leads to very peculiar Bloch rotation dynamics in the fully quantum description, and this feature escapes from the quantum-classical treatment.
For non-quadratic Hamiltonians, the quantum and hybrid dynamics agreed sufficiently well for short times. For long times, the quantum purity tends to decrease in the hybrid description, while its average levels appear unvaried in the fully quantum treatment. This means that we observed a higher level of decoherence in the hybrid description. Also, unlike the predictions of the Ehrenfest model, the QCWE retains the quantum backreaction force on the classical evolution in all the considered cases.
The present results suggest that the Koopman approach represents a valuable step forward in modeling hybrid quantum-classical dynamics. Indeed, while overcoming several outstanding challenges, this approach has recently led to an upgrade model that, unlike the QCWE used here, is successful in retaining the positivity of the classical Liouville density while treating classical phases as unmeasurable quantities [21,22,23]. Current computational efforts are focusing also on this more advanced approach. For now, we emphasize that the present work is a proof-of-principle study of the QCWE mathematical model and we do not dwell upon the computational complexity of our numerical implementation. On the one hand, the presence of characteristic curves for the evolution of the density D = Υ 2 seems to offer some advantages over fully quantum treatments and certain quantum-classical descriptions [35]. On the other hand, the long tail in the initial condition (1.5) may pose new challenges, e.g. involving the statistical sampling. Nevertheless, we believe that the QCWE and its recent nonlinear upgrade provide a valuable platform for the identification of new convenient multiscale tools in computational physics and chemistry. In the context of nonadiabatic molecular dynamics, for example, "the parallel development of rigorous but computationally expensive methods and more approximate but computationally efficient methods" is regarded as absolutely "critical", as discussed in [29]. Other questions also involve the level of applicability of the quantumclassical treatment itself, which may or may not always be a good approximation depending on the fully quantum problem under study. It appears that little is known about this particular problem beyond the standard literature on Born-Oppenheimer molecular dynamics in the adiabatic regime.
Future work will consider more realistic models, beyond the pure-dephasing Hamiltonians studied here. It will also be interesting to consider more complex classical subsystems, e.g., comprising many interacting classical particles. 62210 from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation.

A Numerical solution of the phase-space equations
The transport equations (2.3)-(2.4) for the hybrid model both take the following form: where V (p) = ∂ p H and U (q) = −∂ q H. These are Vlasov equations, well-known in the plasma physics community. They can be solved by projecting the phase-space density f (q, p, t) on a phase-space grid f ij (t) = f (q i , p j , t) and using a split-operator method for the time stepping.
The split-operator technique amounts to solving in sequence the following equations: Although the shifts (A.7) are exact, they will not generally map one grid point onto another, hence the need for an interpolation technique. Several methods have been used in the past, ranging from cubic splines to spectral methods. Here, we have used a finite-volume scheme, which has the advantage of locally conserving the transported quantity exactly [19]. For most simulations reported in this work, we used a computational box [−25, 25] × [−25, 25] (position space × momentum space), with N q = 1000 points in position space and N p = 1000 points in momentum space. The time step was typically ∆t = 0.01.
Grid-based codes generally display good accuracy properties, as they mesh the entire phase space with the same resolution, so that even regions of low phase-space density are well sampled. This is in contrast with trajectory-based methods which sample the phase-space density using pointlike marker "particles": in regions where the density is low the number of markers is also low, leading to poor statistical sampling. In addition, since each of the steps (A.2)-(A.6) has an exact solution in time, grid-based methods do not suffer from limitations due to the Courant-Friedrichs-Lewy (CFL) conditions.
The main limitation of grid-based codes applied to Vlasov-like equations such as (A.1) is their significant computational cost, due to the necessity of meshing the entire phase space, which can become a hindrance for high-dimensional configurations. In the present context, this issue is even more compelling because the initial condition (1.5) decays algebraically instead of exponentially, as is commonly the case for initial conditions described by a Maxwell-Boltzmann density.

B Solution of the Schrödinger equations
If the Hamiltonian is at most quadratic in the phase space variables (q, p), then it is possible to find a semi-analytical solution to the Schrödinger equation: i∂ t ψ ± = H ± ψ ± (B.1) for the spinor ψ = (ψ + , ψ − ) T . For the case of quadratic coupling considered in the main text, the Hamiltonian is: It can be shown that an exact solution to the time-dependent Schrödinger equation takes the Gaussian form provided the initial condition is also in this form. In the above expression, the auxiliary functions a ± (t), b ± (t), c ± (t), d ± (t) and e ± (t) are real functions that depend on the time t, and A ± are normalization constants chosen so that |A + | 2 + |A − | 2 = 1.
By injecting the Gaussian wave function (B.3) into the Schrödinger equation (B.1) and equating the coefficients of the terms in q 0 , q 1 and q 2 , one obtains that the above auxiliary functions should obey the following set of ordinary differential equations: where a dot stands for differentiation with respect to time. As initial condition at t = 0, we take: b ± (0) = 1/2 (minimum uncertainty wave packet), d ± (0) = q 0 and e ± (0) = 0 (Gaussian wave packet centered at q 0 ), and a ± (0) = c ± (0) = 0 (vanishing initial phase).
The equations (B.4) form a system of nonlinear ODEs, which can be solved numerically with little effort, using for instance a standard fourth-order Runge-Kutta scheme. Once the auxiliary functions are known (numerically) for all times, it is easy to reconstruct any observable from the wave function (B.3) at time t. For instance, the mean position is: and the z component of the spin which also shows that σ z is a constant of the motion, as expected. For other variables, the expressions are more convoluted, but can still be easily computed.
For Hamiltonians other than quadratic, the above procedure cannot be used and a fully numerical solution of the Schrödinger equation (B.1) has to be envisaged. For the anharmonic Hamiltonians considered here, we chose a standard Crank-Nicolson scheme, which is secondorder accurate in time and space.