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

HAL
open science

# Evolution of hybrid quantum-classical wavefunctions

François Gay-Balmaz, Cesare Tronci

► To cite this version:

François Gay-Balmaz, Cesare Tronci. Evolution of hybrid quantum-classical wavefunctions. Physica D: Nonlinear Phenomena, 2022, 440, 10.1016/j.physd.2022.133450. insu-03847077

HAL Id: insu-03847077

https://insu.hal.science/insu-03847077

Submitted on 10 Nov 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Evolution of hybrid quantum-classical wavefunctions

François Gay-Balmaz<sup>a,1</sup>, Cesare Tronci<sup>b,c,*1</sup>

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

## ARTICLE INFO

**Article history:**

Received 8 April 2022

Received in revised form 5 July 2022

Accepted 7 July 2022

Available online 16 July 2022

Communicated by Liang Huang

**Keywords:**

Koopman wavefunctions

Exact factorization

Mixed quantum-classical dynamics

Hamiltonian systems

Gauge theories

Geometric mechanics

## ABSTRACT

A gauge-invariant wave equation for the dynamics of hybrid quantum-classical systems is formulated by combining the variational setting of Lagrangian paths in continuum theories with Koopman wavefunctions in classical mechanics. We identify gauge transformations with unobservable phase factors in the classical phase-space and we introduce gauge invariance in the variational principle underlying a hybrid wave equation previously proposed by the authors. While the original construction ensures a positive-definite quantum density matrix, the present model also guarantees the same property for the classical Liouville density. After a suitable wavefunction factorization, gauge invariance is achieved by resorting to the classical Lagrangian paths made available by the Madelung transform of Koopman wavefunctions. Due to the appearance of a phase-space analogue of the Berry connection, the new hybrid wave equation is highly nonlinear and it is proposed here as a platform for further developments in quantum-classical dynamics. Indeed, the associated model is Hamiltonian and appears to be the first to ensure a series of consistency properties beyond positivity of quantum and classical densities. For example, the model possesses a quantum-classical Poincaré integral invariant and its special cases include both the mean-field model and the Ehrenfest model from chemical physics.

© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The coupling between classical and quantum degrees of freedoms represents an outstanding question ever since the rise of quantum mechanics as a physical theory. Following Bohr's work, Landau and Lifshitz [1] defined a quantum measurement as the interaction between a measured quantum system and a classical system commonly referred to as the apparatus. In this way, "quantum mechanics occupies a very unusual place among physical theories: it contains classical mechanics as a limiting case, yet at the same time it requires this limiting case for its own formulation" [1]. However, the search for an extension of quantum theory (and its mathematical footing) to include the coupling to a classical system has been one of the most challenging directions for about a century. The main difficulty is capturing the correlations between the quantum system and the 'apparatus' that are responsible for the decoherence phenomenon, i.e. the loss of pure states in the quantum sector as time goes by. In addition, according to Bohr [2], the classical apparatus itself undergoes some kind of uncontrollable disturbance that is also a consequence of quantum-classical correlations. While Bohr's disturbance has never been fully characterized, a possible interpretation lies in the loss of classical pure states, i.e. the loss of delta-like classical Liouville distributions [3].

The search for an extension of quantum theory to incorporate the interaction with classical systems leads naturally to the idea of a hybrid model in which classical and quantum dynamics are treated on an equal footing: when no interaction is present in the treatment, the theory must reduce to the uncoupled equations of quantum and classical mechanics, that is Schrödinger's (or quantum Liouville) and Newton's (or classical Liouville) equations, respectively. For example, this important property is satisfied by the mean-field model. In standard Poisson bracket and commutator notation, this model is given by

$$ \frac{\partial \rho_c}{\partial t} = \{\text{Tr}(\hat{\rho}_q \hat{H}), \rho_c\}, \qquad i\hbar \frac{\partial \hat{\rho}_q}{\partial t} = \left[ \int \rho_c \hat{H} dqdp, \hat{\rho}_q \right]. \qquad (1.1) $$

<footnote><sup>a</sup> CNRS and École Normale Supérieure, Laboratoire de Météorologie Dynamique, Paris, France</footnote>

<footnote><sup>b</sup> Department of Mathematics, University of Surrey, Guildford, United Kingdom</footnote>

<footnote><sup>c</sup> Department of Physics and Engineering Physics, Tulane University, New Orleans LA, United States</footnote>

<footnote>* Corresponding author at: Department of Mathematics, University of Surrey, Guildford, United Kingdom.</footnote>

<footnote>E-mail address: c.tronci@surrey.ac.uk (C. Tronci).</footnote>

<footnote n="1">1 Both authors contributed equally to the development of this work.</footnote>
Here, $\rho_c$ and $\hat{\rho}_q$ are the classical Liouville density and the quantum von Neumann operator, while the hybrid Hamiltonian $\hat{H}(q, p)$ is a phase-space function taking values in the quantum operators. Unfortunately, the mean-field equations (1.1) neglect quantum-classical correlations [4], thereby ignoring decoherence. Thus, one is motivated to look for a more involved model where the fundamental quantities are functions of both classical and quantum coordinates on phase-space and configuration-space, respectively. This direction was pursued only occasionally in the physics community [5–8] and with little success. For example, in some cases the models fail to reduce to the uncoupled equations of quantum and classical mechanics in the absence of quantum-classical interaction [9]. Lists of consistency criteria for hybrid quantum-classical models were provided in [5] and more recently in [10]. Among the most important criteria, we find the positivity of quantum and classical densities, as well as the decoherence property $\text{Tr}\hat{\rho}_q^2 \neq \text{const.}$

The difficulties encountered in establishing a hybrid quantum-classical theory beyond mean-field led much of the community to move away from this direction, thereby abandoning Bohr's initial view of quantum measurements based on quantum-classical interaction. For example, no-go arguments have been proposed over the years to even exclude the possibility of such hybrid theories [11–13] beyond simple models of mean-field type. A currently common picture involves a fully quantum apparatus in thermodynamic equilibrium, while the link to classical mechanics is provided by the classical limit that is obtained as a result of long-time evolution. Based on Lindblad's master equation [14] this approach completely ignores the motion of the apparatus, whose degrees of freedom are traced out to yield the irreversible dynamics of the quantum density matrix. Due to irreversibility, the fine small-scale structures are eventually smoothened out just as in the classical Fokker-Planck equation. Then, in this picture a classical state emerges as the result of irreversible dissipative dynamics at the quantum scales. Due to the absence of small-scale phenomenology, it has been pointed out that the classical states made available in this way are "arguably among the least classical ones" [15], as it seems to emerge also from recent simulations in the Wigner picture [16]. However, we shall not enter this discussion here.

While quantum-classical coupling first emerges as a foundational question in relation to the measurement process, similar question appear in semiclassical descriptions of gravity [17]. In this context, the absence of a consistent quantum theory of gravity leads to treating the geometric background as a classical system [18]. Then, the question is how an evolving quantum field in an expanding universe affects the expansion rate of the background on which it is defined. Currently available semiclassical theories fail to address this backreaction problem and thus one is motivated to look for more complete models that can include this important effect [5].

Besides the foundational aspects, the development of mixed quantum-classical models is also an active area of research in more applied areas such as chemical physics [19–21], solid state physics [22,23], and more recently spintronics [24]. For example, in the latter case, recent proposals involve the use of classical ferromagnets to enlarge the parameter space in spin control. On the other hand, in quantum chemistry and solid-state physics, hybrid/mixed quantum-classical models emerge as a promising direction to alleviate the curse of dimensionality currently challenging many-particle quantum simulations. One typically identifies the slow motion of some degrees of freedom which are eventually treated as classical in order to mitigate the cost of simulating their evolution. The most celebrated example is provided by Born–Oppenheimer molecular dynamics [25], which arises from methods in slow-manifold reduction [26]. In this model, the nuclei are treated as classical, while quantum electrons are constrained to evolve adiabatically in their ground state. Then, the degeneracies appearing in the associated electronic eigenvalue problem lead to topological singularities, known as “conical intersections”. The numerical treatment of these singularities and their resulting geometric phase pose further difficulties that continue to challenge the community. However, recently the question emerged whether the topological singular nature should be retained at all in order to account for geometric phase effects [27,28].

Despite the paramount importance of Born–Oppenheimer molecular dynamics, the adiabatic approximation underlying electron dynamics is quite limitative and a whole new industry has now emerged to tackle nonadiabatic dynamics. In this context, the equations currently implemented in hybrid quantum-classical methods [19] often lack a sound mathematical footing. The desire for such a fundamental rigorous ground has led some [29,30] to derive mixed quantum-classical algorithms from the Aleksandrov–Gerasimenko (AG) equation [31,32]

$$ \frac{\partial \hat{D}}{\partial t} = -i\hbar^{-1}[\hat{H}, \hat{D}] + \frac{1}{2} (\{\hat{H}, \hat{D}\} - \{\hat{D}, \hat{H}\}). \qquad (1.2) $$

Here, $\hat{D}(q, p)$ is a phase-space density with values in the quantum von Neumann operators, so that the classical and quantum densities are written as $\rho_c = \text{Tr}\hat{D}$ and $\hat{\rho}_q = \int \hat{D} dqdp$. Eq. (1.2) satisfies several important properties [5] and is commonly regarded as a fundamental step beyond the mean-field model. In turn, the AG equation allows for a sign-indefinite density matrix $\hat{\rho}_q$ of the quantum subsystem with the potential drawback of violating the uncertainty principle. In addition, as the AG equation breaks the time-reversal symmetry, one would like to justify this property in terms of an $H$–theorem for the entropy increase. However, such a result is currently unavailable. This lack of entropy arguments also applies to stochastic augmentations of the AG equation [33].

Recently, a new approach has been proposed by the authors [3] upon following a suggestion by George Sudarshan [7,34]. By exploiting Koopman's Hilbert space formulation of classical mechanics [35], one can try to construct a hybrid quantum-classical wavefunction $\Upsilon(q, p, x)$, whose dynamics recovers the mean-field model by writing $\Upsilon(q, p, x) = \chi(q, p)\psi(x)$. Here, $(q, p)$ are classical phase-space coordinates and $x$ is the quantum configuration coordinate. Moreover, $\chi$ is a Koopman wavefunction obeying an equation of the type

$$ i\hbar\partial_t\chi = \{i\hbar H, \chi\} + \varphi\chi, \qquad (1.3) $$

so that the classical Liouville equation $\partial_t\rho_c = \{H, \rho_c\}$ emerges by writing $\rho_c(q, p) = |\chi(q, p)|^2$. However, finding a consistent evolution law for a general hybrid wavefunction $\Upsilon(q, p, x)$ is far from easy and Sudarshan's early attempt led to interpretative issues [13,36,37]. The authors recently addressed this problem by modifying Koopman's original prescription $\rho_c = |\chi|^2$ after selecting the specific phase factor $\varphi = p\partial_pH - H$ given by the usual Lagrangian on phase-space. The resulting hybrid theory leads to the construction of an unsigned operator-valued distribution $\hat{D}(q, p)$, which however yields a positive-definite quantum density matrix $\hat{\rho}_q$. As this feature allows to recover the quantum uncertainty principle, the underlying hybrid model from [3] represents a step forward beyond the AG equation (1.2). However, the classical density $\rho_c$ constructed in this way has a sign-indefinite expression and thus one is led to ask whether an initially positive Liouville density may develop negative values in time. As discussed in [38], a general answer is currently unavailable, although there are infinite families of hybrid systems for which $\rho_c$ was shown to stay positive in time.
The question whether the classical density should be positive-definite in hybrid dynamics is crucial since the development of negative values may lead to interpretative issues. Some of these issues were addressed in [3] by resorting to analogies with Wigner functions, following earlier ideas by Feynman [39]. However, general arguments are still lacking and thus one may ask if the hybrid theory presented in [3] can be modified to accommodate a positive-definite classical Liouville density. A positive answer was provided by the authors in [10] by combining recent mathematical techniques in chemical physics and the geometry of Hamilton's variational principle. The main feature is that, while the original hybrid theory is linear, the new closure variant appears to be nonlinear and the nonlinearity seems to be produced by the emergence of a Berry connection in phase-space. Explicitly, this closure module reads as follows:

$$i\hbar\partial_t\hat{P} + i\hbar\operatorname{div}(\hat{P}(\mathbf{X}_{\hat{H}})) = [\hat{H}, \hat{P}], \qquad (1.4)$$

where $\hat{P}(q, p)$ is a hybrid density playing the role of $\hat{D}$ in (1.2), while we have used the notation $\mathbf{X}_{\hat{A}} = (\partial_p\hat{A}, -\partial_q\hat{A})$ and $(\hat{A}) = \operatorname{Tr}(\hat{A}\hat{P})/\operatorname{Tr}\hat{P}$. Here, the effective Hamiltonian $\hat{H}$ is given by $\hat{H} = \hat{H} + \hbar\hat{\mathcal{F}}(\hat{P})$, where the operator-valued function $\hat{\mathcal{F}}$ is expressed in terms of $\hat{P}$ and $\{\hat{P}, \hat{H}\} + \{\hat{H}, \hat{P}\}$. Neglecting the $\hbar$-contribution to the effective Hamiltonian reduces Eq. (1.4) to the Ehrenfest model, which will be discussed at the end of this paper. In more generality, the full closure model (1.4) has so far been formulated only in terms of hybrid quantum-classical densities. In this paper, we will show that a wavefunction description unfolds several relevant properties, including the existence of a Poincaré integral invariant leading to the construction of a Liouville volume density that can be used to identify a Gibbs entropy functional for the classical sector.

In particular, this paper shows that positivity of the classical density may be achieved by adopting a gauge invariance principle to reflect the fact that classical phases are actually unobservable. This particular remark about the irrelevance of local phases in classical mechanics represents the second fundamental point raised by Sudarshan and it was addressed in his paper by the enforcement of specific superselection rules whose mode of operation has never been clear [13,36,37]. Instead of superselection rules, this paper addresses the point of unobservable classical phases by treating them as a gauge freedom. This is exactly what happens in standard quantum mechanics where global phase factors are irrelevant and wavefunctions are always defined up to a phase. In this context, the Liouville density can be thought of as a Noether charge just like total probability in quantum theory. This paper pursues this analogy by making the original treatment in [3] manifestly gauge invariant under local phase transformations, that is multiplicative phase factors for which the phase is a function on the classical phase-space. The achievement of gauge invariance is made possible by the specific fact that the original quantum-classical wave equation in [3] is already manifestly gauge-covariant. On the other hand, the actual process leading to gauge invariance requires several steps. Starting with the variational principle underlying the original theory, its associated Lagrangian will undergo several exact transformation steps so that it is finally taken into a form where the gauge invariance principle can finally be applied. As we will see, this is a nontrivial process combining hydrodynamic approaches to Liouville phase-space dynamics with specific wavefunction factorization techniques recently developed in chemical physics. The combination of these approaches is made necessary by two facts: (1) the classical Liouville equation is a characteristic equation accompanied by a Lagrangian-path hydrodynamic description, and (2) the extraction of the classical phase in a hybrid context requires writing the quantum-classical wavefunction in a convenient way so that a classical phase can be adequately identified. Fortunately, the process leading to this gauge-invariance principle has a purely classical analogue, which will serve as a pilot example. As a general methodology, we emphasize that this paper will not proceed by operating on a particular set of equations of motion. Instead, we will mostly manipulate its underlying variational principle in order to keep track of basic conservation laws and transformation properties. This is a common method in gauge theory and has a long history in quantum dynamics [40]. For an account of quantum-classical variational principles in the context of the Ehrenfest theorem, see [41].

**Plan of the paper.** Our discussion starts with a discussion of Koopman wavefunctions in classical mechanics. After a preliminary introduction to standard Koopman-von Neumann (KvN) theory in Section 2.1, we show how apparent issues in the original treatment may be addressed by adding a phase term first proposed by van Hove. This is presented in Section 2.2, which covers the Koopman-van Hove (KvH) formulation of classical mechanics. The important role of phase transformations in the Koopman context is discussed in Section 2.3, where the concept of gauge transformation in classical mechanics is also introduced. As the KvH construction appears to differ substantially from KvN theory, one is led to ask whether these two formulations are actually related. As shown in Section 3, standard KvN dynamics may be obtained by simply applying a gauge-invariance principle to the variational formulation of KvH theory. As shown in Section 3.1, the application of this gauge-invariance principle requires techniques from continuum theories, such as the Lagrangian flow paths widely used in hydrodynamics. Then, the gauge-invariance principle leads to an alternative variational formulation of KvN theory, which represents the immediate extension of the well-known variational setting of the classical Liouville equation [42].

After discussing the various Koopman formulations of classical mechanics, the paper proceeds in Section 4 to illustrate the present status of the Koopman hybrid theory for the description of mixed quantum-classical systems [3]. Once hybrid quantum-classical wavefunctions are introduced in Section 4.1, Section 4.2 discusses how classical and quantum densities may be obtained in this context. The rest of the paper proceeds to apply the same gauge-invariance principle from Koopman classical mechanics in the hybrid setting. This requires several steps, the first of which consists in a suitable factorization of the hybrid wavefunction, as discussed in Section 5. The latter is entirely devoted to transforming the variational principle of the original hybrid formulation from the preceding sections in such a way that the gauge principle can finally be applied. After applying the Madelung transform in Section 5.1, once again we observe the emergence of Lagrangian flow paths in the Koopman context, as discussed in Section 5.2.

At this point, the variational principle of the original hybrid formulation is amenable to the application of the gauge principle in such a way to make classical phases unobservable. This is the topic of Section 6, which starts by constructing the gauge-invariant Lagrangian of quantum-classical dynamics. The latter is formulated in Section 6.1 while Section 6.2 presents the full gauge-invariant evolution of hybrid wavefunctions.

