# Hybrid quantum-classical dynamics of pure-dephasing systems *

Giovanni Manfredi, Antoine Rittaud, Cesare Tronci

► To cite this version:

Giovanni Manfredi, Antoine Rittaud, Cesare Tronci. Hybrid quantum-classical dynamics of pure-dephasing systems *. Journal of Physics A: Mathematical and Theoretical, 2023, 56 (15), pp.154002.
10.1088/1751-8121/acc21e hal-04252435

HAL Id: hal-04252435

https://hal.science/hal-04252435v1

Submitted on 20 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
PAPER • OPEN ACCESS

# Hybrid quantum-classical dynamics of pure-dephasing systems*

To cite this article: Giovanni Manfredi et al 2023 *J. Phys. A: Math. Theor.* **56** 154002

View the [article online](#) for updates and enhancements.

You may also like

- On two recent proposals for witnessing nonclassical gravity

Michael J W Hall and Marcel Reginatto

- Entanglement dynamics in a quantum-classical hybrid of two q-bits and one oscillator

L Fratino, A Lampo and H-T Elze

- Dissipative dynamics of an open quantum battery

M Carrega, A Crescente, D Ferraro et al.
# Hybrid quantum-classical dynamics of pure-dephasing systems*

Giovanni Manfredi¹ ID, Antoine Rittaud¹ and Cesare Tronci²,³,** ID

¹ Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, 67000 Strasbourg, France

² Department of Mathematics, University of Surrey, GU2 7XH Guildford, United Kingdom

³ Department of Physics and Engineering Physics, Tulane University, New Orleans, LA 70118, United States of America

E-mail: c.tronci@surrey.ac.uk

Received 31 October 2022; revised 13 February 2023

Accepted for publication 7 March 2023

Published 20 March 2023

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

We consider the interaction dynamics of a classical oscillator and a quantum two-level system for different pure-dephasing Hamiltonians of the type $\hat{H}(q,p) = H_C(q,p)\mathbf{1} + H_I(q,p)\hat{\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

<figure><img src="image_3.png" /></figure>

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote>* Contribution to the Special Collection 'Koopman methods in classical and quantum-classical mechanics'.</footnote>

<footnote>** Author to whom any correspondence should be addressed.</footnote>
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.

Keywords: Koopman wavefunction, mixed quantum–classical dynamics, pure-dephasing systems

(Some figures may appear in colour only in the online journal)

# 1. 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].

## 1.1. *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 $\hat{\rho}$ 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 $\hat{\rho}^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:

$$ \frac{\partial D}{\partial t} + \left\langle \frac{\partial \hat{H}}{\partial p} \right\rangle \frac{\partial D}{\partial q} - \left\langle \frac{\partial \hat{H}}{\partial q} \right\rangle \frac{\partial D}{\partial p} = 0, \quad \frac{\partial \psi}{\partial t} + \left\langle \frac{\partial \hat{H}}{\partial p} \right\rangle \frac{\partial \psi}{\partial q} - \left\langle \frac{\partial \hat{H}}{\partial q} \right\rangle \frac{\partial \psi}{\partial p} = -\frac{i}{\hbar} \hat{H} \psi. \quad (1.1) $$

Here, D is the classical Liouville distribution, while the quantum density matrix is $\hat{\rho} = \int D\psi\psi^\dagger dqdp$. Also, we have denoted $\langle \hat{A} \rangle = \langle \psi, \hat{A}\psi \rangle$, while the hybrid Hamiltonian $\hat{H} = \hat{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–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 (QCWE) previously proposed in [9]. This equation governs the unitary evolution of a square-integrable hybrid wavefunction $\Upsilon(q,p,x)$, where $x$ is the quantum position coordinate. In its simplest form, the QCWE reads $\partial_t\Upsilon + \{\Upsilon, \hat{H}\} = i\hbar^{-1}(p\partial_p\hat{H} - \hat{H})\Upsilon$. 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 $\Upsilon \mapsto e^{i\varphi(q,p)}\Upsilon$, 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–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 QCWE 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 QCWE 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 spin-momentum 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 $\hat{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.

## 1.2. Quantum-classical wave equation (QCWE)

The formulation of the 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–24, 54]. At this stage we simply present its hybrid wave equation in the general form:

$$ \frac{\partial \Upsilon}{\partial t} - \{\hat{H}, \Upsilon\} = \frac{i}{\hbar} \left( \mathcal{A}_q \frac{\partial \hat{H}}{\partial p} - \mathcal{A}_p \frac{\partial \hat{H}}{\partial q} - \hat{H} \right) \Upsilon. \quad (1.2) $$

As mentioned previously, $\Upsilon(q,p,x)$ is a wavefunction depending on the classical phase-space coordinates $z \equiv (q,p)$ as well as the quantum coordinate $x$. The right-hand side of equation (1.2) deserves some comments. First, the vector potential $\mathcal{A} = (\mathcal{A}_q, \mathcal{A}_p)$ is called symplectic potential and its defining relation is

$$ \nabla \mathcal{A} - (\nabla \mathcal{A})^T = -\mathcal{J}, \quad \text{with} \quad \mathcal{J} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}, $$

where the gradient $\nabla = (\partial_q, \partial_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 $\mathcal{A} = (p, 0)$ so that in the purely classical case $\hat{H} = H1$ the parenthesis on the right-hand side of equation (1.2) reduces to the phase-space Lagrangian $\mathcal{L} = p\partial_p H - H$. Indeed, this term on the right-hand side of equation (1.2) has the role of incorporating the phase evolution in the classical sector [16, 24]. Another relevant gauge is given by

$$ \mathcal{A} = (p, 0) - \frac{1}{2} \nabla(pq) = \frac{1}{2}(p, -q). $$

In this gauge, the QCWE reads

$$ \frac{\partial \Upsilon}{\partial t} - \{\hat{H}, \Upsilon\} = \frac{i}{\hbar} \left( \frac{p}{2} \frac{\partial \hat{H}}{\partial p} + \frac{q}{2} \frac{\partial \hat{H}}{\partial q} - \hat{H} \right) \Upsilon, $$
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 equation (1.2) vanish identically.

In the QCWE setting, the quantum density matrix is given by

$$ \hat{\rho} = \int \Upsilon(q,p) \Upsilon^{\dagger}(q,p) dqdp, \qquad (1.3) $$

which is consistently positive definite thereby ensuring that the uncertainty principle is satisfied at all times. Here, the notation is such that $\Upsilon^{\dagger}(q,p)$ is the quantum adjoint, so that both $\Upsilon_1^{\dagger}\Upsilon_2 = \langle\Upsilon_1|\Upsilon_2\rangle$ and $\Upsilon^{\dagger}\Upsilon = \|\Upsilon\|^2$ identify distributions on the phase space. In addition, the classical distribution reads

$$ \rho_c(q,p) = \|\Upsilon(q,p)\|^2 - \mathrm{div}(\|\Upsilon(q,p)\|^2 \mathcal{J} \mathcal{A}(q,p)) + \hbar \mathrm{Im}\{\Upsilon^{\dagger}(q,p), \Upsilon(q,p)\}. \qquad (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 $\Upsilon(q,p,x) = \sqrt{D(q,p)}e^{iS(q,p)/\hbar}\psi(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 equation (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 $\hat{D}(q,p)$ so that $\rho_c = \mathrm{Tr}\hat{D}$ and $\hat{\rho} = \int \hat{D} dqdp$. Explicitly, one has

$$ \hat{D} = \Upsilon\Upsilon^{\dagger} - \mathrm{div}(\mathcal{J}\mathcal{A}\Upsilon\Upsilon^{\dagger}) + i\hbar\{\Upsilon, \Upsilon^{\dagger}\}. $$

This quantity is needed to compute quantum–classical expectation values. For example, the total quantum–classical energy is given by $\mathrm{Tr}(\hat{D}\hat{H}) = \langle\hat{H}\rangle$ and analogously $\langle\hat{A}\rangle = \mathrm{Tr}(\hat{D}\hat{A})$ for an arbitrary hybrid observable $\hat{A}(q,p)$.

A typical initial condition for the QCWE is given by a factorized hybrid wavefunction of the type $\Upsilon_0(q,p,x) = \chi(q,p)\Psi(x)$, where $\chi$ and $\Psi$ are quantum and classical wavefunctions, respectively. The Koopman wavefunction $\chi$ is chosen in such a way that at the initial time $\rho_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]

$$ \Upsilon_0(q,p,x) = \Psi(x) e^{-i \frac{\rho_0}{2\hbar} \sqrt{\frac{1 - (1 + \beta H_0) e^{-\beta H_0}}{\beta H_0^2}}}, \quad \text{with} \quad H_0 = \frac{1}{2} p^2 + \frac{1}{2} (q - q_0)^2. \qquad (1.5) $$

Here, $\beta$ represents an inverse classical temperature. In the fully quantum description, this initial condition corresponds to the tensor product of a quantum state $\Psi(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.
## 2. 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 $\Upsilon = (\Upsilon_+, \Upsilon_-)$. In particular, we focus on pure-dephasing Hamiltonians of the general type [27, 48, 49]

$$ \hat{H}(q,p) = H_C(q,p)\mathbf{1} + (H_I(q,p) + B_0)\hat{\sigma}_z. \quad (2.1) $$

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, $\hat{\sigma}_z = \text{diag}(1, -1)$ is the third Pauli matrix and $\mathbf{1}$ is the identity matrix.

### 2.1. Hybrid wave equation for pure-dephasing Hamiltonians

For generic nonlinear Hamiltonians in this class, the QCWE reads

$$ i\hbar \frac{\partial \Upsilon_{\pm}}{\partial t} = i\hbar\{H_{\pm}, \Upsilon_{\pm}\} - \mathcal{L}_{\pm}\Upsilon_{\pm}, \quad (2.2) $$

with $H_{\pm} = H_C \pm (H_I + B_0)$ and $\mathcal{L}_{\pm}(\mathbf{z}) = \mathbf{z} \cdot \nabla H_{\pm}(\mathbf{z}) / 2 - H_{\pm}(\mathbf{z})$. Our numerical solution takes advantage of the Madelung representation $\Upsilon_{\pm} = \sqrt{D_{\pm}} \exp(iS_{\pm}/\hbar)$. Then, the densities $D_{\pm}$ and phases $S_{\pm}$ obey the following equations in the phase space

$$ \partial_t D_{\pm} = \{H_{\pm}, D_{\pm}\}, \quad (2.3) $$

$$ \partial_t S_{\pm} = \{H_{\pm}, S_{\pm}\} + \mathcal{L}_{\pm}. \quad (2.4) $$

These are the transport equations that will be solved numerically to obtain the quantum-classical hybrid results. The chosen numerical algorithm is a grid-based split-operator technique that advances the functions $D_{\pm}$ and $S_{\pm}$ 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 $\Upsilon_{0\pm} = \sqrt{D_{0\pm}} \exp(iS_{0\pm}/\hbar)$ reads

$$ D_{0\pm} = \frac{1}{2\pi} \frac{1 - (1 + \beta H_0) e^{-\beta H_0}}{\beta H_0^2}, \quad S_{0\pm} = -\frac{pq_0}{2}, \quad \text{with} \quad H_0 = \frac{1}{2} p^2 + \frac{1}{2} (q - q_0)^2 $$

which corresponds to equation (1.5) with $\Psi(x)$ replaced by the spinor $\Psi = (1, 1)/\sqrt{2}$. In this Madelung formulation, the classical density is written as:

$$ \rho_c(\mathbf{z}) = \sum_{\pm} \left[ D_{\pm} + \text{div} \left( D_{\pm} \mathcal{J} \nabla S_{\pm} + \frac{\mathbf{z}}{2} D_{\pm} \right) \right], \quad (2.5) $$

with $\int \rho_c(\mathbf{z})d\mathbf{z} = \int [D_+(\mathbf{z}) + D_-(\mathbf{z})]d\mathbf{z}$, while the quantum density matrix is:

$$ \begin{aligned} \hat{\rho} &= \int \Upsilon(\mathbf{z}) \Upsilon^{\dagger}(\mathbf{z}) d\mathbf{z} = \int d\mathbf{z} \begin{pmatrix} |\Upsilon_+|^2 & \Upsilon_+ \Upsilon_-^* \\ \Upsilon_+^* \Upsilon_- & |\Upsilon_-|^2 \end{pmatrix} \\ &= \int d\mathbf{z} \begin{pmatrix} \frac{D_+}{\sqrt{D_+ D_-} e^{-i\Delta S/\hbar}} & \sqrt{D_+ D_-} e^{i\Delta S/\hbar} \\ \sqrt{D_+ D_-} e^{-i\Delta S/\hbar} & D_-, \end{pmatrix} \end{aligned} \quad (2.6) $$

where $\Delta S = S_+ - S_-$ is the phase difference.
## 2.2. 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 equation (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 $\langle \psi | \hat{\sigma}_z \psi \rangle = 0$ is preserved in time, so that the classical evolution in equation (1.1) simply becomes

$$ \partial_t D + \{D, H_C\} = 0. $$

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 $\Upsilon_{\pm}$ of the hybrid wavefunction decouple entirely, one shows that the two densities

$$ \rho_{\pm} = |\Upsilon_{\pm}|^2 + \frac{1}{2} \mathrm{div}(\mathbf{z} | \Upsilon_{\pm}|^2) + \hbar \mathrm{Im}\{\Upsilon_{\pm}^*, \Upsilon_{\pm}\} $$

obey a classical Liouville equation: $\partial_t \rho_{\pm} = \{H_{\pm}, \rho_{\pm}\}$. Notice that the densities $\rho_{\pm}$ are transported along the characteristics of the time-independent vector field $(\partial_p H_{\pm}, -\partial_q H_{\pm})$ and thus they remain strictly positive at all times. Then, the sign of the classical density $\rho_c = \rho_+ + \rho_-$ 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

$$ i\hbar\partial_t\psi = \hat{H}_C\psi + (\hat{H}_I + B_0)\hat{\sigma}_z\psi \quad (2.7) $$

is given by $W = W_+ + W_-$. Here, $W_{\pm}$ is the Wigner transform of $\psi_{\pm}$ and obeys $\partial_t W_{\pm} = \{\{H_{\pm}, W_{\pm}\}\}$, where $\{\cdot, \cdot\}$ is the Moyal bracket. Then, if $H_{\pm}$ 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 (section 3) and two types of anharmonic potentials and coupling (section 4). Regimes with both constant and time-dependent parameters will be investigated.
### 3. 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

$$H_C = \frac{1}{2} (p^2 + \omega^2 q^2) \qquad (3.1)$$

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 $\omega$ 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.

#### 3.1. Linear coupling

Here, we consider in (2.1) a linear interaction Hamiltonian of the type

$$H_I = \gamma q \qquad (3.2)$$

and the parameter values:

$$\omega = 1, \quad \gamma = 0.5, \quad B_0 = 0.2.$$

The initial condition is given by equation (1.5) with initial shift $q_0 = 0$ and inverse temperature $\beta = 2$. Units in which $\omega = \hbar = m = 1$ are used throughout this paper. A fully quantum system corresponding to equations (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 \otimes \beta$ Jahn-Teller model [41].

Figure 1 shows the Bloch vector dynamics along with purity evolution, where purity is given by $P(t) = \text{Tr} \hat{\rho}^2$. Since in the present pure-dephasing case $\langle \sigma_z \rangle$ is conserved, the motion occurs in the plane $\langle \sigma_z \rangle = 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 $\langle p \hat{\sigma}_x \rangle$. Here, we are particularly interested in the quantum-classical spin-momentum correlations $\langle p \hat{\sigma}_x \rangle - \langle p \rangle \langle \hat{\sigma}_x \rangle$ 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.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Dynamics of the quantum subsystem (linear coupling). The left panel shows the evolution of the Bloch vector in the equatorial disk $\langle\sigma_z\rangle = 0$: the purple color corresponds to earlier times, while yellow indicates later times $t = 20$. The dashed line corresponds to the hybrid evolution, while the thick solid line represents the fully quantum dynamics. The right panel shows the evolution of purity with the same convention. The same graphical conventions and color code will be used throughout the remaining figures.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Spin-momentum correlations (linear coupling). The pale blue curves indicate the evolution of the spin current correlation $\langle p\hat{\sigma}_x\rangle - \langle p\rangle\langle\hat{\sigma}_x\rangle$, for the quantum evolution (solid curve) and the hybrid evolution (dashed curve). Note that in this case $\langle p\rangle = 0$, because no initial shift was present ($q_0 = 0$).</figcaption></figure>

## 3.2. 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

$$H_I = \frac{\eta}{2}(q^2 - p^2), \qquad (3.3)$$
with parameters

$$ \omega = 1, \quad \eta = 0.5, \quad B_0 = 0.3. $$

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 quantum-classical 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 equation (2.1) is absent. While the inverse temperature of the initial condition is fixed ($\beta = 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 $\Upsilon_{\pm} = \sqrt{D_{\pm}}e^{iS_{\pm}/\hbar}$ in the corresponding wave equations. Then, the QCWE leads to the expectation value expressions

$$ \langle\hat{\sigma}_x\rangle, \langle\hat{\sigma}_y\rangle = 2 \int \sqrt{D_+D_-} (\cos(\Delta S/\hbar), -\sin(\Delta S/\hbar)) dqdp, $$

where $\Delta S = S_+ - S_-$. The latter quantity vanishes at all times due to the initial condition (1.5) and the particular form $\partial_t\Upsilon - \{H_C\mathbf{1} + H_I\hat{\sigma}_z, \Upsilon\} = 0$ of the QCWE when $B_0 = q_0 = 0$. Thus, $\langle\hat{\sigma}_y\rangle = 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.

### 3.3. 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 $\omega$ or the coupling constant $\eta$ 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.
<figure><img src="image_6.png" /><figcaption>Figure 3. Dynamics of the quantum subsystem (quadratic coupling). Bloch and purity dynamics corresponding to the quadratic coupling in equation (3.3). The upper panels refer to the case $q_0 = 0$, while the lower panels represent the case $q_0 = 1$.</figcaption></figure>

We start by considering the following temporal expression for $\omega(t)$:

$$ \omega(t) = 1 + \frac{1}{2} \Delta\omega [\tanh(t - t_1) - \tanh(t - t_2)], \quad (3.4) $$

with $\Delta\omega = 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 $\beta = 2$ and the coupling strength is $\eta = 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 = \|\Upsilon\|^2$ resulting from the QCWE reproduces the quantum dynamics of the time-dependent oscillator exactly, as discussed in section 2.2.
<figure><img src="image_7.png" /><figcaption>Figure 4. Spin-momentum correlations (quadratic coupling). Same as figure 2, although now the curves are obtained in the case of quadratic coupling with an initial position displacement $q_0 = 1$. Notice that $\langle p \rangle$ does not vanish in this case.</figcaption></figure>

In the second simulation, we take a constant oscillator frequency ($\omega = 1$) and a time-dependent coupling coefficient $\eta(t)$ given by

$$ \eta(t) = \eta_0 + \frac{1}{2} \Delta\eta [\tanh(t - t_1) - \tanh(t - t_2)], $$

with $\eta_0 = 0.5$ and $\Delta\eta = -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.

## 4. 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 section 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 $\alpha q^4/(1 + \alpha 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.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Dynamics of the quantum subsystem. (a) Bloch dynamics in the equatorial disk $\langle\sigma_z\rangle = 0$; (b) dynamics of the single components $\langle\sigma_x\rangle$ and $\langle\sigma_y\rangle$; and (c) purity dynamics corresponding to the Hamiltonians (3.1) and (3.3). The time-dependent frequency $\omega(t)$ in equation (3.4) is presented in panel (c) as a black solid line.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Dynamics of the quantum subsystem. Bloch (left panel) and purity (right panel) dynamics corresponding to the Hamiltonians (3.1) and (3.3), with a time-dependent coupling constant $\eta(t)$ whose evolution is presented in the right panel as a black solid line (see main text for the analytical formula). The results from the hybrid (dashed lines) and the fully quantum (solid lines) dynamics exhibit overall qualitative agreement.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>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.</figcaption></figure>

## 4.1. Quartic potential correction

In this section, we modify the hybrid dynamics in section 3.2 by adding a quartic potential correction, that is

$$H_C(p,q) = \frac{p^2 + \omega^2 q^2}{2} + \epsilon \frac{q^4}{4}, \quad (4.1)$$

while the coupling terms remain quadratic: $H_I(p,q) = \frac{1}{2}\eta(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 $\omega = 1$, quadratic coupling $\eta = 0.3$, external field $B_0 = 0.1$, and correction parameter $\epsilon = 0.01$. The initial condition is still the one of equation (1.5), with inverse temperature $\beta = 2$.

The results for the Bloch and purity dynamics of the quantum subsystem are displayed in figure 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 figure 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 $\sim 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 $\sim q^2$. This is discussed in the next section.

## 4.2. Rational correction terms

In this case, we modify the dynamics from section 3.2 by adding a rational correction term [13, 39, 46] $\alpha q^4/(1 + 2\alpha q^2)$ to either the classical Hamiltonian $H_C$ or the coupling $H_I$. Hence, we consider hybrid Hamiltonians (2.1) with either
<figure><img src="image_11.png" /><figcaption>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.</figcaption></figure>

$$H_C(p,q) = \frac{p^2 + \omega^2 q^2}{2} + \epsilon \frac{\alpha q^4}{1 + 2\alpha q^2}, \qquad H_I(p,q) = \frac{\eta}{2}(q^2 - p^2), \qquad (4.2)$$

or

$$H_C(p,q) = \frac{p^2 + \omega^2 q^2}{2}, \qquad H_I(p,q) = \frac{\eta}{2}(q^2 - p^2) + \epsilon_c \frac{\alpha q^4}{1 + 2\alpha q^2}. \qquad (4.3)$$

For this set of simulations, we use the following parameters: oscillator frequency $\omega = 1$, quadratic coupling $\eta = 0.5$, external field $B_0 = 0.1$, and inverse temperature $\beta = 2$. In addition, we set $\alpha = 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 \to \infty$.

In the first case, we consider a correction of the classical potential ($\epsilon = 0.1$ and $\epsilon_c = 0$), while in the second case we insert a correction in the coupling ($\epsilon = 0$ and $\epsilon_c = 0.05$). The results are plotted respectively in figures 8 and 9, for short times ($t \in [0, 20]$, top panels) and long times ($t \in [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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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 figure 10. For short times, the hybrid predictions are in excellent agreement with the fully quantum results. For long times, the behavior 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 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.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Uncertainty evolution (rational potential/coupling correction). Evolution of the uncertainty $\langle q^2 \rangle \langle p^2 \rangle$ 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.</figcaption></figure>

**5. 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 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 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 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–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 = ||\Upsilon||^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 quantum-classical 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.

**Data availability statement**

No new data were created or analyzed in this study.

**Acknowledgment**

This work arises from our discussions during the conference *Koopman Methods in Classical and Classical-Quantum Mechanics* held online on 18–23 April 2021, during the times of
COVID-19. We are grateful to the Organizing Committee and in particular to the Wilhelm und Else Heraeus Foundation at Bad Honnef, Germany, for the support provided and for hosting the meeting on their online platform. Also, we are grateful to Paul Bergold, Denys Bondar, Francesco Di Maiolo, Darryl Holm, and Raymond Kapral for their keen remarks on these and related results. G M and C T acknowledge financial support from the Royal Society Grant IES R3 203005. These results were also made possible through the support of Grant 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.

## Appendix A. Numerical solution of the phase-space equations

The transport equations (2.3) and (2.4) for the hybrid model both take the following form:

$$ \partial_t f + V(p) \partial_q f + U(q) \partial_p f - \mathcal{L}(q,p) = 0, \quad (A.1) $$

where $V(p) = \partial_p H$ and $U(q) = -\partial_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:

$$ \partial_t f + V(p) \partial_q f = 0, \quad \text{for } \Delta t/2 \quad (A.2) $$

$$ \partial_t f + U(q) \partial_p f = 0, \quad \text{for } \Delta t/2 \quad (A.3) $$

$$ \partial_t f - \mathcal{L}(q,p) = 0, \quad \text{for } \Delta t \quad (A.4) $$

$$ \partial_t f + U(q) \partial_p f = 0, \quad \text{for } \Delta t/2 \quad (A.5) $$

$$ \partial_t f + V(p) \partial_q f = 0, \quad \text{for } \Delta t/2. \quad (A.6) $$

Note that the above scheme is time-symmetric, and hence exact to second order in the time step $\Delta t$. Each of the above equations possesses an exact solution. For instance, for equation (A.2), one has:

$$ f(q,p,t + \Delta t/2) = f(q - V(p)\Delta t/2, p, t), \quad (A.7) $$

i.e. a shift of $V(p)\Delta t/2$ in the $q$ direction. Similarly, for equation (A.3) one gets a shift of $U(q)\Delta t/2$ in the $p$ direction. The solution of equation (A.4) is trivial.

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] \times [-25, 25]$ (position space $\times$ momentum space), with $N_q = 1000$ points in position space and $N_p = 1000$ points in momentum space. The time step was typically $\Delta 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 conditions.
The main limitation of grid-based codes applied to Vlasov-like equations such as equation (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.

## Appendix 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 \partial_t \psi_{\pm} = H_{\pm} \psi_{\pm} \qquad (B.1)$$

for the spinor $\psi = (\psi_+, \psi_-)^T$. For the case of quadratic coupling considered in the main text, the Hamiltonian is:

$$H_{\pm} = \frac{\omega^2(t)q^2 + p^2}{2} \pm \frac{\eta(t)}{2} (q^2 - p^2) \pm B_0(t). \qquad (B.2)$$

It can be shown that an exact solution to the time-dependent Schrödinger equation takes the Gaussian form

$$\psi_{\pm}(q,t) = A_{\pm} \left(\frac{2b_{\pm}}{\pi}\right)^{\frac{1}{4}} \exp\left(-b_{\pm}(q-d_{\pm})^2 + ia_{\pm}(q-d_{\pm})^2 + ie_{\pm}(q-d_{\pm}) + ic_{\pm}\right), \quad (B.3)$$

provided the initial condition is also in this form. In the above expression, the auxiliary functions $a_{\pm}(t)$, $b_{\pm}(t)$, $c_{\pm}(t)$, $d_{\pm}(t)$ and $e_{\pm}(t)$ are real functions that depend on the time $t$, and $A_{\pm}$ 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:

$$\left\{
\begin{aligned}
\dot{a}_{\pm} &= -\frac{\omega^2(t) \pm \eta(t)}{2} - 2(a_{\pm}^2 - b_{\pm}^2)[1 \mp \eta(t)], \\
\dot{b}_{\pm} &= -4a_{\pm}b_{\pm}[1 \mp \eta(t)], \\
\dot{c}_{\pm} &= -\frac{\omega^2(t) \pm \eta(t)}{2}d_{\pm}^2 + \left(\frac{e_{\pm}^2}{2} - b_{\pm}\right)[1 \mp \eta(t)] \mp B_0(t), \\
\dot{d}_{\pm} &= e_{\pm}[1 \mp \eta(t)], \\
\dot{e}_{\pm} &= -d_{\pm}[\omega^2(t) \pm \eta(t)],
\end{aligned}
\right.
\qquad (B.4)$$

where a dot stands for differentiation with respect to time. As initial condition at $t=0$, we take: $b_{\pm}(0) = 1/2$ (minimum uncertainty wave packet), $d_{\pm}(0) = q_0$ and $e_{\pm}(0) = 0$ (Gaussian wave packet centered at $q_0$), and $a_{\pm}(0) = c_{\pm}(0) = 0$ (vanishing initial phase).

The equations (B.4) form a system of nonlinear ordinary differential equations, 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:

$$\langle q \rangle = d_+ |A_+|^2 + d_- |A_-|^2,$$

and the z component of the spin

$$\langle \sigma_z \rangle = \int |\psi_+|^2 \, dq - \int |\psi_-|^2 \, dq = |A_+|^2 - |A_-|^2,$$
which also shows that $\langle \sigma_z \rangle$ 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 second-order accurate in time and space.

## ORCID iDs

Giovanni Manfredi https://orcid.org/0000-0002-5214-8707

Cesare Tronci https://orcid.org/0000-0002-8868-8027

## References

[1] Abedi A, Maitra N T and Gross E K U 2012 Correlated electron-nuclear dynamics: exact factorization of the molecular wavefunction *J. Chem. Phys.* **137** 22A530

[2] Agostini F, Caprara S and Ciccotti G 2007 Do we have a consistent non-adiabatic quantum–classical mechanics? *Eur. Phys. Lett.* **78** 30001

[3] Akimov A V, Long R and Prezhdo O V 2014 Coherence penalty functional: a simple method for adding decoherence in Ehrenfest dynamics *J. Chem. Phys.* **140** 194107

[4] Albert V V, Scholes G D and Brumer P 2011 Symmetric rotating-wave approximation for the generalized single-mode spin-boson system *Phys. Rev. A* **84** 042110

[5] Aleksandrov I V 1981 The statistical dynamics of a system consisting of a classical and a quantum subsystem *Z. Naturforsch. A* **36** 902–8

[6] Alonso J L, Clemente-Gallardo J, Cuchí J C, Echenique P and Falceto F 2012 Ehrenfest dynamics is purity non-preserving: a necessary ingredient for decoherence *J. Chem. Phys.* **137** 054106

[7] Bermúdez Manjarres A D 2021 Projective representation of the Galilei group for classical and quantum–classical systems *J. Phys. A: Math. Theor.* **54** 444001

[8] Blais A, Huang R-S, Wallraff A, Girvin S M and Schoelkopf R J 2004 Cavity quantum electrodynamics for superconducting electrical circuits: an architecture for quantum computation *Phys. Rev. A* **69** 062320

[9] Bondar D I, Gay-Balmaz F and Tronci C 2019 Koopman wavefunctions and classical-quantum correlation dynamics *Proc. R. Soc. A* **475** 20180879

[10] Boucher W and Traschen J 1988 Semiclassical physics and quantum fluctuations *Phys. Rev. D* **37** 3522–32

[11] Breuer H-P and Petruccione F 2002 *The Theory of Open Quantum Systems* (Oxford: Oxford University Press)

[12] Carlip S 2008 Is quantum gravity necessary? *Class. Quantum Grav.* **25** 154010

[13] Cariñena J F, Rañada M F and Santander M 2007 A quantum exactly solvable non-linear oscillator with quasi-harmonic behaviour *Ann. Phys.* **322** 434–59

[14] Chruściński D, Kossakowski A, Marmo G and Sudarshan E C G 2011 Dynamics of interacting classical and quantum systems *Open. Syst. Inf. Dyn.* **18** 339–51

[15] Crespo-Otero R and Barbatti M 2018 Recent advances and perspectives on nonadiabatic mixed quantum-classical dynamics *Chem. Rev.* **118** 7026–68

[16] de Gosson M A 2004 On the notion of phase in mechanics *J. Phys. A: Math. Gen.* **37** 7297–314

[17] Diósi L 2014 Hybrid quantum-classical master equations *Phys. Scr. T163* 014004

[18] Emary C and Bishop R F 2002 Exact isolated solutions for the two-photon Rabi Hamiltonian *J. Phys. A* **35** 8231–41

[19] Filbet F, Sonnendrücker E and Bertrand P 2001 Conservative numerical schemes for the Vlasov equation *J. Comput. Phys.* **172** 166

[20] Foskett M S, Holm D D and Tronci C 2019 Geometry of nonadiabatic quantum hydrodynamics *Acta Appl. Math.* **162** 1–41

[21] Gay-Balmaz F and Tronci C 2023 Dynamics of mixed quantum-classical spin systems *J. Phys. A: Math. Gen.* accepted (https://doi.org/10.1088/1751-8121/acc145)
[22] Gay-Balmaz F and Tronci C 2022 Evolution of hybrid quantum–classical wavefunctions *Physica D* **440** 133450

[23] Gay-Balmaz F and Tronci C 2022 Koopman wavefunctions and classical states in hybrid quantum-classical dynamics *J. Geom. Mech.* **14** 559–96

[24] Gay-Balmaz F and Tronci C 2019 Madelung transform and probability densities in hybrid quantum–classical dynamics *Nonlinearity* **33** 5383–424

[25] Gerasimenko V 1982 Dynamical equations of quantum–classical systems *Theor. Math. Phys.* **50** 49–55

[26] Ghose P 2015 The unfinished search for wave-particle and classical-quantum harmony *J. Adv. Phys.* **4** 236–51

[27] Gu B and Franco I 2018 Electronic interactions do not affect electronic decoherence in the pure-dephasing limit *J. Chem. Phys.* **149** 174115

[28] Hall M J W and Reginatto M 2016 *Ensembles on Configuration Space* (Berlin: Springer)

[29] Hammes-Schiffer S 2022 Theoretical perspectives on non-Born–Oppenheimer effects in chemistry *Phil. Trans. R. Soc. A* **380** 20200377

[30] Hu W, Gu B and Franco I 2018 Lessons on electronic decoherence in molecules from exact modeling *J. Chem. Phys.* **148** 134304

[31] Huang J-F, Liao J-Q, Tian L and Kuang L-M 2017 Manipulating counter-rotating interactions in the quantum Rabi model via modulation of the transition frequency of the two-level system *Phys. Rev. A* **96** 043849

[32] Hurst J, Hervieux P-A and Manfredi G 2017 Phase-space methods for the spin dynamics in condensed matter systems *Phil. Trans. R. Soc. A* **375** 20160199

[33] Jauslin H R and Sugny D 2010 Dynamics of mixed quantum–classical systems, geometric quantization and coherent states *Mathematical Horizons for Quantum Physics* (Lecture Notes Series, Institute for Mathematical Sciences, National University of Singapore vol 20) (Singapore: World Scientific) pp 65–96

[34] Joseph I 2020 Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics *Phys. Rev. Res.* **2** 043102

[35] Kapral R 2006 Progress in the theory of mixed quantum–classical dynamics *Annu. Rev. Phys. Chem.* **57** 129–57

[36] Koopman B O 1931 Hamiltonian systems and transformations in Hilbert space *Proc. Natl Acad. Sci.* **17** 315

[37] Larson J and Mavrogordatos T 2021 *The Jaynes-Cummings Model and its Descendants* (Bristol: IOP Publishing)

[38] Manfredi G and Feix M R 1996 Theory and simulation of classical and quantum echoes *Phys. Rev. E* **53** 6460

[39] Mitra A K 1978 On the interaction of the type $\lambda x^2/(1+gx^2)$ *J. Math. Phys.* **19** 2018

[40] Ng K M, Lo C F and Liu K L 1999 Exact eigenstates of the two-photon Jaynes-Cummings model with the counter-rotating term *Eur. Phys. J. D* **6** 119–26

[41] O’Brien M C M and Chancey C C 1993 The Jahn-Teller effect: an introduction and current review *Am. J. Phys.* **61** 688–97

[42] Peres A 2000 Classical interventions in quantum systems. I. The measuring process *Phys. Rev. A* **61** 022116

[43] Peres A and Terno D R 2001 Hybrid quantum–classical dynamics *Phys. Rev. A* **63** 022101

[44] Prezhdo O V and Kisil V V 1997 Mixing quantum and classical mechanics *Phys. Rev. A* **56** 162–75

[45] Reichman D, Silbey R J and Suárez A 1996 On the nonperturbative theory of pure dephasing in condensed phases at low temperatures *J. Chem. Phys.* **105** 10500–6

[46] Roy A K, Jalbout A F and Proynov E I 2008 Bound state spectra of the 3D rational potential *Int. J. Quant. Chem.* **108** 827–36

[47] Salcedo L L 1996 Absence of classical and quantum mixing *Phys. Rev. A* **54** 3657–60

[48] Skinner J L 1988 Theory of pure dephasing in crystals *Ann. Rev. Phys. Chem.* **39** 463–78

[49] Skinner J L and Hsu D 1986 Pure dephasing of a two-level system *J. Phys. Chem.* **90** 4931–8

[50] Sudarshan E C G 1976 Interaction between classical and quantum systems and the measurement of quantum observables *Prāmana* **6** 117–26

[51] Rustagi A, Kajale S and Upadhyaya P 2022 Manipulating quantum impurity spins via dynamical modes of nanomagnets (arXiv:2203.03652)

[52] Terno D R 2006 Inconsistency of quantum–classical dynamics and what it implies *Found. Phys.* **36** 102–11
[53] Travenec I 2012 Solvability of the two-photon Rabi Hamiltonian Phys. Rev. A 85 043805

[54] Tronci C and Joseph I 2021 Koopman wavefunctions and Clebsch variables in Vlasov-Maxwell kinetic theory J. Plasma Phys. 87 835870402

[55] Tully J C 1998 Mixed quantum-classical dynamics Faraday Discuss. 110 407–19

[56] von Neumann J 1932 Zur Operatorenmethode in der klassischen Mechanik Ann. Math. 33 587–642

[57] Widom A and Srivastava Y 1992 Lagrangian formulation of Bohr's measurement theory Nuovo Cimento B 107 71–75

[58] Xie Q, Zhong H, Batchelor M T and Lee C 2017 The quantum Rabi model: solution and dynamics J. Phys. A: Math. Theor. 50 113001