Complex Fluid Models of Mixed Quantum–Classical Dynamics

Several methods in nonadiabatic molecular dynamics are based on Madelung’s hydrodynamic description of nuclear motion, while the electronic component is treated as a finite-dimensional quantum system. In this context, the quantum potential leads to severe computational challenges and one often seeks to neglect its contribution, thereby approximating nuclear motion as classical. The resulting model couples classical hydrodynamics for the nuclei to the quantum motion of the electronic component, leading to the structure of a complex fluid system. This type of mixed quantum–classical fluid models has also appeared in solvation dynamics to describe the coupling between liquid solvents and the quantum solute molecule. While these approaches represent a promising direction, their mathematical structure requires a certain care. In some cases, challenging higher-order gradients make these equations hardly tractable. In other cases, these models are based on phase-space formulations that suffer from well-known consistency issues. Here, we present a new complex fluid system that resolves these difficulties. Unlike common approaches, the current system is obtained by applying the fluid closure at the level of the action principle of the original phase-space model. As a result, the system inherits a Hamiltonian structure and retains energy/momentum balance. After discussing some of its structural properties and dynamical invariants, we illustrate the model in the case of pure-dephasing dynamics. We conclude by presenting some invariant planar models.


Madelung hydrodynamics in molecular motion
The search for cost-effective computational approaches to many-body quantum systems has stimulated a great deal of studies, often exploiting Madelung's hydrodynamic formulation of quantum mechanics [58].For example, the hydrodynamic approach has led to the development of several trajectory-based algorithms in quantum molecular dynamics [21,36,53,56,77], beyond the well-established Born-Oppenheimer theory.In this setting, one deals with the coupled motion of both nuclei and electrons in a molecular system for which the full Schrödinger equation represents a formidable challenge.The latter is normally tackled by convenient approaches that alleviate the computational complexity.Among these approaches, Madelung hydrodynamics has the great advantage of restoring the concept of trajectory, which can then be exploited to model the nuclear dynamics while the electronic motion is usually approximated as a finite-dimensional quantum system [3,24,51].
In this hydrodynamic picture, the Lagrangian trajectories associated to the nuclear motion sweep the electronic unitary evolution.In many important cases [3,51], this description is made possible by the observation that any two-particle wavefunction Ψ(r, x, t) may be expressed in terms of the exact factorization [1], i.e.Ψ(r, x, t) = χ(r, t)ψ(x, t; r), with |χ| 2 = 0 and ´|ψ| 2 d 3 x = 1.Here, the semicolon indicates that ψ is a conditional wavefunction so that the coordinate r appears merely as a parameter.Notice that this factorization is directly generalizable to the case of two families of several particles, i.e., nuclei with coordinates {r j } and electrons with coordinates {x k }, although here we will deal the simplest possible case of only one nucleus.In addition, to avoid possible complications from functional analysis, we will consider finite-dimensional truncations of the conditional wavefunction.This means that ψ(x, t; r) is replaced by a field ψ(t; r) taking values in the electronic Hilbert space H = C n .Replacing the exact factorization ansatz in the two-particle Schrödinger equation i ∂ t Ψ = − 2 ∆ r Ψ/M + H(r)Ψ and writing χ = √ De −iS/ , as well as ̺ = ψψ † , yields the system [32,51] Here, and throughout this paper, A|B = Tr(A † B) is the Frobenius inner product on C n×n and the Hamiltonian operator H(r) is an r-dependent Hermitian matrix.The latter is usually expressed as H(r) = V (r)1 + V I (r) + H Q , where V (r) is the scalar potential energy, H Q is an r-independent matrix, and V I (r) is a general interaction term.Also, we have denoted The density matrix of the electronic system is expressed as ρ = ´Ψ(r)Ψ † (r) d 3 r = ´D ̺ d 3 r.Equations (1) describe a system in which the quantum hydrodynamic flow sweeps a finitedimensional quantum system, which also undergoes its own unitary evolution.The latter comprises both a standard term [ H, ̺] of von Neumann type and an 2 -term carrying secondorder gradients.In turn, these terms also feed back in the hydrodynamic momentum equation.Together with the quantum potential V Q , these terms make the level of complexity of this system rather challenging.In particular, the quantum potential in Madelung hydrodynamics has motivated relevant works in PDE analysis [34,35,26,27].We observe that (1) has the structure of a complex fluid model, in which a macroscopic fluid equation is coupled to the evolution of a microscopic component.While in this case the latter is given by the density matrix of the quantum state, in liquid crystal theory this is given by the Q-tensor of the molecular rotational state [25].In physical terms, the quantum potential V Q is known to be responsible for the emergence of high-frequency and short-wavelength patterns that may lead to the formation of regions of space where the density D vanishes (quantum vortices) [10,29,33].Even if one attempted to discard V Q , the other nonlinear terms involving second-order gradients would still lead to analogous complications.Overall, the full hydrodynamic system (1) is accompanied by computational difficulties that still challenge the chemistry community [37,44,68,78].