A discussion of some of the various properties of the new gauge-invariant model is presented in Section 7. After discussing the presence of a Poincaré integral invariant and some of its consequences in Section 7.1, the noncanonical Hamiltonian structure of hybrid dynamics is presented in Section 7.2. Also, Section 7.3 presents a prototype for a quantum-classical density operator: while this appears
to be generally unsigned, it gives rise to positive quantum and classical densities and possesses important covariance properties finally ensuring basic conservation laws such as momentum conservation. Finally, Section 7.4 presents relevant special cases, including the Ehrenfest model from chemical physics. As discussed therein, the present model allows to identify explicitly the source of backreaction forces.

## 2. Koopman wavefunctions in classical mechanics

This section reviews different formulations of classical mechanics based on Koopman wavefunctions. As we will see, the standard Koopman-von Neumann formulation presents important issues which will be addressed by exploiting the variational structure of Koopman dynamics.

### 2.1. Koopman-von Neumann construction

The simplest example of Koopman classical dynamics is given by the Koopman-von Neumann (KvN) equation in the following common form:

$$i\hbar\partial_t\chi = \{i\hbar H, \chi\} =: \hat{L}_H\chi. \qquad (2.1)$$

Here, the factor $i\hbar$ is inserted to produce the Hermitian operator $\hat{L}_H = \{i\hbar H, \}$, known as Liouvillian. In this way, classical dynamics is regarded as a unitary flow on the Hilbert space $\mathcal{H}_C = L^2(\mathbb{R}^2)$ of square-integrable functions on phase-space. Here, we focus on one-dimensional degrees of freedom as the extension to higher dimensions is straightforward. In polar form $\chi = \sqrt{D}e^{iS/\hbar}$, one finds

$$\partial_t D = \{H, D\}, \qquad \partial_t S = \{H, S\}, \qquad (2.2)$$

thereby recovering the usual Liouville equation for the density $D = |\chi|^2$. We observe that the phase can be simply set to zero so that the KvN wavefunction may be identified with a real-valued amplitude $|\chi|$.

Since the KvN equation (2.1) identifies classical dynamics as a type of unitary flow, one may look at its underlying variational principle. Mimicking the Dirac-Frenkel Lagrangian for the standard Schrödinger equation [43], we write [44,45]

$$\delta \int_{t_1}^{t_2} L_{KvN}(\chi, \partial_t \chi) dt = 0, \qquad L_{KvN}(\chi, \partial_t \chi) = \text{Re} \int i\hbar \chi^*(\partial_t \chi - \{H, \chi\}) dqdp. \qquad (2.3)$$

In [3], we showed how this variational principle for KvN dynamics may lead to relevant issues. For example, integration by parts leads to rewriting the conserved Hamiltonian functional $h_{KvN}(\chi) = \text{Re} \int i\hbar\chi^*\{H, \chi\} dqdp$ as

$$h_{KvN}(\chi) = \hbar \int H \operatorname{Im}\{\chi^*, \chi\} dqdp, \qquad (2.4)$$

which differs from the physical energy $\int H|\chi|^2 dqdp$, unless one sets $|\chi|^2 = \hbar \operatorname{Im}\{\chi^*, \chi\}$. Indeed, both quantities $|\chi|^2$ and $\hbar \operatorname{Im}\{\chi^*, \chi\}$ satisfy the classical Liouville equation, although they can be set equal only by allowing for topological singularities in the KvN wavefunction [45]. More importantly, if one initializes a zero phase $S = 0$ in (2.2), then the Hamiltonian functional (2.4) collapses entirely, thereby indicating that this variational formulation may need special care. This difficulty was recently overcome in [42], where an alternative variational formulation was presented. This variant will be discussed later on. Instead, the next section shows how this and other apparent issues may be addressed by modifying the original KvN formulation.

### 2.2. Koopman-van Hove construction

One possibility to overcome the difficulties arising from the conventional variational formulation of KvN theory was presented in [3] and follows from previous work by van Hove [46] and Kostant [47,48] on prequantum theory in geometric quantization [49,50]. Without going into the details, the immediate modification consists in changing the phase dynamics of KvN theory by replacing (2.2) by

$$\partial_t D = \{H, D\}, \qquad \partial_t S + \{S, H\} = p\partial_p H - H =: \mathcal{L}, \qquad (2.5)$$

where $\mathcal{L}$ identifies the usual phase-space expression of the Lagrangian which already appeared in the Introduction. Notice that the same equations may be written in the characteristic form

$$\frac{d}{dt} D(q, p, t) = 0, \qquad \frac{d}{dt} S(q, p, t) = \mathcal{L}(q, p), \qquad \text{along} \qquad (\dot{q}, \dot{p}) = (\partial_p H, -\partial_q H) =: \mathbf{X}_H(q, p),$$

where we have used standard notation for the Hamiltonian vector field $\mathbf{X}_H(q, p)$. We notice that the phase equation is reminiscent of the usual phase dynamics occurring in Hamilton-Jacobi theory; see [51] for more details on this point.

Eqs. (2.5) are naturally written in terms of the wavefunction $\chi = \sqrt{D}e^{iS/\hbar}$, thereby leading to the Koopman-van Hove (KvH) equation

$$i\hbar\partial_t\chi = i\hbar\{H, \chi\} - \mathcal{L}\chi =: \hat{\mathcal{L}}_H\chi. \qquad (2.6)$$

Following [3], here we have introduced the covariant Liouvillian operator $\hat{\mathcal{L}}_H := i\hbar\{H, \} - \mathcal{L}$, which is also known in prequantum theory as prequantum operator [47,50]. This Hermitian operator replaces the standard Liouvillian $\hat{\mathcal{L}}_H := i\hbar\{H, \}$ appearing in the KvN equation (2.1). The sense in which $\hat{\mathcal{L}}_H$ is covariant was explained in [3] in terms of minimal coupling arguments in phase-space. See Section 2.3 for further comments on gauge covariance in KvH theory. As Eq. (2.6) again identifies a unitary flow on $\mathcal{H}_C = L^2(\mathbb{R}^2)$ analogous to the Schrödinger equation, one can again look at the underlying variational principle of Dirac-Frenkel type, that is

$$\delta \int_{t_1}^{t_2} L_{KvH}(\chi, \partial_t \chi) dt = 0, \qquad L_{KvH}(\chi, \partial_t \chi) = \text{Re} \int \chi^*(i\hbar\partial_t\chi - i\hbar\{H, \chi\} + \chi\mathcal{L}) dqdp. \qquad (2.7)$$
We observe that, upon expanding $\mathcal{L} = p\partial_p H - H$ and integrating by parts, this variational principle identifies the conserved Hamiltonian functional

$$h_{KvH}(\chi) = \int H(|\chi|^2 + \partial_p(p|\chi|^2) + \hbar \operatorname{Im}\{\chi^*, \chi\}) dqdp. \qquad (2.8)$$

As pointed out in [3], one notes that in the parenthesis both the first term and the sum of the last two obey the classical Liouville equation. This observation leads to defining the Liouville density in KvH theory as

$$\begin{aligned} \rho_c &= |\chi|^2 + \partial_p(p|\chi|^2) + \hbar \operatorname{Im}\{\chi^*, \chi\} \\ &= D + \partial_p(pD) + \{D, S\}, \end{aligned} \qquad (2.9)$$

which eliminates all the ambiguities otherwise occurring in the KvN construction. For example, as shown in (2.6), the KvH phase undergoes nontrivial dynamics and therefore the Hamiltonian functional (2.8) cannot be made to vanish by simply setting an initial zero phase. In addition, we notice that the distribution (2.9) appearing in (2.8) integrates to 1, unlike the expression $\hbar \operatorname{Im}\{\chi^*, \chi\}$ appearing in (2.4).

While the present discussion omits the mathematical details, we point out that the above expression has a deep geometric structure in terms of momentum maps, as discussed in [38]. Importantly, here we notice that the expression (2.9) of the classical density is not positive-definite. This point should not be regarded as an issue in classical dynamics: since the Liouville equation is a characteristic equation, the sign of the initial condition is preserved in time. For example, following [3], one can initialize the KvH wavefunction $\chi$ corresponding to a Gaussian density.

## 2.3. Phase transformations in Koopman mechanics

A distinctive difference between the KvN and KvH formalisms lies in the different transformation properties underlying the classical Liouville density $\rho_c$ under local $U(1)$ transformations. Indeed, while the KvN prescription $\rho_c = |\chi|^2$ is evidently invariant under $\chi \mapsto \chi e^{i\phi/\hbar}$, the KvH expression (2.9) involves a phase shift $S \mapsto S + \phi$. Nevertheless, as pointed out in [3], the KvH prescription still enjoys a covariance property. To see this, let us introduce the symplectic potential $\mathcal{A} := (p, 0)$ (equivalently, $\mathcal{A}(z) \cdot dz = pdq$) to rewrite (2.9) as

$$\rho_c = D + \operatorname{div}(DJ(\nabla S - \mathcal{A})), \quad \text{with} \quad J = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}, \qquad (2.10)$$

where $J = (\nabla \mathcal{A})^T - \nabla \mathcal{A}$ is the canonical symplectic form. Then, it is clear that combining the phase transformation $\chi \mapsto \chi e^{i\phi/\hbar}$ with the shift $\mathcal{A} \to \mathcal{A} - \nabla \phi$ leaves $\rho_c$ invariant, thereby recovering the usual notion of gauge-covariance. In this sense, the KvH construction is gauge-covariant. More specifically, we observe that the KvH Lagrangian in (2.7) is manifestly covariant in the sense that it satisfies the relation

$$L(e^{i\phi/\hbar}\chi, e^{i\phi/\hbar}\partial_t\chi; \mathcal{A} - \nabla\phi) = L(\chi, \partial_t\chi; \mathcal{A}), \qquad (2.11)$$

where the semicolon indicates that the vector potential $\mathcal{A}$ appears in the Lagrangian merely as a parameter. In more generality, we recall that a covariant Lagrangian allows for an additive total time derivative on the right-hand side of (2.11). Such an additive time derivative is also allowed in the case of gauge-invariant Lagrangians, while manifest gauge-invariance simply means $L(e^{i\phi/\hbar}\chi, e^{i\phi/\hbar}\partial_t\chi) = L(\chi, \partial_t\chi)$. We notice that, while the KvN expression of the Liouville density is indeed manifestly gauge-invariant, the KvN Lagrangian in (2.3) is not even covariant. This is a point which we will come back to in Section 3.2, where KvN theory will be endowed with a manifestly gauge-invariant Lagrangian.

Besides these technical details, we emphasize that this form of gauge transformations has been discussed in several prominent works in hybrid quantum-classical dynamics. In particular, here we refer to [5,52]. Therein, the authors point out that the gauge invariance under local $U(1)$ transformations is absolutely essential in reflecting the correct dynamics in the classical sector of hybrid systems. This is due to the fact that, as a general rule, phases are normally considered unobservable in classical mechanics. For example, as discussed in [52], Sudarshan's early approach to hybrid dynamics aimed at incorporating gauge invariance by enforcing appropriate superselection rules. However, the role of these superselection rules has attracted a certain amount of criticism [13,36,37].

In this work superselection rules are replaced by local phase invariance, much in the same spirit as global phase invariance in standard quantum mechanics. In particular, gauge invariance is achieved by starting with a covariant formulation based on KvH theory, so that the full invariance can be enforced at the level of the underlying variational principle. The next section illustrates this process for the special case of purely classical dynamics.

## 3. From KvH back to KvN

So far, we have observed that KvH theory overcomes the issues emerging in the standard variational formulation of KvN dynamics. However, one may insist that the KvN prescription $\rho = |\chi|^2$ is perfectly sensible so that the question becomes to identify an alternative variational formulation overcoming the issues emerging from the canonical structure for the KvN theory presented in Section 2.1. While this alternative construction was presented in [42], here we will show how this can be obtained as a closure model emerging from KvH theory. This new result will pave the way to the construction of a hybrid quantum-classical model later on.

We start our discussion by presenting a remark that motivates the steps taken later on. We realize that, if we could set $\partial_p S = 0$ and $p = \partial_q S$ in (2.10), the KvH expression of the classical density would recover the original KvN prescription $\rho_c = |\chi|^2$. Importantly, upon taking the gradient of the second in (2.5), we notice that $\nabla S - \mathcal{A}$ satisfies the evolution equation

$$\partial_t(\nabla S - \mathcal{A}) + \mathbf{X}_H \cdot \nabla(\nabla S - \mathcal{A}) + \nabla \mathbf{X}_H \cdot (\nabla S - \mathcal{A}) = 0 \qquad (3.1)$$

so that the initial condition $\nabla S = \mathcal{A}$ is preserved in time. However, this type of initial condition requires introducing topological singularities that are difficult to treat. Nevertheless, it is rather suggestive that replacing $\nabla S \to \mathcal{A}$ takes the classical density $\rho_c$ back to the KvN prescription $\rho_c = |\chi|^2$. In what follows we will take full advantage of this observation, yet avoiding topological singularities.
### 3.1. KvH variational principles and hydrodynamic paths

Instead of allowing for a singular phase, here we will pursue an alternative direction which consists in replacing $\nabla S \to \mathcal{A}$ within a suitable form of the KvH variational principle. Indeed, in order to make this replacement, we first need to bring the variational principle (2.7) into a form that is amenable to appropriate manipulations. For this purpose, we first express the Lagrangian $L_{KvH}$ in terms of phase and amplitude, that is

$$L_{KvH} = \int D(\partial_t S - \{H, S\} - p \partial_p H + H) dqdp. \qquad (3.2)$$

Inspired by Madelung's hydrodynamic method in quantum mechanics [53], we will now introduce the characteristics of the vector field $\mathbf{X}_H$. More specifically, borrowing the terminology from continuum dynamics, we introduce the Lagrangian coordinate $\eta(q_0, p_0, t)$ such that $\dot{\eta}(q_0, p_0, t) = \mathbf{X}_H(\eta(q_0, p_0, t))$. In this hydrodynamic picture, $\mathbf{X}_H$ is regarded as the Eulerian velocity corresponding to its Lagrangian counterpart $\dot{\eta}$. Then, as a consequence of the first equation in (2.5), the density $D$ is transported along $\eta$. Upon introducing the notation $\mathbf{z} = (q, p)$, the density transport can be written as the Lagrange-to-Euler map

$$D(\mathbf{z}, t) = \int D_0(\mathbf{z}_0) \delta(\mathbf{z} - \eta(\mathbf{z}_0, t)) d^2\mathbf{z}_0 = \frac{D_0(\mathbf{z}_0)}{\mathcal{J}_{\eta}(\mathbf{z}_0, t)} \Big|_{\mathbf{z}_0=\eta(\mathbf{z}, t)^{-1}}. \qquad (3.3)$$

Notice that here we are retaining the Jacobian $\mathcal{J}_{\eta} = \det \nabla \eta$ in the denominator for later convenience, although the flow of a Hamiltonian vector field identifies a canonical transformation with unit Jacobian.

Let us now replace (3.3) in the Lagrangian $L_{KvH}$. In particular, let us apply the product rule $D\partial_t S = \partial_t(DS) - S\partial_t D$ to the first term of (3.2) and ignore the irrelevant total derivative $\partial_t(DS)$. Here, the time derivative $\partial_t D$ may be found from (3.3) so that we have

$$\partial_t D + \mathrm{div}(\mathcal{X}D) = 0, \qquad (3.4)$$

with $\mathcal{X}(\mathbf{z}) = \dot{\eta}(\mathbf{z}_0, t)|_{\mathbf{z}_0=\eta^{-1}(\mathbf{z}, t)}$. Notice that, while the original relation $\dot{\eta} = \mathbf{X}_H(\eta)$ leads to $\mathbf{X}_H(\mathbf{z}) = \dot{\eta}(\mathbf{z}_0, t)|_{\mathbf{z}_0=\eta^{-1}(\mathbf{z}, t)}$, the prescribed vector field $\mathbf{X}_H$ has been replaced by the unknown vector field $\mathcal{X}$. Indeed, in the present treatment the Lagrangian coordinate $\eta$ is left as a dynamical variable and analogously for its velocity $\dot{\eta}$. Then, the associated Eulerian vector field $\mathcal{X}$ must be treated in the same way. With these steps in mind, we observe that the Lagrangian $L_{KvH}$ in (3.2) assumes the new form

$$L_{EPKvH}(\mathcal{X}, D, S) = \int D(\nabla S \cdot \mathcal{X} + (\mathcal{A} - \nabla S) \cdot J \nabla H - H) d^2\mathbf{z}, \qquad (3.5)$$

up to an ignorable total time derivative. For the sake of clarity, here we have indicated the functional dependence of $L_{EPKvH}$ on the dynamical variables. We notice that the variational principle associated to this Lagrangian is not of standard type in that, unlike $\delta S$, the variations $\delta D$ and $\delta \mathcal{X}$ are not arbitrary. Instead, these variations must be found from their definitions in terms of $\eta$. This procedure is a result of a reduction process from Lagrangian to Eulerian variables that is known as Euler-Poincaré reduction [54], which justifies the letters EP in the subscript 'EPKvH' above. Once the KvH Lagrangian is taken into the Euler-Poincaré form (3.5), the equations of motion are obtained by simply adopting standard methods in Euler-Poincaré variational principles [55,56]. In this approach, Eq. (3.4) appears as an auxiliary equation following from the definition (3.3) and accompanying the variational principle $\delta \int_{t_1}^{t_2} L_{EPKvH} dt = 0$.

$$\delta D = -\mathrm{div}(D\mathcal{Y}) \qquad (3.6)$$

