Ehrenfest statistical dynamics in chemistry: study of decoherence effects

In previous works, we introduced a geometric route to define our Ehrenfest Statistical Dynamics (ESD) and we proved that, for a simple toy-model, the resulting ESD does not preserve purity. We now take a step further: we investigate decoherence and pointer basis in the ESD model by considering some uncertainty in the degrees of freedom of a simple but realistic molecular model, consisting of two classical cores and one quantum electron. The Ehrenfest model is sometimes discarded as a valid approximation to non-adiabatic coupled quantum-classical dynamics because it does not describe the decoherence in the quantum subsystem. However, any rigorous statistical analysis of the Ehrenfest dynamics, such as the described ESD formalism, proves that decoherence exists. In this article, decoherence in ESD is studied by measuring the change in the quantum subsystem purity and by analysing the appearance of the pointer basis to which the system decoheres, which for our example is composed by the eigenstates of the electronic Hamiltonian.


Introduction
The Schrödinger equation for a combined system of electrons and nuclei is generally too complex and involves too many degrees of freedom to be solvable, neither analytically nor by numerical methods. Approximations need to be made, one of the most important and successful being the classical approximation for a subset of the particles. Hybrid quantum-classical dynamical (HQCD) models are therefore necessary and widely used [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. In a previous work [2] we discussed how these HQCD models are built. Most approaches can be described in two steps: first, a partial deconstruction of the quantum mechanics (QM) of the total system (electrons and nuclei) which simplifies the model, and then a reconstruction that aims to recover the essential properties of the total Schrödinger equation lost in the deconstruction process.
HQCD models in the literature present at least two levels of deconstruction. The first one, called Born-Oppenheimer molecular dynamics (BOMD), is far away from the total Schrödinger equation for electrons and nuclei, as electrons are assumed to remain in the ground state for all times. The second one, closer to the total Schrödinger equation, is called Ehrenfest Dynamics (ED). In ED the nuclei are still classical (as in BOMD) but the electrons are allowed to populate excited states. A recent review on the topic [15] discusses these two approaches. The deconstruction simplifies the model by forcing the separability of the quantum states of the nuclei (which are later considered as classical) and the electrons [17], even if this separability cannot be preserved exactly by the evolution of an interacting quantum system. Therefore the deconstruction ignores entangled states of the quantum nuclear and electronic degrees of freedom. One could naively conclude that this implies a unitary evolution for the electronic subsystem, thus being the cause for the preserved purity in the ED evolution. This reasoning is erroneous, as generic ED evolution is not unitary for the electronic subsystem [2]. Nevertheless, as proved in Theorem 1 of the cited work, ED does preserve purity, and no decoherent effects can be found.
The second step in the definition of HQCD models, the reconstruction of the dynamics, is much more complicated. Many different proposals tackle the difficulty of re-incorporating into HQCD models the essential properties of the total Schrödinger equation that have been lost. One of these properties is the decoherence phenomenon in the electronic subsystem. In the context of HQCD, decoherence is understood as the fact that the neglected wave functions in the classical limit (i.e. those of the nuclei) rapidly lose overlap in time, leading to the destruction of superpositions in the quantum subsystem, i.e. forcing the electronic wave functions into a mixture of pure states [3,5]. Such pure states form the so-called pointer basis for the quantum sybsystem [18][19][20].
Different approaches aim to address the reconstruction of quantum dynamics with quite different tools. In J. C. Tully's Trajectory Surface Hopping (TSH) algorithms [12], for example, the deconstruction goes to BOMD and the reconstruction proceeds by allowing the system to perform certain specially designed stochastic jumps between adiabatic states. These jumps cannot however be well justified from first principles. Another relevant algorithm, widely used in Molecular Dynamics (MD), is the decay of mixing formalism by Truhlar and coworkers [11,16]. In this method, the deconstruction stops at the ED and the reconstruction is developed by adding terms in the dynamics which introduce decoherence. In this formalism, one considers an ensemble of hybrid quantum-classical systems and computes the dynamics of the quantum subsystem using two components: one arising as the fully coherent solution to the Liouville-von Neumann equation and one, ad-hoc, that incorporates electronic decoherence (see expression (34)

below).
There are other studies and proposals to tackle decoherence effects. Among them, we would like to mention Bittner and Rossky [21], Neria and Nitzan [22], Schwartz and coworkers [23] and Subotnik [7,8]. Particularly in the last two approaches, one can find the idea of considering statistical mixtures of hybrid quantum-classical systems in order to study the problem of decoherence. For example, Subotnik used in his works the formalism of the partial Wigner transform, introduced in the context of MD by Nielsen and coworkers [24][25][26], to represent the hybrid quantum-classical system; by adding some extra variables into the picture, it is possible to describe in an efficient way the decoherence effects of some systems.
A completely different route is the one followed by Abedi and coworkers [27][28][29][30], which prescribes an exact factorisation of the total wave function. Then, the classical limit for the nuclear part can be taken, and recently the appearance of decoherence in the resulting dynamics has been analyzed [31].
In this work, we take a step back from these approaches, and examine how uncertainty in the initial conditions may affect Ehrenfest dynamics. In order to adress this issue, a formulation of dynamics in terms of statistical distributions is necessary. In previous works [1,2] we introduced a geometric route to define our Ehrenfest Statistical Dynamics (ESD). We now take a step further by investigating decoherence in ESD. It is important to notice that, due to the interaction between classical cores and quantum electrons in a molecular model, cores act as an environment to the quantum subsystem of the molecule. In this setting, it is thus possible to consider the decoherence hypothesis [18][19][20], by which the environment of a non-isolated quantum system selects a set of orthogonal projectors onto the Hilbert space of the system. In other words, the interaction described in Ehrenfest dynamics causes the existence of a pointer basis for the quantum subsystem.
In our previous work [2] we proved that, for a simple toy-model,the resulting dynamics does not preserve purity. Now, and after carefully defining the meaning of decoherence in the context of HQCD models, we apply ESD to a realistic model: a diatomic, isolated molecular system. We will see that decoherence does occur in ESD and we will compute the changes in purity and the appearance of pointer basis associated to the decoherence phenomenon.
In Section 2 we summarise the main properties of ED and the corresponding statistical model. We also detail in Section 2.3.2 how ESD defines an evolution for the quantum subsystem which does not preserve its purity; we define the corresponding decoherence-time and the pointer basis. Section 3 studies the dynamics of a ionised dimer subject to uncertainty in its initial conditions. By making use of ESD, we prove the appearance of decoherence and we determine the pointer basis for this example, which turns out to be the set of eigenstates of the electronic Hamiltonian. Conclusions and an outlook of the work are presented in Section 4.