Mixed quantum-classical fluid models
We observe that the difficulties mentioned above are removed if we deliberately eliminate the higher gradient terms in the first two equations of (1).This step is usually accompanied by arguments resorting to well-established WKB methods and corresponds to a classical treatment of the hydrodynamic flow: indeed, one verifies that dropping the higher gradient terms in (1) takes the equation for the phase function S(r, t) into a classical Hamilton-Jacobi equation.In the context of molecular dynamics, for example, this step corresponds to approximating nuclear motion as classical while leaving the electronic dynamics as fully quantum.
In the case of equations (1), neglecting the second-order 2 -terms yields the Ehrenfest fluid model for the interaction of a classical system and a quantum system: This complex fluid model is the simplest case of mixed quantum-classical dynamics, which generally comprises a vast set of models that have been devised to avoid the complications arising from fully quantum descriptions.While popular in molecular dynamics, similar approaches have also been proposed in solid state physics [7,54], spintronics [63], and quantum cosmology [13].Despite its mathematical appeal, the Ehrenfest fluid model ( 2) is insufficient to provide a reliable description of the coupled dynamics [4,73].Thus, much of the computational chemistry community has been developing alternative quantum-classical models beyond the Ehrenfest system [22].Among the most popular alternatives to the Ehrenfest model, the surface hopping algorithm was recently shown to violate Heisenberg's uncertainty principle, since the quantum density matrix ρ may become sign-indefinite [15].The same violation of the Heisenberg principle is also allowed by the quantum-classical Liouville equation on phase-space [5,16,43,55].In more generality, certain well-known consistency issues accompany current quantum-classical approaches beyond the Ehrenfest equations [2].On the one hand, since several studies do not seem to be affected by these issues, the existing approaches continue to be regarded as a useful tool.On the other hand, the systematic formulation of quantum-classical models beyond Ehrenfest dynamics remains a subject of current research [16,22,23,28,32,45,55].
While the above models arise as hydrodynamic formulations of multi-particle systems, similar approaches have also appeared in different contexts to study systems involving a genuine complex fluid component.For example, mixed quantum-classical fluid models were proposed in solvation theory to describe the interaction dynamics of macroscopic liquid solvents with quantum solute molecules [17,52].Since solvation phenomena occur in a wide variety of natural and technological processes, this solute-solvent interaction acquires crucial importance in different areas.In particular, the fluid-molecule coupling involves challenging nonlinear effects at different scales, which are revealed in unprecedented detail by time-resolved ultrafast experiments [65].Current solvation approaches are mostly based on different combinations of computationally expensive atomistic approaches with continuum descriptions involving many substantial approximations [66].A direction that remains little explored is the development of quantum-classical multiscale approaches combining microscopic solute dynamics with a mesoscopic solvent description.A remarkable effort in this direction is found in [17,18], where the authors formulate a quantum-classical hydrodynamic model by applying the moment method from kinetic theory to the quantum-classical Liouville equation.Depending on the closure adopted for the operator-valued tensor-density Π, the system presented in [17] reads Here, ρ and ĝ are operator-valued density and momentum-density, respectively, so that the fluid density and momentum are D = Tr ρ and m = Tr ĝ.Also, H eff is an effective Hamiltonian carrying nonlocal terms.Benchmarked in [52] for certain closure schemes, the system (3) provides a new promising perspective in solvation models.
Nonetheless, the necessity to address the consistency requirements of quantum-classical dynamics motivates us to consider also complementary modeling strategies to tackle problems in both molecular dynamics, solvation theory, and other fields.In particular, this paper presents a new class of mixed quantum-classical models lying beyond Ehrenfest dynamics and still satisfying its consistency properties.Similarly to the work in [17,18], the present complex fluid model is formulated as a closure scheme of a more fundamental mixed quantum-classical approach on phase-space.However, instead of applying the closure at the level of the equations of motion, we will operate on the variational principle underlying a recent phase-space formulation that, unlike other available variants, appears to ensure a series of important consistency properties.As a result, the new hydrodynamic system inherits desirable variational and Hamiltonian structures which are analogous to those appearing in the theory of liquid crystals [39,47,70].

Quantum-classical dynamics on phase-space
Several mixed quantum-classical models currently used in molecular dynamics are based on the phase-space formulation of classical mechanics [19,30,55,76].While doubling the dimensionality of classical motion, this approach offers the opportunity for a wider range of possible approximations leading to computational schemes.Over the years, however, these phase-space formulations have not helped to overcome the consistency issues affecting quantum-classical models beyond Ehrenfest dynamics.Until very recently.Blending Koopman's Hilbert-space formulation of classical mechanics with van Hove's unitary representations in symplectic geometry [14,57,75], the present authors proposed a new phase-space model of quantum-classical mechanics that satisfies all the following criteria [16] and still retains the geometric underpinnings of both classical and quantum mechanics [14,40,41]: 1) the classical system is identified by a probability density at all times; 2) the quantum system is identified by a positive-semidefinite density operator ρ at all times; 3) the model is equivariant under both quantum unitary transformations and classical canonical transformations; 4) in the absence of an interaction potential, the model reduces to uncoupled quantum and classical dynamics; 5) in the presence of an interaction potential, the quantum purity ρ 2 has nontrivial dynamics (decoherence property).
If f (q, p, t) denotes the classical Liouville density on phase-space and ψ(t; q, p) ∈ H is the conditional wavefunction representing the quantum dynamics (now depending on the classical phase-space coordinates), the model proposed in [40,41] reads as follows with P = ψψ † , and Here, {•, •} is the canonical Poisson bracket, H(q, p) is the Hermitian Hamiltonian operator parameterized by the classical phase-space coordinates, and we have denoted The density matrix of the quantum subsystem is expressed as ρ = ´f P d 3 qd 3 p.If we discard the -corrections in both the vector field X and the Hermitian generator H , the model reduces to a phase-space variant [79] of the Ehrenfest fluid equations (2).Indeed, after replacing X = P |X H and H = H in (4), with H(q, p) = M −1 |p| 2 /2 + H(q), and taking the hydrodynamic moments (D, MDu) = ´(1, p)f d 3 p and D ̺ = ´f P d 3 p, one recovers equations (2) by neglecting pressure effects.
While equations ( 4)-( 6) appear hardly tractable, a direct calculation of divX reveals that no gradients of order higher than two appear in the equations (4).In addition, the variational and Hamiltonian structures of this system [40,41] reveal the features occurring in quantumclassical coupling.For example, as we will see, equations (4)-( 6) lead to the identification of a quantum-classical Poincaré integral invariant which extends the classical quantity ¸p • dq.In more generality, the system (4)- (6) appears to be the first quantum-classical Hamiltonian model beyond Ehrenfest dynamics that satisfies the consistency criteria listed above.
Given the combination of consistency properties and underlying mathematical structure, we consider the equations ( 4)-( 6) as a platform for the formulation of simplified closure models that can be used in physically relevant cases.Indeed, the characteristic nature of the equations (4) offers the possibility of devising trajectory methods, e.g.following [32,51,71].For example, in [8] we recently showed how particle schemes based on (4)-( 6) succeed in reproducing crucial dynamical features in both quantum and classical sectors, well beyond the Ehrenfest model.However, the model ( 4)-( 6) involves the full classical phase-space, which makes the curse of dimensions particularly challenging.In an attempt to alleviate this problem, here we consider the formulation of a fluid closure on the configuration space.This closure represents a quantumclassical alternative to the quantum hydrodynamic approach in (1) and yet retains correlation effects beyond Ehrenfest dynamics.

Plan of the paper
The paper proceeds by illustrating the original quantum-classical phase-space model in §2.After an overview of its derivation, the underlying action principle is presented in §2.1, along with a discussion of its geometric structure.The Hamiltonian structure and its associated Casimir invariants are found in §2.2, while §2.3 deals with the equivariance properties of the model.
The complex fluid model is formulated in §3.After introducing the relevant fluid variables, the first closure steps are performed in §3.1 via a restriction to the configuration space of the original paths on phase-space.The crucial advance is presented in §3.2, where we show that the last term in (8) may be projected to the configuration space by resorting to a Poisson bracket structure known as Nambu bracket.This process leads to the appearance of two new scalar functions that are transported along the flow, one of which -the backreaction fieldincorporates the forces excerpted by the quantum system in the classical fluid flow.Then, the explicit form of the fluid equations is presented in §3.3, where we emphasize that the new model contains only first-order gradients, in contrast to the fully quantum treatment.
A discussion of the relevant properties of the new model is found in §4.In §4.1, we discuss the appearance of a non-Abelian gauge connection previously found in the context of molecular dynamics.In particular, this quantity seems to play an important role in the expression of the fluid stress tensor.The new complex fluid model possesses infinite families of dynamical invariants which are presented in §4.2, along with a discussion of the underlying Hamiltonian structure.Importantly, the Berry connection appears prominently when expressing the quantum dynamics in terms of wavefunctions, thereby allowing the characterization of a cross-helicity hydrodynamic invariant.An interesting specialization is presented in §4.3, where we show that the proposed model succeeds in retaining backreaction forces otherwise lost in the Ehrenfest treatment.Finally, a class of invariant planar subsystems is illustrated in §4.4.In this case, we show that restricting to two spatial dimensions allows for larger families of dynamical invariants that may be used for Lyapunov stability studies.