follows from taking variations of the Lagrange-to-Euler map (3.3), upon defining the arbitrary displacement vector field $\mathcal{Y}(\mathbf{z}, t) = \delta\eta(\mathbf{z}_0, t)|_{\mathbf{z}_0=\eta^{-1}(\mathbf{z}, t)}$. By proceeding analogously, a more cumbersome vector calculus exercise leads to

$$\delta \mathcal{X} = \partial_t \mathcal{Y} + \mathcal{X} \cdot \nabla \mathcal{Y} - \mathcal{Y} \cdot \nabla \mathcal{X}. \qquad (3.7)$$

We verify that the variations $\delta \mathcal{X}$ and $\delta D$ produce (3.1), which is equivalent to the second equation in (2.5) up to an irrelevant number. In addition, arbitrary variations $\delta S$ lead to $\mathrm{div}(D\mathcal{X}) = \mathrm{div}(D\mathbf{X}_H)$ so that from (3.4) one recovers the first equation in (2.5). Thus, Hamilton's principle $\delta \int_{t_1}^{t_2} L_{EPKvH} dt = 0$ with the Lagrangian (3.5) does indeed recover KvH dynamics in the form (2.5).

### 3.2. Alternative variational approach to KvN theory

Following the observations preceding Section 3.1, let us now perform the replacement

$$\nabla S \to \mathcal{A}$$

in the Euler-Poincaré Lagrangian (3.5). We obtain

$$L_{EPL} = \int D(\mathcal{A} \cdot \mathcal{X} - H) d^2\mathbf{z}. \qquad (3.8)$$

Here, the subscript 'EPL' anticipates the next result: the replacement $\nabla S \to \mathcal{A}$ recovers the Euler-Poincaré variational formulation of the classical Liouville equation. Indeed, upon recalling that $\mathcal{A} = (p, 0)$ is constant in time, the variations (3.6) and (3.7) produce $\mathcal{X} = \mathbf{X}_H$ so that the auxiliary equation $\partial_t D + \mathrm{div}(D\mathcal{X}) = 0$ becomes the classical Liouville equation. Notice that, upon writing $\mathcal{X}$ in terms of its velocity and force components, that is $\mathcal{X}(\mathbf{z}, t) = (u(\mathbf{z}, t), f(\mathbf{z}, t))$, explicitly one has $\mathcal{A} \cdot \mathcal{X} = p u(q, p, t)$. However, it may be convenient to retain the full $\mathcal{A}$-notation since the symplectic potential may acquire different suitable forms depending on the problem under consideration. For example, as noticed in [3], the symplectic potential $\mathcal{A}(\mathbf{z}) = J\mathbf{z}/2$ is particularly useful in the case of quadratic Hamiltonians.
At this point, as presented extensively in [42], one can perform the replacement $D = |\chi|^2$. Taking the square root of (3.3) yields the following Lagrange-to-Euler map for Koopman-von Neumann wavefunctions:

$$ \chi(\mathbf{z}, t) = \frac{\chi_0(\mathbf{z}_0)}{\sqrt{\mathcal{J}_\eta(\mathbf{z}_0, t)}} \Big|_{\mathbf{z}_0=\eta^{-1}(\mathbf{z}, t)}. \qquad (3.9) $$

Consequently, an alternative KvN Lagrangian to the one in (2.3) may be written as

$$ L_{EPK vN} = \int |\chi|^2 (\mathcal{A} \cdot \mathcal{X} - H) d^2 z, \qquad (3.10) $$

which is manifestly gauge invariant under phase transformations $\chi \mapsto e^{i\varphi/\hbar}\chi$. Here, the time derivative $\partial_t\chi$ and the variations $\delta\chi$ follow from (3.9) and read

$$ \partial_t \chi = -\mathcal{X} \cdot \nabla \chi - \frac{1}{2} \chi \operatorname{div} \mathcal{X}, \qquad \delta\chi = -\mathcal{Y} \cdot \nabla \chi - \frac{1}{2}(\operatorname{div} \mathcal{Y})\chi, \qquad (3.11) $$

respectively. The variations $\delta\mathcal{X}$ are still given by (3.7) and they produce the relation $\mathcal{X} = \mathcal{X}_H$, so that the first equation in (3.11) recovers the KvN equation in the form $\partial_t\mathcal{X} = \{H, \mathcal{X}\}$. In the general case, the first equation in (3.11) identifies the transport equation for a wavefunction, so that the well-known density transport in (3.4) is indeed recovered by setting $D = |\chi|^2$. In this construction, no ambiguity occurs when the Koopman phase is set to zero and the total energy coincides with the KvN prescription $\int H|\chi|^2 d^2z$.

Before concluding this section, we present the noncanonical Hamiltonian structure of the KvN equation, as it arises from the variational principle (3.10). In general, a Hamiltonian structure comprises a Hamiltonian functional $h = h(\chi)$ and a Poisson bracket $\{, \}_{\chi}$, so that any functional $f(\chi)$ evolves according to $df/dt = \{f, h\}_{\chi}$. A Poisson bracket is a Lie bracket satisfying the Leibniz product rule $\{f, g\}_{\chi} = \{f, g\}_{\chi}k + \{f, g\}_{\chi}g$. The main difficulty in constructing a Poisson bracket consists in ensuring the Jacobi identity for cyclic permutations. There are lots of examples in the literature where candidate Poisson brackets fail to satisfy the Jacobi identity; see [57] for a discussion in the context of hybrid quantum-classical dynamics. Thus, the safest way to construct a Poisson bracket is usually to derive it from another bracket which is already known to be Poisson. As discussed in [10,42], the Poisson bracket and Hamiltonian functional associated to the KvN equation read, respectively,

$$ \{f, k\}_{\chi} = \int \frac{1}{|\chi|^2} \left( \chi \diamond \frac{\delta f}{\delta \chi} \right) \cdot J \left( \chi \diamond \frac{\delta k}{\delta \chi} \right) d^2 z, \qquad h(\chi) = \int H |\chi|^2 d^2 z, \qquad (3.12) $$

where we introduce

$$ \chi \diamond \frac{\delta f}{\delta \chi} := \frac{1}{2} \operatorname{Re} \left( \frac{\delta h^*}{\delta \chi} \nabla \chi - \chi^* \nabla \frac{\delta f}{\delta \chi} \right) \qquad (3.13) $$

for compactness of notation. As we will see, an analogue of this Poisson structure will reappear in the context of hybrid quantum-classical dynamics.

While this section has presented a prescription for obtaining an alternative variational structure for KvN theory that is immune of the ambiguities discussed in Section 2.1, little was said about how this prescription can be justified within KvH theory. As already mentioned, the equality $\nabla S = \mathcal{A}$ requires a singular phase which we want to avoid in the present work. Thus, the present wavefunction treatment is insufficient by itself to obtain the KvN construction as a type of reduced model of KvH dynamics. An extension of this treatment that overcomes this difficulty is presented in Appendix A, which shows how KvN can be realized as an exact solution of the KvH equation for density matrices.

## 4. Hybrid quantum-classical dynamics

Having characterized the different Koopman formulations of classical mechanics, we move on to discussing the Koopman approach to mixed quantum-classical systems. In particular, here we will review a recent approach proposed by the authors and based on the KvH construction [3,38]. This approach is heavily inspired by Sudarshan's early work [7], which however was based on KvN theory. As we have seen in Section 2.2, the KvH approach differs substantially from KvN by the fact that the former has a nontrivial phase dynamics. In the present context this phase dynamics is crucial in realizing the quantum-classical coupling.

While the approach discussed in this section succeeds in capturing essential properties such as quantum uncertainty, it is not clear whether the associated phase-space distribution of the classical subsystem remains positive in time. So far, this particular point was partly addressed in [38], where we identified an infinite family of hybrid systems for which the classical subsystem has a positive Liouville density at all times. However, a general statement is still lacking. Later on in this paper we will overcome this potential issue by presenting a closure model ensuring a positive classical density. At this stage, we simply review the original approach which our closure model is based on.

### 4.1. Hybrid wavefunctions

As explained in [3], the interaction of a quantum and a classical particle can be formulated by starting with the KvH equation for the two-particle wavefunction $\Upsilon(\mathbf{z}_1, \mathbf{z}_2)$ and then applying canonical quantization to quantize one of them. This approach was also partly followed in [58]. The quantization process follows the usual rules: in particular, we make the replacements $q_2 \to \hat{x} = x$ and $p_2 \to \hat{p} = -i\hbar\partial/\partial x$. Then, upon setting $\partial\Upsilon/\partial p_2 = 0$ and dropping the subscripts, we obtain the quantum-classical wave equation:

$$ i\hbar\partial_t\Upsilon = \{i\hbar\hat{H}, \Upsilon\} + (\hat{H} - p\partial_p\hat{H})\Upsilon. \qquad (4.1) $$
Here, $\Upsilon(q, p, x)$ is a hybrid wavefunction, that is a square-integrable function of both the classical and quantum coordinates $z = (q, p)$ and $x$, respectively. Also, $\hat{H}(q, p, \dot{x}, \dot{p})$ is an operator-valued function on the classical phase-space. Following the discussion in Sections 2.1 and 2.2, here we can write the equation above as

$$i\hbar\partial_t\Upsilon = \hat{\mathcal{L}}_{\hat{H}}\Upsilon, \quad \text{with} \quad \hat{\mathcal{L}}_{\hat{H}} := \{i\hbar\hat{H}, \} - p\partial_p\hat{H} + \hat{H}. \qquad (4.2)$$

The quantum-classical Liouvillian operator $\hat{\mathcal{L}}_{\hat{H}}$ is again Hermitian so that, similarly to quantum and classical mechanics, hybrid dynamics is described as a unitary flow on the space of quantum-classical wavefunctions.

The hybrid wave equation (4.1) possesses a profound geometric structure in symplectic geometry as discussed in [38]. While here we will not enter the details, we will however emphasize that Eq. (4.1) possesses a Dirac-Frenkel variational principle of the type (2.3) and (2.7). In particular, Eq. (4.1) arises as the Euler-Lagrange equation associated to the Lagrangian

$$L_{QC}(\Upsilon, \partial_t \Upsilon) = \mathrm{Re} \int (i\hbar \Upsilon^* \partial_t \Upsilon - \Upsilon^* \hat{\mathcal{L}}_{\hat{H}} \Upsilon) d^2z \, dx, \qquad (4.3)$$

where the last term identifies the total energy $h(\Upsilon) = \mathrm{Re} \int \Upsilon^* \hat{\mathcal{L}}_{\hat{H}} \Upsilon \, d^2z \, dx$.

## 4.2. Quantum and classical densities

The Lagrangian (4.3) is important in that its associated total energy naturally identifies a distribution-valued von Neumann operator playing the role of a hybrid quantum-classical density. Indeed, integration by parts yields

$$\mathrm{Re} \int \Upsilon^* \hat{\mathcal{L}}_{\hat{H}} \Upsilon \, d^2z \, dx = \mathrm{Tr} \int \hat{H} \hat{D} \, d^2z \qquad (4.4)$$

$$\text{with} \qquad \hat{D}(\mathbf{z}) := \Upsilon(\mathbf{z})\Upsilon^{\dagger}(\mathbf{z}) + \partial_p(p\Upsilon(\mathbf{z})\Upsilon^{\dagger}(\mathbf{z})) + i\hbar\{\Upsilon(\mathbf{z}), \Upsilon^{\dagger}(\mathbf{z})\}. \qquad (4.5)$$