The case of molecular systems: Ehrenfest model
In Ehrenfest Dynamics (ED), the particles in a molecule are split in two sets [32]: • N Q quantum particles, typically the most external (valence) electrons, • which are coupled to the N C cores that will be considered classical.
In the following, the coordinates of electrons and cores will be denoted by 3-dimensional vectors r 1 , . . . , r NQ and R 1 , . . . , R NC , respectively. The Hamiltonian operator for this molecular system is given by (atomic units will be used hereafter) with R = ( R 1 , . . . , R NC ) ∈ R 3NC , and where M J is the mass of the J-th core and ∇ 2 J is the Laplacian operator with respect to R J . The operator H e (R), called the electronic Hamiltonian of the molecule, depends parametrically on the core positions: where Z J is the charge of the J-th core and ∇ 2 j is the Laplacian operator with respect to electronic coordinates r j .
The first approximation to the solution of the otherwise unassailable Schrödinger equation for this Hamiltonian is the separability of electrons and cores: with |ξ representing the state of the cores and |ψ the electronic state. This separability assumption ignores the possibility of electron-nuclear entanglement [33]. From these separable states, the next step is the substitution of the wave function for the cores by single classical trajectories, characterised by positions R = ( R 1 , . . . , R NC ) ∈ R 3NC and momenta P = ( P 1 , . . . , P NC ) ∈ R 3NC variables. The states in the hybrid quantum-classical model are then represented by elements in the following spaces: • the phase space of the N C classical cores corresponds to • the quantum Hilbert space H describing the states of the most external electrons.
The original Schrödinger equation of the full quantum model is approximated by a set of coupled differential equations when written on the hybrid quantum-classical variables introduced above [13,14,17,[34][35][36]. These are called the Ehrenfest equations. For a system composed of a set of N C classical particles, described by points ξ = (R, P ) ∈ M C , and a quantum subsystem, described by a wavefunction |ψ ∈ H, the Ehrenfest equations are: 2.2 Hybrid quantum-classical Ehrenfest statistical model

Hybrid mechanical systems
In previous works [1,2], we introduced a geometric formulation of ED for hybrid quantum-classical models and its extension to Statistical Mechanics. The first step was to adapt the geometrical formulation of Quantum Mechanics [34,[37][38][39][40][41][42] for the description of hybrid quantum-classical systems.
The main idea of the construction is to combine the classical and the quantum degrees of freedom in a form similar to the composition of classical systems. This is achieved by providing a description of the quantum subsystem formally equivalent to the description of the classical one. Let us take any basis {|e j } j of the Hilbert space H and consider the real and imaginary parts of the corresponding set of coordinates: The real numbers (q 1 , p 1 , q 2 , p 2 , · · · ) can be understood as coordinates on a real differentiable manifold 2 M Q , whose points are in one-to-one correspondence with the vectors in H. This is a Kähler manifold [34,[37][38][39][40][41], with a Poisson bracket {·, ·} Q and a symmetric product of functions (·, ·) Q determined by the properties of the Hermitian product in H. Their coordinate expressions are Schrödinger's equation can be written as a system of Hamilton equations with respect to the quantum Poisson bracket {·, ·} Q [34,[37][38][39][40][41]. Therefore, a purely quantum system can be considered as a "classical" Hamiltonian system with respect to the symplectic structure on the manifold M Q .
This geometric framework allows to combine the quantum and the classical parts of Ehrenfest's equations in a Hamiltonian system defined on the product manifold M C ×M Q , where M C is the phase space of the classical degrees of freedom defined above. The composition is completely analogous to the composition of two classical systems: the phase space of the whole system is the Cartesian product of the phase spaces of the subsystems, and a global Poisson bracket (or equivalently a global symplectic structure) is obtained by adding those of the subsystems: where {·, ·} C is the canonical classical Poisson bracket on M C . In a geometrical formalism, the physical observables are represented by functions on the global phase space of the system. From the geometric formalism of Quantum Mechanics, we know that a purely quantum operator, identified as a Hermitian operator A on the Hilbert space of the quantum system, is represented geometrically by the function Hamiltonian dynamics associated to these functions preserve the Kähler structure on the manifold M Q [34,[37][38][39][40][41], represented by the products of functions in (7). This property should also be true for generic observables on quantum-classical systems. As a consequence, for each ξ ∈ M C any observable has to be represented by a Hermitian operator A(ξ) acting on the Hilbert space of the quantum subsystem [40]. Thus, a function f A (ξ, ψ) representing a hybrid observable can always be written as When introducing mixed quantum-classical observables, it is relevant to analyse their algebraic structures. In Classical Mechanics, observables are represented by smooth functions on the phase space, which close a Lie algebra with respect to the classical Poisson bracket. For observables on quantum systems, as first shown by Heslot [40], it is enough to consider functions of the form of (9), which also close a Lie algebra, in this case with respect to the quantum Poisson bracket {·, ·} Q defined in (7). In the case of mixed quantum-classical systems, the natural extension of quantum functions, represented by (10), no longer close a Lie algebra with respect to the new Poisson bracket (8). Indeed, its action over two such functions, which are quadratic on the quantum degrees of freedom, would in general produce a non-quadratic function. Thus, from an algebraic point of view, a Lie algebra with respect to the quantum-classical Poisson bracket is obtained by considering all the smooth functions on M QC , as already proposed by some of us [1]. Observe that only the subset of functions of the form (10) corresponds to physical observables, while all the remaining functions have in principle no physical meaning. The function f H associated to the molecular Hamiltonian represents the energy of the hybrid quantum-classical system [1,2]. It is given by the expression: The electronic Hamiltonian H e (R) for finite-dimensional systems is obtained in an analogous way to (2). By using function f H and the Poisson bracket (8), it is possible to define a Hamiltonian vector field whose integral curves are the solutions to the Ehrenfest equations [1].

Ehrenfest statistical systems
It can be concluded that Ehrenfest equations, written in the appropriate way, are formally analogue to a classical Hamiltonian system. Therefore, from Liouville theorem we can define a volume on M C × M Q which is invariant under the dynamics. This allows us to define a statistical mechanical model where Ehrenfest equations define the dynamics of the microstates. The statistical description is defined by a distribution density function F QC on the manifold M C × M Q satisfying the following normalisation condition: where dµ QC (ξ, ψ) represents the volume element on M C × M Q defined by the symplectic volumes dµ C (ξ) and dµ Q (ψ) on the classical and quantum phase spaces, respectively. From an empirical perspective, the description of molecular systems by means of probability distributions seems natural, since any macroscopic system will have a certain distribution of states, associated to thermal motion or to a simple uncertainty of the exact state of the particles. Notice that, as we have a distribution depending on two types of degrees of freedom, we can consider the corresponding marginal distributions defined as: Marginal distributions F C and F Q encode respectively the corresponding uncertainties of the classical and quantum degrees of freedom.
By using these probability distribution we can measure the average value of any magnitude, classical, quantum or hybrid. The statistical average of a hybrid observable A, represented by a function f A (ξ, ψ), is obtained as We can also define a density matrix to represent the quantum subsystem in a more standard way [2] by averaging the projectors on the quantum states by the distribution F QC : The last equality shows that the density matrix ρ is the representation of the marginal distribution F Q as a density operator. It is immediate to verify that the density matrix ρ allows us to write the average value of a pure quantum observable B as Analogously, we can consider a ξ-dependent operator in order to represent in the same language the total probability distribution F QC : Observe that this is operator is not a density matrix, as in general Trρ C (ξ) = 1. It is possible to relate (15) and (17) as

Simple examples
The expressions obtained above are valid for any generic probability distribution F QC . A "pure state", in which the state of both the classical and quantum subsystem are completely determined to be ξ 0 and ψ 0 , corresponds to: The corresponding ρ C (ξ) operator, defined by (17), is To consider some uncertainty on the state of the molecule, we may study a situation such as that describes a molecule with uncertainty in the classical degrees of freedom, described by a probability distribution f (ξ) on M C . This is the case studied in a previous work [2] for a simplified model where the distribution was chosen as With the definition in (17), the ρ C (ξ) operator associated to (21) is

General systems
The choice of discrete distributions for the quantum degrees of freedom may seem a very special case, but as we are going to see it covers all possible situations. The reason for that arises from Gleason's theorem [43] and the given definition for ρ C (ξ) in (17). Gleason proved that the state of any quantum system can be encoded under the form of a density matrix. Extending this theorem to the case of hybrid quantum-classical systems, for a fixed value of the classical degrees of freedom ξ * , the density F QC (ξ * , ψ) must determine an operator ρ C (ξ * ) with the following properties: • it is Hermitian: Normalisation (12) of the full probability density F QC implies that the family of operators It is immediate to verify that expression (14) for the average value A of a hybrid observable A can be written as It is also immediate to understand that there are multiple probability densities which produce the same result for all possible observables on the system. Indeed, consider a generic probability density F QC and its associated operator ρ C (ξ). It is possible to consider the spectral decomposition of this operator, which allows us to rewrite it as where λ k (ξ) is the k-th eigenvalue of the operator and χ k (ξ) represents the corresponding eigenvector. We can now define a new probability density as The probability densities F QC and F ′ QC define the same operator ρ C (ξ), and therefore they produce the same expectation values for any hybrid observable of the system. In other words, these probability densities are indistinguishable and physically equivalent.
We can conclude that, for hybrid systems with finite dimensional quantum subsystems, it suffices to consider a family of discrete quantum distributions indexed by the classical degrees of freedom as in (27) to cover all physically meaningful situations. In the analysis of statistical systems, however, it is useful to consider both approaches, as each one has its advantages: the operators ρ C (ξ) are biunivocally related to the state of hybrid quantum-classical systems, whereas the probability density F QC can be more convenient for some mathematical manipulations, e.g. using Liouville's theorem.

Dynamical evolution
As seen above, the geometric formalism shows us that ED is of Hamiltonian type [1], which implies that the volume element dµ QC is preserved by the evolution. As a consequence, in the case of statistical systems, the evolution can be translated onto the space of probability distributions by means of the Liouville equation [2,44,45]: where f H is given for the Ehrenfest model by (11) and {·, ·} QC represents the Poisson bracket on M C × M Q is defined by (8). While the Ehrenfest equation for pure systems obviously cannot account for the phenomenon of decoherence, this statistical Ehrenfest equation does allow for decoherence to take place, as we will now show.
The time-dependent density matrix of the system is given by Liouville's equation (28) determines the evolution of the time-dependent probability density F QC (ξ, ψ; t) and from it we can immediately recover the density matrix for any time. In the case of a pure system: the evolution of the corresponding density matrix is given by: which resembles the expression of von Neumann's equation for quantum systems. The dependence on the classical positions R(t), however, causes this evolution to be non-linear, and therefore non-unitary. Nevertheless, the purity is still preserved, as this is simply a rewriting of the equations of the standard Ehrenfest model (5).
However, more complex distributions will evolve following more complex equations, decoherence will appear, von Neumann's equation will no longer hold, and purity will no longer be preserved. Consider that we add an extra term to our previous pure state, in the form: for a generic functionF C (ξ), with the only constraint of having zero integral, in order to satisfy the normalisation condition (12) for F QC (ξ, ψ). The dynamics for the density matrix of the quantum subsystem iṡ where ξ(ξ 0 , ψ 0 , t) and ψ(ξ 0 , ψ 0 , t) are the classical and quantum parts, respectively, of the Ehrenfest trajectory with initial conditions (ξ 0 , ψ 0 ). The dynamical equation for the density matrix has a structure similar to the equation introduced by Truhlar and co-workers for density matrices [11,16]: whereρ C (t) models the coherent evolution (i.e evolution according to the Liouville-von Neumann equation) andρ D (t) represents the dissipation term. By comparing this decomposition with the dynamical equation (33), we can conclude that the description of hybrid systems with uncertainty evolving under ESD allows us to obtain an a priori formulation for dissipative terms. In Sections 3.1 and 3.3, the particular characteristics of decoherence in ESD for a simple molecular system will be described.

Purity change, decoherence time, and stable projectors
Consider an isolated quantum system. The evolution of its density matrix ρ is governed by von Neumann's equation. In this case, its purity, is a constant of motion. For a (purely quantum) composite system, however, the situation is different. The evolution of the reduced density matrix of a subsystem, obtained through a partial trace on the total one, does not verify in general von Neumann's equation, and its purity is no longer conserved. An initial pure state, therefore, evolves into an statistical mixture of pure states, i.e. a mixed state. In some circumstances, an equation of motion can be derived or assumed for the reduced density matrix (e.g. of Lindblad type), describing the coupling of the reduced system to an external bath that may cause not only changes in purity, but also decoherence. The statistical hybrid model described above is also a composite system, although one of the subsystems is no longer quantum but classical. We should expect to find changes in purity and presence of decoherence, as in the case of purely quantum composite systems. In the context of environment-induced decoherence [46], decoherence in an open quantum system occurs when the interaction with an environment causes a non-reversible evolution of the quantum state of the system, represented by a density matrix, whose non-diagonal terms vanish for long times. Physically, the system typically evolves from a single state to a statistical mixture of states, whose projectors form the spectral decomposition of the final density matrix. One of the key points in the description of decoherence is that these states are not arbitrary: the decoherence hypothesis [18][19][20] states that these are elements from a fixed set selected by the environment, called the pointer basis. Thus, the presence of decoherence in the evolution of a quantum system can be deduced from the identification and description of a pointer basis for the system.
More specifically, in the context of HQCD, decoherence may appear by analising the evolution of the quantum subsystem; in this case, the classical part of the system acts as the environment that may select the pointer basis for the quantum subsystem. Due to the particular mathematical description of these models, decoherence in HQCD is due to the fact that neglected wave functions in the classical limit rapidly lose overlap in time. The pointer basis of the system is formed because the evolution of the classical subsystem destroys the superposition of states in the quantum subsystem [3,5,[18][19][20]. In conclusion, the interaction described in Ehrenfest dynamics causes the existence of a pointer basis for the quantum subsystem, selected by the classical subsystem, consisting of projectors onto subspaces of the Hilbert space of the quantum system which are both stable and attractors of the dynamics of the density matrices.
It is noticeable that our approach, which does not depend on any particular choice of basis for the system, allows us to consider a definition of decoherence time ∆ which is completely intrinsic. As the decrease in purity is a characteristic of decoherent evolution, we will consider that, for a system loosing coherence it time, such decoherence time is reached when its purity P stabilises.
In order to determine the pointer basis of a system, two aspects will be considered: stability and attractiveness. Consider a hybrid system evolving from a pure state and with some uncertainty on the initial values of the classical degrees of freedom. If decoherence occurs, then after a period determined by the decoherence time, the system will have evolved into an approximately stable mixture of states. Thus, in each particular example, the pointer basis can be identified as the set of those stable states, if any, whose projectors conform the density matrix of the system. In the next section, a simple diatomic molecule will be analysed both analytically and numerically, and we will determine the pointer basis of this molecular system.

Definition of the system
We have chosen as an example an ionised dimer, a simple yet non trivial model of molecule. For simplicity, we consider two atoms with only one valence electron, and all the inner electrons are identified together with the corresponding nucleus as a single classical particle called a 'core'. Summing up, the whole ionised molecule is described by two classical cores of charge +1 and a single quantum valence electron. Each core is given a mass of 23 a.m.u.; in this way, the system could be seen as a simple model for a ionised dimer of sodium, Na + 2 . We used the Octopus software [47,48], that uses an intuitively simple grid-based basis to represent the quantum degrees of freedom. The eigenstates and eigenvalues of the electronic Hamiltonian H e (R), as well as the dynamics, are computed in regular rectangular 3-dimensional grid contained in a spherical simulation box. The convergence tests showed that a grid with radius 13 a.u. and spacing 1.2 a.u is dense enough for a precise description of the classical and quantum degrees of freedom of our molecule. The dimension of the Hilbert space is thus given by the number of points in the real space grid. We did not impose absorbing boundary conditions at the walls of the chamber, since the simulations did not lead the electronic cloud to approach them.
The ion-electron interaction in Octopus is handled with pseudo-potentials, which avoid the divergence of the real Coulomb interaction. In this case, we opted for a simple soft-Coulomb potential, whose expression for a particle with charge Z: We take Z = 1 for the simulation of the ionised dimers. We have chosen the value α = 3 a.u., which reproduces approximately the experimental properties of the Na 2 neutral molecule. With this new potential, it is possible to consider the eigenvalue equation for the electronic Hamiltonian: with E 0 (R) ≤ E 1 (R) ≤ . . . The eigenstates of the electronic Hamiltonian are assumed to be normalised. Figure 1a shows the first eigenvalues of this operator. For each set R, the eigenvectors of the electronic Hamiltonian form a suitable basis for the Hilbert space of electronic states. It is possible to simplify the problem by restricting the positions of the cores. In the following, the set of positions of the cores, hereafter denoted by D, is restricted to positions R along the x axis, and such that the intercore distance varies only between 11 and 17 a.u. Figure 1b shows the behaviour of the eigenstates of H e (R): by taking a reference position R 0 , the change in the j-th eigenstate is given by the projection | φ j (R)|φ j (R 0 ) |. As the figure suggests, it can be proved that It is useful for the problem at hand to consider a basis for the electronic states formed by eigenstates of the electronic Hamiltonian at a reference position R 0 of the cores: For any R ∈ D, the electronic Hamiltonian H e (R) is approximately diagonal in this basis, which simplifies the computations. The next step in the description of the problem is the analysis of the dynamics determined by the Ehrenfest model (5). For given initial conditions (R 0 , P 0 , ψ 0 ), integration of the Ehrenfest model gives the evolution (R(t), P (t), ψ(t)) of the quantum and classical degrees of freedom. Initial parameters R 0 and P 0 are taken along the x axis, so that R(t) ∈ D for all t ≥ 0.
With these tools, it is interesting to analysis the adiabaticity of the Ehrenfest model when applied to the described example. For this purpose, let us choose as initial quantum state any of the eigenstates of the electronic Hamiltonian at R 0 , i.e. ψ (j) (0) = φ j (R 0 ). We have computed numerically the projections φ k (R(t))|ψ (j) (t) for different values of t Figure 1c represents these projections for j = k = 0, 1, 2. We have found that these projections take approximately value 1 when j = k. As for each R ∈ D elements in the basis {φ j (R)} are orthogonal to each other, the rest of projections are approximately zero. This leads to the following approximation We can conclude that, for the considered example, the system behaves in an approximately adiabatic way. Summing up, the two relevant approximations found for this model are (38) and (40). Both can be combined to approximate the evolution of generic initial conditions by ED (5) (as long as positions R(t) of the cores stay in the allowed set D): It is important to notice that adiabaticity (neither exact nor approximate) is not a general property of ED, as the Ehrenfest equations introduces non-adiabatic couplings [1]. If the state is written at each t as then from the Ehrenfest equation (5) it is immediate to compute that with d jk J (R) the non-adiabatic couplings. However, in our examples, the second term in the right hand side turns out to be negligible.

Statistical uncertainty for a single dimer: Decay of purity and pointer basis
We are interested in the appearance of decoherence effects: purity changes and pointer basis. For this reason, let us consider that the uncertainty in the initial conditions of our dimer is described by the following probability density function: with R = (R 1 , R 2 ) and P = (P 1 , P 2 ) the position and momenta, respectively, of the two cores in our model, and ψ the quantum state of the valence electron. Observe that the only uncertainty is assumed in the classical degrees of freedom, while the quantum subsystem is taken in the initial state ψ 0 ∈ M Q . For simplicity, the initial ionic positions are also fixed, with R 0 the equilibrium position of the cores along the x axis for the given quantum state M Q . As explained above, initial momenta P j,0 should be taken along the x axis. Numerical results show that a value of N = 41 is high enough to provide good statistical results. The evolution of the probability distribution can be written as From this expression we can write the corresponding density matrix as: Let us now analyse some simple examples, which will become useful in order to predict the changes in purity and the existence of a pointer basis for the molecule. As a first approach to the problem, consider in the initial probability distribution (44) that the initial quantum state is an eigenstate of the electronic Hamiltonian, i.e. ψ 0 = φ k ∈ B. An approximation to the density matrix of the system can be obtained by substituting (41) in (46) for every set of initial conditions (R 0 , P j,0 , ψ 0 ). Observe that, due to the small variation of the eigenstates of H e (R) and to the approximate adiabaticity of the dynamics, eigenstates of the electronic Hamiltonian are approximately constant (except for a global phase). As a consequence, the purity of the quantum subsystem is close to 1 for all the evolution. Observe that, although no decoherence seems to occur in this example, this is the first hint of the existence of a pointer basis for the molecule.
Consider now a second example. In order to extract some conclusion on the behaviour of the system, we choose N = 2 and we take the following expression for the initial probability distribution: where the initial quantum state is chosen as According to (41), the quantum trajectories ψ 1 (t), ψ 2 (t), determined by the Ehrenfest equations with initial conditions (R 0 , P 1,0 , ψ 0 ), (R 0 , P 2,0 , ψ 0 ) respectively, are the following for j = 1, 2. Substituting in (46) it is possible to estimate the density matrix of the quantum subsystem: where we have defined the quantities Observe that g j (t) are the gaps between the ground state and the first excited state energies for each initial conditions at each time t. Notice that the density matrix ρ(t) corresponds to the sum of two projectors, each one projecting onto the subspace generated by the vector Even if the difference |g 1 (t) − g 2 (t)| between the gaps is very small, integration in time makes the exponent ∆ j (t) non-negligible. Thus, the difference between the vectors |ψ 1 (t) and |ψ 2 (t) will become periodically relevant for the description of the quantum state. We can verify this assumption by computing the spectral decomposition of the density matrix ρ(t), defined in (15). The dependence of its eigenvalues with time is presented in Figure 2a. We notice that the number of eigenvalues different from zero is indeed a function of time, depending on the orthogonality of vectors χ 1 (t) and χ 2 (t).
A more general analysis can be done by extending the above computation to generic values of N in the initial distribution (47). If the quantum state is approximated by (41), the density matrix of the quantum subsystem has the following form: which is a sum of projectors onto the subspaces generated by Notice that the corresponding sum of projectors tend to zero when N grows, since we are adding N time-dependent vectors of norm one, moving with different velocities, in the linear space generated by φ 0 and φ 1 . The case for N = 3 is represented in Figure 2b. The case N = 41 is represented in Figure 2c. Notice how we obtain two eigenvalues approximately equal and a series of remaining eigenvalues representing the effect of the change of the position of the cores in the spectrum of the Hamiltonian H e (R) and the non-vanishing part from the projector sum, with a negligible weight. Indeed, for large values of N , coefficients in the sum in (53) are an approximately random set of complex numbers with modulus one, their sum being zero. Thus, for large enough values of t, the density matrix tends to The asymptotic value of the purity for a generic case can be computed based on the identification of the pointer basis with eigenstates of the electronic Hamiltonian, as detailed above. If the initial state is a linear superposition of n such states, with equal coefficients, then the purity tends to a value of 1/n. For generic linear combinations, the asymptotic value of purity can be computed as: Observe that the complex phase of the coefficients is irrelevant for the computations, as the eigenstates of the electronic Hamiltonian can always be redefined as c j φ j = |c j |φ ′ j . Additionally, the above computations show the existence of pointer basis in ESD. According to the meaning of decoherence in the context of HQCD [3,5], described above, the pointer basis appears as the set of stable states under the non-unitary evolution of the system. From (41), it can be concluded that eigenstates of the electronic Hamiltonian, i.e. elements in the basis B, are for ESD the elements in the pointer basis of the quantum system. Summarising, we have observed that the classical subsystem (the cores), acting as an environment onto the valence electron, produce a dynamics that, for times longer than the decoherence time, selects a small set of the possible quantum states. In this particular example, these are the eigenstates of the electronic Hamiltonian. These conclusions can be a starting point for a full analysis of pointer basis in molecular systems and in ESD, which could be performed by considering more general initial parameters and larger molecular systems. Observe that this is in agreement with the decoherence hypothesis [19], as the pointer basis does not depend on the initial state of the quantum subsystem. It is also important to notice that the results here obtained are characteristic of ESD. Other decoherent mechanisms will present different features and the pointer basis might appear in a completely different fashion.

Numerical simulations of an ionised dimer
In order to illustrate the results in Section 3.2, we have computed explicitly the changes in the purity of a simple model of an ionised dimer, as described above. As initial conditions, we take a certain initial state ψ 0 of the quantum subsystem. For this state, the equilibrium position R 0 of the cores are taken as their initial positions. The nuclei are allowed to move only along the axis of the molecule and keeping their center of mass stationary. Some uncertainty is allowed in the values of the initial momenta, ranging from 0 to 40 × 10 −5 a.u., in both directions.
With these values, the kinetic energy of the cores is never large enough to dissociate the dimer. Thus, different regimes of oscillation are considered. In conclusion, the uncertainty in the initial conditions is given in the form of (44).
Three different choices for the quantum subsystems are considered. In each case, the initial position of the cores is always the equilibrium position for the given quantum state: • Case A: The initial quantum state is an eigenstate of the electronic Hamiltonian: • Case B: The initial state is a linear superposition of two eigenstates of the electronic Hamiltonian: • Case C: The initial state is a linear superposition of three eigenstates of the electronic Hamiltonian: The results plotted in Figure 3a are in agreement with the behaviour given in (56). ESD leads to changes in the purity of the quantum states [1,2], unlike standard ED, which was purity-preserving. Let us analyse in detail each example. Case A reproduces the case in which the initial quantum state is an eigenstate of the electronic Hamiltonian at R 0 (see the beginning of Section 3.2). The purity is approximately constant and equal to 1 for all the evolution. On the contrary, cases B and C show a sharp drop on the purity (on the left side of Figure 3a), and reach approximately stationary values of 1/2 and 1/3, respectively.
The decoherence time can be read from the plots as the time it takes for the system to reach a stationary value. This definition is of course not precise, but in any case this time is approximately 1200 and 2500 a.u. for cases B and C, respectively. This time is very short in comparison with the typical oscillation time of our molecule, which has been computed to be approximately 30000 a.u.
Observe that the numbers N of different trajectories considered in the mixtures are essential in order to understand the behaviour of the system. Figures 3b and 3c represents the purity of cases B and C, respectively, for different numbers of trajectories, chosen in increasing order of initial speeds of the cores. It can be inferred that a low number of trajectories causes large oscillations in the purity, while adding trajectories causes the stabilisation of the value of the  (57), (58) and (59), respectively. The cores are initially at the equilibrium position for the corresponding electronic state, and with different initial speeds along the molecule axis, as described in the text. The system evolves according to ESD. (a) Evolution of the purity of the density matrix ρ(t) for the proposed cases. The purity reaches an approximately stable asymptotic value of 1/n, with n the number of eigenstates of the electronic Hamiltonian whose linear superposition determine the initial state of the quantum subsystem. (b-c) Evolution of the purity of the density matrix ρ(t) with respect to time and the value of N in (44). The initial quantum state corresponds to (58) in (b) and to (59) in (c). It can be observed that the decoherence time and the asymptotic value of the purity stabilise for larger values of N . purity after the decoherence time. If much large cases were considered, these oscillations are expected to disappear. Decoherence time would then determine when the system has become a mixture of states from the pointer basis, and the decoherence hypothesis introduced in Section 1 is satisfied [19].

Conclusions
Non-adiabatic transitions and decoherence are two of the most important phenomena in chemical dynamic reactions. To which extent does a given theoretical model include each effect is a complicated question. Several methods have been developed to deal with both notions [15].
In previous works [1,2], we introduced a geometric route to define the Ehrenfest Statistical Dynamics (ESD). We proved that, for a simple toy-model, the resulting ESD is purity nonpreserving. Now, after having carefully defined the meaning of decoherence in the context of hybrid quantum-classical dynamics (HQCD), we have tested out that, in ESD, decoherence does appear.
In fact, in this article we have proved that, when applied to a realistic molecular model, ESD defines a noncoherent evolution for the electronic degrees of freedom of the molecule. In our computations, the state space encoding the quantum degrees of freedom is considered as the finite dimensional vector space obtained by considering the values of the electronic wave-function on a grid, instead of the more usual approach for Ehrenfest dynamics that considers the truncated expansion in the electronic Hamiltonian eigenbasis. It is important to notice that this choice and the intrinsic (i.e., tensorial) nature of our formalism ensures that our results are completely independent of any choice of basis for the system. Thus, in this intrinsic way, we have proved how the statistical nature of the model is able to encode, in the simple language of Ehrenfest dynamics, non-coherent phenomena which are observed in realistic situations.
Some features of the decoherence phenomenon have been investigated. In particular, the purity of the electronic state has been observed to decrease in a short time, reaching an asymptotic value. Also, we have observed the existence of a pointer basis, which in our example turns out to be composed of the eigenstates of the electronic Hamiltonian of the molecular system. Notice that this fact cannot be predicted for other systems and situations (i.e. situations in which the evolution is not approximately adiabatic). In conclusion, the present paper shows that decoherence and pointer basis can be observed in simple models evolving under ESD. Further studies are required in order to extend this analysis to more complex systems.