Geometric structure of the phase-space model
As anticipated, we will now review the formulation of the quantum-classical model ( 4)- (6), with special emphasis on its geometric structure [40].Before embarking on this discussion, we will start by presenting the theoretical background.In first place, one adopts Koopman's Hilbertspace formulation of classical mechanics [57]: writing f = |χ| 2 in the classical Liouville equation ∂ t f = {H, f } yields i ∂ t χ = i {H, χ} − ϕχ, where ϕ(q, p) is an arbitrary function.Then, since the operator i {H, }+ϕ is self-adjoint, the Koopman wavefunction χ(q, p, t) evolves unitarily, just in the same way as quantum Schrödinger wavefunctions.Guided by a series of arguments in representation theory, van Hove fixed the particular phase term ϕ = p • ∂ p H − H [75], that is the phase-space form of the Lagrangian usually written as p • q − H.This step makes Koopman's original formulation compatible with Hamilton-Jacobi theory [42].
A first quantum-classical model was obtained by writing the Koopman-van Hove equation i ∂ t χ = i {H, χ} − (p • ∂ p H − H)χ for two interacting systems and then quantizing one of them [14,42].While ensuring a series of important consistency criteria, the resulting quantumclassical wave equation i ∂ t Υ = i { H, Υ} − (p • ∂ p H − H)Υ generally allows for the classical Liouville density to become unsigned.Here, Υ ∈ L 2 (R 6 ) ⊗ H, where we recall that H is the quantum Hilbert space.The issue of classical positivity was addressed by resorting to a physical principle proposed by George Sudarshan [67]: any model of mixed quantum-classical dynamics should prevent classical phases from having observable effects.In other words, Hamilton-Jacobi functions are not measurable.This principle led the authors to treat the infinite-dimensional group of S 1 -valued functions on phase-space as a gauge group.Indeed, as shown in [41], applying a S 1 -symmetry by a suitable closure relation leads to treating classical phases as a gauge freedom.This is in direct analogy with the role of unit complex numbers in standard quantum mechanics.While applying the S 1 -symmetry to the action principle for the quantum-classical wave equation was found to restore positivity of the classical density, this procedure requires resorting to the Lagrangian trajectories underlying classical phase-space motion.This requirement is accommodated by combining the exact factorization method [1] with the Euler-Poincaré variational theory of continuum media [49], which lies as the basis for the new treatment.However, the resulting model i (∂ t + X • ∇)Υ + i (div X )Υ/2 = H Υ (up to phase factors) is highly nonlinear and a real intuition on its structure can only come from investigating its variational formulation.Here, X and H are given in ( 5)-( 6) upon writing f P = ΥΥ † .Also, if we insist on considering only finite-dimensional quantum systems, then Υ(q, p) is a Koopman wavefunction taking values in H = C n and satisfying ´Υ † Υd 6 z = 1.
To avoid dealing with unnecessary phase factors, the remainder of this section considers this model as given in ( 4)-( 6), that is in terms of the quantities f and P .

Action principle and dynamics on group orbits
The basis for the model equations ( 4)-( 6) is given by their underlying action principle δ ´t2 We refer to §5.1 in [41] for the derivation of this Lagrangian from ( 7) and the corresponding Euler-Poincaré variations.Here, we have denoted z = (q, p) and , = Re | , while A = (p, 0) is the coordinate representation of the Liouville one-form A = A • dz = p • dq, so that ω = −dA is the canonical symplectic form.Also, ξ(z, t) is a skew-Hermitian operator and P (z, t) := ψ(z, t)ψ(z, t) † .As usual in Euler-Poincaré variational theories, the variations δf and δX are constrained so that where Y is arbitrary.Analogously, we have where ζ is an arbitrary skew-Hermitian operator.These variations arise by standard Euler-Poincaré reduction from Lagrangian to Eulerian variables [49].In particular, if η(z 0 , t) is the diffeomorphic Lagrangian path on phase-space and U(z, t) is a unitary operator, we have the following relations between Lagrangian and Eulerian variables: where f 0 and P 0 are initial conditions.Here, η * f 0 = (f 0 / det ∇η) • η −1 denotes the pushforward of the initial density by η.We will denote the space of densities on the phase-space T * Q by Den(T * Q), which is identified with the distributional dual F (T * Q) * of scalar functions F (T * Q) via the L 2 -pairing.Also, we have P ∈ F (T * Q, Her(H)), where Her(H) is the space of Hermitian operators on H and F (T * Q, Her(H)) denotes the space of mappings T * Q → Her(H).
While f is positive-definite, P is only positive-semidefinite, and we will regard these properties as being preserved in time from suitable initial conditions.The first two relations in (11) indicate that the evolution of occurs on orbits of the semidirect-product group where Diff(T * Q) is the group of phase-space diffeomorphisms, and F (T * Q, U(H)) denotes the space of functions on T * Q taking values into the group U(H) of unitary operators on the quantum Hilbert space H.In particular, these orbits are determined by the group action given by the composition of the standard conjugation representation of F (T * Q, U(H)) and the pushforward action of Diff(T * Q).This evolution on group orbits will be crucial in later sections.In this paper, we will consider Q = R 3 for simplicity, although the treatment extends naturally to a smooth manifold.This geometric structure reflects directly in the equations ( 4)-( 6).Indeed, taking the time derivative of the first two relations in (11) respectively, and we observe the appearance of Lie-transport operators.Furthermore, upon taking variations of ( 8), the action principle δ ´t2 and we have used A := P , A .Then, after various manipulations we recover the system (4)- (6).The purely quantum and classical cases are found by restricting to the cases X H = 0 and H = H1, respectively [40].We readily realize that it would be rather unwise to try to obtain some information on this model by looking at the explicit equations ( 4)-( 6) obtained by expanding the functional derivatives δh/δf and δh/δ P .Instead, we are motivated to accept the guidance of their underlying geometric structure.While we have seen how the latter manifests in the variational formulation ( 8)- (10), we now consider the corresponding Hamiltonian setting.