Here, '$\dagger$' denotes the quantum adjoint so that, for example, the measure-valued operator $\Upsilon(\mathbf{z})\Upsilon^{\dagger}(\mathbf{z})$ has matrix elements $\Upsilon(\mathbf{z}, x)\Upsilon^*(\mathbf{z}, x')$, while $\Upsilon^{\dagger}(\mathbf{z})\Upsilon(\mathbf{z}) = \int |\Upsilon(\mathbf{z}, x)|^2 dx$ is a phase-space density. The right-hand side of the relation (4.4) confers the distribution-valued operator $\hat{D}(\mathbf{z})$ the role of a hybrid density for the calculation of expectation values. Given a hybrid observable represented by an operator-valued function $\hat{A}(\mathbf{z})$, we prescribe the following definition of its expectation value: $\langle \hat{A} \rangle = \mathrm{Tr} \int \hat{D}(\mathbf{z})\hat{A}(\mathbf{z}) \, d^2z$. The hybrid Ehrenfest theorem for expectation value dynamics was presented in [3], where we showed that the total quantum-classical momentum $\langle p + \hat{p} \rangle$ remains conserved in time for translation-invariant Hamiltonians. Recent work [59] has investigated this result in the context of the Galilean covariance of the quantum-classical wave equation (4.2).

We notice that the hybrid density operator (4.5) is sign-indefinite. Despite this potential issue, we observe that the quantum density matrix

$$\hat{\rho} = \int \hat{D}(\mathbf{z}) \, d^2z = \int \Upsilon(\mathbf{z})\Upsilon^{\dagger}(\mathbf{z}) \, d^2z \qquad (4.6)$$

is actually positive-definite by construction and this property ensures that the Heisenberg uncertainty principle is satisfied. This is a crucial ingredient in the study of quantum-classical coupling. For example, the AG equation (1.2) fails to satisfy this property. To see how decoherence is captured in the present hybrid context, we write the dynamics of the quantum density matrix as [3]

$$i\hbar \frac{\partial \hat{\rho}}{\partial t} = \int [\hat{H}, \hat{D}] \, d^2z. \qquad (4.7)$$

In general, this equation does not allow for an initial pure quantum state $\hat{\rho}_0 = \psi_0 \psi_0^\dagger$ to remain pure over time. Indeed, upon defining the quantum purity as $\mathrm{Tr}(\hat{\rho}^2)$, we see that $d\mathrm{Tr}(\hat{\rho}^2)/dt \neq 0$ thereby leading to quantum decoherence. As we will see below, purity non-preservation also affects the classical evolution.

Let us look closely at the expression of the classical Liouville distribution, that is

$$\rho_c(\mathbf{z}) = \mathrm{Tr} \hat{D}(\mathbf{z}) = \int (|\Upsilon(\mathbf{z})|^2 + \partial_p(p|\Upsilon(\mathbf{z})|^2) + \hbar \mathrm{Im}\{\Upsilon^*(\mathbf{z}), \Upsilon(\mathbf{z})\}) \, dx. \qquad (4.8)$$

Importantly, this quantity is generally sign-indefinite. At present, we do not know whether an initially positive $\rho_c$ stays always positive in time or may develop negative values. The evolution on the classical distribution $\rho_c$ is governed by the equation

$$\frac{\partial \rho_c}{\partial t} = \mathrm{Tr}(\hat{H}, \hat{D}). \qquad (4.9)$$

So far, we could only prove that this equation does not change the initial sign of $\rho_c$ whenever the hybrid Hamiltonian $\hat{H}$ depends only on a set of mutually commuting operators; see [38] for further details. This property was confirmed in [3] for the case of a pure-dephasing problem with quadratic coupling. In this particular case, we observed that both the quantum density (4.6) and the classical density (4.8) remain positive in time while the hybrid density (4.5) develops negative eigenvalues.

As the question whether (4.8) remains always positive is currently open, one may use analogies with Wigner functions to justify the possible emergence of negative values in the classical distribution in the presence of quantum-classical interaction [3]. Alternatively, one may want to enforce a positive Liouville density in some way. We will pursue this second direction later on in this paper.

Before closing this section, we observe that Eq. (4.9) does not allow for Klimontovich particle solutions. Specifically, a $\delta$-like initial condition $\rho_{c0}(\mathbf{z}) = \delta(\mathbf{z} - \zeta_0)$ gets spread across phase-space as time goes by. While this may sound somewhat reminiscent of diffusion processes, we emphasize that no diffusion mechanism takes place in this case; at least, not in any standard sense. Since Klimontovich states are classical pure states [60–62], we conclude that this classical form of purity non-preservation is analogous to the decoherence
appearing in the quantum sector. This point raises the foundational question whether this phenomenon in the classical sector is related to the uncontrollable disturbance invoked by Bohr in his original interpretation of quantum measurements.

Also, we notice that while the loss of classical pure states leads to a form of statistical uncertainty, the latter is very different from Heisenberg's relation $\Delta\langle\hat{x}\rangle\Delta\langle\hat{p}\rangle \ge \hbar/2$ which establishes a well-characterized lower bound. Thus, the form of classical decoherence experienced by the classical system cannot be considered as a loss of its 'classicality'; instead, it simply means that this system must be described by a statistical phase-space distribution in the standard Liouville picture.

## 5. Exact factorization of hybrid wavefunctions

In this section, we formulate a closure scheme of the quantum-classical wave equation (4.2) so that both the quantum and classical densities are positive at all times. While the quantum density operator (4.6) is already positive-definite by construction, we will ensure a positive classical density by exploiting a gauge-invariance principle, as outlined in the Introduction. In particular, we will modify the original theory from Section 4 to enforce a gauge symmetry with respect to local phase factors for which the phase is a function of the classical phase space. For the purely classical case, a similar procedure was followed in Sections 3.1 and 3.2. However, before this procedure can be applied in the hybrid context, we face the difficulty that the classical degrees of freedom cannot be isolated in a simple way since the current form of the theory involves a hybrid wavefunction $\Upsilon(\mathbf{z}, x)$ in which quantum and classical coordinates, respectively $x$ and $\mathbf{z} = (q, p)$, are treated on an equal footing.

In order to circumvent this difficulty, we resort to a method from chemical physics [63]. Known under the name exact factorization, in our context this method simply consists in rewriting the hybrid wavefunction as follows:

$$ \Upsilon(\mathbf{z}, x, t) = \chi(\mathbf{z}, t)\psi(x, t; \mathbf{z}), \quad \text{with} \quad \int |\psi(x, t; \mathbf{z})|^2 dx = 1, \qquad (5.1) $$

where the semicolon indicates that the classical coordinates appear in the expression of the quantum wavefunction $\psi$ merely as parameters. In other words, one has a Koopman wavefunction $\chi(\mathbf{z}, t)$ and a Schrödinger wavefunction $\psi(x, t; \mathbf{z})$, where the latter is parameterized by the classical coordinates $\mathbf{z} = (q, p)$. While so far the quantum state space was identified with the infinite-dimensional space $L^2(\mathbb{R})$ of square-integrable Schrödinger wavefunctions, in what follows we shall consider an arbitrary quantum Hilbert space $\mathcal{H}_Q$ with inner product $\langle \cdot | \cdot \rangle$, induced norm $\| \cdot \|$, and real-valued pairing

$$ \langle \cdot, \cdot \rangle = \mathrm{Re} \langle \cdot | \cdot \rangle. \qquad (5.2) $$

For later purpose, it is convenient to introduce the shorthand notation

$$ \langle \hat{A} \rangle = \langle \psi, \hat{A} \psi \rangle \qquad (5.3) $$

for any Hermitian operator-valued function on phase-space $\hat{A} = \hat{A}(\mathbf{z})$. We hope that no confusion arises with the expectation value notation introduced in Section 4.2. Notice that both wavefunctions $\chi$ and $\psi$ in the factorization (5.1) are defined up to an arbitrary phase factor $e^{i\psi(t, \mathbf{z})}$.

The name exact factorization arises from the fact that the relation (5.1) generally identifies an exact solution of (4.1), as long as $\chi$ is nowhere vanishing. Over the years, factorizations of the type (5.1) appeared in the context of standard quantum mechanics. For example, this is the common approach to the hydrodynamic formulation of the Pauli equation [64]. However, it was only in [63] that this type of wavefunction factorization was recognized to have more general validity. The geometric underpinning of (5.1) was studied in [65–67] and the results therein provide the basis for the present work.

Importantly, we observe that, upon writing $\chi = \sqrt{D}e^{iS/\hbar}$, the exact factorization (5.1) transforms the original expression (4.8) into the form

$$ \rho_c = D + \mathrm{div}(DJ(\nabla S + \mathcal{A}_B - \mathcal{A})), \qquad (5.4) $$

where we have introduced the Berry connection

$$ \mathcal{A}_B := \langle \psi, -i\hbar\nabla\psi \rangle. \qquad (5.5) $$

Thus, similarly to the arguments in Section 3, we realize that setting

$$ \nabla S + \mathcal{A}_B = \mathcal{A} \qquad (5.6) $$

would lead to a positive definite classical density $\rho_c = |\chi|^2$. Motivated by this observation, we will perform various steps in order to apply the method in Sections 3.1 and 3.2 within the present hybrid setting. As we will see, the relation (5.6) leads to a gauge-invariant quantum-classical theory for which classical phases merely represent a gauge freedom and are thus unobservable.

The discussion in the following sections provides an alternative treatment to that in [10], where the same closure model was obtained in terms of density matrices by exploiting the mathematical methods from geometric mechanics [55,56]. Instead, the present treatment follows a more direct approach mostly based on wavefunctions.

### 5.1. Madelung transform and variational approach

As we aim to develop a closure model by using the tools illustrated in Section 3, it is important to combine the variational approach to the Madelung transform with the Euler-Poincaré construction, as presented in Section 3.1. Here we restrict to focus only on the classical degrees of freedom and thus we choose to apply a partial Madelung transform only to the classical part of the wavefunction in (5.1). Thus, upon writing $\chi = \sqrt{D}e^{iS/\hbar}$, we replace the ansatz $\Upsilon(\mathbf{z}, x) = \sqrt{D(\mathbf{z})}e^{iS(\mathbf{z})/\hbar}\psi(x; \mathbf{z})$ into the hybrid Lagrangian (4.3), thereby obtaining

$$ L_{EF} = \int D(\partial_t S - \langle \psi, i\hbar\partial_t \psi \rangle) d^2z + h(D, S, \psi), \qquad (5.7) $$
where the subscript 'EF' stands for Exact Factorization and the Hamiltonian $h(D, S, \psi)$ arising from (4.3) is written as

$$ 
\begin{aligned}
h(D, S, \psi) &= \int D(\mathbf{X}_{\hat{H}}) \cdot (\nabla S + \mathcal{A}_B) + (\hat{\mathcal{L}}_{\hat{H}} - \mathcal{A}_B \cdot \mathbf{X}_{\hat{H}}) \, d^2z \\
&= \int D(\mathbf{X}_{\hat{H}}) \cdot (\nabla S + \mathcal{A}_B - \mathcal{A}) + (\psi, \hat{H}\psi - i\hbar\tilde{\mathbf{X}}_{\hat{H}} \cdot \nabla\psi) \, d^2z.
\end{aligned} 
\quad (5.8) $$

Note that (5.7) is the quantum-classical extension of the Lagrangian (3.2) for the classical case. Here, we recall (5.3) and we have denoted

$$ \mathbf{X}_{\hat{H}} = (\partial_t \hat{H}, -\partial_q \hat{H}), \qquad \tilde{\mathbf{X}}_{\hat{H}} := \mathbf{X}_{\hat{H}} - \langle \mathbf{X}_{\hat{H}} \rangle, \quad (5.9) $$

while $\hat{\mathcal{L}}_{\hat{H}}$ was defined in (4.2). Also, all differential operators act on the classical phase-space and we recall the symplectic potential $\mathcal{A} = (p, 0)$ from Section 3. More importantly, in (5.8) we have added and subtracted the Berry connection $\mathcal{A}_B$. In addition, we notice that the condition $\|\psi\|^2 = 1$ is not enforced here as a constraint. As we will see, this condition is preserved in time by the final equations of motion.

Arbitrary variations $\delta S$ and $\delta D$ in the action functional associated to (5.7) yield

$$ \partial_t D + \mathrm{div}(D(\mathbf{X}_{\hat{H}})) = 0, \quad (5.10) $$

$$ \partial_t S + \langle \mathbf{X}_{\hat{H}} \rangle \cdot \nabla S = \langle \psi, i\hbar \partial_t \psi \rangle - \langle \psi, \hat{\mathcal{L}}_{\hat{H}} \psi \rangle \quad (5.11) $$

and we will next focus on the equation for the quantum wavefunction $\psi$.

## 5.2. Quantum wavefunction and classical flow paths

In order to write the quantum evolution equation, it is convenient to introduce the functional

$$ f(D, \psi) := \int D\langle\psi, (\hat{\mathcal{L}}_{\hat{H}} - \mathcal{A}_B \cdot \mathbf{X}_{\hat{H}})\psi\rangle \, d^2z. \quad (5.12) $$

Then, arbitrary variations $\delta \psi$ in $\int_{t_1}^{t_2} L_{EF} dt$ yield [10]

$$ i\hbar(\partial_t + \langle \mathbf{X}_{\hat{H}} \rangle \cdot \nabla)\psi = (\nabla S + \mathcal{A}_B) \cdot \mathbf{X}_{\hat{H}}\psi + \frac{1}{2D} \frac{\delta f}{\delta \psi}, \quad (5.13) $$

where we compute

$$ \frac{\delta f}{\delta \psi} = 2D\hat{H}\psi - 2D(\mathcal{A} + \mathcal{A}_B) \cdot \mathbf{X}_{\hat{H}}\psi - 2i\hbar D\tilde{\mathbf{X}}_{\hat{H}} \cdot \nabla\psi - i\hbar\psi \mathrm{div}(D\tilde{\mathbf{X}}_{\hat{H}}). $$

While the current explicit form of Eq. (5.13) is not particularly insightful, we will be able to make relevant considerations by recognizing that the functional (5.12) can be entirely expressed in terms of the local density matrix $\hat{\rho}(z) = \psi(z)\psi(z)^\dagger$ as follows:

$$ f = \frac{1}{2} \int D\left(\langle \mathbf{X}_{\hat{H}}, i\hbar[\hat{\rho}, \nabla\hat{\rho}] \rangle - 2\langle \mathcal{A} \cdot \mathbf{X}_{\hat{H}} - \hat{H}, \hat{\rho} \rangle\right) \, d^2z. \quad (5.14) $$

Here, we use the real-valued pairing $\langle \hat{A}, \hat{B} \rangle = \mathrm{Re}(\hat{A}|\hat{B}) = \mathrm{Re} \mathrm{Tr}(\hat{A}^\dagger \hat{B})$. Notice the appearance of the non-Abelian gauge connection $[\hat{\rho}, \nabla\hat{\rho}]$ already emerged in Mead's work on molecular geometric phases [68]. More importantly, the chain-rule relation $\delta f / \delta \psi = 2(\delta f / \delta \hat{\rho})\psi$ allows us to rewrite Eq. (5.13) as

$$ i\hbar(\partial_t + \langle \mathbf{X}_{\hat{H}} \rangle \cdot \nabla)\psi = \left( (\nabla S + \mathcal{A}_B) \cdot \mathbf{X}_{\hat{H}} + \frac{1}{D} \frac{\delta f}{\delta \hat{\rho}} \right) \psi, \quad (5.15) $$

where we compute

$$ \frac{\delta f}{\delta \hat{\rho}} = D\hat{H} - D\mathcal{A} \cdot \mathbf{X}_{\hat{H}} + \frac{i\hbar}{2} (D\langle\hat{\rho}, \hat{H}\rangle + D\langle\hat{H}, \hat{\rho}\rangle + [\hat{\rho}, \{D, \hat{H}\}]). $$

At this point, the form of Eq. (5.15) reveals the structure underlying the quantum evolution. Indeed, we realize that, since $\delta f / \delta \hat{\rho}$ acts on $\psi(z, t)$ as a Hermitian operator and so does $(\nabla S + \mathcal{A}_B) \cdot \mathbf{X}_{\hat{H}}$, the parenthesis on the right-hand side of (5.15) generates a unitary propagator parameterized by the phase-space coordinates. This implies $(\partial_t + \langle \mathbf{X}_{\hat{H}} \rangle \cdot \nabla)\|\psi\|^2 = 0$, so that the partial normalization condition $\|\psi(z, t)\|^2 = 1$ is preserved in time. More importantly, we notice the appearance of the material derivative $\partial_t + \langle \mathbf{X}_{\hat{H}} \rangle \cdot \nabla$ already present in (5.11). This indicates that the unitary quantum evolution occurs in a phase-space frame moving with the Lagrangian flow path $\eta(z_0, t)$ generated by the vector field $\langle \mathbf{X}_{\hat{H}} \rangle$, which transports the density $D$ as in (5.10). This important remark was recently exploited in [10]. In the fully quantum case, a similar approach was previously followed in [65].

Having uncovered the hydrodynamic Lie-transport features underlying the exact factorization system, we are motivated to cast the Lagrangian (5.7) in Euler-Poincaré form, which is a more natural setting for hydrodynamic continuum theories. Following the arguments in Section 3.1, we consider the Lagrange-to-Euler map (3.3) and write the Lagrangian (5.7) in Euler-Poincaré form as

$$ L_{EPEF} = \int D(\nabla S \cdot \mathcal{X} + \langle \psi, i\hbar \partial_t \psi \rangle) \, d^2z - h(D, S, \psi), \quad (5.16) $$

where $h(D, S, \psi)$ is given in (5.8). Note that (5.16) is the quantum-classical extension of the Lagrangian (3.5) for KvH classical mechanics. We recall that the vector field $\mathcal{X}$ is defined in terms of the Lagrangian phase-space path $\eta(z_0, t)$ as $\dot{\eta}(z_0, t) = \mathcal{X}(\eta(z_0, t), t)$, so that one
has the constrained variations (3.6)–(3.7) along with arbitrary variations $\delta S$. One can verify that, up to an irrelevant time-dependent number, the Lagrangian (5.16) leads to Eq. (5.11), while the auxiliary equation $\partial_t D + \text{div}(D\mathcal{X}) = 0$ recovers (5.10). In addition, the quantum evolution equation (5.13) follows from taking arbitrary variations $\delta \psi$.

So far the exact factorization (5.1) has not provided much insight into the nature of quantum-classical coupling, which still looks very intricate. However, as anticipated, the variational setting made available by the Lagrangian (5.16) allows us to apply the method in Section 3 thereby leading to a closure model ensuring that the classical density (5.4) is positive definite at all times. This is done in the next section.

# 6. A nonlinear hybrid wave equation

Let us now extend the method from Section 3 to the hybrid setting and use that as a tool for ensuring a positive-definite classical density.

## 6.1. The quantum-classical Lagrangian

In first place, we recall that the exact factorization (5.1) transforms the original expression (4.8) into (5.4). Thus, here we want to be able to use the relation (5.6) by following the method in Section 3. In this way, the classical density becomes $\rho_c = |\chi|^2$, which is positive-definite. Notice that, unlike the classical case treated in Section 3, the relation (5.6) may not be preserved by the hybrid system (5.10), (5.11), and (5.13). As we show in Appendix C, the wavefunction treatment can be extended to von Neumann operators, in which case relation (5.6) can be realized as a dynamical closure ansatz. In the present discussion, however, we will simply follow the procedure from Sections 3.1 and 3.2, which is entirely based on wavefunctions. In particular, we will make use of the relation (5.6) by replacing $\nabla S \to \mathcal{A} - \mathcal{A}_B$ in the Lagrangian (5.16). Then, the latter becomes gauge-independent and acquires the form

$$ \ell(\mathcal{X}, D, \psi, \partial_t \psi) = \int D((\mathcal{A} - \mathcal{A}_B) \cdot \mathcal{X} + \langle \psi, i\hbar \partial_t \psi - \hat{H}\psi + i\hbar \tilde{\mathbf{X}}_{\hat{H}} \cdot \nabla \psi \rangle) d^2z, \quad (6.1) $$

where the last two terms are found by replacing (5.6) in (5.8). At this point, the relation (5.4) yields $\rho_c = D = |\chi|^2$, thereby recovering a positive classical density at all times. Indeed, we observe that the initial sign of $D$ is preserved by its dynamics, which is given by transport equation $\partial_t D + \text{div}(D\mathcal{X}) = 0$.

Upon using $D = |\chi|^2$, one may restore the hybrid quantum-classical wavefunction $\Upsilon = \chi \psi$. For this purpose, we notice that $\langle \Upsilon, -i\hbar \nabla \Upsilon \rangle = D\nabla S + D\langle \psi, -i\hbar \nabla \psi \rangle$, where we recognize the appearance of the Berry connection (5.5) on the right hand-side. Since the addition of the pure differential $\nabla S$ is irrelevant, we are led to writing

$$ \mathcal{A}_B = \frac{\langle \Upsilon, -i\hbar \nabla \Upsilon \rangle}{\|\Upsilon\|^2}. \quad (6.2) $$

Here, the angle brackets denote the quantum pairing, that is $\langle \Upsilon_1, \Upsilon_2 \rangle = \text{Re}(\Upsilon_1^\dagger(z)\Upsilon_2(z))$, so that the notation (5.3) leads to writing $\langle \hat{A} \rangle = \langle \Upsilon, \hat{A}\Upsilon \rangle / \|\Upsilon\|^2$. With this notation, the Lagrangian (6.1) becomes

$$ \ell_{QC}(\mathcal{X}, \Upsilon, \partial_t \Upsilon) = \int \langle \Upsilon, i\hbar \partial_t \Upsilon + (\mathcal{A} - \mathcal{A}_B) \cdot \mathcal{X} \Upsilon \rangle d^2z - h_{QC}(\Upsilon), \quad (6.3) $$

where the quantum-classical Hamiltonian functional now reads

$$ \begin{aligned} h_{QC}(\Upsilon) &= \int \langle \Upsilon, \hat{\mathcal{L}}_{\hat{H}} \Upsilon + (\mathcal{A} - \mathcal{A}_B) \cdot \mathbf{X}_{\hat{H}} \Upsilon \rangle d^2z \\ &= \int \langle \Upsilon, \hat{H}\Upsilon - i\hbar \tilde{\mathbf{X}}_{\hat{H}} \cdot \nabla \Upsilon \rangle d^2z. \end{aligned} \quad (6.4) $$

Here, the variation $\delta\Upsilon$ is arbitrary and $\delta\mathcal{X}$ is given by (3.7). Notice that, while the original hybrid Hamiltonian $h_{QC}(\Upsilon) = \int \langle \Upsilon, \hat{\mathcal{L}}_{\hat{H}} \Upsilon \rangle d^2z dx$ appearing in (4.3) is manifestly covariant with respect to gauge transformations $\Upsilon(z) \mapsto e^{i\psi(z)/\hbar}\Upsilon(z)$, the new reduced Hamiltonian (6.4) is manifestly gauge invariant. This gauge invariance under local phase factors is a long-sought property in the theory of quantum-classical wavefunctions and its absence in previous hybrid models has been referred to as a “severe problem” [5]. The fact that the present theory is immune from these issues, thereby ensuring positive quantum and classical densities, represents an important step forward in the theory of quantum-classical coupling.

In addition, we observe that the total energy identified by the Hamiltonian functional (6.4) clearly differs from the usual expression of the expectation value of the Hamiltonian operator, that is the first term $\int \langle \Upsilon, \hat{H}\Upsilon \rangle d^2z$. A possible interpretation of this fact is that the energy balance of quantum-classical interactions involves an extra work – here given by the second term in (6.4) – that is produced by correlation effects. This seems to be the possibility suggested by the authors of [69]. Alternatively, we may simply use integration by parts to isolate a form of the hybrid quantum-classical density that is alternative to that presented in (4.5). This second approach will be discussed in Section 7.3. Instead, the next section proceeds by writing the equation for $\Upsilon$.

## 6.2. Hybrid wave equation

At this point, we are ready to take variations and write the hybrid wave equation for $\Upsilon$. Indeed, taking arbitrary variations $\delta\Upsilon$ in (6.3) yields a nonlinear wave equation of the form

$$ i\hbar \frac{\partial \Upsilon}{\partial t} + i\hbar \mathcal{X} \cdot \nabla \Upsilon + \frac{i\hbar}{2} \Upsilon \operatorname{div} \mathcal{X} = \frac{1}{2} \frac{\delta h_{QC}}{\delta \Upsilon} - \mathcal{X} \cdot \mathcal{A} \Upsilon, \quad (6.5) $$
where we compute

$$ \frac{1}{2} \frac{\delta h_{QC}}{\delta \Upsilon} = (\hat{H} - \mathcal{A}_B \cdot \tilde{\mathbf{X}}_{\hat{H}})\Upsilon - i\hbar \tilde{\mathbf{X}}_{\hat{H}} \cdot \nabla \Upsilon - \frac{1}{2} i\hbar (\text{div} \, \tilde{\mathbf{X}}_{\hat{H}})\Upsilon. \qquad (6.6) $$

Here, we have used (6.4) and noticed that $\delta\langle\mathbf{X}_{\hat{H}}\rangle/\delta\Upsilon = 2\|\Upsilon\|^{-2}\tilde{\mathbf{X}}_{\hat{H}}\Upsilon$. Then, Eq. (6.5) becomes

$$ i\hbar \frac{\partial \Upsilon}{\partial t} + i\hbar(\mathbf{X} + \tilde{\mathbf{X}}_{\hat{H}}) \cdot \nabla \Upsilon + \frac{i\hbar}{2} \text{div}(\mathbf{X} + \tilde{\mathbf{X}}_{\hat{H}}) \Upsilon = (\hat{H} - \mathcal{A}_B \cdot \tilde{\mathbf{X}}_{\hat{H}}) \Upsilon - \mathbf{X} \cdot \mathcal{A} \Upsilon. \qquad (6.7) $$

As a consequence, we observe that Eq. (6.5) leads to the following continuity equation for the classical Liouville density $\rho_c = \|\Upsilon\|^2$:

$$ \partial_t \|\Upsilon\|^2 + \text{div}(\|\Upsilon\|^2 \mathbf{X}) = 0. \qquad (6.8) $$

This shows that the vector field $\mathbf{X}$ indeed identifies the Lagrangian trajectories of the classical flow, which advances the classical density according to (3.3), that is $\rho_c(\mathbf{z}, t) = \int \rho_{c0}(\mathbf{z}_0)\delta(\mathbf{z} - \eta(\mathbf{z}_0, t)) d^2\mathbf{z}_0$.

At this stage, however, the vector field $\mathbf{X}$ still needs to be determined. In order to find its expression, we first see that taking variations $\delta\mathbf{X}$ in (6.3) according to (3.7) and making use of (6.8) yields

$$ \begin{aligned} \mathbf{X} \cdot \nabla \mathcal{A} + \nabla \mathbf{X} \cdot \mathcal{A} &= \partial_t \mathcal{A}_B + \mathbf{X} \cdot \nabla \mathcal{A}_B + \nabla \mathbf{X} \cdot \mathcal{A}_B \\ &= \frac{1}{2\|\Upsilon\|^2} \left( \left( \frac{\delta h_{QC}}{\delta \Upsilon}, \nabla \Upsilon \right) - \left( \Upsilon, \nabla \frac{\delta h_{QC}}{\delta \Upsilon} \right) \right) + \nabla(\mathbf{X} \cdot \mathcal{A}). \end{aligned} \qquad (6.9) $$

The second equality is easily verified. Indeed, upon introducing $\Upsilon = \Upsilon/\|\Upsilon\|$, we use $\mathcal{A}_B = \langle\Upsilon, -i\hbar\nabla\Upsilon\rangle$ to notice that (6.5) leads to $i\hbar\|\Upsilon\|(\partial_t + \mathbf{X} \cdot \nabla)\Upsilon = (1/2)\delta h_{QC}/\delta\Upsilon - \mathbf{X} \cdot \mathcal{A}\Upsilon$. The latter can then be used to expand the terms after the first equality in (6.9). Therefore, upon using the defining relation $(\nabla\mathcal{A})^T - \nabla\mathcal{A} = J$, we obtain

$$ \begin{aligned} \mathbf{X} &= \frac{1}{2\|\Upsilon\|^2} J \left( \langle \Upsilon, \nabla \frac{\delta h_{QC}}{\delta \Upsilon} \rangle - \left\langle \frac{\delta h_{QC}}{\delta \Upsilon}, \nabla \Upsilon \right\rangle \right) \\ &= \frac{1}{\|\Upsilon\|^2} J \left[ \|\Upsilon\|^2 \langle \nabla(\hat{H} - \tilde{\mathbf{X}}_{\hat{H}} \cdot \mathcal{A}_B) \rangle + \mathcal{A}_B \langle \Upsilon, (\text{div} \, \tilde{\mathbf{X}}_{\hat{H}})\Upsilon \rangle + 2 \langle \nabla \Upsilon, \tilde{\mathbf{X}}_{\hat{H}} \cdot (i\hbar\nabla)\Upsilon \rangle - \nabla \langle \Upsilon, \tilde{\mathbf{X}}_{\hat{H}} \cdot (i\hbar\nabla)\Upsilon \rangle \right], \end{aligned} $$

where the second equality is proved in Appendix B using (6.6). Alternatively, a further vector calculus exercise shows that

$$ \mathbf{X} = \langle \mathbf{X}_{\hat{H}} \rangle + \frac{1}{2\|\Upsilon\|^2} \left( \langle \|\Upsilon\|^2 J \hat{\Gamma} \cdot \nabla, \mathbf{X}_{\hat{H}} \rangle - \langle \mathbf{X}_{\hat{H}} \cdot \nabla, \|\Upsilon\|^2 J \hat{\Gamma} \rangle \right). \qquad (6.10) $$

where $\hat{\Gamma} = i\hbar\|\Upsilon\|^{-2}[\Upsilon\Upsilon^\dagger, \nabla(\Upsilon\Upsilon^\dagger)] = 2\mathcal{A}_B\Upsilon\Upsilon^\dagger + i\hbar\|\Upsilon\|^{-2}(\Upsilon\nabla\Upsilon^\dagger - (\nabla\Upsilon)\Upsilon^\dagger)$ is a non-Abelian gauge potential already appeared in Eq. (5.14).

Then, together with (6.10), Eq. (6.7) gives a nonlinear hybrid wave equation (NHWE) obtained by a suitable modification of the variational structure underlying the original quantum-classical wave equation (4.2). The NHWE is very complicated in several ways, due to the cumbersome expression of the phase-space vector field $\mathbf{X}$. Thus, suitable numerical methods will need to be developed to approach the associated initial value problem and we are currently pursuing this direction. Here, we will simply notice that, in both the purely quantum and the purely classical cases, one has $\tilde{\mathbf{X}}_{\hat{H}} = 0$. Then, in the purely quantum case one recovers the standard Schrödinger equation $i\hbar\partial_t\Upsilon = \hat{H}\Upsilon$, while in the purely classical case, one has $\mathbf{X} = \mathbf{X}_H$ and (6.7) returns the KvH equation $i\hbar\partial_t\Upsilon = \hat{\mathcal{L}}_H\Upsilon$. Notice that in the present context the phase term $-(\mathbf{A} \cdot \mathbf{X}_H - H)$ appearing in the prequantum operator $\hat{\mathcal{L}}_H$ is an irrelevant phase factor. Indeed, we need to recall that the classical density is now given by $\rho_c = \|\Upsilon\|^2$, rather than the expression in (4.8).

The fluctuation vector field $\tilde{\mathbf{X}}_{\hat{H}} = \mathbf{X}_{\hat{H}} - \langle\mathbf{X}_{\hat{H}}\rangle$ deserves some words. Notice that $\langle\tilde{\mathbf{X}}_{\hat{H}}\rangle = 0$ and if

$$ \hat{H}(q, p) = \hat{H}_Q(\hat{x}, \hat{p}) + H_C(q, p) + \hat{V}_I(q, \hat{x}), \qquad (6.11) $$

then the expression

$$ \tilde{\mathbf{X}}_{\hat{H}} = (0, \hat{F} - \langle \hat{F} \rangle), \quad \text{with} \quad \hat{F} = -\partial_q \hat{V}_I(q, \hat{x}), \qquad (6.12) $$

identifies the force fluctuation $\tilde{F} := \hat{F} - \langle \hat{F} \rangle$ from the Hellmann-Feynman average $\langle \hat{F} \rangle$ [70]. As we will see later on, neglecting the fluctuation force $\hat{F}$ leads to a particularly simple model.

In the remainder of this paper, we will discuss various interesting features arising from the NHWE. Before moving on to this discussion, we point out that the NHWE may be regarded as a closure scheme based on von Neumann operators and their Wigner transform, similarly to the discussion in Appendix A. This is the topic of Appendix C.

## 7. Discussion

This section presents a series of implications of the NHWE model comprised by Eqs. (6.7) and (6.10). As we will see, this model enjoys several properties including the presence of invariants associated to a Hamiltonian structure combining those of quantum and classical mechanics alone. Later on in this section, we will consider how the NHWE model recovers important special cases occurring in the chemical physics literature and the comparison between the full model and its special cases will lead us to make considerations on the presence of the quantum backreaction on the classical trajectories. At first, we will present how Poincaré's integral invariant transfers to the quantum-classical setting.
## 7.1. Poincaré integral invariant

The NHWE model discloses an interesting feature: as a consequence of the first line in (6.9), we obtain the hybrid Poincaré integral invariant

$$ \frac{d}{dt} \oint_{c(t)} (pdq - \mathcal{A}_B(\mathbf{z}, t) \cdot d\mathbf{z}) = 0, \qquad (7.1) $$

where $c(t) = \eta(c_0, t)$ and $c_0$ is an arbitrary loop in the classical phase-space. In terms of differential forms, one may introduce the Berry curvature

$$ \mathcal{B} = \nabla \mathcal{A}_B - (\nabla \mathcal{A}_B)^T = \frac{2\hbar}{\|\nabla\|^2} (\mathrm{Im}\langle\nabla\Upsilon|\nabla\Upsilon\rangle - \mathrm{Im}\langle\Upsilon|\nabla\Upsilon\rangle \wedge \mathrm{Re}\langle\Upsilon|\nabla\Upsilon\rangle), $$

where we have used the wedge product notation $\mathbf{v} \wedge \mathbf{w} = \mathbf{vw} - \mathbf{wv}$. Then, applying Stokes theorem to (7.1) yields the following evolution for the exact two-form $\Omega(t) = J - \mathcal{B}(t)$:

$$ \Omega_{jk}(\eta(\mathbf{z}_0, t), t) d\eta^j(\mathbf{z}_0, t) \wedge d\eta^k(\mathbf{z}_0, t) = \Omega_{0jk}(\mathbf{z}_0) dz_0^j \wedge dz_0^k, \qquad (7.2) $$

which is analogous to the relation for Hamiltonian flows in symplectic geometry. However, in this case the differential 2-form $\Omega(t)$ is time-dependent and it remains symplectic at all times if it is so initially.

Notice that, by the usual identification between maximal forms and densities, the Liouville theorem ensures that the density $\det\Omega$ satisfies the continuity equation $\partial_t(\det\Omega) + \mathrm{div}(\mathcal{X}\det\Omega) = 0$. Consequently, any one-variable function $\Phi(x)$ generates an invariant functional

$$ \mathcal{I}(\Upsilon) = \int \|\Upsilon\|^2 \Phi \left( \frac{\det \Omega}{\|\Upsilon\|^2} \right) d^2 z. $$

For example, upon recalling (6.2) and its preceding lines, we observe that in the purely classical case the Berry connection is a pure differential so that the invariant $\int \|\Upsilon\|^2 \ln(\det\Omega/\|\Upsilon\|^2) d^2 z$ recovers the usual expression of Gibbs' entropy. The identification of entropy functionals in the context of quantum-classical dynamics is an interesting question that is currently under investigation.

The existence of a Poincaré integral invariant unfolds some of the symplectic properties of the NHWE (6.7). As recently shown in [10] by using the density matrix formalism, the closure model provided by the NHWE possesses a very rich geometric content and here we continue this discussion by presenting the Hamiltonian structure underlying Eqs. (6.7) and (6.10).

## 7.2. Density matrix and Hamiltonian structure

As the NHWE model is derived from the variational principle associated to the Lagrangian (6.3), it is expected that its equations of motion are Hamiltonian with respect to the Hamiltonian functional (6.4). In order to present the Hamiltonian structure, it is convenient to rewrite the NHWE (6.5) in terms of the density matrix $\hat{\mathcal{P}}(\mathbf{z}) = \Upsilon(\mathbf{z})\Upsilon^\dagger(\mathbf{z})$. This is an easy step upon recognizing that the Hamiltonian functional (6.4) may be rewritten as

$$ h_{QC}(\hat{\mathcal{P}}) = \int \left\langle \hat{\mathcal{P}}, \hat{H} + \frac{i\hbar}{2\mathrm{Tr}\hat{\mathcal{P}}} \left[ \nabla\hat{\mathcal{P}}, \mathbf{X}_{\hat{H}} \right] \right\rangle d^2z \qquad (7.3) $$

and by resorting to the chain rule, that is $\delta h_{QC}/\delta \Upsilon = 2(\delta h_{QC}/\delta \hat{\mathcal{P}})\Upsilon$. These relations allow rewriting (6.5) as $i\hbar\partial_t \hat{\mathcal{P}} + i\hbar \mathrm{div}(\mathcal{X}\hat{\mathcal{P}}) = [\delta h_{QC}/\delta \hat{\mathcal{P}}, \hat{\mathcal{P}}]$. In addition, the vector field $\mathcal{X}$ may be written in terms of $\hat{\mathcal{P}}$ by using the chain rule in the first line of (6.10). One finds $\mathcal{X} = J(\hat{\mathcal{P}}, \nabla(\delta h_{QC}/\delta \hat{\mathcal{P}}))/\mathrm{Tr}\hat{\mathcal{P}} = \langle \mathbf{X}_{\delta h_{QC}/\delta \hat{\mathcal{P}}} \rangle$, where we have used $(\hat{A}) = (\hat{\mathcal{P}}, \hat{A})/\mathrm{Tr}\hat{\mathcal{P}}$. In conclusion, one obtains Eq. (1.4) with $\hat{H} = \delta h/\delta \hat{\mathcal{P}}$ and

$$ \frac{\delta h_{QC}}{\delta \hat{\mathcal{P}}} = \hat{H} + \frac{1}{2\mathrm{Tr}\hat{\mathcal{P}}} \left( 2i\hbar\{\langle\hat{\mathcal{P}}, \hat{H}\rangle + \langle\hat{H}, \hat{\mathcal{P}}\rangle\} + i\hbar\{[\ln \mathrm{Tr}\hat{\mathcal{P}}, \hat{H}], \hat{\mathcal{P}}\} - \langle i\hbar\{\hat{\mathcal{P}}, \hat{H}\} + i\hbar\{\hat{H}, \hat{\mathcal{P}}\}\rangle \mathbf{1} \right). \qquad (7.4) $$

Notice that Eq. (1.4) is gauge-independent and thus the pure phase term $-\mathcal{A} \cdot \mathcal{X}\Upsilon$ in (6.5) plays no role in this setting. In addition, upon using $\langle \mathbf{X}_{\hat{H}} \rangle = \langle \rho, \mathbf{X}_{\hat{H}} \rangle = \mathbf{X}_{(\rho, \hat{H})} - \langle \mathbf{X}_{\rho}, \hat{H} \rangle$, with $\rho = \hat{\mathcal{P}}/\mathrm{Tr}\hat{\mathcal{P}}$, the Leibniz product rule yields the expression of $\mathcal{X}$ in (6.10) where $\hat{H} = i\hbar[\rho, \nabla\rho]$.

Eq. (1.4) was shown to be Hamiltonian in [10]. In particular, upon introducing the notation $A : B = \mathrm{Tr}(AB)$, its Poisson bracket reads

$$ \{k, h\}_{\hat{\mathcal{P}}} = - \int \left\langle \hat{\mathcal{P}}, \frac{i}{\hbar} \left[ \frac{\delta k}{\delta \hat{\mathcal{P}}}, \frac{\delta h}{\delta \hat{\mathcal{P}}} \right] \right\rangle d^2 z + \int \frac{1}{\mathrm{Tr} \hat{\mathcal{P}}} \left( \hat{\mathcal{P}} : \left\{ \frac{\delta k}{\delta \hat{\mathcal{P}}}, \frac{\delta h}{\delta \hat{\mathcal{P}}} \right\} : \hat{\mathcal{P}} \right) d^2 z, \qquad (7.5) $$

and this is accompanied by the Hamiltonian functional (7.3). A thorough discussion of the Poisson bracket (7.5) was given in [10]. The bracket structure for the NHWE is found from (7.5) by simply using the chain rule. One obtains

$$ \{k, h\}_{\Upsilon} = \int \underbrace{\left( \frac{1}{2\hbar} \mathrm{Im} \left\langle \frac{\delta k}{\delta \Upsilon} \middle| \frac{\delta h}{\delta \Upsilon} \right\rangle \right)}_{\text{Quantum Schrödinger bracket}} + \underbrace{\frac{1}{\|\Upsilon\|^2} \left( \Upsilon \diamond \frac{\delta k}{\delta \Upsilon} \right) \cdot J \left( \Upsilon \diamond \frac{\delta h}{\delta \Upsilon} \right)}_{\text{Classical Koopman-von Neumann bracket}} d^2 z. \qquad (7.6) $$

where

$$ \Upsilon \diamond \frac{\delta h}{\delta \Upsilon} := \frac{1}{2 \|\Upsilon\|^2} \left( \left\langle \frac{\delta h}{\delta \Upsilon}, \nabla \Upsilon \right\rangle - \left\langle \Upsilon, \nabla \frac{\delta h}{\delta \Upsilon} \right\rangle \right). $$

Importantly, due to its construction, this bracket structure becomes Poisson when restricted to gauge-invariant functionals, which is indeed the case for (6.4). Notice that the bracket (7.6) returns Eq. (6.7) in the gauge obtained by dropping the pure phase term $-\mathcal{A} \cdot \mathcal{X}\Upsilon$.
The bracket (7.6) has a particularly suggestive form in that it is the sum of the canonical quantum bracket for the evolution of Schrödinger wavefunctions and the classical KvN bracket (3.12) from (3.12). For earlier attempts to construct Poisson brackets in quantum-classical dynamics, we address the reader to [6,71].

### 7.3. Hybrid density operator

In this section we are interested in the analogue of the hybrid density operator (4.5) in the context of the NHWE model. Upon performing a similar integration by parts to the second equality in (4.4), we observe that the NHWE Hamiltonian functional (6.4) may be rewritten as

$$ \int \langle \gamma, \hat{\mathcal{L}}_{\hat{H}} \gamma + (\mathcal{A} - \mathcal{A}_B) \cdot \mathbf{X}_{\hat{H}} \gamma \rangle d^2z = \mathrm{Tr} \int \hat{D}\hat{H} d^2z, \qquad (7.7) $$

where

$$ \hat{D} := \gamma \gamma^{\dagger} - \mathrm{div}(J \mathcal{A}_B \gamma \gamma^{\dagger}) + i\hbar\{\gamma, \gamma^{\dagger}\} \qquad (7.8) $$

and all quantities are evaluated at the same phase-space point $z = (q, p)$. Notice that the original expression in (4.5) is obtained by simply replacing $\mathcal{A}_B \to \mathcal{A}$ in (7.8). As we showed in [10], the hybrid operator (7.8) may also be written in terms of the variable $\hat{P} = \gamma\gamma^{\dagger}$. Similarly to (4.5), the new density operator is sign-indefinite while it consistently recovers a positive quantum density matrix $\hat{\rho}_q = \int \hat{D}d^2z$ and a positive classical density $\rho_c = \mathrm{Tr}\hat{D}$. We recall from the original treatment in Section 4 that the expression (4.8) was not positive-definite.

More importantly, in [10] we proved that the hybrid operator (7.8) is covariant with respect to both canonical transformations $\eta$ in phase-space and unitary transformations $\mathcal{U}$ of the quantum state space. Upon denoting the functional dependence of $\hat{D}(z)$ on $\gamma$ by $\hat{D}[\gamma(z)]$, we have the properties

$$ \hat{D}[\gamma(\eta(z))] = \hat{D}[\gamma](\eta(z)) \quad \text{and} \quad \hat{D}[\mathcal{U}\gamma] = \mathcal{U}\hat{D}[\gamma]\mathcal{U}^{\dagger}. $$

Here, $\mathcal{U}$ is a quantum unitary operator while $\eta$ is a classical canonical transformation, so that $J_{jk} d\eta^j(z) \wedge d\eta^k(z) = J_{jk} dz^j \wedge dz^k$. As showed in [10], these covariance properties lead to casting the equations for the classical density $\rho_c = \|\gamma\|^2$ and the quantum density matrix $\hat{\rho}_q = \int \gamma \gamma^{\dagger} d^2z$ as

$$ \frac{\partial \rho_c}{\partial t} = \mathrm{Tr}\{\hat{H}, \hat{D}\}, \qquad i\hbar \frac{d\hat{\rho}_q}{dt} = \int [\hat{H}, \hat{D}] d^2z, \qquad (7.9) $$

which are formally identical to the analogous relations (4.9) and (4.7) in the original treatment. The relations (7.9) govern the Ehrenfest equations for quantum and classical expectation values. In general, we define the expectation value of a hybrid observable $\hat{A} = \hat{A}(z)$ as $\langle \hat{A} \rangle_{\gamma} := \mathrm{Tr} \int \hat{D}\hat{A} d^2z$. Then, if $\hat{A} = A_C(z)$ is a classical observable and $\hat{A}_Q$ is a quantum observable (independent of $z$), their Ehrenfest equations are

$$ \frac{d}{dt} \langle A_C \rangle_{\gamma} = \langle \langle A_C, \hat{H} \rangle \rangle_{\gamma}, \quad \text{and} \quad i\hbar \frac{d}{dt} \langle \hat{A}_Q \rangle_{\gamma} = \langle \langle \hat{A}_Q, \hat{H} \rangle \rangle_{\gamma}. $$

While no simple Ehrenfest equation is available for general hybrid observables, the relations above provide already much insight. For example, it is easy to see that hybrid Hamiltonians of the type $\hat{H} = m^{-1}\hat{p}^2/2 + M^{-1}p^2/2 + \hat{V}(q - \hat{x})$ yield conservation of that total momentum, that is $d\langle p + \hat{p} \rangle_{\gamma}/dt = \langle \langle p, \hat{V} \rangle \rangle - i\hbar^{-1}[\hat{p}, \hat{V}]_{\gamma} = 0$.

### 7.4. Specializations and further comments

This section is devoted to presenting relevant special cases of the NHWE model in quantum-classical dynamics. As we already discussed after Eq. (6.10), the purely quantum and purely classical cases are immediately recovered by the present treatment. Here, we will show how simple quantum-classical models are also naturally recovered. Eventually, we will conclude this section by making considerations on the appearance of the quantum backreaction on the classical Lagrangian paths.

**Mean-field model.** If $\gamma(z, x) = \psi(x)\sqrt{D(z)}e^{iS(z)/\hbar}$ (with $\psi$ independent of $z$), then the NHWE Hamiltonian (6.4) reduces to $h = \int D(\hat{H}) d^2z$ and $\mathcal{A}_B = \nabla S$. Thus, the term $\langle \gamma, \mathcal{A}_B \cdot \mathcal{X} \gamma \rangle$ reduces to a pure time derivative and the Lagrangian (6.3) drops to $\ell_{MF} = \int D(\mathcal{A} \cdot \mathcal{X} + \langle \psi, i\hbar\dot{\psi} - \hat{H} \rangle) d^2z$, whose associated variational principle returns the mean-field equations (1.1) for $\rho_c = D$ and $\hat{\rho}_q = \psi\psi^{\dagger}$.

**Ehrenfest model.** If we assume $\tilde{\mathbf{X}}_{\hat{H}} \simeq 0$ and discard the corresponding term in the NHWE Hamiltonian (6.4), then the latter reduces to $h(\gamma) = \langle \hat{H} \rangle_{\gamma}$. In this case, we have $\mathcal{X} = \langle \mathbf{X}_{\hat{H}} \rangle$ and the NHWE (6.7) reduces to

$$ i\hbar \frac{\partial \gamma}{\partial t} + i\hbar \langle \mathbf{X}_{\hat{H}} \rangle \cdot \nabla \gamma + \frac{i\hbar}{2} \gamma \mathrm{div} \langle \mathbf{X}_{\hat{H}} \rangle = \hat{H}\gamma - \langle p \partial_p \hat{H} \rangle \gamma, \qquad (7.10) $$

In terms of the measure-valued density matrix $\hat{P} = \gamma\gamma^{\dagger}$, the hybrid wave equation (7.10) reads $i\hbar\partial_t\hat{P} + i\hbar\mathrm{div}(\hat{P}\langle\mathbf{X}_{\hat{H}}\rangle) = [\hat{H}, \hat{P}]$. Recently, this model was also obtained as a closure model directly from the original theory in Section 4; see [72]. The equations for the classical density $\rho_c = \mathrm{Tr}\hat{P}$ and the quantum density matrix $\hat{\rho}_q = \int \hat{P}d^2z$ are

$$ \frac{\partial \rho_c}{\partial t} + \mathrm{div}(\rho_c \langle \mathbf{X}_{\hat{H}} \rangle) = 0, \qquad i\hbar \frac{d\hat{\rho}_q}{dt} = \int [\hat{H}, \hat{P}] d^2z, \qquad (7.11) $$

respectively. We notice that, unlike the mean-field model (1.1), the present construction allows to capture some decoherence effects. Indeed, the second in (7.11) implies that the quantum purity $\|\hat{\rho}_q\|^2 = \mathrm{Tr}(\hat{\rho}_q^2)$ possesses nontrivial dynamics. Similar conclusions were
reached previously in [73], although in most cases the decoherence levels obtained by this model – usually known as the Ehrenfest model – are considerably lower than those obtained from fully quantum predictions. Another limitation of Ehrenfest dynamics arises in the context of pure-dephasing systems, that is hybrid systems involving a quantum two-level subsystem and whose hybrid Hamiltonian is $\hat{H} = H_c + H_I\hat{\sigma}_k$, where $\hat{\sigma}_k$ is any of the three Pauli matrices. In this case, Ehrenfest dynamics gives $\partial_t\langle\hat{\sigma}_k\rangle + \langle\mathbf{X}_{\hat{H}}\rangle \cdot \nabla\langle\hat{\sigma}_k\rangle = 0$ with $\langle\mathbf{X}_{\hat{H}}\rangle = \mathbf{X}_{H_c} + \mathbf{X}_{H_I}\langle\hat{\sigma}_k\rangle$. Then, if $\langle\hat{\sigma}_k\rangle = 0$ initially, this relation holds indefinitely and the classical evolution in (7.11) returns simply $\partial_t\rho_c = \{H_c, \rho_c\}$. Thus, we observe that the classical dynamics decouples completely from the quantum motion. This quantum-classical decoupling occurring for pure-dephasing systems is among the main drawbacks of the Ehrenfest model. Instead, this issue is absent in the NHWE, or its density-matrix formulation given by (1.4), with $\hat{H} = \delta\hbar/\delta\hat{P}$ and (7.4). Indeed, in this case one obtains $(\partial_t + \mathbf{X} \cdot \nabla)\langle\hat{\sigma}_k\rangle = D^{-1}\{H_I, D(\langle\hat{\sigma}_k\rangle^2 - 1)\}$, with $\mathbf{X} = \mathbf{X}_{H_c} + \mathbf{X}_{H_I}\langle\hat{\sigma}_k\rangle + \rho_c^{-1}(\langle\rho_c\hat{\mathbf{I}}\rangle, \hat{\sigma}_k) \cdot \nabla\mathbf{X}_{H_I} - \mathbf{X}_{H_c} \cdot \nabla\langle\rho_c\hat{\mathbf{I}}\rangle, \hat{\sigma}_k)$ /2 and $\rho_c\hat{\mathbf{I}} = i\hbar[\hat{\mathbf{P}}, \nabla\hat{\mathbf{P}}]$. See [10] for more discussions on how the closure model proposed here adapts to the case of quantum two-level systems.

