Geometry of nonadiabatic quantum hydrodynamics

The Hamiltonian action of a Lie group on a symplectic manifold induces a momentum map generalizing Noether's conserved quantity occurring in the case of a symmetry group. Then, when a Hamiltonian function can be written in terms of this momentum map, the Hamiltonian is called `collective'. Here, we derive collective Hamiltonians for a series of models in quantum molecular dynamics for which the Lie group is the composition of smooth invertible maps and unitary transformations. In this process, different fluid descriptions emerge from different factorization schemes for either the wavefunction or the density operator. After deriving this series of quantum fluid models, we regularize their Hamiltonians for finite $\hbar$ by introducing local spatial smoothing. In the case of standard quantum hydrodynamics, the $\hbar\ne0$ dynamics of the Lagrangian path can be derived as a finite-dimensional canonical Hamiltonian system for the evolution of singular solutions called `Bohmions', which follow Bohmian trajectories in configuration space. For molecular dynamics models, application of the smoothing process to a new factorization of the density operator leads to a finite-dimensional Hamiltonian system for the interaction of multiple (nuclear) Bohmions and a sequence of electronic quantum states.


Factorized wave functions in quantum molecular dynamics
Quantum molecular dynamics deals with the problem of solving the molecular Schrödinger equation i ∂ t Ψ = T n + T e + V n + V e + V I Ψ =: HΨ , (1.1) which governs the quantum evolution for a set of nuclei interacting with a set of electrons. In the equation above, Ψ({r k }, {x l }, t) is called the molecular wave function. The notation is such that {r k } = {r k : k = 1, · · · , N n } and {x l } = {x l : l = 1, · · · , N e } denote, respectively, N n nuclear and N e electronic coordinates. Each r k corresponds to a nucleus of mass M k whilst all electrons have the same mass m. The notation T and V in (1.1) refers to the kinetic energy and potential energy operators, while the subscripts n and e denote nuclear and electronic energies, respectively. The subscript I refers to the interaction potential between nuclei and electrons. More explicitly, one has where the potentials are of Coulomb type [63]. Without loss of generality, in this paper we shall consider only two particles, one nucleus and one electron, with coordinates denoted by r and x, respectively. As the molecular Schrödinger equation is practically intractable for standard computational methods, a series of different closures and approximations for extracting its dynamics have been developed over almost a century. Since the work of Born and Oppenheimer [16] in 1927, many efforts have been devoted to going beyond the adiabatic approximation in molecular dynamics, e.g., in the Jahn-Teller transition [84]. In the standard approach, one separates out the nuclear kinetic energy term by writing the molecular Hamiltonian operator as the sum, (1. 3) Here, H e is called electronic Hamiltonian. Since H e also includes the potential energy operators V n and V I (both depending on r), one may write H e = H e (r). The Hamiltonian operator H e (r) acts on the electronic Hilbert space H e , identified with the space of L 2 −functions of the electronic coordinates, x, which depend parametrically on the nuclear coordinates, r. Thus, its eigenvalue equation reads H e (r)φ n (x; r) = E n (r)φ n (x; r) .
At every point r, the eigenvectors φ n (x; r) provide an orthonormal frame in H e and the level sets of the eigenvalues E n (r) comprise hypersurfaces in the nuclear coordinate space, called potential energy surfaces (PES) [63,85]. As customary in the chemical physics literature, for simplicity, here we assume a discrete electronic spectrum. Formally, one can solve the molecular Schrödinger equation by writing the so-called Born-Huang expansion [17] Ψ(r, x, t) = ∞ n=0 Ω n (r, t)φ n (x; r) .
At this point, one can proceed by writing the equations for the coefficients Ω n (r, t), which in chemical physics are interpreted as nuclear wave functions.
Historically, various additional approximations have been made that have treated the nuclei as classical particles. In particular, the lowest-order truncation of the Born-Huang expansion is the Born-Oppenheimer (BO) approximation, [16] Ψ(r, x, t) = Ω 0 (r, t)φ 0 (x; r) .
In the BO approximation, the electrons remain in the ground state identified with the lowest nuclear eigenvalue E 0 (r) (adiabatic hypothesis). Within the physical chemistry community, it is widely accepted that stable molecular configurations correspond to the minima of the lowest energy PES, E 0 (r), see, e.g., [63].
Despite the several successes of the BO approximation, the adiabatic hypothesis appears to be too restrictive for realistic computer simulations. Consequently, a great deal of work has been devoted to formulating mathematical models for nonadiabatic dynamics [8,21]. Some of these approaches still exploit the Born-Huang expansion, while others introduce a fully timedependent ansatz for the molecular wave function Ψ. Among the most acknowledged models for capturing nonadiabatic effects, the mean-field model is probably the simplest. The standard mean-field factorization ansatz is given by Ψ(r, x, t) = χ(r, t)ψ(x, t) , (1.4) where both χ and ψ are wave functions (in their corresponding Hilbert spaces; respectively, H n and H e ). After finding the wave equations for χ and ψ, semiclassical methods are typically applied to the nuclear wave function, χ, thereby treating the nuclei as classical particles. Two different approaches are commonly used in dealing with factorized wave functions for molecular dynamics: 1) Frozen Gaussian wavepackets, [38,56], whose underlying geometric structure is based on coherent states [69]; and 2) Bohm's hydrodynamic approach [15] which reduces the nuclear dynamics to a Hamilton-Jacobi equation. While the geometry of the first approach has been illustrated in [67,12,66] (see also [24] for related discussions), the present paper deals only with Bohm's hydrodynamic approach [15]. While capturing nonadiabatic effects, the mean-field model does not adequately reproduce particle correlations between nuclei and electrons, [63]. Thus, the mean-field model is apparently too simple to apply accurately to realistic situations in which correlations are important. Therefore, alternative methods have been developing during the past few decades. For example, in Tully's surface hopping algorithm [79,80], probability amplitudes are used to design a stochastic algorithm that enforces "hopping" between different energy levels E n (r). In recent years, an augmented factorization scheme has also been proposed [1,2], following earlier works by Hunter [46] and going back to von Neumann's book [81]. This approach is currently known as exact factorization (EF), which reads as follows Ψ(r, x, t) = χ(r, t)ψ(x, t; r) .
(1. 5) In the exact factorization approach, the electronic wave function ψ depends parametrically on the nuclear coordinates r; so, it can be regarded as a mapping from the nuclear coordinate space into the electronic Hilbert space H e . In this sense, the exact factorization provides a time dependent generalization of the BO approximation. Although the classical limit of the nuclear wave function χ(r, t) could also have been taken by exploiting Gaussian wavepackets, the present work will investigate exact factorization by employing Bohm's hydrodynamic approach.
For an excellent review of the Gaussian wavepacket approach to nonadiabatic electronic effects, see [49]. This paper aims to investigate and compare the hydrodynamic approaches for both the mean-field models and the exact factorization ansatz in the context of geometric mechanics.
Within this framework, separating out the nuclear kinetic energy in the molecular Hamiltonian (1.3) corresponds in the hydrodynamic approach to transforming into a Lagrangian coordinate frame moving with the nuclei. Then, the momentum map associated to the evolution of the nuclear wave function collectivizes in the sense of Guillemin and Sternberg, [35,36,55,60,61]. This means that equivariant momentum maps transform canonical Hamiltonian dynamics into motion on coadjoint orbits generated by the action of a Lie group on the dual of its Lie algebra. Eventually, Lie-Poisson reduction leads to a new hydrodynamic formulation of nonadiabatic dynamics in which hyperbolicity is retained, rather than setting → 0 and passing to the Hamilton-Jacobi equation.
1. The remainder of Section 1 introduces background material which links standard elements of quantum mechanics with familiar objects in the setting of geometric mechanics. The fundamental variational principles and symplectic Hamiltonian structure in nonrelativistic quantum mechanics appear in Section 1.2.
2. The transformation to quantum hydrodynamics is discussed in Section 2.1. In Section 2.2, Bohmian trajectories [15] are reinterpreted as Lagrangian paths associated with the quantum hydrodynamic flow. Section 2.4 regularizes the → 0 limit of standard quantum hydrodynamics by suitably applying a spatial smoothing operator to the fluid variables of both the collectivized Hamiltonian and the corresponding Lagrangian before taking → 0. The resulting smoothed quantum fluid equations are found to admit singular solutions supported on delta functions. We call these singular solutions 'Bohmions', because the delta functions on which they are supported move along Lagrangian paths of the regularized quantum fluid Hamiltonians. Section 2.5 shows how the cold-fluid closure of Wigner distributions corresponds to a classical closure of mixed state dynamics arising from the Liouville-von Neumann equation.
3. In the mean-field approximation of coupled nuclear and electronic systems, the wave function is separated into a product of two independent factors, as in equation (1.4), above. Thus, the mean-field factorization of the wave function neglects the classicalquantum correlations between nuclei and electrons. Section 3 reviews the mean-field model and derives its quantum fluid representation in the geometric mechanics setting.
4. The exact factorization (EF) model [1,2,6,3] captures some of the nuclear and electronic correlation effects which are neglected in the mean field approximation, by letting the electron wave function depend on the nuclear spatial parameters, as in equation (1.5), above. Section 4 discusses the EF model in both the wave function and density matrix representations, then derives its quantum fluid representation in the geometric mechanics setting.
5. In Section 5 a new model is introduced by invoking a factorization ansatz at the level of the molecular density operator. Then, combining the classical closure of nuclear mixed states with the smoothing process presented in Section 2.4 leads to an entirely fintedimensional Hamiltonian system for the interaction of nuclear Bohmion solutions with an ensemble of quantum electronic states. Two different finite-dimensional schemes are presented depending on whether the smoothing process is applied in the Hamiltonian or in the variational formalism.

Quantum Lagrangians, Hamiltonians, and momentum maps
This section introduces the standard setting of the Hamiltonian approach to quantum dynamics by focusing on the evolution of pure quantum states. Later sections of this work will introduce von Neumann's density operators and their evolution for mixed states. However, the present section considers only Schrödinger-type equations. Among the most commonly used tools in chemical physics, the Dirac-Frenkel (DF) variational principle [28] for the evolution of the wave function is expressed in terms of the symplectic Hamiltonian structure of Schrödinger's equation. For a time-dependent quantum state ψ(t) in the Hilbert space H , the DF variational principle is expressed as a phase space Lagrangian, which reads Here, the bracket operation, · , · , defines the real-valued pairing which is induced by the natural inner product ψ 1 |ψ 2 on H , given by x, in which, e.g., ψ * 1 denotes the complex conjugate of ψ 1 , and ψ † 1 carries an implied integration.
The Schrödinger equation i ψ = Hψ follows as the Euler-Lagrange equation for the DF variational principle in (1.6), in which H is the self-adjoint Hamiltonian operator constructed from the canonical operators Q and P (the so called canonical observables). Thus, H = H( Q, P ) and [ Q, P ] = i ½. Notice that, since H is self-adjoint, the DF Lagrangian in (1.6) is U(1)−invariant so that the condition ψ 2 L 2 = 1 is naturally preserved. This amounts to conservation of the total probability. As presented in [12], the Euler-Poincaré formulation [42] of pure state dynamics is derived from the DF variational principle upon letting ψ(t) = U(t)ψ 0 with U(t) ∈ U(H ). Here, U(H ) denotes the group of unitary operators on H . In earlier years, this strategy was also exploited in [53,73] upon restricting U(t) to be the unitary representation of a finite-dimensional Lie group. For example, if U(t) is a representation of the Heisenberg group, substituting the ansatz ψ(t) = U(t)ψ 0 into the DF variational principle yields canonical Hamiltonian motion on phase space [72].
Notice that in (1.6), the functional h(ψ) = ψ, Hψ identifies the total energy of the system and, thus, it is deemed the Hamiltonian functional. The functional h(ψ) is sometimes called Dirac Hamiltonian, to distinguish it from the Hamiltonian operator, H. Depending on the context, the operator H and the functional h(ψ) may both be called the 'Hamiltonian'. More general systems (such as the nonlinear Schrödinger equation) can be obtained by replacing ψ, Hψ by a suitable functional h(ψ). In this case, the normalization condition ψ 2 L 2 = 1 must be incorporated as a constraint, [68], as where λ(t) is a real-valued Lagrange multiplier. For such constrained systems, the Euler-Lagrange equations yield the projective Schrödinger equation [52] ( along with the condition which implies that h(ψ) is U(1)−invariant, since for a phase shift δh = iψ, δh/δψ . Then, Noether's theorem for the U(1) symmetry of the constrained Lagrangian in (1.8) again implies conservation of ψ 2 = ψ, ψ , since the Lagrangian is invariant under infinitesimal phase shifts. Consequently, the constraint ψ 2 − 1 = 0 is satisfied and we may write, simply, which is the Hamiltonian form of the class of Schrödinger equations. The Hamiltonian structure for the class of Schrödinger equations is encoded in the following symplectic form on H : ω(ψ 1 , ψ 2 ) = 2 Im ψ 1 |ψ 2 = 2 iψ 1 , ψ 2 . In turn, this symplectic structure leads to the Poisson bracket given by This Poisson bracket then yields the corresponding Hamiltonian equation (1.11) via the expected relationḟ = {f, h}.
Both the DF variational principle and the Hamiltonian structure presented above will be used again and again throughout this paper to illuminate the geometric features of current models in nonadiabatic molecular dynamics. The next section will review the geometric setting of Bohm's quantum hydrodynamics in terms of its Hamiltonian structure. In particular, the next section will show that the Hamiltonian functional ψ| Hψ collectivizes, in the sense of Guillemin and Sternberg [35,36,55,60,61], through the momentum maps leading from Schrödinger's equation to quantum hydrodynamics. A (left) Hamiltonian action of a Lie group G on a symplectic manifold (M, ω) induces the momentum map J : M → g * , where g * is the dual space to the Lie algebra g of G. In the special case when M is a symplectic vector space (so that M = V ), then the momentum map is defined by where ξ(x) denotes the infinitesimal generator associated to the linear G−action on V . For example, if G = SO(3), the momentum map J(q, p) evaluates the angular momentum at each point (q, p) ∈ R 6 . In more generality, if M is replaced by a Poisson manifold (so that M = P ) with Poisson bracket {·, ·} P , the momentum map (if it exists) is defined as This relation defines the Lie-Poisson bracket on g * , given in terms of the adjoint action of the Lie algebra on itself, ad : g × g → g, denoted as ad ξ ζ = [ξ, ζ] for any Lie algebra elements ξ, ζ ∈ g. Upon denoting the pairing by · , · g : g * × g → R, the Lie-Poisson bracket on g * reads as [62] f (1.14) Momentum maps are ubiquitous in quantum mechanics. For example, the projection operator J(ψ) = −i ψψ † is a momentum map for the left action of the unitary group U(H ) on the quantum state space H . Other important examples are given by quantum expectation values [13] and covariance matrices of Gaussian wavepackets [66].

Quantum hydrodynamics
This section illustrates the geometry of the hydrodynamic setting of quantum mechanics, which has its foundations in the Madelung transform [58,59]. After reviewing the geometry of halfdensities and their momentum maps, the latter are exploited to derive the quantum hydrodynamic (QHD) equations.
Consequently, the space Den 1/2 (R 3 ) inherits the standard symplectic form on L 2 (R 3 ). As a result, we have the 1-form density, where the last equality follows from writing the wave function ψ in polar form, ψ = √ De iθ . The momentum map J (ψ) coincides (up to the mass factor, m) with the well-known probability current from quantum mechanics. It is also well known that the physical Hamiltonian operator H = − 2 ∆/2m + V (x) transforms the total energy into the form as can be verified by a direct calculation. Here, the quantity D = |ψ| 2 arises as another momentum map which is associated to the action of (local) phase transformations ψ(x) → e iϕ(x) ψ(x). Indeed, this action on the phase of the wave function has the momentum map |ψ| 2 = D.
Relation (2.4) shows that the Hamiltonian functional ψ| Hψ collectivizes, in the sense of Guillemin and Sternberg [35,36], through the momentum maps J (ψ) and |ψ| 2 . That is, the Hamiltonian in (2.4) may be expressed solely in terms of the collective variables µ and D, given by ψ| Hψ = h(µ, D) with µ = J (ψ) and (µ, D) ∈ (X(R 3 ) C ∞ (R 3 )) * . Here, X(R 3 ) C ∞ (R 3 ) denotes the semidirect-product Lie algebra of the semidirect-product Lie group Diff(R 3 ) C ∞ (R 3 , S 1 ), whose elements (η, ϕ) act from the left on the space This formula extends the action in (2.1) to include a local phase shift. The important feature here is that, under the collectivization (J (ψ), |ψ| 2 ) → (µ, D), the Hamiltonian h(µ, D) given by (2.4) belongs to a widely studied class of Hamiltonians possessing the Lie-Poisson bracket structure. This structure has been derived from the Euler-Poincaré formulation of ideal classical continuum dynamics with advected quantities in [42]. In particular, upon defining the velocity vector field via the reduced Legendre transform, u = δh/δµ = m −1 µ/D, the Euler-Poincaré Lagrangian associated to the Hamiltonian (2.4) reads Throughout this paper, we shall denote Euler-Poincaré Lagrangians by ℓ and ordinary Lagrangians by L. In (2.6), the velocity vector field is given by u =η • η −1 ∈ X(R 3 ), while the Eulerian density D is defined as for a reference density, D 0 = D 0 (x) d 3 x. In the last definition, the symbol η * denotes the operation of push-forward by the map η ∈ Diff(R 3 ); so, η * D 0 denotes the push-forward of the reference density, D 0 by the map η. Push-forward by the smooth flow η is called advection in hydrodynamics. In this context, the Lagrangian particle path of a fluid parcel is given by the smooth, invertible, time-dependent map, η t : R 3 → R 3 , as follows, η t x 0 = η(x 0 , t) ∈ R 3 for initial reference position η 0 x 0 = η(x 0 , 0) = x 0 . After this definition, there should be no confusion between η t ∈ Diff(R 3 ) and η t x 0 = η(x 0 , t) ∈ R 3 . The subscript t is omitted in most of this paper, for simplicity of notation. Now, taking variations in Hamilton's principle δ´t 2 t 1 ℓ(u, D) dt = 0 for the reduced Lagrangian (2.6) yields the following quantum hydrodynamics (QHD) equations, Here, the quantum potential arises from taking variations of the middle term of the reduced Lagrangian in (2.6), which can be rearranged as Remark 2.1 (Effects of the quantum potential) The appearance of the amplitude of the wave function in the denominator of the quantum potential in (2.9) implies that its effects do not necessarily fall off with distance. That is, the effects of the quantum potential need not decrease, as the amplitude of the wave function decreases. Moreover, the middle term in (2.6) is known as the Fisher-Rao norm, which is well-known in information theory. For further discussion of the information geometry in quantum mechanics, see, e.g., [18].
Equations (2.8) follow from Hamilton's principle for the reduced (collective) Lagrangian (2.6), upon using the following constrained variations from the Euler-Poincaré theory of ideal fluids with advected quantities, derived in [42], Here, the arbitrary vector field w = δη • η −1 ∈ X(R 3 ) vanishes at the endpoints in time. The density D is an advected quantity, satisfying the mass transport equation in (2.8).

Bohmian trajectories, Lagrangian paths & Newton's Law
In quantum hydrodynamics, the role of the Lagrangian path η ∈ Diff(R 3 ) is of paramount importance. Namely, it plays the role of a hidden variable in the Bohmian interpretation of quantum dynamics [15]. Indeed, in this framework the path η is the fundamental dynamical variable, while the wave function is simply transported in time along the Lagrangian motion of η(x 0 , t), which in turn satisfiesη This relation defines the so-called Bohmian trajectory, which is precisely the Lagrangian fluid path of the hydrodynamic picture! At this point, it is important to emphasize that the (infinite-dimensional) Bohmian trajectories η(x 0 , t) are completely different from the point particle trajectories q(t) (finitedimensional), which arise when the quantum dispersion is neglected. In order to clarify this point, it is convenient to rewrite the Lagrangian (2.6) in terms of the Bohmian trajectory by using (2.7) and (2.11). We have where the quantum potential is written in terms of η as The dynamics of the Bohmian trajectory η is given by the Euler-Lagrange equation [85] . We emphasize that the dynamics of the Bohmian trajectory η is not equivalent to point particle dynamics. In principle, the latter could be obtained by setting a point-like initial density of the type D 0 (x 0 ) = δ(x 0 − q 0 ) and then integrating the Euler-Lagrange equation over D 0 . However, this type of initial condition is not allowed by the structure of the quantum potential. For this reason, asymptotic semiclassical methods are required to properly derive the effects of the quantum potential in a weak limit as 2 → 0. For more details, see, e.g., [48]. Nonetheless, the Newtonian limit neglects the order O( 2 ) quantum dispersion term in the Lagrangian (2.6) (or, equivalently, (2.12)) and varies the remainder. The resulting equation for the Bohmian trajectory becomes D 0 mη + ∇ η V (η) = 0. It is clear that the point-particle initial condition D 0 (x 0 ) = δ(x 0 − q 0 ) is now allowed and thus denoting q(t) = η(q 0 , t) and integrating over space yields Newton's Law mq + ∇V (q) = 0.
In the Eulerian picture, one proceeds analogously by discarding the O( 2 ) in the Lagrangian (2.6), so that the equations of motion (2.8) restrict to Then, one considers the relations (2.7) and (2.11) between Lagrangian and Eulerian quantities. We observe that the initial particle-type initial density Integrating (2.13) over space again recovers Newton's Law for q(t).
A common alternative method to derive Newton's Law exploits the analogy with the Hamilton-Jacobi equation of classical mechanics. Since u = m −1 µ/D = ∇θ according to the momentum map relation for the collective variable J(ψ) in (2.3), the first of these restricted QHD equations happens to recover the Hamilton-Jacobi equation for S = θ, as for geometrical optics with the classical Hamiltonian H(q, p) = |p| 2 /2M + V (q). This is not necessarily convenient for a fluids interpretation, though, because solutions of Hamilton-Jacobi equations may become singular (e.g., form caustics) even for smooth initial data.
Remark 2.2 (Regularization of an "ultraviolet catastrophe" for 2 → 0) The Newtonian limit of QHD (2.8), obtained by simply neglecting the contribution to the Euler-Poincaré equations from the quantum potential in (2.9) turns out to be problematic. In particular, because the potential (which plays the role of a pressure term) is assumed to be independent of time, the Newtonian limit system (2.13) is not strictly hyperbolic. This observation is a well known signal, [54], that the solution behaviour in the classical limit 2 → 0 can become singular, as 2 multiplies the highest spatial derivative. This is especially clear when the wave function is written in the usual WKB form, as ψ = √ D exp(iS/ ), where S is the action integral for the Schrödinger equation. Indeed, in the 2 → 0 limit, the gradient of the quantum potential produces highly oscillatory spatial behaviour. See, e.g., [30,48] and references therein for discussions of the weak convergence of the rapidly oscillatory solutions obtained in passing the WKB description of the Schrödinger equation to the classical limit as 2 → 0. As indicated in [48], one convenient way to carry out the limit 2 → 0 is to apply the Wigner transform of the wave function [82,86]. The real-valued Wigner function then acts as the quantum-equivalent of the phase-space distribution function in classical mechanics; although the quantum version introduces mathematically technical features, such as Moyal operator brackets, instead of Poisson brackets. One could also treat the rapid oscillations as being stochastic and use probability theory to obtain the expected solution as a classical limit, [74]. Thus, in retrospect, one can appreciate the role of non-zero 2 in the quantum potential (2.9) in equations (2.8) as a dispersive regularization of what would otherwise have led to a type of ultraviolet catastrophe [22] for the restricted (Hamilton-Jacobi) QHD in (2.13), as the solutions of the restricted QHD equations form caustics.

Lie-Poisson structure of quantum hydrodynamics
In terms of the variables (µ, D), the collective Hamiltonian in (2.4) for the QHD equations in (2.8) reads (2.14) In these variables, the QHD equations can be written in Hamiltonian form, with a Lie-Poisson bracket written symbolically as, see, e.g., [42] ∂ ∂t in which each box in (2.15) indicates where to substitute elements of the last column of variational derivatives of the Hamiltonian in the matrix multiplication. Here ad * denotes the coadjoint action of the Lie algebra X(R 3 ) on its dual, the 1-form densities X( . The coadjoint action ad * : g × g * → g * in (2.15) is dual to the adjoint action ad : g × g → g, via the pairing · , · g : g * × g → R in which the Lie-Poisson bracket in equation (1.14) was defined, see, e.g., [44,62], ad * ∂h/∂µ µ , ∂f /∂µ g := µ , ad ∂h/∂µ (∂f /∂µ) g . The symbol £ u in (2.15) denotes Lie derivative with respect to the vector field u =η • η −1 ∈ X(R 3 ). For example, the corresponding Lie derivative of the density D(x, t)d 3 x is given by for elements of the tensor space a ∈ V and on its dual δh/δa ∈ V * . In the example of the advected density, we have D ∈ Den(R 3 ). The corresponding notation is defined explicitly for QHD in components by (2.18) The Hamiltonian operator C in (2.15) is linear in the dynamical variables (µ, D). The corresponding Lie-Poisson bracket is given explicitly by This is the Lie-Poisson bracket dual to the semidirect product Lie algebra X( For further discussions of Lie-Poisson brackets, see, e.g., [62,44] and references therein.

Regularized QHD and Bohmion solutions
In remark 2.2, we have emphasized that the order O( 2 ) term in the QHD Hamiltonian (2.14) may be regarded as a dispersive regularization of the "ultraviolet catastrophe" which occurs in the quantum fluid Hamiltonians as 2 → 0. The order O( 2 ) term is an energy penalty for high gradients |∇D| 2 /D = 4|∇ √ D| 2 that yields only a weak classical limit as 2 → 0 [48]. In this section we proceed by regularizing the Lagrangian or Hamiltonian to allow for singleparticle solutions. As we have observed in Section 2.2, the O( 2 )−terms in QHD prevent the existence of particle-like solutions so that Bohmian trajectories can only be identified with Lagrangian paths following the characteristic curves of the Eulerian fluid velocity. Thus, the O( 2 )−terms in QHD must be treated with particular attention. Instead of adopting semiclassical methods to take the limit 2 → 0, this section presents an alternative strategy consisting in regularizing the O( 2 )−terms by a smoothening process. More particularly, we shall discuss Bohmian trajectories which can be computed from regularized QHD Hamiltonians and Lagrangians, whose fluid variables have been spatially smoothed; so that their 2 → 0 limit is no longer singular.
Depending on which terms are regularized, different particle motions may emerge. We present two regularization strategies. The first simply smoothens all the terms in the Hamiltonian, while the second only smoothens the O( 2 )−terms. Although all the equations of motion derived here are Hamiltonian equations on a canonical phase space, they may or may not be in the usual Newtonian form, depending on which regularization scheme is adopted. In particular, the first regularization scheme is adopted in the Hamiltonian framework to regularize the hydrodynamic momentum and density, while the second scheme is based on the variational approach following the standard Bohmian method.
Hamiltonian regularization. The first regularization introduces the mollified collective Hamiltonian for regularized quantum hydrodynamics (RQHD). This is simply obtained by replacing the variables µ and D in (2.14) by the corresponding spatially smoothed variables, where K(x, s) is a positive definite, symmetric smoothing kernel which falls off at least exponentially in |x − s|. For example, we may take the kernel K(x, s) to be the Green's function of the Helmholtz operator (1 − α 2 ∆), where α is a length scale and the limit α → 0 returns the original hydrodynamic variables. This choice of smoothing kernel is also an energy penalty for high gradients. It promotes functions (µ, D) that are bounded in the L 2 norm to functions (μ,D) that are bounded in the H 1 Sobolev norm. Similarly, choosing K(x, s) to be the Greens function for an integer power p of the Helmholtz operator would smooth gradients to promote functions of (μ,D) from being bounded in L 2 , to being bounded in the Sobolev space, H p . Yet another choice would be to take K(x, s) to be a Gaussian.
Upon replacing the variables (µ, D) in (2.14) for Hamiltonian h(µ, D) by the spatially smoothed variables (μ,D) as h(μ,D) = h(K * µ, K * D), we find the following Hamiltonian equations for the original variables, (µ, D), Of course, the regularized quantum hydrodynamics (RQHD) equations arising after replacing the Hamiltonian in (2.14) by h RQHD = h(μ,D) must take the same Lie-Poisson form as in equations (2.15), modified now to read The variations of h RQHD (µ, D) are given by the following L 2 functions for an appropriate choice of the kernel K, rewritten in the same form as for the unsmoothed variables, At this point, the advantage of having regularized by simply smoothing the variables in the Hamiltonian by K * will emerge. Namely, while the physical meaning of the various expressions in the Hamiltonian has been preserved, the solutions in the original variables (µ, D) can now be singular and finite dimensional along the Lagrangian paths for the diffeomorphism, η in Section 2.2. Specifically, equations (2.21) for h RQHD = h(μ,D) admit Lagrangian paths as particle-like singular solutions for (µ, D), which we propose to call 'Bohmions'. These are given by the singular momentum map with a w a = 1. The momentum map (2.25) was presented in [45] as the immediate extension of the singular momentum map represented by the first relation in (2.25) and first discovered in [41]. This singular momentum map underlies the 'peakon' singular solutions for the nonlinear wave/fluid equations in [20,43]. These previously investigated 'peakon' singular solutions were obtained, respectively, by smoothing the momentum in the kinetic energy for either geodesic motion on Diff(R 3 ) in the case of [20], and by smoothing both the momentum and the depth for shallow water dynamics in [43]. Before proceeding further, let us emphasize that the formal process leading to the singular solutions (2.25) has provided a Lagrangian-particle representation of the dynamics. However, this does not imply that the variables D = |ψ| 2 and µ = Im(ψ * ∇ψ) will become delta functions. Nevertheless, as we shall see below, the dynamics emerging from the singular expressions (2.25) does reveal the Lagrangian-particle content of quantum hydrodynamics. Namely, the Bohmian 'particles' are merely fluid labels following Lagrangian flow trajectories. As such, the Bohmion singular solutions in (2.25) represent flow lines in quantum hydrodynamics as virtual particles. Of course, this is not a radically new idea, since particle methods have a long and successful history in computational fluid dynamics. Upon restricting to consider smoothing kernels of the type K(x, s) = K(x−s), substitution of the Bohmion singular solutions (2.25) Then, according to equivariance of the momentum map (2.25) discovered in [41,45], the dynamics of (q, p) satisfy the canonically conjugate Hamiltonian equations in phase space, (2.28) Both the momentum term and the quantum term proportional to 2 in the canonical equations (2.28) provide extensive, potentially long-range coupling among the singular particle-like Bohmion solutions, because of the presence ofD in the denominators of these terms in the Hamiltonian h RQHD in (2.22). However, the limit 2 → 0 in the canonical Bohmion equations (2.28) is no longer singular.
Lagrangian regularization. So far, the discussion has been entirely based on the Hamiltonian structure of QHD. In this section, we shall take an alternative route: instead of regularizing all terms in the QHD Hamiltonian, we shall regularize only the O( 2 )−terms in the QHD Lagrangian (2.6). As a first step, we perform the substitution D →D in the O( 2 )−term of the original QHD Lagrangian (2.6), which then becomes Then, upon following the arguments in the previous section, we set the initial condition Then, upon recalling the Lagrange-to-Euler map in (2.7), the Eulerian density becomes D(x, t) = a w a δ(x − q a (t)). In turn, this expression can be inserted into the regularized QHD Lagrangian (2.29), which becomes for which Hamilton's principle produces the Euler-Lagrange equation In analogy to the arguments in the previous section, we emphasize that the formal process outlined here reveals a form of Newtonian dynamics, which is suitably modified by the regularized expression of the quantum potential. While the formal relation (2.30) provides a particle description, it does not imply, for example, that a smooth initial probability distribution D = |ψ| 2 would evolve to concentrate into a delta function.

Density operators and classical closures
So far, the discussion has focused uniquely on wave functions for pure quantum states. However, mixed quantum states are a more general class of states that can be represented by a Hermitian, unit-trace, and positive-definite integral operator ρ satisfying the Liouville-von Neumann equation Pure states are regarded as a special case of mixed states under the identification ρ = ψψ † , or equivalently ρ(x, x ′ ) = ψ(x)ψ * (x ′ ). Here, the notation ρ(x, x ′ ) is used for the kernel (matrix element, in the physics literature) of the integral operator ρ. In more generality, for an arbitrary sequence {ψ n (x)} of N square-integrable functions, the density operator is given by with n w n = 1. For the momentum map aspects of quantum mixed states, see [65,77]. Equation . (2.34) Given the above left action, the verification detailed in the Appendix A shows that its infinitesimal generator may be written as where P k = −i ∂ k and thus, by using (1.13), one can prove that the corresponding momentum map is given in matrix element notation as For the special case of pure states, one verifies that ρ(x, x ′ ) = ψ(x)ψ * (x ′ ) recovers the momentum map (2.3). However, in the general case of mixed quantum states, the dynamics of J (ρ ) cannot be expressed only in terms of µ(x) itself and D(x) := ρ(x, x) as in the case of pure states [85]. Rather, mixed states lead to a multi-fluid system that is obtained by combining the arguments in Section 2.1 with the relation (2.33). Nevertheless, here we show that the classical limit of mixed state dynamics (as given by the Liouville-von Neumann equation (2.32)) can be obtained from an exact closure by allowing for the operator ρ to be sign-indefinite. This should come as no surprise, since classical states violate the uncertainty principle in such a way that the density operator can no longer be positive-definite. The proposed classical closure for the Liouville-von Neumann equation (2.32) is expressed as where mD(x)u(x) = µ(x), as one can show by a direct verification. With the ansatz above, the total energy ρ, H becomes ρ, H =´(|µ| 2 /(2mD) + DV ) d 3 x, which coincides with the QHD Hamiltonian (2.14) after dropping the 2 −term. Thus, the corresponding equations of motion naturally coincide with the classical hydrodynamic limit (2.13) in terms of Newton's Law, as discussed in Section 2.2. However, it is important to emphasize that, unlike pure states, here the fluid velocity u is no longer an exact differential and thus the corresponding hydrodynamic flow preserves the nontrivial circulation¸γ u(x) · dx, for an arbitrary loop γ moving with velocity u. Then, this produces the vorticity dynamics ∂ t ω = curl(u × ω), with ω = curl u. In addition, since the probability density is no longer defined in terms of a square-integrable function, the discussion in Section 2.2 can be extended by allowing for the multi-particle initial condition (2.30). Then, upon defining q a (t) = η(q a (0), t), the relations (2.7) lead to the Eulerian density D(x, t) = N a=1 w a δ(x − q a (t)), where each q a (t) satisfies Newton's Law for N non-interacting particles.
We conclude this section by discussing the nature of the closure (2.37) in terms of the Wigner function W (x, p) = (2π ) −3´ρ (x + y/2, x − y/2) e −i p·y d 3 y. Without entering discussions about the Wigner-Moyal formulation of quantum mechanics, we shall simply address the reader to [86] and present the closure (2.37) as the operator ρ associated to the following Wigner function: W (x, p) = D(x)δ(p − mu(x)) . (2.38) Similar considerations of the cold-fluid closure (2.38) have already appeared [48] in the context of the semiclassical limit for pure state dynamics. See also [19] and references therein for the use of delta-function closures in hybrid quantum-classical dynamics. Finally, we emphasize again that the Wigner function in (2.38) does not identify a quantum state. This is analogous to what happens for the quantum harmonic oscillator: in this case, the Wigner-Moyal equation coincides with the classical Liouville equation thereby allowing for delta-function solutions returning classical motion. However, delta-function Wigner distributions do not correspond to quantum states, as their associated density operator is signindefinite.

The mean-field ansatz
Here we apply the Euler-Poincaré method [42] to formulate the mean-field model in terms of reduction by symmetry. Then, the molecular wave function for a nucleus and an electron is written as Ψ := Ψ(r, x) ∈ L 2 (R 3 × R 3 ). We think of the coordinate r as that corresponding to the nucleus, while x corresponds to the electron. In physics, a mean-field ansatz introduces a factorization of the type Ψ(r, where χ(r, t) is regarded as the wave function describing nuclear dynamics (corresponding to the subsystem we want to treat classically) and both ψ and χ are normalized with respect to the coordinate upon which they depend. The Hamiltonian operator now takes the form H = H( Γ, Z), where we introduce the notation Γ = (r, −i ∇ r ) and Z = (x, −i ∇ x ). Upon recalling the pairing Ψ 1 , Ψ 2 = Re´Ψ * 1 (r, x)Ψ 2 (r, x) d 3 r d 3 x, insertion of the ansatz (3.1) into the action principle (1.6) yields a Lagrangian L : T H n × T H e → R, given by where the second equality uses the natural pairings on the respective Hilbert spaces H n and H e , and we have introduced the effective Hamiltonian At this point, a further approximation is often introduced; namely, one assumes that the nuclear dynamics can be treated as classical. This assumption produces a mixed classicalquantum system. In what follows, we will introduce a geometric approach which restricts the nuclear evolution to classical particle trajectories.

Quantum hydrodynamics and nuclear motion
In this section, we derive the quantum fluid picture for the mean-field model. We assume the effective Hamiltonian operator (3.3) may be computed from a Hamiltonian operator of the form where H e := − 2 ∆ x /2m + V e (x), while M and m denote the nuclear and electronic masses, respectively. Consequently, upon recalling (2.3), we have where we have suppressed the r−dependence for convenience of notation. Thus, upon denoting D = |χ| 2 and µ = J (χ), the mean-field Hamiltonian functional h(µ, D, ψ) := ψ|H ′ ψ reads This Hamiltonian functional is a mapping h : is understood to be the space of 1-form densities on R 3 . At this point, according to the procedure outlined in Section 2.1, we may perform the partial Legendre transform u = δh/δµ = M −1 µ/D, and write the mean-field Lagrangian in the following collective, or reduced, form: The reduced mean-field Lagrangian in (3.6) defines a map ℓ : T H e × (X(R 3 ) × Den(R 3 )) → R which can be regarded as a mixed hydrodynamic/phase-space Lagrangian. If the term corresponding to the quantum potential were simply discarded in taking the classical restriction of neglecting 2 in the nuclear dynamics, the reduced mean-field Lagrangian (3.6) would become, An analogous result could be obtained by following an alternative procedure which would exploit the density operator formalism for the nuclear dynamics, as indicated in Section 2.5.
In this case, one would obtain the same Lagrangian (3.7), although with ∇ × u = 0. At this point, one can apply Hamilton's variational principle by taking arbitrary variations δψ and constrained variations (2.10) for u and D. In general, a Lagrangian of this type yields the following Euler-Poincaré equations of motion [42] (∂ t + £ u ) δℓ δu = D∇ δℓ δD , (3.8) where £ u denotes the Lie derivative along the vector field u. These equations take the following forms, upon specializing to the mean-field Lagrangian (3.7): i ψ = H e +ˆDV I (x) d 3 r ψ . Again, as explained in Section 2.2, setting D(r, t) = δ(r − q(t)) and integrating (3.11) over space yields classical trajectories. Eventually, the corresponding classical system reads This classical restriction preserves the variational structure, whose Lagrangian is now given by L(q,q, ψ,ψ) = M|q| 2 /2 − V n (q) + ψ, i ψ − ( H e + V I (q, x))ψ .
Equations (3.14) represent the standard mean-field model as it is usually implemented in molecular dynamics simulations [63] (although here we have focused on the simplest case of one nucleus and one electron). As we can see in the previous equation, the classical-quantum coupling in this model occurs solely through the interaction potential V I .
Unfortunately, this quantum fluid picture of the mean-field model is not satisfactory in many cases, because the mean-field factorization (3.1) disregards the classical-quantum correlations between nuclei and electrons. A more advanced model capturing part of these correlation effects will be presented in the next section.

Exact factorization
First appearing in von Neumann's book [81] and later developed by Hunter [46], the following factorization ansatz has also been called exact factorization in recent work [1, 2, 6, 3]: Ψ(r, x, t) = χ(r, t)ψ(x, t; r) . (4.1) Here, the electron degree of freedom ψ depends parametrically on the nuclear coordinate r. This means that ψ is a smooth map ψ ∈ C ∞ (R 3 , H e ) from physical space to the Hilbert space H e = L 2 (R 3 ) of electronic wave functions. Furthermore, the factorization (4.1) invokes the partial normalization condition (PNC) which as a result of (4.1) ensures that´|Ψ(r, x, t)| 2 d 3 x = |χ(r, t)| 2 , so that D(r, t) := |χ(r, t)| 2 may be interpreted as the nuclear probability density.

Wave functions vs. density operators
The assumption of exact factorization (4.1) in the wave function transfers to the density operator, which is then written as where we have dropped the time-dependence for convenience of notation. Then, the physical electron density operator is given by the partial trace, written in matrix element notation as The corresponding nuclear density operator is in which we notice that the PNC (4.2) does not apply. This means the quantities χ and ψ are not true wave functions for the nuclei and electrons (which may not even exist in the presence of decoherence, i.e. quantum mixing). However, we shall continue to refer to them as such, because they retain certain mnemonic relationships. We remark that expectation values of nuclear observables involve integration over the r-parameters of the electron "wave functions". More specifically, the expectation value A n for a nuclear observable A n (r, r ′ ) is given by (again, in matrix element notation) A n :=˜d 3 r d 3 r ′´d3 x χ * (r ′ )ψ * (x; r ′ )A n (r ′ , r)ψ(x; r)χ(r). As we shall see, this structure of the nuclear density operator leads to important consequences in the development of the exact factorization theory. At this stage, we shall only emphasize that all the relations above also apply naturally in the context of the Born-Oppenheimer approximation [47], thereby indicating again that the interpretation of nuclear and electronic motion in terms of genuine wave functions needs to be revisited. For example, backreaction effects generated by the presence of ψ in (4.5) can lead to nuclear decoherence effects since indeed one has ρ 2 n = ρ n . This is a general feature of classical-quantum coupling [11], which in fact erodes purity in both the classical and the quantum subsystems.

General equations of motion
In this section, we generalize the approach to exact factorization used in [1,2,6,3] by allowing for an arbitrary Dirac Hamiltonian functional h(χ, ψ) and thereby extending the treatment in Section 1.2. Thus, inserting the ansatz (4.1) in the Dirac-Frenkel Lagrangian and then enforcing the PNC in (4.2) via a Lagrange multiplier λ(r, t) gives where we have introduced the notation: ψ 1 |ψ 2 (r) = ψ † 1 ψ 2 (r) =´ψ * 1 (x; r)ψ 2 (x; r) d 3 x, to denote the natural L 2 inner product on H e . Naturally, the λ equation enforces the PNC, and computing the χ Euler-Lagrange equation yields Consequently, upon using ∂ t |χ| 2 = Im(χ * δh/δχ), we can compute the ψ equation, as Upon taking the real part of the inner product of this equation with ψ, one obtains that λ = ψ, (δh/δψ) /2 − |χ| 2 ψ, i ∂ t ψ . Analogously, taking the imaginary part of the inner product of equation (4.8) with ψ yields the following compatibility condition, Im ψ δh δψ = Im χ * δh δχ . (4.9) In conclusion, this produces the final form of the ψ equation as in terms of the notation ψ † · := ψ| · . Before closing this section, we remark that the difference between equation (4.9) and the compatibility condition in (1.10) is its non-zero right-hand side, which arises because the inner product is taken only over the electronic degrees of freedom. That is, equation (4.9) depends on the nuclear coordinate r.

Local phases and gauge freedom
One may observe that the exact factorization (4.1) is defined only up to compensating local phase shifts of the nuclear and electronic wave functions. Namely, the replacements ψ(x, t; r) → ψ ′ (x, t; r) = e −iθ(r,t) ψ(x, t; r) χ(r, t) → χ ′ (r, t) = e iθ(r,t) χ(r, t) (4.11) leave the EF product wave function Ψ(r, x, t) = χ(r, t)ψ(x, t; r) invariant for an arbitrary local phase θ(r, t). This is a typical example of gauge freedom in a field theory.
To specify the evolution completely, one may either transform to gauge invariant functions, such as the electric and magnetic fields in electromagnetism, or one may fix the gauge by imposing one condition per degree of gauge freedom. Gauge fixing is not always the best option, though, because it may obscure physical effects arising due to local breaking of gauge symmetry. An example is the Berry phase, which arises from locally breaking the gauge symmetry of phase shifts in nonrelativistic quantum mechanics [10]. See also [83] for broadly ranging discussions of geometric phases in physics.
Another convenient choice consists in fixing ψ|i ∂ t ψ = 0, so that the ψ equation in (4.10) becomes 2i |χ| 2 ∂ t ψ = δh/δψ − ψ|(δh/δψ) ψ. This gauge is called the temporal gauge (or Weyl gauge) in electromagnetism and it has been adopted recently in [1,4,75]. We remark that gauge theory is also important in other aspects of chemical physics; for example, see [57] for applications of gauge theory in molecular mechanics.

The Hamiltonian functional
We now return to the Hamiltonian operator H in (1.1), written as the sum in (1.3), where H e is the electron Hamiltonian operator defined in the Introduction. Henceforth, we will suppress the subscript r, and write ∇ r simply as ∇. In this simplified notation, we define the Berry connection [10] as where the notation A(r, t) suggests that this quantity plays a role as a gauge field, analogous to the magnetic vector potential in electromagnetism. Here, we recall that · | · denotes the natural L 2 inner product on H e and · denotes the corresponding norm (whose values again depend on the nuclear coordinate r).
We also define the effective electronic potential, ǫ(ψ, ∇ψ), given by (4.13) The last term in (4.13) is the trace of the real part of the complex quantum geometric tensor where we denote A j := ψ| − i ∂ j ψ . The imaginary part of Q ij is proportional to the Berry curvature B ij := ∂ i A j − ∂ j A i ; namely, 2 Im(Q ij ) = B ij [51]. The emergence of the trace of its real part, T ij = Re(Q ij ), (4.15) in the electron energy in (4.13) indicates the geometry underlying the present formulation.
Notice that the interpretation of ǫ(ψ, ∇ψ) in (4.13) as an effective electronic potential departs slightly from that found in the literature, where this quantity is called the gauge invariant part of the time-dependent potential energy surface [4,5,75]. After some manipulation involving completing the square and integration by parts, the Hamiltonian for exact factorization in [1,2] can be expressed as This formula for the EF Hamiltonian agrees with the the chemical physics literature; see, e.g., [75,4,5]. For the Hamiltonian (4.16), one computes so that the Euler-Lagrange equations (4.7)-(4.10) specialize to, These equations also agree with the results found in the recent chemical physics literature [1,2]. In this case, one can show that λ vanishes identically (in agreement with [3]) and that the compatibility condition (4.9) is indeed satisfied.

Hydrodynamic approach
In this section, we again transform into the hydrodynamic picture for the nuclear wave function χ. Upon applying the same procedure as in the mean-field case, now denoting collective fluid variables as D = |χ| 2 and µ := J (χ) = Im(χ * ∇χ), the Hamiltonian functional (4.16) reads At this stage, we will treat the quantity ψ in (4.20) as a parametric variable, whose variations will be taken independently of those for µ and D. Applying the partial Legendre transform in this case yields u = δh/δµ = M −1 (µ + DA)/D, with A defined in terms of ψ in (4.12), and we write the EF Lagrangian in the following form: We now apply Hamilton's variational principle by taking arbitrary variations in δψ and δλ, and Euler-Poincaré variations (2.10) for u and D. As before, a Lagrangian of this type yields general equations of motion (3.8)-(3.10), this time along with a standard Euler-Lagrange equation in λ, which naturally recovers the PNC. Then, the equation for the electronic wave function ψ reads where we have introduced the functional which we will call the electronic Hamiltonian functional. Upon making use of the effective electronic potential ǫ(ψ, ∇ψ) defined in (4.13), one obtains the functional derivative (4.24) whose insertion into (4.22) yields Again, we obtain a compatibility condition from varying the Lagrange multiplier λ in (4.21). In this case, the relation (4.9) becomes Im ψ δF δψ = 0, (4.26) as one can see by replacing the Hamiltonian (4.20) in (4.9). In addition, one can show that F is invariant under local phase transformations. As we will see in Section 5, this local U(1) invariance will ultimately lead to a density matrix formulation of the electronic dynamics. The Lagrangian (4.21) yields the electron dynamics (4.22), as well as the following Euler-Poincaré equations of nuclear hydrodynamics. Specifically, the Euler-Poincaré variations (2.10) yield (4.28) Here, V Q denotes the nuclear quantum potential (2.9) for the nuclei, while B := ∇ × A is the Berry curvature, effectively a magnetic field generated by the electrons, and E := − ∂ t A − ∇ ψ|i ∂ t ψ plays the role of an electric field generated by the electrons. In this analogy, in the Berry connection defined in equation (4.12) plays the role of the coupling constant (charge) in the electromagnetic force on a charged particle. The quantity ψ|i ∂ t ψ in the definition of E can be fixed by selecting a particular gauge; for example, possible options are the temporal gauge, ψ|i ∂ t ψ = 0, and the hydrodynamic gauge ψ|i ∂ t ψ = A · u. A more explicit expression of E can be found by using (4.22), thereby leading to the equation where we recall the notation · , · = Re · | · in (1.7) for the real-valued pairing. Then, we have This may be written in components, in terms of the real part of the quantum geometric tensor Then, the equations (4.27)-(4.28) and (4.22) take the form, We emphasize that the equations of motion (4.31)-(4.33) could also be derived from the coupled Schrödinger equations (4.18)- (4.19). Indeed, (4.31)-(4.32) can be derived from (4.18) by writing χ = √ De iS/ and by finding the evolution for u = (∇S + A)/M. In addition, (4.33) is equivalent to (4.19), as it can be verified by replacing (4.17) in (4.10).
Notice that equation (4.31) does not conserve the spatial integral of the nuclear momentum density, Here, S is the local phase of the wave function χ = √ D e iS/ and D = |χ| 2 , while DA is part of the nuclear momentum density. This non-conservation of the hydrodynamic momentum should come as no surprise. In fact, this is already apparent in the original system (4.18)-(4.19) which, instead, conserves the total momentum where P e = ψ| − i ∇ x |ψ . Thus, the total motion has two momentum contributions: one from the r-gradient and the other from the x-gradient. In particular, the momentum density of the nuclei is given by ´Ψ * (r, x)(−i∇ r )Ψ(r, x) d 3 x = µ+DA where A is the Berry connection defined in equation (4.12), and µ := J(χ) = Im(χ * ∇χ).

Newtonian limit and Lorentz force
The Newtonian limit performed in the chemical physics literature neglects the order quantum potential term in the Lagrangian (4.21) and varies the remainder, treating the nuclei as classical particles. The corresponding dynamics may be obtained by neglecting V Q in (4.27) and by replacing Mu = ∇S + A. One then verifies that this process leads to the following Hamilton-Jacobi equation: ∂ t S + M −1 |∇S + A| 2 /2 = ψ|i ∂ t ψ −ǫ(ψ, ∇ψ), which corresponds to charged particle motion in a Maxwell field, governed by the equation [5] Mq = −E −q × B − ∇ǫ(ψ, ∇ψ) . (4.36) The same result can be obtained by setting D(r, t) = δ(r − q(t)) in (4.28) to produceq(t) = u(q(t), t). Then, after multiplying (4.27) by D(r, t), integration of the delta function over physical space returns the classical equation (4.36) above. An important point here is that the customary operation in chemical physics of neglecting the quantum potential term in the Lagrangian (4.21) can be problematic. Normally, this step would invoke the limit 2 → 0. However, here this process would also lead to discarding the terms M −1 ( 2 ∇ψ 2 − A 2 )/2 in the effective electronic potential (4.13), thereby taking the exact-factorization model into the standard mean-field theory. This crucial issue will be resolved in Section 5 by performing the exact factorization at the level of the molecular density operator.
We should also comment on the Lorentz force appearing in (4.36). We notice that the combination of this electromagnetic-type force with the potential energy contribution ∇ǫ suggests that the conventional picture of nuclei evolving on potential energy surfaces fixed in space may be oversimplified. As we shall see, despite the claims made in [76], the force E + u × B cannot vanish without requiring major modifications of the electron energy function ǫ(ψ, ∇ψ). Such modifications would result in singular solution behaviour for the Berry curvature that is unexpected for the exact factorization model.
In the present context, one may regard the assumption of E+u×B = 0 as an incompressible magnetohydrodynamic (MHD) approximation of the quantum fluid-plasma equations in (4.27) and (4.28). Setting D = 1, for which (4.28) implies divu = 0, and taking the curl of (4.29) yields Vanishing of the second term on the right side of equation (4.37) would require a functional relation between ψ and δF/δψ. To see the implications of this requirement, we specialize to the 2D case for which the vector A(r, t) := ψ| − i ∇ψ lies in the plane, so that B = curlA = Bẑ points normal to the plane. Hence, equation (4.37) becomes whereẑ · ∇f × ∇h = J(f, h), is the Jacobian between functions f and h on the (x, y) plane. As expected, vanishing of the right hand side of equation (4.38) requires a functional relation between ψ and δF/δψ. This occurs, for example, in the Ginzburg-Pitaevskii description of Bose-Einstein condensation (BEC) [70], in which case the B-equation (4.38) admits singular Berry-curvature solutions of the form B(x, y, t) = k Γ k δ(x − x k (t))δ(y − y k (t)), for constants Γ k [27]. The emergence of this type of singular behaviour in the solutions of equation (4.38) is precluded by the dependence on both ψ and ∇ψ of the electron energy function ǫ(ψ, ∇ψ) in equation (4.23) for the exact factorization model.

Circulation dynamics for the Berry connection
The dynamics of the Berry connection A(r, t) := ψ| − i ∇ψ in (4.12) governs the circulation of the nuclear fluid around closed loops which move with the velocity u = (µ/D + A)/M. To see this, we write the motion equation (4.31) as the Lie derivative of a circulation 1-form, (4.39) where we have used the relation (4.34) in the second line and d denotes spatial differential in r. Equating the right hand sides of equation (4.39) and integrating around an arbitrary closed loop γ(t) moving with the nuclear flow velocity u(r, t) annihilates the exact differential terms and produces the following circulation dynamics for the Berry connection: This means that the nuclear circulation integral¸γ (t) A · dr, interpreted as the Berry phase obtained by integration around a loop moving with the nuclear fluid, is generated dynamically by an interplay between nuclear and electronic properties. Likewise, the evolution of the Berry curvature B = curlA follows by applying the Stokes theorem to relation (4.40). Thus, the flux of the Berry curvature through a surface S whose boundary ∂Σ = γ(t) is a closed loop moving with the nuclear fluid (simply known as the Berry phase) satisfies In terms of the real and imaginary parts of the quantum geometric tensor, Q ij in (4.14), this becomes Equation (4.42) expresses the quantum geometric mechanics of the correlated nuclear and electronic degrees of freedom in the EF model. Namely, the nuclear probability density D and expectation of the gradient ∇ H e are coupled dynamically with the real and imaginary parts of the quantum geometric tensor, Q ij .

Electron dynamics in the nuclear frame
So far, we have presented the geometric aspects of the exact factorization model which are currently available in the literature. This section takes a step further to consider the evolution of the electron density matrix. Upon following the arguments in [5], we shall write the electron dynamics in the Lagrangian frame moving with the nuclear hydrodynamic flow.
Recalling the notation H e for the electronic Hilbert space, we begin by introducing the group, C ∞ (R 3 , U(H e )), of smooth mappings from the physical space into the unitary group of the electronic Hilbert space. Then, we make the following evolution ansatz for the electronic wave function: ψ(t) = (Uψ 0 ) • η −1 , or more explicitly ψ(x, t; r) = U(η −1 (r, t), t)ψ 0 (x; η −1 (r, t)) , (4.43) where η is the nuclear hydrodynamic path obeyingη = u(η, t) and U(r, t) ∈ C ∞ (R 3 , U(H e )) is a local unitary operator on H e . The evolution ansatz (4.43) results in the following equation for the time evolution of ψ: where we have defined ξ := (UU −1 )•η −1 . Upon substituting these relations into the Lagrangian (4.21), one finds At this point, we shall prove that the electron energy function ǫ(ψ, ∇ψ) can be written uniquely in terms of the density operator ρ = ψψ † and proceed for the rest of this Section by employing a density matrix description of the electronic dynamics. Indeed, we find by direct calculation that where the tensor T with components T ij = Re(Q ij ) is the real part of the quantum geometric tensor in equation (4.14). Regardless of the minus sign in equation (4.14) the term in the parentheses on the right side of equation (4.46) is positive, since the left side is positive. Relations resembling (4.46) also appear in standard quantum mechanics when writing the Fubini-Study metric on the projective Hilbert space PH , see e.g. [26]. Hence, we conclude that the electron energy function (4.13) may be expressed as, where one defines A|B := Tr(A † B) by using the generalized trace. In matrix element notation, one has A|B : The key formula (4.47) enables us to write the previous Euler-Poincaré Lagrangian equivalently as At this point, along with (2.10), one finds the variational relations with u ∈ X(R 3 ), D ∈ Den(R 3 ), ξ ∈ C ∞ (R 3 , u(H e )) and iρ ∈ C ∞ (R 3 , u(H e )). Namely, cf. [33,31,34,39] (4.50) Remark 4.1 (Analogies with complex fluids) We take this opportunity to make the connection between the hydrodynamic exact factorization system and previous investigations of the geometry of liquid crystal flows, as found in [32,78,33,31,34,39]. In this comparison, the electronic wave function ψ(r, t) ∈ C ∞ (R 3 , H e ) is replaced by the director, an orientation parameter field n(r, t) ∈ C ∞ (R 3 , S 2 ); the unitary evolution operator U(r, t) ∈ C ∞ (R 3 , U(H e )) becomes a rotation matrix R(r, t) ∈ C ∞ (R 3 , SO (3)); and one still considers the coupling to the fluid velocity u(r, t) given by the action of diffeomorphisms η ∈ Diff(R 3 ). Indeed, with these replacements, one has a reduced Lagrangian of the same type, ℓ(u, ξ, D, n), and the resulting Euler-Poincaré equations are equivalent to those in (4.50).
For convenience, we rewrite the electronic Hamiltonian functional in (4.23) as F (D, ρ) = D( ρ| H e + M −1 2 ∇ρ 2 /4) d 3 r. Consequently, the fluid velocity equation in (4.50) becomes: which indeed reduces to equation (4.31) upon specializing to pure states ρ = ψψ † . Also, we notice the following analogue of relation (4.30): div(D∇ρ), ρ . This result implies that spatially uniform pure initial states (such that ρ 2 e = ρ e ) become mixed states as time proceeds. Thus, in agreement with, e.g., [64], the exact factorization model captures electronic decoherence effects (that is, quantum state mixing) from pure initial states; since the density matrix evolution is no longer unitary.
Upon collecting equations, now we may specialize the general system (4.50) to our case as ∂ j D∇ρ, ∂ j ρ , (4.53)

Hamiltonian structure
The Hamiltonian structure of equations (4.53) can be conveniently rewritten in terms of the quantities m := MDu ,ρ := Dρ , (4.54) so that the total energy (4.20) reads Upon restricting to pure quantum states, for which ρ, ∇ρ = 0 and ρ, ∇ρ = D∇D, we find Here, the minus sign arises as in the calculation in (4.46). The term in the parentheses on the right side of this equation is positive, though, since the left side is positive. Consequently, for pure quantum states, the Hamiltonian (4.55) may be expressed as The appearance of the amplitude of the probability density D in the denominators of the 2 terms in the Hamiltonian in (4.57) implies that the dynamical effects of the quantum terms in (4.57) need not decrease, as the squared amplitude of the wave function decreases. In terms of the variables (m, D,ρ) defined in (4.54), the system of equations in (4.53) may be rewritten equivalently, as follows: where h(m, D,ρ) is taken to be the Hamiltonian in (4.57), and homogeneous boundary conditions are assumed, under integration by parts. The change of variableρ → i ρ shows that this bracket is Lie-Poisson on the dual of the following Lie algebra L comprising a direct sum of semidirect product actions: where in (4.61) the angle brackets · , · r denote L 2 pairing in the r coordinates, and the square brackets denote the components of the adjoint action of the semidirect-product Lie algebra L, whose r-coordinate pairings are given explicitly in equation (4.60).

Density operator factorization and singular solutions
Our discussion in Section 4.7 shows that the only contribution to the circulation of the hydrodynamic flow arises from the Berry connection associated to the electronic function ψ. This is due to the fact that the hydrodynamic velocity M −1 µ/D = Im(χ * ∇χ)/|χ| 2 is an exact differential and therefore has zero vorticity. However, we showed in Section 2.5 that this restriction can be relaxed by considering density operators. In the same section we also showed that the classical closure of mixed state dynamics allows for multi-particle trajectories arising from the initial condition (2.30), which in turn is not compatible with the standard QHD definition D 0 = |χ 0 | 2 . In this section we shall include mixed state dynamics by extending the exact factorization model to a density operator formulation.

Factorization of the molecular density operator
In order to generalize the exact factorization (4.1) to density operators, we recall the relation (4.3) and extend it to consider a molecular density operator of the form In order to formulate the dynamical equations for such a factorization ansatz, it is convenient to consider a sequence {Ψ k (r, x)} and exploit (2.33) to write where k w k = 1 and each Ψ k satisfies a separate (uncoupled) Schrödinger equation with Hamiltonian (1.1). We remark that (5.2) is the equivariant momentum map for unitary transformations of the {Ψ k } recently studied in [77]. Correspondingly, the overall dynamics of {Ψ k } is produced by the variational principle At this stage, we consider a ψ(x; r) satisfying (4.2), and we restrict to the case Ψ k (r, x) = χ k (r)ψ(x; r). Consequently, the ansatz (5.1) is recovered by setting ρ n (r, r ′ ) = k w k χ k (r)χ * k (r ′ ) and the above Lagrangian becomes Here, the Lagrange multiplier enforces ψ 2 = 1 and the Hamiltonian reads as h({χ k }, ψ) = k w k˜χ * k (r)ψ * (x; r)Ĥ(χ k ψ)(r, x) d 3 r d 3 x. Now, we restrict ρ n to undergo unitary evolution by writing, for each k, χ k (r, t) = U n (t)χ (0) k (r), where U n ∈ U(H n ) is a time-dependent unitary operator on the nuclear quantum space H n = L 2 (R 3 ). Then, the Lagrangian above in (5.4) can be rewritten in terms of the nuclear density matrix ρ n (r, r ′ ) and its diagonal elements D(r) := ρ n (r, r) as ℓ( ξ , ρ n , ψ,ψ) = Reˆ i ρ n (r, r ′ )ξ(r ′ , r) + D(r) ψ|i ψ + λ ψ 2 − 1 dr − h(ρ n , ψ) , (5.5) where ξ :=U n U −1 n , ρ n = U n ρ (0) n U −1 n , and the Hamiltonian is: in which ǫ(ψ, ∇ψ) is the same as in (4.13) and P denotes the nuclear momentum operator. The Lagrangian (5.5) produces dynamics for ρ n which is identical to the dynamics for χ(r)χ * (r ′ ) emerging from the equation (4.7). In this generalized case, the compatibility condition (4.9) is replaced by Im ψ(r) δh δψ (r) = 2Imˆρ n (r, r ′ ) δh δρ n (r ′ , r) d 3 r ′ . (5.7) Now that the variational principle and the Hamiltonian functional are completely characterized, we may proceed by restricting nuclear dynamics to undergo classical motion. To this purpose, we combine the two approaches described in Sections 2.4 and 2.5, by applying the regularization technique after performing the classical closure for density operators.

Classical closure and singular solutions
In following the discussion in Section 2.5, we wish to collectivize the Hamiltonian in terms of the hydrodynamic quantities (µ, D). This can be achieved by applying the closure (2.37) to ρ n in Section (5.1), that is ρ n (r, r ′ ) = D(r/2 + r ′ /2) exp(iM(r − r ′ ) · v(r/2 + r ′ /2)/ ), with MD(r)v(r) = µ(r). Lengthy but straightforward computations using matrix elements (or, by Wigner transforming, direct applications of Weyl calculus) eventually take the Hamiltonian (5.6) into the hydrodynamic form h(µ, D, ψ) = 1 2Mˆ| which coincides with the Hamiltonian (4.20), except for the quantum potential term. Now however, note that the nuclear hydrodynamic variables (µ, D) are no longer given in terms of a unique wavefunction, so that curl(M −1 µ/D) = 0. Thus, equation (4.40) becomes This means that the dynamics of the nuclear circulation integral¸γ (t) u · dr is now interpreted as a genuine hydrodynamic Kelvin theorem for the circulation around a loop moving with the nuclear fluid with Eulerian velocity u := M −1 (µ/D + A).
At this point, one can Legendre transform the Hamiltonian (5.8) and follow the treatment in Section 4.8 to express the electron function ψ in the nuclear frame. Then, rearranging yields the new Hamiltonian which again coincides with (4.55) except for the quantum potential, although a similar potential (the last term) with opposite sign is produced according to the relation (4.56). The equations of motion associated to the Hamiltonian (5.10) can now be easily formulated by applying the Lie-Poisson bracket structure (4.60). The last two equations in (4.58) do not change, although the force arising from the quantum potential is modified accordingly in the momentum equation.
Remarkably, Hamiltonian systems with Lie-Poisson brackets of the type (4.60) for the exact factorization (EF) model have already been studied and their geodesic motions have been shown to admit singular solutions for certain choices of quadratic Hamiltonians, in [45]. Following this previous work, in this Section we will extend the momentum map (2.25) which led to the Bohmion singular solutions in Section 2.4, to a singular solution momentum map for a regularized version of the Lie-Poisson system associated to the Hamiltonian (5.10). Again, in analogy with Section 2.4, here we shall present two different treatments: while the first is uniquely based on the Hamiltonian structure, the second exploits the associated variational principle.

Hamiltonian regularization in 1D
For further simplification, we consider a one dimensional nuclear coordinate. Upon considering the Lie-Poisson structure given by the Hamiltonian (5.10) and the bracket (4.60), we can write the dynamical variables according to the momentum map given as follows [45]: ̺ a (t)δ(r − q a (t)) .

(5.11)
Here, the D−equation in (4.58) implies that the w a with a = 1, 2, . . . , N in (5.11) are all constant. In contrast, as we shall see, theρ−equation for the regularized Hamiltonian obtained from introducing smoothed variables in (5.10) will imply an evolution equation for the coefficients ̺ a (t) = ϕ a (t)ϕ † a (t) in (5.11). The corresponding singular solutions represent Bohmions in the density matrix formulation.
In addition to the last term in the Hamiltonian (5.10), which was already regularized in Section 2.4, a further barrier to singular solutions is represented by the term involving the gradient ofρ. Hence, we again smoothen our variables by replacing h(m, D,ρ) → h(m,D,ρ) form = K * m,D = K * D, andρ = K * ρ, in the collectivized Hamiltonian (5.10). Then, the corresponding regularized Hamiltonian h REG for the singular solutions is given by Thus, substitution of the singular solutions (5.11) into the regularized Hamiltonian (5.10) yields the Hamiltonian (5.12) in terms of its canonical phase space variables (q a , p a ) with a = 1, 2, . . . , N , augmented as expected by the equation for the matrix ρ a , given by substitution of the last equation in (5.11) into the middle equation of (4.58) as where we have introducedH e (q a ) =´H e (r) K(r − q a ) dr and we do not sum on the index a.
As for the QHD Bohmions in Section 2.4, the equivariance of the momentum map (5.11) discovered in [41,45] implies that the dynamics of ({q}, {p}) satisfy the canonically conjugate Hamiltonian equationṡ q a = δh REG δp a = u(q a (t), t) andṗ a = − δh REG δq a , for a = 1, 2, . . . , N, (5.14) where K ′ denotes the derivative of the smoothing kernel K with respect to its argument. As remarked in Section 2.4 for QHD, these singular solutions for the regularized dynamics extend the 'peakon' singular solutions for nonlinear wave equations in [20,43]. As such, the peakon-like equations in (5.11) do not possess the usual Newtonian form. This is because of the smoothing which was introduced in the kinetic energy term. Again, the term in the canonical equations (5.14) arising from summands in h REG in (5.12) proportional to 2 provides extensive, potentially long-range coupling among the singular particle-like solutions, because of the presence ofD in the denominators of these terms in the Hamiltonian h REG in (2.22). However, as in the previous section, the limit 2 → 0 in the canonical Bohmion equations (5.14) in the density matrix formulation is regular.

Lagrangian regularization in 1D
At this point, we perform the analogous procedure to that in the last part of Section 2.4 now for the above model obtained by factorizing the molecular density matrix. As before, instead of regularizing all terms as in the Hamiltonian approach, we shall regularize only the O( 2 )−terms. We now consider the cold-fluid closure of the Lagrangian (4. Now, the evolution ofρ in terms of the Lagrangian path η and the electronic propagator U(r 0 ) is given asρ = η * ρ =´ρ(r 0 , t) δ(r − η(r 0 , t)) dr 0 , withρ(r 0 , t) = U(r 0 , t)ρ 0 (r 0 )U † (r 0 , t), so that (5.17) yields ρ(r, t) = N a=1 ̺ a (t)δ(r − q a (t)) , with ̺ a (t) := U(q (0) a , t)̺ (0) a U † (q (0) a , t) .
Here, we set ̺ (0) a = ϕ (0) a ϕ (0) a † so that ̺ a (t) = ϕ a (t)ϕ † a (t) is a projection at all times. Furthermore, we evaluatê where ξ a (t) := ∂ t U(q (0) a , t) U † (q (0) a , t) , and we have recalled q (0) a = η −1 (q(t), t) as well as ξ(r) = ∂ t U(η −1 (r), t) U † (η −1 (r), t). Then, insertion of (2.30) and (5.17)  While the Euler-Lagrange equations for the trajectories q a are obvious, the electronic dynamics is obtained as an Euler-Poincaré equation by using the variational relation δξ a = ∂ t ν a − [ξ a , ν a ], with ν a arbitrary and vanishing at the endpoints. As usual, this is obtained by an explicit calculation using the definition of ξ a (t). Eventually, we obtain the following quantum equation which coincides with (5.13) except thatH e is now replaced by the unfiltered operator H e . The comparison of solution behaviour between these two regularized Bohmion models in the density matrix formulation will be discussed elsewhere by using computer simulations.

Conclusions
In this paper, we have exploited momentum maps to collectivize a sequence of molecular quantum chemistry models for factorized nuclear and electronic wave functions, thereby obtaining a sequence of quantum fluid models with shared semidirect-product Lie-Poisson structures. After reviewing the Born-Oppenheimer product of nuclear and electronic wave functions, we started with mean-field theory, and then passed to a recent development called 'exact factorization' (EF) for nonadiabatic correlated electron-nuclear dynamics, which has been reported to describe decoherence of pure electron quantum states into mixed states. In the last part, we extended the exact factorization approach to apply for density operators. In Section 2.4, we mollified the weakly convergent WKB → 0 limit by applying a smoothing operator to the quantum variables in the collectivized Hamiltonian for regularized quantum hydrodynamics. This smoothing operation preserved the Hamiltonian structure of the quantum fluid model and it resulted in the discovery of singular delta function solutions called 'Bohmions' for smooth quantum fluid Hamiltonians. Depending on which terms are regularized in the Hamiltonian, different sets of Bohmion equations are available.
In the development of the paper, we showed that the Hamiltonian formulation of the collectivized quantum fluid equations for EF possesses the same Lie-Poisson bracket structure as in earlier work on perfect complex fluids (PCF), such as liquid crystals, [32,78,33,31,34,39]. The parallel between EF and PCF is that the nuclear fluid velocity vector field Lie-transports both the nuclear probability density and the electron density matrix, while the latter also has its own unitary dynamics in the moving frame of the nuclear fluid. This picture was also extended to present a new PCF dynamical model based on the factorization of the molecular density operator.
In the PCF formulation of the nonadiabatic electron problem, smoothing all terms in (5.10) yields singular momentum maps corresponding to the 'peakon' solutions of the well known EPDiff equation [41]. In one spatial dimension, this more general class of Bohmions is governed by a countably infinite set of canonical Hamiltonian equations in phase space, in analogy to the solitons for the Camassa-Holm equation [20]. The countably infinite phase space system can be truncated to a multi-particle phase space system at any finite number of Bohmions, because the Hamiltonian dynamics does not create new Bohmions. In fact, the Bohmion collectivized solutions discussed in Sections 2.4 and 5.2 comprise a semidirect-product version of the class of 'peakon' solutions for the CH equation [20] which arise from the well-known singular momentum map for the entire class of EPDiff equations [41].
A second approach to Bohmion dynamics was presented in Section 5.2.2. This approach was developed in the variational framework by smoothening only the O( 2 )−terms in the Lagrangian (5.15). Although the analogy to the peakon solutions of the Camassa-Holm equation no longer holds entirely, the resulting dynamical system still consists of a countable set of finite dimensional Hamiltonian equations.
Future work will take further advantage of the analogy between continuum dynamics and the collectivization of quantum dynamics via momentum maps. For example, products of delta functions in different spaces can be introduced, corresponding to Bohmion dynamics for the different factorized wave functions of many interacting molecules. This approach is reminiscent of the closure models arising in time dependent Hartree (TDH) theory [37] for quantum dynamics in nuclear physics. Approaches such as these have long been applied in several fields of science, including molecular chemistry, nuclear physics, and condensed matter physics, as well as in celestial mechanics, in hopes of lifting the "curse of dimensions" which tends to be ubiquitous in many-body problems [14].
The third step uses Jacobi's formula for the derivative of the determinant. Next, we show that ξ = i −1 { u k , P k }/2 thereby recovering (2.35). This is verified as follows: where we have used the matrix elements Thus we have proved that the infinitesimal generator for the right action (A.1) is indeed given by (A.2), withξ = i −1 { u k , P k }/2. Correspondingly, the infinitesimal generator for the left action (2.34) is given by ξ(ρ) = −[ξ, ρ], which then proves (2.35).