Hamiltonian structure and Poincaré invariant
As we will see below, the Hamiltonian structure of the model ( 4)-( 6) is rather peculiar and we are not aware of similar structures occurring elsewhere in continuum mechanics.First, it is convenient to introduce the weighted variable P = f P so that f = Tr P and A = P, A / Tr P.Then, the Euler-Poincaré Lagrangian (8) may be expressed entirely in terms of the variables (X , ξ, P).Going through the same steps as above leads to rewriting (12) as X = X δh/δ P and [i ξ − δh/δ P, P] = 0, where h = ´ H Tr P + i { P, H} d 6 z.Then, the Hamiltonian equation where we have introduced the convenient notation A : B = Tr(AB).The proof that the bracket ( 14) is Poisson involves a combination of results in Lagrangian and Poisson reduction [41].On the one hand, we observe that the second term is naturally inherited from the usual Lie-Poisson structure underlying the quantum Liouville equation.On the other hand, the first term is a new type of Poisson bracket which reduces to the Lie-Poisson structure for the classical Liouville equation for f = Tr P if both functional derivatives are a multiple of the identity.Equation ( 13) easily leads to characterizing the Casimir invariant C 1 = Tr ´f Φ( P/f ) d 6 z for any matrix analytic function Φ.Also, upon writing P = f ψψ † , one finds the quantumclassical Poincaré integral invariant ¸c(t) (p•dq + ψ, i dψ ) = const.for any loop c(t) = η(c 0 , t) in phase-space, where we recall that η(t) is the flow of X .Then, by Stokes theorem, one identifies a Lie-transported quantum-classical two-form on T * Q, that is Ω(t) = η * Ω(0) , with Ω(t) := ω + Im dψ(t)| ∧ dψ(t) , so that Ω(t) is symplectic if it is so initially.As a result, one finds the additional class of Casimir invariants where Ω ∧ Ω ∧ Ω is a volume form and Θ is any scalar function of one variable.These Casimirs may be used to construct quantum-classical extensions of Gibbs/von Neumann entropies [40].

Quantum-classical von Neumann operator
We observe that the Hamiltonian energy functional h in ( 12) is not simply given by the usual average of the Hamiltonian operator H. Indeed, the −term seems to play a crucial role in taking the model ( 4)-( 6) beyond simple Ehrenfest dynamics.As discussed in [40], this suggests that the quantum-classical correlations trigger extra energy terms that are not usually considered.A similar situation occurs in standard quantum mechanics with spin-orbit coupling, that is the coupling of the quantum spin to the force acting on the orbital degrees of freedom.See Remark 4.1 below.Obtained from the semirelativistic limit of the Dirac equation, this effect is often discarded and yet it is crucial in a variety of contexts.
One may still rearrange the last relation in (12) in such a way that the total energy is formally given by an average of H. Following this route leads to rewriting the last relation in (12) as h = Tr ´ D H d 6 z, where is a measure-valued von Neumann operator.Then, classical and quantum densities are simply given by taking the trace and integral of D, respectively.Unlike the quantum density operator, we observe that the hybrid operator D is not sign-definite.Remarkably, however, D enjoys the equivariance properties where η is a symplectic diffeomorphism on T * Q and U ∈ U(H).These two properties ensure the following dynamics in the classical and quantum sector, respectively [41]: For example, upon denoting Σ = i[ P , X P ], the first property is verified directly by writing H identifies a bivector.With this expression of X one recognizes that divX = Tr{ P + div Σ, H}/2 − div Tr({log f, H} Σ)/2 involves only first-and second-order gradients.In more generality, due to the equivariance properties ( 16), the hybrid von Neumann operator acquires an important role in different aspects, such as the quantum-classical momentum balance [40].

Formulation of the fluid closure
Given the level of complexity of the model ( 4)-( 6), it is desirable to develop suitable closure schemes that can substantially simplify the treatment.In particular, here we present a fluid closure in which the classical system is represented by a fluid with density and momentum respectively, while the measure-valued density P(q, p) on phase-space is replaced by ˆf (q, p) P (q, p) where Tr ̺(q) = 1.In particular, we focus on quantum-classical Hamiltonians of the type While we could implement the moment method directly on the equations of motion, as customary in the kinetic theory of gases, here we perform the fluid closure at the level of the variational structure ( 8)-( 10) underlying the full phase-space model.Before proceeding, we emphasize that the current treatment does not include dissipation and diffusion.Indeed, we rely on the possibility to add standard diffusion and dissipation terms a posteriori, once the conservative equations are written explicitly.Similarly, pressure effects arising from the purely classical kinematic term p • ∂ q f /M in the first equation of (4) will also be added a posteriori by resorting to a conventional equation of state.

Restriction to cotangent lifts
So far the classical motion has been identified with diffeomorphic paths η on phase-space via their generating vector field X defined in (11).Each of these paths acts on the classical density as well as on the unitary operators governing the quantum dynamics as in (11).In order to obtain a fluid description of the classical motion, the latter needs to be restricted to occur on the configuration space.Since we want to operate on the Eulerian action principle associated to the Lagrangian (8), we need to restrict the phase-space vector field X to generate only transformations on the configuration space Q.We also ask for these transformations to possess a natural extension to phase-space.Such a class of transformations is well known under the name of cotangent lifts (or simply lifts) of diffeomorphisms.In particular, any diffeomorphism η ∈ Diff(Q) induces a phase-space diffeomorphism η L ∈ Diff(T * Q), that is the 'lift' of η, which reads η L (q, p) = η(q), ∇η −1 (η(q)) • p .Notice that we have η L * A = A, where A = p • dq.Likewise, in terms of Eulerian variables, any vector field u(q) identifies a lift on phase-space, which reads u L (q, p) = u(q), −∇u(q) • p .Thus, as a first step, we will restrict phase-space paths to cotangent lifts of paths on the configuration space, and this amounts to replacing ˆf (q, p)A(q, p) • X (q, p) d 6 z −→ ˆm(q) • u(q) d 3 q in (8).Here, we have used the second definition in (17) while the vector field u is defined in terms of Lagrangian paths as η = u• η.Notice that, on phase space, one also has ηL = u L • η L .
In pursuing our restriction to subgroup transformations on the configuration space, we also need to deal with the quantum propagators U ∈ F (T * Q, U(H)), which are again mappings defined on the classical phase space.In order to restrict to the configuration space, here we will consider only transformations U(q) in the subgroup F (Q, U(H)).Then, the quantum generator of motion ξ in (11) Consequently, by recalling (18), we make the following replacement in (8): ˆf (q, p) P (q, p), i ξ(q, p) d 6 z −→ ˆD(q) ̺(q), i ξ(q) d 3 q.
Another important step concerns the term ´f P , H d 6 z in (8).At this stage we adopt the cold-fluid closure f (q, p) = D(q)δ(p − m(q)/D(q)), thereby discarding pressure effects.Upon using (19), this expression of the phase-space density leads to the further replacement Notice that the cold-fluid closure is used here only as an intermediate step restricting to the simplest possible case.We will restore pressure effects later on by adding a fluid internal energy term to the right of the relation above.