**Quantum decoherence and backreaction force.** While quantum decoherence is absent in the mean-field model (1.1) and is poorly captured by the Ehrenfest wave equation (7.10), the full NHWE model (6.7) seems to suggest that decoherence and other quantum-classical correlation effects arise from the fluctuation vector field $\mathbf{X}_{\hat{H}}$, which is given by (6.12) for the standard-type Hamiltonian operator (6.11). Then, the fluctuation force $\tilde{\mathbf{F}} = \hat{F} - \langle\hat{F}\rangle$ around the Hellmann-Feynman average $\langle\hat{F}\rangle$ emerges in the present framework as the main responsible for quantum-classical correlations. Besides quantum decoherence, this Hellmann-Feynman fluctuation $\tilde{\mathbf{F}}$ also comprises the so-called *quantum backreaction*, that is the feedback force on the classical trajectories. While the mean-field vector field $\mathbf{X}_{\langle\hat{H}\rangle}$ occurring in the first equation $\partial_t\rho_c + \text{div}(\rho_c\mathbf{X}_{\langle\hat{H}\rangle}) = 0$ of (1.1) ignores this important effect, part of the backreaction appears in the Ehrenfest vector field $\langle\mathbf{X}_{\hat{H}}\rangle = \mathbf{X}_{\langle\hat{H}\rangle} + (\langle\mathbf{X}_{\hat{H}}\rangle - \mathbf{X}_{\langle\hat{H}\rangle})$ from the first equation in (7.11). Unfortunately, the Ehrenfest model is known to lack most backreaction effects [21]. In the more general case of Eq. (6.10), however, we observe that the rest of the backreaction effects is retained by the difference $\mathbf{X} - \langle\mathbf{X}_{\hat{H}}\rangle$, which for the Hamiltonian (6.11) is triggered entirely by the fluctuation force in $\mathbf{X}_{\hat{H}} - \langle\mathbf{X}_{\hat{H}}\rangle = (0, \tilde{\mathbf{F}})$. Then, the quantity $\tilde{\mathbf{F}}$ emerges as a fundamental object in quantum-classical coupling and we shall call it the *backreaction force*.

## 8. Summary and outlook

Based on Sudarshan's early attempt in [7], this paper has formulated a dynamical theory of hybrid quantum-classical systems by following two sequential steps. In the first, we followed our original work by blending the variational structure of Koopman and Schrödinger wavefunctions after exploiting the KvH construction to retain classical phases within the formalism. In the second, we took a step beyond our original model by making classical phases unobservable via a gauge-invariance principle. As we showed, the second step is crucial in ensuring that the classical state is always identified by a positive-definite Liouville density on phase-space, whose transport equation is now directly related to Noether's theorem. While classical phases were made unobservable by enforcing superselection rules following Sudarshan's work, the latter led to substantial issues which are here circumvented by resorting to a simpler gauge principle. The idea of a gauge principle reflects previous comments made by Ghose [52], and much earlier by Boucher and Traschen [5], on the role of classical phases in hybrid dynamics.

As we showed, the introduction of a gauge principle in the original variational principle requires nontrivial sophisticated steps in geometric mechanics combining tools from continuum dynamics, such as Lagrangian paths and Lie transport, with Hilbert-space methods from quantum and Koopman classical mechanics, e.g. phase-space wavefunctions and local density matrices. In particular, a key ingredient is the exact factorization method allowing the identification of the classical phase and ultimately leading to the application of the gauge principle.

The final NHWE model comprised by (6.7) and (6.10) consists of a formidable equation for a hybrid QC wavefunction. Unlike the original formulation in (4.2), the NHWE equation is made highly nonlinear by the presence of the Berry connection in the Lagrangian (6.3). Yet, in this context the Berry connection is crucial in realizing the QC correlations. Given the level of complexity of the NHWE model, trajectory-based computational tools analogue to those in [65,67] are currently under development in order to account for several degrees of freedom and strong subsystem-bath couplings. Also, further closure models are desirable to alleviate the difficulties arising from a full phase-space treatment. At present, fluid closure models are under investigation.

Despite its complexity, the NHWE model was shown to retain several remarkable features that are not easy to find in alternative hybrid theories. For example, the model possesses a hybrid Poincaré integral invariant which is unavailable in the original formulation. As showed, this leads to defining a time-dependent hybrid symplectic form whose associated Liouville volume can be used to construct dynamical invariants. Indeed, the latter are allowed in the form of Casimir functionals by the presence of a noncanonical Hamiltonian structure that is written in (7.6) as the sum of a quantum Schrödinger term and a classical Koopman term. In turn, the existence of Casimir invariants allows the identification of Lyapunov-stable equilibria by following the methods in [74]. In addition, as discussed in the previous section, the NHWE overcomes the quantum-classical decoupling occurring in pure-dephasing Ehrenfest dynamics, thereby ensuring the persistence of a quantum backreaction on the classical evolution.

We conclude by emphasizing two more different directions arising from this work. At a foundational level, we would like to understand how the results in this paper can be used to formulate a dynamical theory of quantum measurement. As pointed out by Peres [75], the latter can be envisioned as the combination of a reversible *pre-measurement* mechanism and an irreversible process. The latter is due to the fact that the classical apparatus is realistically in a thermal equilibrium. This description requires the addition of entropy sources in the classical sector, thereby indicating a possible way forward in the formulation of a measurement theory.

On a more applied level, we would also like to apply the present theory to devise new models for quantum-classical spin hybrids [24]. The general idea is to reach previously unaccessible regions in the control parameter space of spin control protocols by exploiting the coupling to large dipole moments of classical magnets. This picture involves the interaction dynamics of a classical magnetization vector with quantum spins, in which case the former obeys Landau-Lifshitz-type dynamics. We plan to pursue this interesting direction in future work.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Cesare Tronci reports financial support was provided by John Templeton Foundation. Cesare Tronci reports financial support was provided by The Royal Society.
Acknowledgements

We are grateful to Andrés Bermúdez, Francesco Di Maiolo, Darryl Holm, Ilon Joseph, Phil Morrison, and Robert MacKay for their keen remarks during the development of this work. These results were 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. Partial support by the Royal Society, United Kingdom Grant IES\R3\203005 is also greatly acknowledged.

Appendix A. Von Neumann operators in Koopman mechanics

This appendix presents an extension of the treatment in Section 3 that overcomes the difficulties arising from the use of the relation $\nabla S = \mathcal{A}$. As discussed in various instances, this relation requires a singular phase which may pose several problems. Here, we will overcome this difficulty by extending the wavefunction treatment from Section 3 to consider the more general setting made available by von Neumann operators and their Wigner transforms.

Let us start from KvH theory. Instead of considering the KvH equation (2.6) for the wavefunction $\chi(\mathbf{z}, t)$, here we consider the Koopman mixture

$$ \hat{\Theta}(\mathbf{z}, \mathbf{z}') = \sum_{a=1}^{N} w_a \chi_a(\mathbf{z}) \chi_a^*(\mathbf{z}'), \qquad (A.1) $$