The quantum backreaction
So far, the quantities in the Lagrangian (8) were restricted to the configuration space, except for the last term.The latter is responsible for the quantum backreaction on the classical flow beyond the usual Hellman-Feynman force averages [31] already appearing on the right-hand side of the momentum equation in the Ehrenfest model (2).The last term in ( 8) is responsible for a much intricate structure of the equations ( 4)- (6), which involve gradients of the quantum variable P and make the identification of a closure far from obvious.Here, we will proceed in stages and pursue the simplest possible closure method to retain nontrivial backreaction effects.Similarly to the procedure in the previous section, here we will continue to operate on the Lagrangian (8), thereby avoiding the necessity to control the various gradients in ( 4)- (6).We start by rearranging the last term in (8) as where the second divergence involves only position coordinates in configuration space.The first equality above follows from a slight rearrangement involving the projection of i{ P , H} = i∂ k P X k H on its Hermitian part as well as integration by parts.Also, the second equality follows from (19).Upon resorting to the cold-fluid closure, we recall (18) to write ̺ = P | p=m/D , and eventually perform the replacement div At this point, our closure problem amounts to finding a closure for the operator-valued current κ.In particular, we need to find a suitable expression of κ in terms of ̺ and ∇̺.For the sake of simplicity, here we assume a linear relation between κ and ∇̺.
To proceed further, we will be guided by the properties of the original model ( 4)- (6).For example, we notice that the term div[ P , f X P ] = f [∂ k P , X k P ]+[ P , X k P ]∂ k f in the first line of (20) does not produce gradients of order higher than one.Then, in order to avoid the emergence of higher-order gradients in div the assumption of linearity between κ and ∇̺ leaves us with the choice κ = −β × ∇̺ , where β(q, t) • dq is a real-valued differential one-form and the minus sign is chosen for later convenience.Notice that Tr ̺ = 1 =⇒ Tr κ = 0 and, in the case ̺ = ψψ † , we also have ̺| κ = 0.While this possible closure will be explored elsewhere, here we deal with the case β = c∇b for some scalar functions b(q, t) and c(q, t), that is κ = c∇̺ × ∇b.
With this expression of κ, the Lagrangian (8) finally reduces to where defines a Poisson bracket structure of Nambu type [62] on any functions F and G and for any function b.Since the latter is responsible for the backreaction terms in the new Lagrangian ( 21), we will refer to b as the backreaction field.We observe the striking analogy between the last term in (21) and the last term in (8).For later purpose, we notice that ̺, i {̺, Remark 3.1 (Extension to higher dimensions) The present treatment can be extended to consider the entire coordinate space R 3d for a molecule with d nuclei.In this case, each term in the Lagrangian (21) would be replaced by its natural higher-dimensional generalization.In particular, the last term would now involve the the direct sum of the same Nambu bracket d times.More explicitly, the Nambu bracket (22) would be replaced by the more general Poisson bracket {F, G} b (r 1 , . . ., Besides b and c, which are dealt with below, the Lagrangian (21) involves the following quantities: Here, the last two expressions arise from the restrictions performed in §3.1, while the first two are obtained by combining the definitions in ( 17)- (18) with the first two relations in (11), as well as replacing η → η and ξ → ξ.In particular, the first two in (23) yield the auxiliary equations which accompany the variational principle associated to (21).
The scalar field c is inserted in the expression of κ in such a way that the last integral in (21) can be made to converge.In particular, c(q, t) is conveniently prescribed as the ratio between the classical density D(q, t) and the Lie-transported volume element µ = η * µ 0 , where µ 0 is the elementary volume form in physical space.The appearance of the volume form in the denominator arises from the formal definition of Nambu-Poisson structures, as found in [74].Then, upon setting µ 0 = 1 in R 3 , we have Here, we emphasize the difference between the push-forward of the classical density η * D 0 = (D 0 / det ∇η) • η −1 and the simpler composition operation D 0 • η −1 , corresponding to the pushforward of the scalar function D 0 /µ 0 .The equation of motion for c(q, t) is thus So far, nothing has been said about the evolution of the backreaction field b in (22).The simplest choice would involve a time-independent function.However, we will now show that this is not compatible with the structure of the original model ( 4)- (6).A certain insight is provided by the observation that the terms involving H in (21) evidently plays the same role as the hybrid von Neumann operator D in (15).By pursuing the analogy with the original model ( 4)-( 6), here we will insist that D must satisfy analogous equivariance properties to those in (16).In the case that b is fixed and D is regarded as a map of the dynamical variables (D, ̺, c), evidently we would have, in push-forward notation, Instead, if b is time-dependent, so that D depends on the four variables (D, ̺, b, c), then we observe that In particular, here we make the following prescription for the evolution of b: that is b is also Lie transported as a scalar function.As we will see, this relation ensures conservation of the total momentum exchanged between the classical and the quantum systems.

Quantum-classical fluid equations
Having characterized the fluid Lagrangian (21) along with the evolution of the dynamical variables, we are now ready to write the equations of motion explicitly.In particular, Hamilton's variational principle δ ´t2 t 1 l dt = 0 involves the following Euler-Poincaré variations: and where w = δη • η −1 and ζ = δUU −1 • η −1 are both arbitrary and vanishing at the endpoints.These variations follow directly from ( 23) and (27).In addition, δm is arbitrary.Before proceeding further, in this section we will make use of the partial Legendre transform m = DMu in (21) and we will restore pressure by inserting an internal energy term of the type − ´DE(D) d 3 q.While pressure effects are absent in molecular dynamics, they become important for solute-solvent coupling in solvation hydrodynamics [17,52].In that context, the motion of a quantum molecule interacting with the surrounding fluid solvent cannot generally ignore the solvent pressure.Upon including the internal energy, the Lagrangian (21) becomes In analogy with the phase-space construction in §2.1, and as a result of the relations ( 23) and ( 27), the action principle associated to this Lagrangian restricts the evolution of (D, ̺, b, c) to occur on orbits of the semidirect product Diff(Q) F (Q, U(H)) acting on the space Den We notice the presence of the extra variables b and c arising from the closure method and previously absent in the phase-space approach.In addition, in the present case, the Lagrangian fluid trajectories evolve along a velocity vector field u that has its own equation of motion and will characterize the fluid momentum evolution.This situation differs from that in the phase-space model, where the Lagrangian paths were determined by a vector field which is given in the first equation of ( 12) by an algebraic expression involving the other dynamical variables.Despite these differences, the fluid system resulting from the action principle δ ´t2 t 1 ℓ dt = 0 has a geometric structure that is shared with common models of complex fluids [39,47,70].
While further geometric structures appearing in the present closure model will be presented later on, at the moment we are interested in the explicit form of the equations of motion.Upon denoting p = D 2 E ′ (D), the action principle associated to the Lagrangian (30) leads to as well as We refer to Appendix A for more calculational details.These quantum-classical fluid equations can be immediately compared to the fully quantum equations in (1).We briefly summarize our considerations as follows: • with respect to the fully quantum treatment in (1), the quantum-classical system (31)- (32) involves two extra variables, the backreaction field b and the quantity c, which are Lie transported as scalar functions; • the quantum potential V Q is now absent and its effects are now replaced by classical pressure effects from standard barotropic fluids; • unlike the quantum equations ( 1), the quantum-classical system (31)-( 32) recovers uncoupled quantum and classical dynamics if ∇ H = ∇H C 1, for some scalar function H C ; • while the fully quantum equations (1) contain terms of order 2 , only terms of the order appear in the the quantum-classical system ( 31)-( 32); • unlike the quantum equations ( 1), the quantum-classical system ( 31)-( 32) involves only first-order gradients; • unlike the solvation model (3), equations ( 31)- (32) ensure that the operator ̺ is positive semi-definite at all times.This follows from the second equation in ( 24); • the system ( 31)-( 32) may be readily extended to include a nonlocal dependence of the Hamiltonian H on the variables D and ρ = D ̺, as occurring in the effective Hamiltonian H eff of the solvation model ( 3) in [17].
Thus, despite their cumbersome appearance and the presence of two auxiliary frozen-in functions, the quantum-classical fluid equations represent a considerable simplification over the fully quantum hydrodynamics, where second-and third-order gradients are responsible for the appearance of interference patterns that pose severe challenges for current trajectory-based algorithms [78].
We also notice that, upon dropping the internal energy and setting an initially constant backreaction field, the system (31) recovers the quantum-classical Ehrenfest equations (2).In particular, the terms containing the b-field are unambiguously associated to the spatial inhomogeneities in the quantum state variable ̺ and these inhomogeneities are responsible for the quantum backreaction on the classical fluid flow in the current model.Furthermore, we emphasize again that the internal energy has been assumed here to depend only on the fluid density.While this is suitable for barotropic fluids, adiabatic flows may also be realized by letting E = E(D, S) depend also on the specific entropy S, which is then Lie transported as an additional scalar function.On this occasion, however, we will restrict to the barotropic case for the sake of simplicity.

Mead connection and the stress tensor
As we observed in §3.2, the quantum-classical von Neumann operator (25) plays the same role as its phase-space variant from §2.3.Remarkably, this operator is expressed in terms of a non-Abelian connection form.Defined on a U(H)-bundle over the configuration manifold Q = R 3 , this quantity is usually given as i[̺, ∇̺], where the imaginary unit is inserted so that this differential one-form takes values in the space Her(H) of Hermitian operators on H.For later convenience, here we will refer to the connection form as the Mead connection.Indeed, up to a numerical factor, the same quantity first appeared in Mead's work on geometric phases in molecular systems [61], although its role has not been much investigated so far.Within the context of mixed quantum-classical dynamics, a phasespace extension of the Mead connection made its first appearance in the original model ( 4)- (6).Indeed, as shown in [40,71], equation ( 5) may be rewritten as , and we recall z = (q, p).As we will see in this section, the Mead connection appears again in the present fluid closure, where it is defined on the configuration space.
As anticipated, the operator ( 25) can be expressed in terms of the Mead connection as D = D ̺ + div(c Γ × ∇b), and we notice that the classical density and the quantum density matrix are given by suitable projections of D, that is D = Tr D and ρ = ´ D d 3 q, respectively.Notice that, unlike other quantum-classical models widely popular in chemical physics [15,55], here both the classical and the quantum density are positive.In more generality, the model in (31) satisfies the consistency criteria outlined at the beginning of §1.3.In particular, the third property therein involves both quantum unitary transformations and classical transformations on the configuration space (that is, cotangent lifts on phase space).The equivariance properties possessed by the equations (31) follow from the equivariance of their underlying variational principle, which in turn hinges on the properties (26) possessed by the von Neumann operator D. In addition, following Proposition 5.11 in [41], the relations (26) naturally lead to In the case of quantum-classical dynamics involving an infinite-dimensional quantum subsystem, we write H(q) = H(q, x, p), with [x, p] = i 1.Then, the relations above ensure momentum conservation for translation-invariant Hamiltonians such as The equations (34) are verified explicitly by rewriting (31), after several rearrangements, as where we have introduced the stress tensor One verifies that T is symmetric by using the argument that any matrix of the form (a × b)c T + (c × a)b T + (b × c)a T is symmetric by the Jacobi identity.As shown in Appendix B, the symmetry of T follows from the fact that ̺, i {̺, H} b is invariant under rotations of ∇ H, ∇̺, and ∇b.Notice that this property would fail if the evolution equation ( 27) was dropped and the backreaction field b was treated as a constant: in that case, the last relation in ( 29) would be replaced by δb = 0, so that the resulting stress tensor would fail to be symmetric.In addition, we observe that the Mead connection plays a predominant role not only in the construction of the hybrid von Neumann operator, but also in the characterization of the fluid stresses.Before concluding this section, we observe that the momentum equation can be further rearranged as and, using the definition (33), we notice certain similarities between the first line above and the expression of the original phase-space vector field (5).

Hamiltonian structure and cross helicity
Having characterized the quantum-classical fluid system, its stress tensor, and the geometric quantities appearing therein, in this section we will show how more insight can be obtained by looking at the underlying Hamiltonian structure.As usual, the latter is given in terms of a Poisson structure and a Hamiltonian functional which identifies the energy of the system.The Hamiltonian structure of the fluid system (31) may be easily found by considering the action principle δ ´t2 t 1 l dt = 0 associated to (21) and introducing the convenient variable ρ = D ̺.Then, we have δ ´t2 Notice that here we have restored the internal energy in (21) in order to retain pressure effects.
The candidate Poisson structure is found by applying the relation ḟ = {f, h} and using the equations of motion associated to an arbitrary Hamiltonian.This step leads to the bracket which is Lie-Poisson on the dual of the semidirect-product Lie algebra Here, X(R 3 ) denotes the space of vector fields on R 3 , while u(H) identifies the Lie algebra of skew-Hermitian operators on the quantum Hilbert space H.We observe that this bracket structure is similar to that underlying the fully quantum hydrodynamic equations (1), although in that case the second terms in the first and second line are absent.Lie-Poisson structures of the type (38) have appeared in several instances over the years, in the context of complex fluid models [39].
The bracket (38) possesses the Casimir invariants where [46], while Φ and F are an arbitrary scalar function and an arbitrary matrix analytic function, respectively.We remark that the Casimirs C 2 extend a similar class of invariants appearing in the dynamics of ferromagnetic fluids [48].While here we treat C 1 and C 2 separately for convenience, we notice that they may be combined into only one functional of the type Tr ´DF(b, c, Λ n , D −1 ρ)d 3 q.
We will now show that more insight may be obtained upon writing ρ = Dψψ † , in which case C 1 becomes a trivial constant and Λ n vanishes, so that C 2 = ´DΦ(b, c) d 3 q.From the first two in (23), we find that the evolution of ψ(t) is ψ = (Uψ 0 ) • η −1 , so that the time derivative gives ∂ t ψ = ξψ − u • ∇ψ, and thus being the Berry connection [9].With this substitution and the introduction of the canonical fluid momentum M = m − DA, the phase-space Lagrangian (21) becomes Also, upon retaining the internal energy term, the Hamiltonian functional reads In this case, since δψ and δM are arbitrary, and M = D(Mu − A), the equations (31) become where we have denoted Notice that the equations (32) remain unchanged.Besides the usual Lorentz force, well known in molecular dynamics [3], we observe the presence of an extra hydrodynamic force involving the fluctuation force operator F which identifies the quantum fluctuations around the Hellman-Feynman average F = − ψ, ∇ Hψ [31].We remark that the last term in the momentum equation above is exactly the same as the last term in the first equation of (31), although the use of the conditional state vector ψ now unfolds the occurrence of the fluctuation force.This occurrence in the present fluid model does not come as a surprise, since the same quantity already appears in the phase-space model ( 4)-( 6), as shown in [40].Importantly, the momentum equation above can be written in terms of the Lie-derivative as (∂ t + £ u )(Mu − A) = −∇(δh/δD) + D −1 (δh/δc)∇c + D −1 (δh/δb)∇b.On the one hand, upon computing δh/δb = div ψ, i c F × ∇ψ and δh/δc = ∇b • ψ, i F × ∇ψ this leads to the following circulation dynamics: where γ = η • γ 0 is any loop moving with the Lagrangian fluid path η.The Lie-derivative form of the momentum equation reads equivalently as D(∂ t + £ u )(Mu −A) = −D∇(δh/δD) + (δh/δb)∇b + (δh/δc)∇c.Consequently, taking the dot product with ∇c × ∇b, one obtains the cross helicity invariant Here, the name follows from similar invariants appearing in the hydrodynamics of magnetized plasmas [20].An even more similar expression arises in hybrid kinetic-fluid models [72].
We have seen how the use of the conditional state vector ψ, while unfolding the role of the fluctuation force F, also allows for a systematic characterization of circulation and cross helicity.All these Casimirs may be used, for example, for a Lyapunov stability study via the energy-Casimir method [50].

Pure-dephasing systems
As an informative specialization of our fluid model ( 31)- (32), in this section we consider a particular type of quantum-classical systems, recently studied in [59].In the fully quantum formulation, a pure-dephasing system is given by a Hamiltonian operator of the type H = H 0 (x, p) + H I (x, p, A), where A is an operator commuting with the canonical observables (x, p).For the present discussion, it is convenient to make extensive use of the notation A = ̺, A .Also, for the sake of simplicity, here we will consider the case when A is a given Pauli matrix σ k and the fully quantum Hamiltonian is given as We will now present a comparison of results obtained from the study of pure-dephasing dynamics in quantum hydrodynamics, the Ehrenfest fluid model, and quantum-classical hydrodynamics.As we will see, the latter succeeds in retaining the quantum backreaction on the classical flow, which instead is lost in Ehrenfest dynamics.
In first place, we specialize the quantum hydrodynamics equations to the Hamiltonian (44), so that upon replacing H = V 0 1 + V I σ k in (1), we have where we have rearranged the quantum evolution equation.We observe that the overall expectation value σ k = ´D σ k d 3 r remains constant in time, thereby ensuring physical consistency.In particular, the initial condition σ k = 0 is preserved by the dynamics.Nevertheless, the same does not hold for the local expectation σ k , which indeed has nontrivial dynamics.This means that the quantum degrees of freedom feed back in the hydrodynamic flow via both last two force terms in the momentum equation.Thus, despite the simplicity of our pure-dephasing Hamiltonian (44), the full quantum hydrodynamics equations remain rather challenging.
The situation changes drastically in the case of the Ehrenfest fluid model.Indeed, if we replace H = V 0 1 + V I σ k in (2), we have In this case, it is clear that the initial condition σ k = 0 is preserved by the dynamics.This means that the hydrodynamic flow decouples entirely from the quantum motion, so that the Ehrenfest model fails to capture any quantum backreaction in pure-dephasing systems [59].
Finally, let us now consider the quantum-classical fluid model in ( 31)- (32).If we replace H = V 0 1 + V I σ k in (32), we obtain This dynamics is generally nontrivial, so that, unlike the the Ehrenfest model and similarly to the quantum case, the initial value σ k is not generally preserved in time.This means that the quantum evolution keeps feeding back into the classical flow.However, unlike quantum hydrodynamics, the fluid flow decouples entirely in the case when V I is spatially constant, in agreement with the item 4) of the consistency criteria discussed in the second paragraph of §1.3.
The present study has shown that the quantum-classical hydrodynamic model ( 31)-( 32) overcomes the problematic force cancelations inherited from Ehrenfest dynamics.Indeed, in the new model the quantum backreaction persists through the force terms containing the backreaction field b, which then acquires a fundamental role.While the persistence of backreaction is shared with quantum hydrodynamics, the absence of second-and third-order gradients in the model equations ( 31)-( 32) represents a substantial difference, which may lead to important simplifications from the viewpoint of both numerical and functional analysis.
Before concluding this section, we emphasize that, while emerging as drastic simplifications of more realistic problems, pure-dephasing systems have been considered in a variety of different fields, from optical physics [12], to quantum chemistry [64], and solvation dynamics [52].In the context of solvation hydrodynamics, quantum-classical pure-dephasing Hamiltonians were considered in [52], although in that case the purely classical potential V 0 is replaced by a nonlocal convolution of the solvent density D. While this takes the problem to a higher level of difficulty, none of the previous arguments would change in the presence of nonlocal potentials.

Invariant planar subsystems
In the search for a simplified version of the quantum-classical fluid equations ( 31)- (32), it is instructive to look for a lower-dimensional invariant subsystem.While we observe that the quantum backreaction is lost if we specialize our model to one spatial dimension, the same is not true for the two-dimensional case.Indeed, if we restrict to consider a planar flow with u = (v, 0) and v ∈ X(R 2 ), then b(x, y, z, t) = βz identifies an exact solution of the second equation in (32).If all the other variables are restricted to depend only on the planar coordinates, then the equations ( 31)- (32) become where we have denoted c = βc and {A, Here, no confusion should arise with the canonical phase-space bracket used in earlier sections.As it stands, the only substantial difference between the above planar subsystem and the full 3D model ( 31)-( 32) is that the backreaction field has now become merely a numerical parameter β, now incorporated in the field c = βc.Nevertheless, the backreaction terms persist in the momentum equation, with the exception of the vertical forces.We notice that these equations are again Hamiltonian with the Poisson bracket given by (38) without the last terms in the second line, and the Hamiltonian (37) with the replacement c{ , } b → c{ , } in the last term.This means that the functional C 1 in ( 39) is still a dynamical invariant, while the second Casimir drops to C 2 = ´DΦ(c, {Λ n }) d 3 q.In turn, upon considering the case ̺ = ψψ † , the third invariant (43) becomes C 3 = ´ΩΘ(c) d 3 q, where Ω = e 3 •∇×(Mv − A) = ω + Im{ψ † , ψ} is the canonical vorticity and Θ is an arbitrary function.In particular, we have ∇ × A = Be 3 with B = Im{ψ † , ψ}, and the equations of motion (42) resulting in the case ̺ = ψψ † reduce to In this case, we have the following circulation law: where γ(t) is again an arbitrary planar loop moving with the Lagrangian fluid flow.
In the attempt to further simplify our quantum-classical fluid equations, we now consider the case of an incompressible fluid flow.In this case, we have div v = 0 and the volume form d 3 q in physical space is preserved in time.For example, this situation could apply to an incompressible fluid solvent interacting with a quantum solute molecule.Then, since c = βD (recall the discussion in §3.2), the equations (47) become Here, the pressure p is now a Lagrange multiplier enforcing the condition div v = 0. We recognize that, even in the incompressible two-dimensional case, the backreaction forces persist and this is true also for pure-dephasing systems.Similarly to the case treated in the previous paragraph, the system (48) possesses the Casimir invariant C = ´(ΩΘ(D) + Φ(D)) d 3 q.

Conclusions and perspectives
The idea of restoring trajectories in quantum dynamics is extremely inviting, as it offers the possibility of new convenient computational schemes borrowing methods from classical simulation codes.However, despite the considerable work carried out over the last few decades on Madelung's hydrodynamics, the latter posses severe difficulties that continue to challenge the community.In this scenario, mixed quantum-classical models represent a promising perspective.While most of these models continue to be used in several test cases with a certain success, their underlying equations suffer from well-known consistency issues.The possible violation of Heisenberg's uncertainty principle in the most popular method [15] is only one example [2].
Having proposed a new Hamiltonian phase-space model overcoming these issues, here we have dealt with the problem of formulating a Hamiltonian fluid closure which goes beyond Ehrenfest dynamics and yet satisfies its consistency properties.As we showed, the proposed complex fluid model succeeds in capturing quantum-classical correlations in the case of puredephasing systems, which is precisely where the Ehrenfest model fails.Certain pure-dephasing systems were shown to be challenging also for other models alternative to Ehrenfest [14].Importantly, the proposed model has also the advantage of involving only first-order gradients.This point represents a major simplification over Madelung hydrodynamics, whose main difficulties arise from the appearance of higher-order gradients in its equations of motion.This simplification happens at the expense of introducing two extra scalar fields that are transported by the flow.While the field c is naturally linked to the fluid density, the backreaction field b deserves further attention.Indeed, it is not clear how b should be initialized.One has a few possibilities in this regard.
A possible way to initialize b consists in adopting a further closure, that is expressing it as a transported scalar that is constructed from the remaining variables.For example, the quantity Tr F (̺) is transported for any matrix analytic function F and one may set b = Tr F (̺) for some F .For example, one may think of setting b = α Tr(̺ ln ̺), for some parameter α.The use of entropy functionals in hydrodynamic closures of mixed quantum-classical dynamics was first proposed in [17].More simply, one may also set b = ̺ 2 .However, if ̺ is initialized as a projection, these expressions reduce to a constant thereby eliminating the backreaction.A possible alternative is to initialize b by its equilibrium value b e .Simple equilibria of the model equation ( 31)-( 32) may be easily found by setting δ(h + C) = 0, where C is one of the Casimirs treated in §4.2, or any combination thereof.For example, take the second in (39) with the particular choice C = ´DΦ(b) d 3 q, where Φ is an arbitrary function.Since we are interested in the equilibria of b, we can simply set δ(h + C)/δb = 0. Absorbing the constants into Φ, we obtain Φ ′ (b e ) = div c e ̺e , i∇̺ e × ∇ H /D e .Then, assuming that Φ ′ is invertible, we find that b e must be a function of the quantity div c e ̺e , i∇̺ e × ∇ H /D e , thereby characterizing a possible initial profile of the backreaction field.
As we have seen, interesting connections to spin-orbit coupling emerge and this is a completely unexplored direction suggesting that quantum-classical coupling may be modeled by the same algebraic structures appearing in the semirelativistic limit of the Dirac equation.Then, a natural open question concerns the comparison between the fully quantum treatment of spin-orbit coupling in quantum hydrodynamics and the analogue treatment in mixed quantum-classical dynamics.Addressing this question requires extending the present model to the case in which the coupling depends on the momentum.This extension is currently under development.
Finally, one may wonder about possible numerical implementations of the proposed model.Given the simpler level of difficulty, we plan to start our computational efforts by focusing on the planar subystems.As the backreaction field is absorbed into a constant number, not only does this case allow for less computational resources, but also eliminates the necessity of dealing with the question about the initial profile of b.Given the presence of several advection equations and a variational formulation, this case would also offer a testing ground for recent structurepreserving finite-element schemes developed within our groups for the long-time simulation of fluid problems [38].The second equation in (31) follows directly from the second in (49) by using the expression of the functional derivatives and the advection equation for D. Note that the variable ξ has been eliminated, and that the resulting equation for ̺ is compatible with its advection given in the second of (24).

B Stress tensor calculations
We provide here some details concerning the form of the right hand side of the momentum equation in (35), which involves the sum of two terms: the von Neumann operator term D, d H and the stress tensor term div T. We also explain the symmetry of T as emerging from the rotational invariance of the Lagrangian in terms of ∇̺, β = c∇b, ∇ H. so that the stress tensor T as defined in (52) gives the expression (36).

Structure
Symmetry of the stress tensor.We here show how the symmetry of T is related to the rotational invariance of the term ̺, i c{̺, H} b = ̺, i β•∇̺×∇ H under the transformation , for all R ∈ SO(3).This again is more easily seen by using the general Lagrangian (51).We consider the following SO(3) invariance of ǫ: ǫ( ξ, D, ̺, ∂ j ̺R j i , β j R j i , Ĥ, ∂ j HR j i ) = ǫ( ξ, D, ̺, for all R ∈ SO(3).In continuum mechanics, such types of invariance are related to objectivity and material frame indifference [60].Taking the derivative with respect to R at the identity in the direction ξ ∈ so(3), we have for all ξ ∈ so(3).This is equivalent to the symmetry of ∂ǫ ∂β ⊗ β − ∂ǫ ∂∇̺ , ⊗∇̺ − ∂ǫ ∂∇ Ĥ , ⊗∇ Ĥ .Recalling the general definition of T in (52), this is equivalent to the symmetry of the tensor T.