with $w_a \ge 0$ and $\sum_{a=1}^N w_a = 1$. The von Neumann operator $\hat{\Theta}$ mimics the density matrix from standard quantum theory and satisfies the equation

$$ i\hbar \partial_t \hat{\Theta} = [\hat{\mathcal{L}}_H, \hat{\Theta}], \quad \text{where} \quad \hat{\mathcal{L}}_H := i\hbar\{H, \}\} + \mathcal{L} \qquad (A.2) $$

is the covariant Liouvillian from (2.6). Then, the variational structure for the Koopman mixture is given by the natural extension of the KvH Lagrangian $L_{KvH}$ in (2.7), that is

$$ L_{\text{mix}} = \sum_{a=1}^{N} w_a L_{KvH}(\chi_a, \partial_t \chi_a), \qquad (A.3) $$

where the subscript 'mix' refers to the fact that we are dealing with a mixture. This Lagrangian yields a sequence of identical KvH equations for each $\chi_a$, with the same Hamiltonian $H$. Then, we can write $\chi_a = \sqrt{D_a} e^{iS_a/\hbar}$ and perform exactly the same steps in the previous section. For example, the Lagrange-to-Euler map (3.3) now reads $D_a(\mathbf{z}, t) = \int D_{a,0}(\mathbf{z}) \delta(\mathbf{z} - \eta(\mathbf{z}_0, t)) d^2\mathbf{z}_0$, or equivalently in terms of the Jacobian $J_\eta = \det \nabla \eta$. We note that, since the $\chi_a$'s all obey the same KvH equation, the Lagrangian path $\eta$ is the same for each $a = 1, \dots, N$. Upon defining

$$ D = \sum_{a=1}^{N} w_a D_a, \quad \text{and} \quad \sigma = \sum_{a=1}^{N} w_a D_a \nabla S_a, \qquad (A.4) $$

the classical Liouville density reads

$$ \rho_c = D + \text{div}(J(\sigma - D\mathcal{A})) \qquad (A.5) $$

and the Lagrangian (A.3) is taken into the Euler-Poincaré form

$$ L_{EPmix} = \int (\boldsymbol{\sigma} \cdot \boldsymbol{\chi} + (D\mathcal{A} - \boldsymbol{\sigma}) \cdot J\nabla H - DH) d^2\mathbf{z}, \qquad (A.6) $$

which replaces the previous expression in (3.5). It may be worth noticing that the quantities (A.4) may also be written entirely in terms of the operator $\hat{\Theta}$ as [65]

$$ D(\mathbf{z}) = \hat{\Theta}(\mathbf{z}, \mathbf{z}')|_{\mathbf{z}'=\mathbf{z}}, \quad \text{and} \quad \sigma(\mathbf{z}) = \frac{1}{2}[\hat{\lambda}, \hat{\Theta}]_{+}(\mathbf{z}, \mathbf{z}')|_{\mathbf{z}'=\mathbf{z}}, \qquad (A.7) $$

where $[\cdot, \cdot]_+$ denotes the anticommutator and we have introduced the operator $\hat{\lambda} = -i\hbar\nabla$, so that $[\hat{\mathbf{z}}, \hat{\lambda}] = i\hbar\mathbf{1}$. Given the appearance of the above quantities in (A.5), here we remark that the KvH formulation generally allows for a sign-indefinite operator $\hat{\Theta}$, as long as the corresponding classical distribution $\rho_c$ is positive. Indeed, as we will see below, it is convenient to consider von Neumann operators beyond the mixture type (A.1), which however represents a useful tool to obtain the Lagrangian $L_{EPmix}$.

In the variational formulation associated to (A.6), the variation $\delta\sigma$ is arbitrary and one makes use of the variations (3.6)-(3.7), along with (3.4). Equivalently, one may also write $\sigma = D\mathcal{A}$ with arbitrary $\delta\mathcal{A}$. Then, Hamilton's principle $\delta \int_{t_1}^{t_2} L_{EPmix} dt = 0$ yields the analogous of (3.1), that is, upon writing $\mathcal{A} := \sigma/D$,

$$ \partial_t(\mathcal{A} - \mathcal{A}) + \mathbf{X}_H \cdot \nabla(\mathcal{A} - \mathcal{A}) + \nabla \mathbf{X}_H \cdot (\mathcal{A} - \mathcal{A}) = 0. \qquad (A.8) $$

We realize that, since $\mathcal{A}$ is no longer an exact differential, the preserved condition $\mathcal{A} = \mathcal{A}$ may be set at the initial time without incurring in topological singularities. In this way, the classical density (A.5) collapses to $\rho_c = D$, which again satisfies $\partial_t D = \{H, D\}$. But what does the relation $\mathcal{A} = \mathcal{A}$ correspond to in terms of the von Neumann operator $\hat{\Theta}$ in (A.2)? Can we obtain the relation $\mathcal{A} = \mathcal{A}$ as a closure model expressing $\sigma = D\mathcal{A}$ as a function $\sigma = \sigma(D, \mathcal{A})$? In order to address these points, it is convenient to use Wigner functions.

Let us consider the Wigner transform of Eq. (A.2), that is

$$ \partial_t W(\mathbf{z}, \mathbf{\lambda}, t) = \langle [\mathcal{L}_H(\mathbf{z}, \mathbf{\lambda}), W(\mathbf{z}, \mathbf{\lambda}, t)] \rangle, \quad \text{where} \quad \mathcal{L}_H(\mathbf{z}, \mathbf{\lambda}) = \mathbf{X}_H(\mathbf{z}) \cdot \mathbf{\lambda} - \mathcal{L}(\mathbf{z}) \qquad (A.9) $$
is the Weyl symbol of the covariant Liouvillian in (A.2), $W(\mathbf{z}, \lambda, t)$ is the Wigner transform of $\hat{\Theta}(\mathbf{z}, \mathbf{z}', t)$, and $\{\cdot, \cdot\}$ is the Moyal bracket in the double phase-space coordinates $(\mathbf{z}, \lambda)$. We notice that, since $L_H(\mathbf{z}, \lambda)$ is linear in $\lambda$, the moment hierarchy $\int \lambda^n W d^2z$ closes at any order greater than or equal to 1. This can be seen explicitly by taking the $L^2$-pairing of (A.9) with the arbitrary analytic function $\lambda^n \Phi_n(\mathbf{z})$, and then using the properties of the Moyal bracket. Thus, we may consider the closed dynamics for the variables given by $n=0, 1$, that is $D = \int W d^2z$ and $\sigma = \int \lambda W d^2z$, respectively. In particular, we notice that any local Gaussian function of the type $W(\mathbf{z}, \lambda, t) = (2\pi \Sigma^2)^{-1} D(\mathbf{z}, t) \exp(-\Sigma^{-2}|\lambda - \mathcal{A}(\mathbf{z}, t)|^2/2)$ is an exact solution of (A.9) for any $\Sigma$, provided that $\mathcal{A} = \sigma/D$ satisfies (A.8) and $D$ satisfies the classical Liouville equation. Thus, setting $\mathcal{A} = \mathcal{A}$ we obtain that the Wigner function

$$W(\mathbf{z}, \lambda, t) = \frac{D(\mathbf{z}, t)}{2\pi \Sigma^2} e^{-\frac{|\lambda - \mathcal{A}(\mathbf{z})|^2}{2\Sigma^2}} = \frac{D(\mathbf{z}, t)}{2\pi \Sigma^2} e^{-\frac{(\lambda - \mu)^2 + \lambda_p^2}{2\Sigma^2}}, \quad (A.10)$$

is an exact solution of (A.9) provided $D$ satisfies $\partial_t D = \{H, D\}$. The second equality in (A.10) follows immediately by recalling $\mathcal{A}(\mathbf{z}) = (p, 0)$. The von Neumann operator associated to the Wigner function (A.10) is found as

$$\hat{\Theta}(\mathbf{z}, \mathbf{z}', t) = \int W\left(\frac{\mathbf{z} + \mathbf{z}'}{2}, \lambda, t\right) e^{-i\lambda \cdot (\mathbf{z} - \mathbf{z}') d^2 \lambda} = D\left(\frac{\mathbf{z} + \mathbf{z}'}{2}, t\right) e^{-\frac{1}{2\hbar} \left(i(q-q')(p+p') + \frac{\Sigma^2}{\hbar} |\mathbf{z}-\mathbf{z}'|^2\right)}. \quad (A.11)$$

In the limit case $\Sigma \to 0$, the von Neumann operator associated to the singular Wigner function $W = D\delta(\lambda - \mathcal{A})$ was presented explicitly in [76].

Here, the Wigner function (A.10) appears as a closure for the von Neumann equation (A.2), or its Wigner-Moyal correspondent in (A.9). Actually, the function (A.10) is an exact solution of Eq. (A.9), which emerges as the natural extension of KvH theory to Wigner distributions. We notice that, as mentioned above, the operator (A.11) is sign-indefinite and thus is more general than the mixture (A.1). Nonetheless, this unsigned operator yields a positive distribution (A.5), so that the resulting closure model coincides with the classical Liouville equation for $\rho_c = D$ (or KvN equation, by writing $D = |\chi|^2$).

Equivalently, the classical Liouville equation may be obtained by simply replacing $\sigma \to D\mathcal{A}$ in the Euler-Poincaré variational principle (A.6) for Koopman mixtures. The latter observation allows us to avoid dealing with von Neumann operators when looking for a reduction process taking the KvH construction to the simpler KvN dynamics.

## Appendix B. Calculations for the diamond operator

From (6.6), we compute the diamond $\Upsilon \diamond \delta h / \delta \Upsilon$ term by term as follows

$$\begin{align*}
\langle \hat{H} \Upsilon + \mathcal{A}_B \cdot (\mathbf{X}_{\hat{H}}) \Upsilon, \nabla \Upsilon \rangle - \langle \Upsilon, \nabla(\hat{H} \Upsilon + \mathcal{A}_B \cdot (\mathbf{X}_{\hat{H}}) \Upsilon) \rangle &= -\|\Upsilon\|^2 \langle \nabla \hat{H} \rangle - \|\Upsilon\|^2 \nabla (\mathcal{A}_B \cdot (\mathbf{X}_{\hat{H}})) \\
\langle \mathbf{X}_{\hat{H}} \cdot (i\hbar\nabla)\Upsilon, \nabla \Upsilon \rangle - \langle \Upsilon, \nabla(\mathbf{X}_{\hat{H}} \cdot (i\hbar\nabla)\Upsilon) \rangle &= 2\langle \mathbf{X}_{\hat{H}} \cdot (i\hbar\nabla)\Upsilon, \nabla \Upsilon \rangle - \nabla\langle \Upsilon, \mathbf{X}_{\hat{H}} \cdot (i\hbar\nabla)\Upsilon \rangle \\
&= 2\langle \mathbf{X}_{\hat{H}} \cdot (i\hbar\nabla)\Upsilon, \nabla \Upsilon \rangle + \nabla\langle \|\Upsilon\|^2 \mathcal{A}_B \cdot (\mathbf{X}_{\hat{H}}) \rangle - \nabla\langle \Upsilon, (\mathbf{X}_{\hat{H}} - \langle \mathbf{X}_{\hat{H}} \rangle) \cdot (i\hbar\nabla)\Upsilon \rangle \\
\langle \mathbf{X}_{\hat{H}} \cdot \mathcal{A}_B \Upsilon, \nabla \Upsilon \rangle - \langle \Upsilon, \nabla(\mathbf{X}_{\hat{H}} \cdot \mathcal{A}_B \Upsilon) \rangle &= \langle \mathbf{X}_{\hat{H}} \cdot \mathcal{A}_B \Upsilon, \nabla \Upsilon \rangle - \|\Upsilon\|^2 \langle \nabla(\mathbf{X}_{\hat{H}} \cdot \mathcal{A}_B) \rangle - \langle \Upsilon, \mathbf{X}_{\hat{H}} \cdot \mathcal{A}_B \nabla \Upsilon \rangle \\
&= - \|\Upsilon\|^2 \langle \nabla(\mathbf{X}_{\hat{H}} \cdot \mathcal{A}_B) \rangle
\end{align*}$$

Also if we apply

$$\frac{1}{2} \left( \left\langle \frac{\delta h}{\delta \Upsilon}, \nabla \Upsilon \right\rangle - \left\langle \Upsilon, \nabla \frac{\delta h}{\delta \Upsilon} \right\rangle \right) = \left\langle \frac{\delta h}{\delta \Upsilon}, \nabla \Upsilon \right\rangle - \frac{1}{2} \nabla \left\langle \frac{\delta h}{\delta \Upsilon}, \Upsilon \right\rangle,$$

we have

$$\frac{1}{2} \left( \langle i\hbar\Upsilon \operatorname{div}(\mathbf{X}_{\hat{H}}), \nabla\Upsilon \rangle + \langle i\hbar\Upsilon, \operatorname{div}(\mathbf{X}_{\hat{H}}) \nabla\Upsilon + \Upsilon\nabla \operatorname{div}(\mathbf{X}_{\hat{H}}) \rangle \right) = \|\Upsilon\|^2 \mathcal{A}_B \operatorname{div}(\mathbf{X}_{\hat{H}}) \\
= 2\mathcal{A}_B \langle \Upsilon, (\mathbf{X}_{\hat{H}} - \langle \mathbf{X}_{\hat{H}} \rangle) \cdot \nabla\Upsilon \rangle$$

and

$$2\langle i\hbar(\mathbf{X}_{\hat{H}}) \cdot \nabla\Upsilon, \nabla\Upsilon \rangle - \langle\Upsilon, (i\hbar(\mathbf{X}_{\hat{H}}) \cdot \nabla\Upsilon)\rangle = 2\langle(\mathbf{X}_{\hat{H}}) \cdot (i\hbar\nabla)\Upsilon, \otimes\nabla\Upsilon\rangle + \nabla\langle\|\Upsilon\|^2 \mathcal{A}_B \cdot (\mathbf{X}_{\hat{H}})\rangle,$$

which altogether recovers (6.10).

## Appendix C. Closure method and hybrid von Neumann operators

In Section 6, we showed how the relation (5.6) is used to achieve a positive-definite classical density while taking the Lagrangian (5.16) into the manifestly gauge-invariant form (6.1). In this appendix, we will show how the resulting gauge-invariant Lagrangian $\ell_{QC}$ in (6.3) emerges in the context of the original theory from Section 4 by resorting to von Neumann operators and their Wigner functions.

Following the discussion in Appendix A, we prefer to partly abandon the wavefunction picture of Koopman dynamics in favour of the mixture picture, which was mentioned in [10] only briefly. The reason for resorting to mixtures is that $\mathcal{A}_B$ often vanishes at the initial time (for example, when $\psi$ is real) and thus the relation (5.6) would require the presence of topological singularities. In more generality, we wish to understand how the hybrid Lagrangian in (6.3) emerges as a type of closure model by adopting an ansatz, similarly to the arguments in Appendix A. In particular, in analogy to (A.1), we combine mixtures with the exact factorization (5.1) by constructing a von Neumann operator with matrix elements

$$\hat{\Theta}(\mathbf{z}, \mathbf{z}', x, x') = \sum_{a=1}^{N} w_a \gamma_a(\mathbf{z}, x) \gamma_a^*(\mathbf{z}', x'), \quad \text{with} \quad \gamma_a(\mathbf{z}, x) = \chi_a(\mathbf{z}) \psi(x; \mathbf{z}), \quad (C.1)$$
and $\int |\psi|^2 dx = 1$. Mixtures of this type are analogous to those considered in [65] within the context of nonadiabatic quantum hydrodynamics. Then, following the same arguments as in Appendix A, we will proceed by inserting the exact factorization $\gamma_a = \chi_a \psi$ in the mixture Lagrangian $\sum_a w_a L_{QC}(\gamma_a, \partial_t \gamma_a)$, where $L_{QC}$ is given in (4.3). After writing $\chi_a = \sqrt{D_a}e^{iS_a/\hbar}$ and using the definitions (A.4), the steps from Section 5 take us to the following Euler-Poincaré Lagrangian (and thus to its associated dynamics):

$$L = \int \left( \boldsymbol{\sigma} \cdot \boldsymbol{X} + D \langle \psi, i\hbar \partial_t \psi - \hat{H} \psi + i\hbar \tilde{X}_{\hat{H}} \cdot \nabla \psi - (D^{-1}\boldsymbol{\sigma} + \mathcal{A}_B - \mathcal{A}) \cdot \mathbf{X}_{\hat{H}} \psi \rangle \right) d^2z. \quad (C.2)$$

In this case, the original expression (5.4) of the classical density changes to $\rho_c = D + \text{div}(\boldsymbol{J}\boldsymbol{\sigma} + \boldsymbol{D}\boldsymbol{J}\mathcal{A}_B - \boldsymbol{D}\boldsymbol{J}\mathcal{A})$ and one would like to set $D^{-1}\boldsymbol{\sigma} + \mathcal{A}_B = \mathcal{A}$ in such a way to obtain a positive-definite density $\rho_c = D$. Unlike the purely classical case treated in Appendix A, the hybrid relation $D^{-1}\boldsymbol{\sigma} + \mathcal{A}_B = \mathcal{A}$ is not preserved in time by the dynamics arising from (C.2). Here, we will devise a closure scheme so that a particular ansatz at the level of the von Neumann operator is applied to the Lagrangian (C.2) thereby returning (6.1).

We realize that the operator (C.1) emerges as a special case of a wider class of operators, which will be particularly convenient in designing our closure scheme. Specifically, let us consider matrix elements of the type

$$\hat{\Theta}(\mathbf{z}, \mathbf{z}', x, x') = \hat{\theta}(\mathbf{z}, \mathbf{z}')\psi(x; \mathbf{z})\psi^*(x'; \mathbf{z}'), \quad (C.3)$$

where $\hat{\theta}$ is itself a von Neumann operator acting only on the classical sector. We observe that (C.3) leads to the relations (A.7) with $\hat{\Theta}$ replaced by $\hat{\theta}$. In this setting, the Hamiltonian functional $\hat{h}(\hat{\theta}, \psi)$ is given by (5.8) upon replacing $D\nabla S \to \boldsymbol{\sigma}$ and using (A.7) to express $D$ and $\boldsymbol{\sigma}$ in terms of $\hat{\theta}$. Thus, if we follow the procedure in Appendix A and define $\boldsymbol{\Lambda} := D^{-1}\boldsymbol{\sigma}$, finding a closure for which $\boldsymbol{\Lambda} + \mathcal{A}_B = \mathcal{A}$ amounts to finding a von Neumann operator $\hat{\theta}$ so that $[\hat{\boldsymbol{\lambda}}, \hat{\boldsymbol{\theta}}]_+(\mathbf{z}, \mathbf{z}')|_{\mathbf{z}'=\mathbf{z}} = 2\hat{\theta}(\mathbf{z}, \mathbf{z}')|_{\mathbf{z}'=\mathbf{z}}(\mathcal{A} - \mathcal{A}_B)$. We will address this problem by resorting to Wigner transforms, following the treatment in Appendix A. We anticipate that, similarly to the approach in that section, we will now allow for the operator $\hat{\theta}$ to be sign-indefinite, as long as its associated phase-space distribution $\rho_c = D + \text{div}(\boldsymbol{J}(\boldsymbol{\sigma} + \boldsymbol{D}\mathcal{A}_B - \boldsymbol{D}\mathcal{A}))$ remains positive.

Let us introduce the Wigner transform

$$w(\mathbf{z}, \boldsymbol{\lambda}) = \frac{1}{(\hbar\pi)^2} \int \hat{\theta}(\mathbf{z} - \mathbf{z}', \mathbf{z} + \mathbf{z}') e^{2i\boldsymbol{\lambda} \cdot \mathbf{z}' / \hbar} d^2\mathbf{z}' \quad (C.4)$$

of the von Neumann operator $\hat{\theta}$, so that

$$D(\mathbf{z}) = \int w(\mathbf{z}, \boldsymbol{\lambda}) d^2\boldsymbol{\lambda}, \qquad D(\mathbf{z})\mathcal{A}(\mathbf{z}) = \int \boldsymbol{\lambda} w(\mathbf{z}, \boldsymbol{\lambda}) d^2\boldsymbol{\lambda}. \quad (C.5)$$

Then, we see that the relation $\mathcal{A} + \mathcal{A}_B = \mathcal{A}$ is achieved by a slight modification of the closure (A.10), that is

$$w(\mathbf{z}, \boldsymbol{\lambda}) = \frac{D(\mathbf{z})}{2\pi \Sigma^2} \exp\left(-\frac{|\boldsymbol{\lambda} - \mathcal{A}(\mathbf{z}) + \mathcal{A}_B(\mathbf{z})|^2}{2\Sigma^2}\right). \quad (C.6)$$

Here, we recall that the Berry connection $\mathcal{A}_B$ is expressed in terms of $\psi$ by (5.5). Thus, we are left with the conclusion that the closure relation $\mathcal{A} = \mathcal{A} - \mathcal{A}_B$ taking the Lagrangian (5.16) into the closure model Lagrangian (6.1) is achieved by a von Neumann operator of the type (C.3), where $\hat{\theta}$ has a Wigner transform given by (C.6) for some $\Sigma$. Similarly to (A.11), we compute

$$\hat{\theta}(\mathbf{z}, \mathbf{z}') = D(\boldsymbol{\zeta}) \exp\left(-\frac{i}{2\hbar}(\mathbf{z} - \mathbf{z}') \cdot (\mathcal{A}(\boldsymbol{\zeta}) - \mathcal{A}_B(\boldsymbol{\zeta})) - \frac{\Sigma^2}{2\hbar^2}|\mathbf{z} - \mathbf{z}'|^2\right) \Big|_{\boldsymbol{\zeta}=(\mathbf{z}+\mathbf{z}')/2}, \quad (C.7)$$

where we recall $\mathcal{A}(\mathbf{z}) = (p, 0)$. In the limit $\Sigma \to 0$, this type of von Neumann operator was previously found in [65], although the expression (C.7) now involves phase-space coordinates. By following analogous steps to those in [65] (see Section 5.1 therein), one shows that replacing (C.7) in the variational principle underlying the dynamics of the hybrid von Neumann operator in (C.3) ultimately leads to the Lagrangian (6.3), thereby recovering the NHWE system (6.7)–(6.10). The explicit relations involved in this argument are omitted here and we refer to Section 5.1 in [65] for the technical details in the fully quantum context.

We notice that the Wigner function (C.6) extends the usual fluid closure method from kinetic theory to consider Maxwellian distributions on phase-space. Although their von Neumann operator is sign-indefinite, as we pointed out in Appendix A, Wigner distributions of general Maxwellian type are occasionally used in mixed-state quantum hydrodynamics [77,78]. In the present context, however, we will not have to deal with the various aspects concerning Maxwellian distributions. Indeed, as we already showed, our NHWE model is simply obtained by replacing $\nabla S \to \mathcal{A} - \mathcal{A}_B$ in the Lagrangian (5.16) from the original theory.

## References

[1] L.D. Landau, E.M. Lifshitz, Quantum Mechanics, in: Non-Relativistic Theory, Elsevier Butterworth-Heinemann, Oxford, UK, 1977.

[2] N. Bohr, Can quantum-mechanical description of physical reality be considered complete? Phys. Rev. 48 (1935) 696–702.

[3] D.I. Bondar, F. Gay-Balmaz, C. Tronci, Koopman wavefunctions and classical-quantum correlation dynamics, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 475 (2229) (2019) 20180879.

[4] D. Fang, S. Jin, C. Sparber, An efficient time-splitting method for the Ehrenfest dynamics, Multiscale Model. Simul. 16 (2018) 900–921.

[5] W. Boucher, J. Traschen, Semiclassical physics and quantum fluctuations, Phys. Rev. D 37 (1988) 3522–3532.

[6] O.V. Prezhdo, V.V. Kisil, Mixing quantum and classical mechanics, Phys. Rev. A 56 (1997) 162–175.

[7] E.C.G. Sudarshan, Interaction between classical and quantum systems and the measurement of quantum observables, Prāmaṇa 6 (3) (1976) 117–126.

[8] A. Widom, Y. Srivastava, Lagrangian formulation of Bohr’s measurement theory, Nuovo Cimento B 107 (1) (1992) 71–75.

[9] M.J.W. Hall, M. Reginatto, Ensembles on Configuration Space, Springer, 2016.

[10] F. Gay-Balmaz, C. Tronci, Koopman wavefunctions and classical states in hybrid quantum-classical dynamics, J. Geom. Mech. submitted for publication, arXiv:2108.01482.

[11] C. Marletto, V. Vedral, Why we need to quantise everything, including gravity, Npj Quant. Inf. 29 (2017).

[12] L.L. Salcedo, Absence of classical and quantum mixing, Phys. Rev. A 54 (4) (1996) 3657–3660.

[13] D.R. Terno, Inconsistency of quantum-classical dynamics, and what it implies, Found. Phys. 36 (1) (2006) 102–111.
[14] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48 (1976) 119.

[15] M. Schlosshauer, K. Camilleri, What classicality? Decoherence and Bohr's classical concepts, AIP Conf. Proc. 1327 (2011) 26-35.

[16] R. Cabrera, D.I. Bondar, K. Jacobs, H.A. Rabitz, Efficient method to generate time evolution of the Wigner function for open quantum systems, Phys. Rev. A 92 (2015) 042122.

[17] M. Albers, C. Kiefer, M. Reginatto, Measurement analysis and quantum gravity, Phys. Rev. D 78 (2008) 064051.

[18] M. Bojowald, D. Ding, Canonical description of cosmological backreaction, J. Cosmol. Astropart. Phys. 3 (2021) 083.

[19] R. Crespo-Otero, M. Barbatti, Recent advances and perspectives on nonadiabatic mixed quantum-classical dynamics, Chem. Rev. 118 (15) (2018) 7026-7068.

[20] I. Tavernelli, B.F.E. Curchod, U. Rothlisberger, Mixed quantum-classical dynamics with time-dependent external fields: A time-dependent density-functional-theory approach, Phys. Rev. A 81 (2010) 052508.

[21] J.C. Tully, Mixed quantum-classical dynamics, Faraday Discuss. 110 (1998) 407-419.

[22] N. Ben Abdallaha, M.J. Cáceres, J.A. Carrillo, F. Vecil, A deterministic solver for a hybrid quantum-classical transport model in nanoMOSFETs, J. Comput. Phys. 228 (2009) 6553-6571.

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

[24] A. Rustagi, S. Kajale, P. Upadhyaya, Manipulating quantum impurity spins via dynamical modes of nanomagnets. arXiv:2203.03652.

[25] M. Born, R. Oppenheimer, Zur quantentheorie der molekeln, Ann. Phys. 389 (20) (1927) 457-484.

[26] R.S. MacKay, Slow manifolds, in: T. Dauxois, R.S. MacKay (Eds.), Energy Localisation and Transfer, in: Adv. Ser. Nonlin. Dyn., vol. 22, 2004, pp. 149-192.

[27] R. Requist, C.R. Proetto, E.K.U. Gross, Asymptotic analysis of the Berry curvature in the $E \otimes e$ Jahn-Teller model, Phys. Rev. A 96 (6) (2017) 062503.

[28] J.I. Rawlinson, C. Tronci, Regularized Born-Oppenheimer molecular dynamics, Phys. Rev. A 102 (2020) 032811.

[29] R. Kapral, Surface hopping from the perspective of quantum-classical Liouville dynamics, Chem. Phys. 481 (2016) 77-83.

[30] J.A. Subotnik, W. Ouyang, B.R. Landry, Can we derive Tully's surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence, J. Chem. Phys. 139 (2013) 214107.

[31] I.V. Aleksandrov, The statistical dynamics of a system consisting of a classical and a quantum subsystem, Z. Naturforsch. 36a (1981) 902-908.

[32] V. Gerasimenko, Dynamical equations of quantum-classical systems, Theor. Math. Phys. 50 (1982) 49-55.

[33] L. Diósi, J.J. Halliwell, Coupling classical and quantum variables using continuous quantum measurement theory, Phys. Rev. Lett. 81 (1998) 2846.

[34] D. Chruściński, A. Kossakowski, G. Marmo, E.C.G. Sudarshan, Dynamics of interacting classical and quantum systems, Open Syst. Inf. Dyn. 18 (4) (2011) 339-351.

[35] B.O. Koopman, Hamiltonian systems and transformations in Hilbert space, Proc. Natl. Acad. Sci. 17 (1931) 315.

[36] C. Barceló, R. Carballo-Rubio, L.J. Garay, R. Gómez-Escalante, Hybrid quantum-classical formulations ask for hybrid notions, Phys. Rev. A 86 (2012) 042120.

[37] A. Peres, D.R. Terno, Hybrid quantum-classical dynamics, Phys. Rev. A 63 (2001) 022101.

[38] F. Gay-Balmaz, C. Tronci, Madelung transform and probability densities in hybrid quantum-classical dynamics, Nonlinearity 33 (10) (2019) 5383-5424.

[39] R. Feynman, Negative probability, in: H. Hiley, P. Peat (Eds.), Quantum Implications: Essays in Honour of David Bohm, Routledge & Kegan Paul Ltd., 1987, pp. 235-248.

[40] P.H. Kramer, M. Saraceno, Geometry of the Time-Dependent Variational Principle in Quantum Mechanics, Springer-Verlag, Berlin, Germany, 1981.

[41] E. Bonet Luz, C. Tronci, Geometry and symmetry of quantum and quantum-classical variational principles, J. Math. Phys. 56 (8) (2015) 082104.

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

[43] J. Frenkel, Wave Mechanics; Advanced General Theory, Clarendon Press, Oxford, UK, 1934.

[44] I.Y. Dodin, Geometric view on noneikonal waves, Phys. Lett. A 378 (2014) 1598-1621.

[45] I. Joseph, Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics, Phys. Rev. Res. 2 (4) (2020) 043102.

[46] L. van Hove, On Certain Unitary Representations of an Infinite Group of Transformations (Ph.D. thesis), Word Scientific, 1951, 2001.

[47] B. Kostant, Line bundles and the prequantized Schrödinger equation, in: Colloquium on Group Theoretical Methods in Physics, Centre de Physique Théorique, Marseille, 1972, pp. IV.1-IV.22.

[48] B. Kostant, Quantization and unitary representations, in: Lectures in Modern Analysis and Applications III, in: Lecture Notes in Math, vol. 170, Springer, 1970, pp. 87-208.

[49] F. Faure, Prequantum chaos: Resonances of the prequantum cat map, J. Mod. Dyn. 1 (2) (2007) 255-285.

[50] J.M. Souriau, Quantification géométrique, Comm. Math. Phys. 1 (1966) 374-398.

[51] M.A. de Gosson, On the notion of phase in mechanics, J. Phys. A: Math. Gen. 37 (2004) 7297-7314.

[52] P. Ghose, The unfinished search for wave-particle and classical-quantum harmony, J. Adv. Phys. 4 (3) (2015) 236-251.

[53] E. Madelung, Quantentheorie in hydrodynamischer form, Z. Phys. 40 (3-4) (1927) 322-326.

[54] D.D. Holm, J.E. Marsden, T.S. Ratiu, The Euler-Poincaré equations and semidirect products with applications to continuum theories, Adv. Math. 137 (1998) 1-81.

[55] D.D. Holm, T. Schmah, C. Stoica, Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions, Oxford University Press, Oxford, UK, 2009.

[56] J.E. Marsden, T.S. Ratiu, Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems, Springer, New York, USA, 1999.

[57] A. Sergi, Non-Hamiltonian commutators in quantum mechanics, Phys. Rev. E 72 (2005) 066125.

[58] H.R. Jauslin, D. Sugny, Dynamics of mixed quantum-classical systems, geometric quantization and coherent states, in: Mathematical Horizons for Quantum Physics, in: Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap. vol. 20, World Scientific, 2010, pp. 65-96.

[59] A.D. Bermúdez Manjarres, Projective representation of the Galilei group for classical and quantum-classical systems, J. Phys. A 54 (2021) 444001.

[60] P.R. Chernoff, J.E. Marsden, Some remarks on Hamiltonian systems and quantum mechanics, in: Foundations of Probability Theory, Statistical Inference, and Statistical Theories of Science III, Vol. 6c, Univ. Western Ontario Ser. Philos. Sci., 1976, pp. 35-53.

[61] Shirokov Yu. M., Quantum and classical mechanics in the phase space representation, Sov. J. Part. Nucl. 10 (1979) 1-18.

[62] C. Tronci, Momentum maps for mixed states in quantum and classical mechanics, J. Geom. Mech. 11 (4) (2019) 639-656.

[63] A. Abedi, N.T. Maitra, E.K.U. Gross, Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction, J. Chem. Phys. 137 (22) (2012) 22A530.

[64] I. Bialynicki-Birula, M. Cieplak, J. Karminski, A.M. Furdyna, Theory of Quanta, Oxford University Press, 1992.

[65] M.S. Foskett, D.D. Holm, C. Tronci, Geometry of nonadiabatic quantum hydrodynamics, Acta Appl. Math. 162 (2019) 1-41.

[66] M.S. Foskett, C. Tronci, Holonomy and vortex structures in quantum hydrodynamics, in: A. Fathi, P.J. Morrison, T.M. Seara, S. Tabachnikov (Eds.), Hamiltonian Systems: Dynamics, Analysis, Applications, in: Math. Sci. Res. Inst. Pub., vol. 72, Cambridge University Press, 2022.

[67] D.D. Holm, J.I. Rawlinson, C. Tronci, The bohman method in nonadiabatic quantum hydrodynamics, J. Phys. A 54 (2021) 495201.

[68] A.C. Mead, The geometric phase in molecular systems, Rev. Modern Phys. 64 (1) (1992) 51-85.

[69] S.M. Carroll, J. Lodman, Energy non-conservation in quantum mechanics, Found. Phys. 51 (2021) 83.

[70] R. Feynman, Forces in molecules, Phys. Rev. 56 (4) (1939) 340-343.

[71] O. Prezhdo, A quantum-classical bracket that satisfies the Jacobi identity, J. Chem. Phys. 124 (2006) 201104.

[72] C. Tronci, F. Gay-Balmaz, From quantum hydrodynamics to Koopman wavefunctions II, Lecture Notes in Comput. Sci. 12829 (2021) 311-319.

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

[74] D.D. Holm, J.E. Marsden, T.S. Ratiu, A. Weinstein, Nonlinear stability of fluid and plasma equilibria, Phys. Rep. 123 (1985) 1-116.

[75] A. Peres, Classical interventions in quantum systems, I. the measuring process, Phys. Rev. A 61 (2000) 022116.

[76] F. Gay-Balmaz, C. Tronci, From quantum hydrodynamics to Koopman wavefunctions I, Lecture Notes in Comput. Sci. 12829 (2021) 302-310.

[77] K.H. Hughes, I. Burghardt, Maximum-entropy closure of hydrodynamic moment hierarchies including correlations, J. Chem. Phys. 136 (2012) 214109.

[78] J.V. Lill, M.I. Haftel, G.H. Herling, Mixed state quantum mechanics in hydrodynamical form, J. Chern. Phys. 90 (9) (1989) 4940-4950.