Canonical mean-field molecular dynamics derived from quantum mechanics

Canonical quantum correlation observables can be approximated by classical molecular dynamics. In the case of low temperature the ab initio molecular dynamics potential energy is based on the ground state electron eigenvalue problem and the accuracy has been proven to be $\mathcal O(M^{-1})$, provided the first electron eigenvalue gap is sufficiently large compared to the given temperature and $M$ is the ratio of nuclei and electron masses. For higher temperature eigenvalues corresponding to excited electron states are required to obtain $\mathcal O(M^{-1})$ accuracy and the derivations assume that all electron eigenvalues are separated, which for instance excludes conical intersections. This work studies a mean-field molecular dynamics approximation where the mean-field Hamiltonian for the nuclei is the partial trace $h:={\rm Tr}(H e^{-\beta H})/{\rm Tr}(e^{-\beta H})$ with respect to the electron degrees of freedom and $H$ is the Weyl symbol corresponding to a quantum many body Hamiltonian $\widehat{H}$. It is proved that the mean-field molecular dynamics approximates canonical quantum correlation observables with accuracy $\mathcal O (M^{-1}+ t\epsilon^2)$, for correlation time $t$ where $\epsilon^2$ is related to the variance of mean value approximation $h$. Furthermore, the proof derives a precise asymptotic representation of the Weyl symbol of the Gibbs density operator using a path integral formulation. Numerical experiments on a model problem with one nuclei and two electron states show that the mean-field dynamics has similar or better accuracy than standard molecular dynamics based on the ground state electron eigenvalue.

1. Classical approximation of canonical quantum observables 1.1. Introduction to the approximations. We study approximation of quantum timecorrelation observables at the quantum canonical ensemble for a system consisting of nuclei (slow degrees of freedom) and electrons (fast degrees of freedom) at the inverse temperature β = 1/(k B T ), where k B is the Boltzmann constant and T > 0 is the temperature. We work in Hartree atomic units in which the reduced Planck constant = 1, the electron charge e = 1, the Bohr radius a 0 = 1 and the electron mass m e = 1. Thus the semiclassical parameter in the subsequent analysis is given by the ratio of nucleus and electron masses M . For example, in the case of a proton-electron system the ratio is M = m p /m e ≈ 1836. The full quantum system is described by the Hamiltonian operator which includes the kinetic energy of nuclei and the electronic kinetic energy operator together with the operator describing interaction between electrons, with coordinates x e , and nuclei, with coordinates x, In this work, in the spirit of Born-Oppenheimer (adiabatic) approximation, we replace the time evolution of electrons by the Schrödinger electron eigenvalue problem. We represent the electronic kinetic energy operator and the interaction operator, H e = − 1 2 ∆ xe + V e (x, x e ) as a matrix-valued potential V (x) obtained by a representation of the operator − 1 2 ∆ xe + V e (x, x e ) on a finite-dimensional (d-dimensional) subspace of suitable normalized electronic eigenfunctions {φ k } d k=1 , as V (x) kℓ = φ k , H e (x, ·)φ ℓ , described precisely in Section 2. Hence we work with the Hamiltonian operator The first term − 1 2M ∆ x ⊗ I represents the kinetic energy of the nuclei where I is the d × d identity matrix. The second term, V (x), is the matrix-valued potential approximation to H e and does not depend on M . We assume that this finite-dimensional approximation of the electronic operator results in a Hermitian matrix-valued smooth confining potential V : R N → R d×d that depends on the positions x i ∈ R 3 of nuclei i = 1, 2, . . . , N ′ , where we set N = 3N ′ . For the sake of simplicity, we assume that the nuclei have the same mass; in the case with different nuclei masses M becomes a diagonal matrix, which can be transformed to the formulation (1.1) with the same mass by the change of coordinates M 1/2 1x = M 1/2 x. The large nuclei/electron mass ratio M ≫ 1 is the basis of semiclassical analysis and implies a separation of time scales, for which nuclei represent slow and electrons much faster degrees of freedom. The Weyl quantization takes this scale separation into account. In particular for the Hamiltonian operator H the corresponding matrix valued Weyl symbol becomes H(x, p) = 1 2 |p| 2 I + V (x) for the nuclei phase-space points (x, p) ∈ R N × R N , as described more precisely in Section 2.
In order to study correspondence between the quantum time-correlation function and its classical counter part we work in Heisenberg representation for the time-dependent quantum observables given by self-adjoint operators A t and B t . We employ the Weyl quantization to link the quantum dynamics given by the Heisenberg equation to classical Hamiltonian equations of motions on the phase space (x, p) ∈ R N × R N of nuclei positions and momenta and to averaging with respect to a suitable canonical Gibbs distribution on the phase space.
More precisely, given a quantum system defined by the Hamiltonian H acting on wave functions in L 2 (R N , C d ) ≡ [L 2 (R N )] d we denote ρ = e −β H the density operator for a quantum Hamiltonian operator H at the inverse temperature β > 0 and consider quantum correlation observables based on the normalized trace with the initial data z 0 := (x 0 , p 0 ) ∈ R 2N of nuclei positions and momenta. The trace Tr (e −βH(z 0 ) ) is over the space R d×d of d × d matrices which can be also viewed as the trace operator with respect to electron degrees of freedom under the finite dimensional approximation of the electronic Hamiltonian. On the other hand Tr (e −β H ) represents the trace acting on the space of trace operators on L 2 (R N , C d ) which we can view as the trace with respect to both nuclei and electron degrees of freedom.
Two main questions arise: (a) how should the mean field Hamiltonian approximation h : R 2N → R be chosen, and (b) how small is the corresponding estimate for the approximation error T qm − T md ? Assume Ψ : R N → C d×d is a differentiable orthogonal matrix. Based on certain regularity assumptions on A 0 , B 0 , V and Ψ, we prove in Theorem 2.1 that for the mean-field Hamiltonian h : R N × R N → R defined by (1.6) h(z) := Tr (H(z)e −βH(z) ) Tr (e −βH(z) ) , and the symbols A 0 and B 0 , which are independent of the electron coordinates, we have where the parameters ǫ 2 j are the variances Tr e −βH L 1 (R 2N ) , Tr e −βH Ψ L 1 (R 2N ) , (1.8) with the definition H Ψ (x, p) := |p| 2 2 I+Ψ * (x)V (x)Ψ(x) using the Hermitian transpose Ψ * (x) .
We note that the mean-field Hamiltonian can be written h(x, p) = |p| 2 2 + Tr (V (x)e −βV (x) ) Tr (e −βV (x) ) where λ i (x), i = 0, . . . , d − 1, are the eigenvalues of V (x), and λ * : R N → R is the obtained mean-field potential. Therefore the mean-field h is independent of the large mass ratio parameter M , so that the dynamics (1.5) is independent of M and consequently the nuclei move a distance of order one in unit time. We see that the mean-field h = Tr (He −βH )/Tr (e −βH ) is the mean value with respect to the Gibbs density. The error term ǫ 2 1 can be written as the corresponding normed variance (1.10) and at points x where the eigenvalues λ i (x) are separated a suitable choice of Ψ(x) is the matrix of eigenvectors to V (x) which implies (1.11) On the other hand to have sets with coinciding eigenvalues is generic in dimension two and higher, see [24], and there the matrix Ψ(x) is in general not differentiable. Therefore (1.11) will typically not hold everywhere. Section 5 presents numerical experiments where also the size of the error terms are analyzed in different settings. In Section 2 we review necessary background on Weyl calculus and state the main theoretical result, namely the error estimate (1.7), as Theorem 2.1, together with the ideas of its proof. Sections 3 and 4 present the proof of the theorem. Section 5 presents numerical experiments on the approximation error T qm (t) − T md discussed in the next Section 1.2.
1.2. Numerical comparisons. In Section 5 we present some numerical examples with varying settings of t, 1/M, ǫ 2 i , related to different potentials V with the purpose to study the following questions: Is the estimate (1.7) sharp or does the error in practise behave differently with respect to t, 1/M and ǫ 2 i ? Is the main contribution to the error coming from approximation of the the matrix-valued potential by a scalar potential in the quantum setting or from the classical approximation of quantum dynamics based on scalar potentials ? Can the mean-field dynamics improve approximation compared to using molecular dynamics based on the ground state eigenvalue λ 0 instead of λ * ? Theorem 2.1 does not give precise answers to these questions. The aim of this section is to provide some insight from several numerical experiments on a model problem, chosen to avoid the computational difficulties for realistic systems with many particles. Therefore we use one nuclei in dimension one, N = 1, and two electron states, defined by the Hamiltonian where I is the 2 × 2 identity matrix, V : R → R 2×2 given by with the two eigenvalues (1.14) plotted in Figure 1 (Case E in Table 1).  Figure 1. The eigenvalue functions λ 0 (x) and λ 1 (x) of the potential matrix V (x), and the corresponding mean-field potential λ * (x), for the parameters c = 1, δ = 0.1 (Case E). Section 5 presents numerical results comparing quantum mechanics to the three different numerical approximations based on: the ground state potential λ 0 , mean-field potential λ * and excited state dynamics. The excited state molecular dynamics studied in [14] uses several paths related to different electron eigenvalues and is defined by with the initial condition z j 0 = (x 0 , p 0 ) = z 0 .
Numerical results on equilibrium observables show that mean-field and excited state molecular dynamics are more accurate than molecular dynamics based only on the ground state.
In the case of correlation observables with τ > 0, mean-field and excited state molecular dynamics give in general different approximations.

Correlation observables.
Observations on quantum dynamics for a system in interaction with a heat bath at thermal equilibrium can be approximated by correlations (1.2) in the canonical ensemble, cf. [6,20,8,13]. For instance, the classical observable for the diffusion constant includes the time-correlation x(τ )·x(0). Hence the corresponding quantum correlation (1.2) would for this case use A τ = x τ I and B 0 = x 0 I, and The numerical results in Section 5 for time-dependent observables are mainly based on the momentum auto-correlation which is related to the diffusion constant D by the Green-Kubo formula [12] since the velocity is equal to the momentum in our case with unit particle mass.
The different numerical experiments in Table 1 are chosen by varying the parameters such that all three, two, one or no molecular dynamics approximate well. In Case A, with low temperature and large eigenvalue gap, all three molecular dynamics approximate the quantum observable with similar small error and also the error terms 1/M , ǫ 2 1 and ǫ 2 2 are small. In Case B, with small difference of the eigenvalues (i.e. c is small), the mean-field and excited states dynamics is more accurate than ground state dynamics and the error terms 1/M , ǫ 2 1 and ǫ 2 2 are still small. The result is similar in Case C, with an avoided crossing (i.e. δ is small) and small difference of the eigenvalues. In Case D, with high temperature and larger difference of the eigenvalues, only the excited state dynamics provides accurate approximations to the quantum observables and the error terms ǫ 2 1 and ǫ 2 2 are large. Finally in Case E, when the difference of the eigenvalues are sufficiently large and we have an avoided crossing, molecular dynamics is accurate only for short correlation time τ and the error terms ǫ 2 1 and ǫ 2 2 are large.  (5.4). Low value of β means a high system temperature, the parameter c measures the difference between the two eigenvalues γ λ := R |λ 1 − λ 0 | e −βλ 0 dx/ R e −βλ 0 dx and the parameters c and δ determine the eigenvalue gap. The value of ǫ 2 2 is computed following (1.11) with Ψ(x) the matrix of eigenvectors to V (x). Figures 11b and 12b show that in Case D and Case E, where mean-field and ground state approximations are not accurate, the approximation error is dominated by the part corresponding to replacing the matrix valued potential by a scalar potential in the quantum formulation and not the part of the error resulting from classical approximation of quantum mechanics for scalar potentials. More precise conclusions relating the error terms 1/M , tǫ 2 1 and t 2 ǫ 2 2 in (1.7) to the numerical experiments are in Section 5. The numerical experiments here also show that the mean-field dynamics has similar or better accuracy compared to ground state molecular dynamics. It would be interesting to do this comparison for realistic problems.
1.3. Relation to previous work. Classical approximation of canonical quantum correlation observables have been derived with O(M −1 ) accuracy for any temperature, see, e.g., [23,14]. Computationally this accuracy requires to solve classical molecular dynamics paths related to several electron eigenvalues, while the mean-field dynamics has the advantage to use only one classical path at the price of loosing accuracy over long time.
Classical limits of canonical quantum observables were first studied by Wigner, [25]. His proof introduces the "Wigner" function for scalar Schrödinger equations and uses an expansion in the Planck constant to relate equilibrium quantum observables to the corresponding classical Gibbs phase space averages.
To derive classical limits in the case of matrix or operator-valued Schrödinger equations previous works, see [23], diagonalize the electron eigenvalue problem, which then excludes settings where the eigenvalues coincide at certain points due to the inherent loss of regularity at such points. The mean-field formulation presented here avoids diagonalization of the electron eigenvalue problem at points with low eigenvalue regularity.
The classical limit with a scalar potential V , e.g. the electron ground state eigenvalue, has been studied by three different methods: (1) Solutions to the von Neumann quantum Liouville equation, for the density operator, are shown to converge to the solution of the classical Liouville equation, using the Wigner function [17], also under low regularity of the potential, cf. [7]. These results use the Wigner function and compactness arguments, which do not provide a convergence rate. (2) In the second method, used in our work, the two main mathematical tools are a generalized version of Weyl's law and quantization properties, as described by semiclassical analysis, e.g., in [26] and [18]. The generalized Weyl's law links the trace for canonical quantum observables to classical phase space integrals, related to Wigner's work. The quantization properties compare the quantum and classical dynamics and provide convergence rates. Our study differs from previous works that also used similar tools, e.g., [23,2]. The standard method to bound remainder terms in semiclassical expansions, based on the Planck constant , use the Calderón-Vaillancourt theorem to estimate operator norms. Such approach yields error bounds with constants depending on the L 1 (R 2N )-norms of Fourier transforms of symbols and the potential. The L 1 (R 2N )-norm of a Fourier transformed function can be bounded by the L 1 (R 2N ) norm of derivatives of order N of the function. Therefore the obtained constants in O( α ) error estimates are large in high dimension N . In our work here we apply dominated convergence to obtain error estimates based on low regularity of symbols and the potential, while Fourier transforms in L 1 (R 2N ) are only required to be finite and do not enter in the final error estimates. (3) The third alternative in [9] provides a new method that also avoids the large constants in the Calderón-Vaillancourt theorem in high dimensions by using convergence with respect to a generalized Wasserstein distance and different weak topologies.
A computational bottleneck in ab initio molecular dynamics simulations of canonical correlation observables is in solving electron eigenvalue problems at each time step. An alternative to approximate quantum observables is to use path integral Monte Carlo formulations in order to evaluate Hamiltonian exponentials. The Hamiltonian exponentials come in two forms: oscillatory integrals in time t ∈ R, based on e it H for the dynamics, and integrals for Gibbs function e −β H that decay with increasing inverse temperature β ∈ (0, ∞). The high variance related to the oscillatory integrand e it H means that standard computational path integral formulations for molecular dynamics are applied only to the statistics based on the partition function Tr (e −β H ) while the dynamics is approximated classically.
Two popular path integral methods are centroid molecular dynamics and ring-polymer molecular dynamics, see, e.g., [11] or [22]. In these methods the discretized path integral is interpreted as a classical Hamiltonian with a particle/bead for each degree of freedom in the discretized path integral. For the centroid version the dynamics is based on the average of the particle/bead positions, i.e., the centroid, with forces related to a free energy potential for the partition function thereby forming a mean-field approximation. It is related to the mean-field approximation (1.4) and (1.6) but differs in that in our work the forces are based on the mean Hamiltonian, for the partial trace over electron degrees of freedom, instead of on the partition function for the discretized path integral with respect to both nuclei and electron degrees of freedom, centered at the centroid. In ring-polymer molecular dynamics classical kinetic energy is added for each bead forming a Hamiltonian with harmonic oscillators in addition to the original potential energy. Consequently the phase-space is related to coupled ring polymers, one for each original particle.
There is so far no convergence proof for centroid nor ring-polymer molecular dynamics. Therefore it would be interesting to further study their relation to the mean-field model we analyse here. The mean-field formulation (1.4) and (1.6) can also offer an alternative to standard eigenvalue solutions by using a path integral formulation of the partial trace over the electron degrees of freedom, in the case of sufficiently large temperature avoiding the fermion sign problem, see, e.g., [4]. Another difference to previous work is that the convergence proof here derives a precise asymptotic representation of the Weyl symbol for the Gibbs density operator using a path integral formulation, providing an example that using path integrals for the Gibbs density can also result in simplification of the theory.

The main result and background from Weyl calculus
In this section we state the main theorem and review necessary tools from semiclassical analysis and functional integration.
To relate the quantum and classical observables we employ Weyl calculus for matrixvalued symbols. First, we introduce functional spaces that we use in the sequel: (i) the Schwartz space of matrix-valued functions on the phase space where we denote a point in the phase space z = (x, p) and for a multi-index of non-negative integers α = (α 1 , . . . , α 2N ), we have the partial derivatives ∂ α z = ∂ α 1 z 1 . . . ∂ α 2N z 2N of the order |α| = i α i and similarly we have z γ = z γ 1 1 . . . z γ 2N 2N for the multi-index γ. For a matrixvalued symbol A we also use the notation A ′ (z) for the tensor (A ′ (z)) m ij := ∂ zm A ij (z) and A ′′ (z) for the 4th-order tensor (A ′′ (z)) mn ij := ∂ 2 zmzn A ij (z). The dual space of tempered distributions is denoted S ′ .
(ii) the space of L 2 vector-valued wave functions We define the Weyl quantization operator of a matrix-valued symbol A ∈ S as the mapping A → A that assigns to the symbol A the linear operator A : H → H defined for all Schwartz functions φ(x) by and extended to all wave functions φ ∈ H by density. The expression (2.1) shows that the kernel K A on H is the Fourier transform in the second argument of the symbol A(x, p) and consequently the Weyl quantization is well defined for symbols in S ′ , the space of tempered distributions.
For example, the symbol H(x, p) := 1 2 |p| 2 I + V (x) yields the Hamiltonian operator We formulate the main result as the following theorem estimating the mean-field approximation.
Theorem 2.1. Let Ψ : R N → C d×d be a differentiable mapping into orthogonal matrices and define V Ψ := Ψ * V Ψ, for the Hermitian potential V : R N → C d×d . Assume that the components of the Hessian V ′′ Ψ are in the Schwartz class and the scalar symbols a 0 : R 2N → C and b 0 : R 2N → C are infinitely differentiable and compactly supported. Furthermore, suppose that there is a constant k such that V + kI is positive definite everywhere, Tr (e −β H ) is finite, and there is a constant C such that then there are constants c,M ,t, depending on C and β, such that the quantum canonical observable (1.3), with A 0 = a 0 I and B 0 = b 0 I, can be approximated by the mean-field molecular dynamics (1.4)-(1.5) with the error

2.1.
Overview and background to the proof. This section provides background and motivation to the proof of the theorem in three subsections. The first subsection reviews application of Weyl calculus for the dynamics, the second one is on a generalized form of Weyl's law in order to relate the quantum trace to phase space integrals and the third subsection introduces path integrals and their application in the context of our result.
2.1.1. Weyl calculus and dynamics. This section first introduces the central relation between commutators and corresponding Poisson brackets for the classical limit of dynamics. Given two smooth functions v(x, p), w(x, p) on the phase space we define the Poisson bracket We denote the gradient operator in the variable z = (x, p) as ∇ z = (∇ x , ∇ p ) and ∇ ′ z = (∇ p , −∇ x ), hence the Poisson bracket is expressed as To relate the quantum and classical dynamics for particular observables with symbols of the type a 0 I treated in Theorem 2.1, we assume that the classical Hamiltonian flow with the initial data (x 0 , p 0 ) = z 0 ∈ R 2N . Given a scalar Schwartz function a 0 we define a smooth function on the flow z t (z 0 ) and we have that a t (z 0 ) satisfies The corresponding quantum evolution of the observable A, for t ∈ R, is defined by the Heisenberg-von Neumann equation which implies the representation The basic property in Weyl calculus that links the quantum evolution (2.7) to the classical (2.6) is the relation where the remainder symbol r a t is small. In Lemma 3.2 we show that The main obstacle to establish the classical limit for dynamics based on the matrix-valued operator Hamiltonian H is that matrix symbols do not commute, i.e. [H, A t ] = 0, which implies the additional larger remainder iM 1/2 [H, A t ] in (2.8). Therefore the usual semiclassical analysis perform approximate diagonalization of [ H, A t ], see [23,14]. Diagonalization of V introduces eigenvectors that are not smooth everywhere unless the eigenvalues are separated, due to the inherent loss of regularity for eigenvectors corresponding to coinciding eigenvalues, see [24]. To have a conical intersection point with coinciding eigenvalues is generic in dimension two, and in higher dimensions the intersection is typically a codimension two set, see [24]. The non smooth diagonalization has been so far a difficult obstacle to handle with the tools of Weyl calculus. The aim in our work is therefore to avoid diagonalization everywhere by analyzing a mean-field approximation differently.
In order to apply (2.8) we use Duhamel's principle, see [5]: the inhomogeneous linear equation in the variable A t − a t ≡ A t − a t I, where A t satisfies the evolution (2.7) and a t is defined by (2.5), can be solved by integrating solutions to the homogeneous problem with respect to the inhomogeneity The quantum statistics has a similar remainder term to (2.9), namely the difference ρ − e −βH of the Weyl symbol ρ for the quantum Gibbs density operator ρ = e −β H and the classical Gibbs density e −βH . To characterize asymptotic behaviour as M → ∞ of this difference we employ representation of the symbol ρ based on Feynman-Kac path integral formulation, as presented in Section 2.1.3.

Generalized Weyl law.
To link the quantum trace to a classical phase space integral we will use a generalized form of Weyl's law, see [26,23]. The semiclassical analysis is based on the fact that the H-trace of a Weyl operator, with a d × d matrix-valued symbol, is equal to the phase-space average of its symbol trace. Indeed we have by the definition of the integral kernel in (2.1) for A ∈ S In fact also the composition of two Weyl operators is determined by the phase-space average as follows, see [23].
Lemma 2.2. The composition of two Weyl operators A and B, with A ∈ S and B ∈ S satisfies Proof. The well-known proof is a straight forward evaluation of the integrals involved in the composition of two kernels and it is given here for completeness. The kernel of the composition is so that the trace of the composition becomes Tr A(y, p)B(y, p) dp dy , using the change of variables (y, The composition of three operators does not have a corresponding phase-space representation. We will instead use the composition operator # (Moyal product), defined by A B = A#B, to reduce the number of Weyl quantizations to two, e.g., as A B C D = A#B C#D. More precisely, the Moyal product of two symbols has the representation For general background we refer the reader to [26] or [18]. The isometry between Weyl operators with the Hilbert-Schmidt inner product, Tr ( A * B), and the corresponding L 2 (R N ×R N , C d×d ) symbols obtained by Lemma 2.2 shows how to extend from symbols in S to L 2 (R N ×R N , C d×d ) by density of S in L 2 (R N ×R N , C d×d ), see [23]. We will use the Hilbert-Schmidt norms A 2 HS = Tr ( A * A) = Tr ( A 2 ) and Tr A 2 to estimate Weyl operators and Weyl symbols, respectively. We show in Lemma 3.2 that having the Weyl symbols and V ′′ in the Schwartz class imply that dominated convergence can be applied in the phase space integrals obtained from the generalized form of Weyl's law.

Feynman-Kac path integrals.
In order to analyse the symbol ρ for the Gibbs density operator ρ = e −β H we will use path integrals (in so called imaginary time), as in [6], [1] and [21], based on Feynman-Kac formula applied to the kernel of the Gibbs density operator and its corresponding Weyl quantization. We start with the kernel representation To motivate the construction of the path integral representation we first identify the Weyl symbol of the density operator in the case of a scalar potential V , i.e., the case d = 1. In order to emphasize that we consider the scalar case, we denote this Weyl symbol ρ s (x, p), and thereby the associated kernel is denoted K ρs . Direct application of the Feynman-Kac formula, see Theorem 7.6 in [15], implies that the kernel can be written as the expected value We recall the definition (2.1) of Weyl quantization for the (scalar) symbol ρ s (x, p) from which we obtain the expression for the symbol of an operator associated with the kernel K ρs (x, y) Using the substitution x ′ = x + y 2 , and y ′ = x − y 2 , i.e., and letting imply by combining (2.14), (2.13) and the transformed path (2.15) We have obtained the path integral representation of the symbol corresponding to the Gibbs density operator e −β H in the case of the scalar potential V We can proceed along similar lines in the case of the matrix-valued potential V studied here. The Feynman-Kac formula has been derived for operator-valued potentials V , see [21], and in the case where the potential V (x) is a general matrix, the exponential (2.16) for the time evolution, is replaced by the corresponding matrix-valued process t ∈ (0, β) and Υ + 0 = I. We will use the notation W β t := 1 2 W β − W t . The steps in the derivation of (2.16) then imply that the symbol in the case of a matrix-valued potential can again be expressed by . To estimate ρ − e −βH we use the symmetry property that the Weyl symbol, ρ, for the Gibbs density operator ρ = e −β H is a Hermitian matrix. Indeed, we have e −β H represented as an L 2 -integral operator with the kernel K ρ and since is real and Hermitian also e −β H = ∞ n=0 (−β H) n /n! is real and Hermitian. Therefore the Weyl symbol corresponding to the Gibbs density operator ρ satisfies Either of these integral representations show that ρ is Hermitian. The same steps as above leading to (2.16) can by (2.18) be applied to the change of variables x ′ = x − y/2 and y ′ = x + y/2, which implies that we also have Therefore we have the symmetrized representation . Using the path-integral representation of the symbol ρ(x, p) we prove in Section 4.2 Lemma 2.3. Assume that the bounds in Theorem 2.1 hold. Then

21)
and Tr The section proves Theorem 2.1 based on three steps: Step 1. use Duhamel's principle recursively to analyse the dynamics based on H, Step 2. use estimates of remainders for Weyl compositions and the Weyl Gibbs density to analyse the statistics at t = 0, Step 3. repeat Step 1 and Step 2 with H replaced byH := Ψ * #H#Ψ to approximately diagonalize H.
Proof of the theorem.
Step 1. Lemma 3.2 shows that the commutator has the representation and the cyclic invariance of the trace together The right hand side can again be estimated by applying Duhamel's principle (3.3), now to B s−t and b s−t , as follows We have by Cauchy's inequality and the cyclic invariance of the trace (3.5) The following two lemmas, proved in Section 4, estimate the remainder terms T 0 and T 3 Lemma 3.1 (Mean-field approximation). Assume that the bounds in Theorem 2.1 hold, then and if c and d are scalar valued and if c is scalar valued where the limits hold in L 1 (R 2N ) and L ∞ (R 2N ). Furthermore the function a t : R 2N → C, defined in (2.5), is in the Schwartz class and there holds The combination of (3.2)-(3.5), Lemmas 3.1 and 3.2 imply Similarly the symmetrized difference has the same bound obtained by interchanging the role of A and B in (3.7).
Step 2. Here we estimate the second term in the left hand side of (3.8), which can be split into The first term in the right hand side in (3.9) has by Lemma 2.2 the classical molecular dynamics approximation Tr a 0 z t (z 0 ) b 0 (z 0 )e −βH(z 0 ) + Tr r ab (z 0 )e −βH(z 0 ) dz 0 , Tr r ab (z 0 )e −βH(z 0 ) dz 0 It remains to estimate the second term in the right hand side of (3.9). We have where the second last equality follows by interchanging the order of the trace and the integration with respect to β and using the cyclic invariance of the trace. The first equality is obtained by splitting the integral as : h(z) > m} and using that ρ is uniformly bounded in L ∞ (R 2N ): the second integral is zero for m sufficiently large, as verified in (3.14) below, and in the compact set we apply dominated convergence.
Verification of Lm |a t (z 0 )|dz 0 = 0. The integration with respect to the initial data measure dz 0 can be replaced by integration with respect to dz t since the phase space volume is preserved, i.e. the Jacobian determinant is constant for all time by Liouville's theorem. We first verify that By assumption R N Tr e −βV (x) dx < ∞. Therefore the smallest eigenvalue λ 0 (x) of V (x) also satisfies R N e −βλ 0 (x) dx < ∞ and as h(x, p) ≥ |p| 2 /2 + λ 0 (x) we have R 2N e −βh(x,p) dx dp < ∞, which combined with the assumption ∇ x h(x, p) L ∞ (R 2N ) ≤ C establishes (3.13). By using the two properties h z t (z 0 ) = h(z 0 ) and h(z) → ∞ as |z| → ∞ together with the compact support of a 0 we obtain (3.14) In conclusion we have which combined with (3.8) proves the theorem for Ψ = I.
Step 3.1. Let α be any complex number and define for t ∈ R the exponential ȳ t := Ψ * e tαĤ Ψ.
Therefore the transformed variable and the cyclic property of the trace implies that the quantum observable satisfies (3.16) where the initial symbols are given bȳ Step 3.2. We have as derived in [14] from the composition in (2.12): where by the orthogonality Ψ * Ψ = I Let Ψ(x) be the orthogonal matrix composed of the eigenvectors to V (x), then the matrix is diagonal, with the eigenvectors λ i (x) of V (x) forming the diagonal d×d matrix Λ(x). The non diagonal part 1 4M ∇Ψ * (x)·∇Ψ(x) ofH(x, p) is small if Ψ(x) is differentiable everywhere. If the eigenvectors to V are not differentiable in a point x * , we may use a regularized version of Ψ in a neighbourhood of x * to form an approximate diagonalization of H.
Step 3.3. The derivation in Step 1 can now be repeated with H replaced byH and A, B byĀ,B. Duhamel's principle (2.10) implies and we obtain by (3.16) as in (3.7) where H is replaced byH in Da, Db, T 0 and T 3 , so that The two terms T 0 and T 3 have the bounds in Lemmas 3.1 and 3.2 with H replaced bȳ H. It remains to show that the initial errors satisfy (3.18) Tr To prove (3.18) let ρ = e −β H , then by the composition (2.12) and Lemma 2.2 we obtain To estimate the initial errorĀ 0 − a 0 we use Lemma 3.2 which by the composition (2.12) implies Here the orthogonality Ψ * Ψ = I implies ∇(Ψ * Ψ) = 0. We obtain as in (3.11) the limit which proves (3.18).

Proof of Lemmas
This section estimates the remainder terms T 0 and T 3 and the statistical error ρ − e −βH . The error term T 3 in Lemma 3.2 is due to remainders from classical approximation while T 0 in Lemma 3.1 is the main term regarding the mean-field approximation. The statistical error is estimated in Lemma 2.3.

Proof Lemma 3.1.
Proof. We consider first the case Ψ = I, i.e.H = H. We have Tr (r a s #b s−t )ρ dz ds Tr (r a s #b s−t )ρ) dz The third equality uses that (H − h) commutes with e −βH and the forth and the last equality is obtained by interchanging the order of the trace and the integral with respect to τ in combination with the cyclic invariance of the trace. Cauchy's inequality, the positive definiteness of e −βH (H − h) 2 and H − h depending only on x imply that the right hand side has the bound

This partition implies
and by assumption there is a constant k such that V + kI is positive definite everywhere. Therefore we have d dt which establishes Tr (Υ + β ) 2 ≤ e kβ and shows that for independent Wiener processes W and W ′ we obtain by Cauchy's inequality We also observe that the generalized Weyl's law implies that ρ is in L 2 (R 2N , C d×d ), namely which proves (2.22).

Proof Lemma 3.2.
Proof. To estimate remainder terms we will use the composition operator for Weyl symbols defined by c#d = c d. The composition operator has the representation which can be written as an expansion using the Fourier transform F, defined for f : R 4N → R by (4.10)

Its inverse transform implies
and Taylor expansion of the exponential function yields . (4.11) The pointwise limit of the remainder term can be estimated by dominated convergence In addition we need convergence in L 1 (R 2N ), as a function of z = (x, p), to apply dominated convergence in the phase-space integrals. We have Therefore we obtain remainder terms that are uniformly bounded in L 1 (R 2N ) provided the Fourier transform (iξ) α Fc(ξ) of ∂ α z c(z) satisfies (4.12) and similarly for d We will apply the composition expansion to functions in the Schwartz class so that (4.12) holds.
We conclude that Schwartz functions c and d satisfy using splitting of the phase space integral as in (3.12), which together with (4.13), (4.14), (4.15) and Lemma 4.1 below proves the lemma.

By summation and maximization over indices we obtain the integral inequalities
(4.17) The functions max ij |α|≤2 ∂ α z 0 ∂ z j z i (t, z 0 ) can therefore be estimated as in [10] by Gronwall's inequality, which states: if there is a positive constant K and continuous positive functions Gronwall's inequality applied to (4.17) implies on the ground state, on the mean-field energy surface, and on a weighted average of all eigenstates (denoted the excited state dynamics below). In particular, we study whether the mean-field approximation can be more accurate than using only the ground state.
In order to demonstrate the proposed mean-field molecular dynamics approximation, we devise the model problem as described by equations (1.12) to (1.14) in Subsection 1.2, where the difference between two eigenvalues λ 1 (x) − λ 0 (x) = 2c √ x 2 + δ 2 can be tuned by the two parameters c and δ. For δ small this model relates to the avoided crossing phenomenon in quantum chemistry where the two potential surfaces get almost intersected at a certain point, see [24]. The assumptions in Theorem 2.1 are satisfied for positive δ but not for δ = 0 and therefore the approximation error is expected to vary with different δ. Tr(e −βĤÂ ) We apply a fourth-order finite difference scheme for the Laplacian operator in the Hamiltonian (1.12) to approximate the equilibrium density µ qm (x) in (5.2). The numerical implementation is explained with more details in Appendix A.1.
As an approximation to the quantum canonical ensemble average (5.1), we consider the normalized trace Tr ( e −βH A) / Tr ( e −βH ) and apply the Lemma 2.2 and equation (1.16) to write the mean-field observable as Tr ( e −βH ) = R 2N Tr(e −βH(x,p) a(x, p)I) dx dp R 2N Tr(e −βH(x ′ ,p ′ ) ) dx ′ dp ′ = R 2N a(x, p)(e −βλ 0 (x) + e −βλ 1 (x) ) e − β|p| 2 2 dx dp where H(x, p) and A(x, p) = a(x, p)I with a : R 2N → C are the Weyl symbols corresponding to the operatorsĤ andÂ, respectively. Specifically for an observable depending only on the position x, we obtain (5.3) Tr ( e −βH A) The classical mean-field density µ mf in (5.3) can also be rewritten as a weighted average The weights q 0 and q 1 can be interpreted as the probability for the system to be in the corresponding electron eigenstate λ 0 and λ 1 , respectively obtained by integration with the corresponding Gibbs density. We first plot the equilibrium quantum mechanics density µ qm using the formula (5.2), and compare it with the classical mean-field density µ mf in (5.3). They are also compared with the density based only on the ground state µ gs , with the formula which is obtained from the classical density formula (5.4) by taking the probability for the excited state as q 1 = 0, and the probability for the ground state as q 0 = 1. The mean-field density (the dashed yellow line) is quite close to the quantum mechanics density curve (the solid blue line), implying a better accuracy than the ground state density (the solid violet line).
In Figure 2 the first reference density curve with quantum mechanics formula (5.2) is plotted in blue colour with a solid line. The density curve µ mf (x), obtained from the classical mean-field formula (5.3), is plotted as the yellow dashed line and agrees well with the quantum mechanics density µ qm (x). Besides, the mean-field density µ mf (x) incurs much smaller error than the ground-state density µ gs (x) (the violet solid curve) in approximating the µ qm (x) density. For Figure 2, we use the parameters M = 1000, c = 1, and β = 1 such that with the eigenvalue gap δ = 0.1, the system has a probability q 1 = 0.16 to be in the excited state.  Figure 3. The point-wise difference between the quantum mechanics density µ qm and the mean-field density µ mf with inverse temperature β = 1. The dashed violet curve with M = 100 has so small an error that it is almost indiscernible from the solid yellow curve with M = 1000.
In Figure 3 we depict a point-wise difference between the classical mean-field density µ mf (x) and the quantum mechanics density µ qm (x), for different values of the mass ratio M . The inverse temperature is still taken as β = 1 for the eigenvalue gap δ = 0.1 and c = 1, so that the probability for the excited state is kept as q 1 = 0.16. It is observed from Figure 3 that as M increases the error in the classical mean-field density approximation decreases.
In order to study the dependence of the approximation error µ qm − µ mf L 1 on the mass ratio M , we vary M for three different inverse temperatures, with the corresponding eigenvalue gaps δ such that the probability to be in the excited state remains to be q 1 = 0.16. As seen from Figure 4, the O(M −1 ) dependence of the error in the equilibrium density using the classical mean-field approximation is in accordance with the theoretical result of Theorem 2.1.
Besides the M -dependence of the classical approximation, we also experiment with a relatively large inverse temperature β = 10 for mass ratio M = 100, with parameters c = 1, δ = 0.1. The quantum density µ qm together with its classical mean-field and ground state approximations µ mf and µ gs are plotted in the Figure 5. The large value of β implies a rather low temperature, which leads to a tiny probability for the electron excited state as q 1 = 7 × 10 −7 . Consequently the density functions concentrate near the minimum of the ground state eigenvalue, and there is almost no difference between the mean-field and the ground state density curves.

The model problem.
We apply the mean-field molecular dynamics to approximate the auto-correlation function between the momentum observables p 0 (at time 0) and p τ (at time τ ). In the Heisenberg representation the time evolution of the momentum observable is given by   Figure 4. Dependence of the L 1 -error between the quantum density µ qm and the classical mean-field density approximation µ mf , shown in log-log scale.  Figure 5. Equilibrium density µ qm with the classical mean-field and ground state approximations µ mf and µ gs , with inverse temperature β = 10, mass ratio M = 100. The probability of electron excited state q 1 = 7 × 10 −7 is tiny.
We study the two-eigenvalue model with the potential matrix V (x) as defined by (1.13) in Subsection 1.2. For computing the quantum correlation function, we approximate the initial position observable x 0 and the initial momentum observable p 0 by the matrices respectively, where we discretize a sufficiently large computational domain Ω = [x 0 , x K ] with uniform grid points x k = x 0 + k∆x, for k = 0, 1, · · · , K and ∆x = |x K −x 0 | K . In the definition of the matrix P , the real symmetric matrix H d is of size 2(K + 1) × 2(K + 1), corresponding to a fourth order finite difference approximation of the Hamiltonian operator H. More details about this approximation are provided in Section A.1, and the definition of matrix H d is given in (A.2). The matrix H d generates the approximations We apply the eig function of Matlab to obtain the eigenpairs (e n , φ n ) of the H d matrix, and rearrange them to obtain the eigendecomposition Thus the right-hand side of (1.3) with A τ = p τ and B 0 = p 0 can be approximated by where in the second equality we use the cyclic property of the trace. By applying the mean-field molecular dynamics formula (1.4) with momentum observables p 0 and p τ in the model problem, we have the approximation for the time-correlation function as where z t := (x t , p t ) solves the Hamiltonian system with an initial state z 0 = (x 0 , p 0 ) ∈ R 2 . The Hamiltonian system (5.8) is solved numerically with the second-order velocity Verlet scheme, see [16]. More details about this numerical implementation is in Appendix A.2.
We also apply the classical molecular dynamics formula for correlation functions introduced in [14,Section 2.3.2], which considers the contribution from the ground state and the excited states. For our specific example the excited state dynamics approximation of momentum correlation observable is given by (5.9) T es (τ ) : with the weights q 0 and q 1 as defined in (5.4), where z j τ = (x j τ , p j τ ), j = 0, 1 solves the Hamiltonian dynamicsẋ j τ = p j τ , p j τ = −∇λ j (x j τ ) , with the initial condition z j 0 = (x 0 , p 0 ) = z 0 and λ 0 (x), λ 1 (x) as defined in (1.14). In addition to these three expressions for time-correlation, T qm (the quantum mechanics correlation), T mf (the mean-field approximation), and T es (the classical excited state approximation), we also compute the approximation based only on the ground state contribution, T gs , obtained by setting the probability q 1 for the excited state equal to be zero, i.e., (5.10) T gs (τ ) = where z τ = (x τ , p τ ) solves the Hamiltonian dynamics with the potential λ 0 (x) with initial state z 0 = (x 0 , p 0 ). As discussed in (1.16), the approximations with mean-field dynamics or excited state dynamics at an initial time τ = 0 are always identical, i.e., T mf (0) = T es (0). For positionrelated observables, the classical ground state dynamics approximation T gs (0) will be, in general, different from T mf (0) and T es (0), which is consistent with our preceding observations on equilibrium density function in Section 5.1, and is confirmed in the upcoming Figure 13.
In practise ground state molecular dynamics in the canonical ensemble is relatively well developed for realistic molecular systems and successful mean-field approximations appear in centroid and ring-polymer molecular dynamics [19]. Direct computations of excited state dynamics for realistic systems seem less attractive, due to the challenge to efficiently compute excited electron states, cf. [3]. Here we compare these three alternatives for a simple model problem in the hope of giving some information also on realistic systems.

5.2.2.
Numerical results. Following the discussion on the variances ǫ 2 1 and ǫ 2 2 with equations (1.8), (1.10) and (1.11), we survey five different cases with varying parameter settings of c, δ, and inverse temperature β, and make a comparison between the performances of different molecular dynamics approximations in each case. A summary of the parameters in each case is given in the Table 1. Case A: Low temperature with large eigenvalue gap, β = 3.3, c = 1, δ = 1, ǫ 2 1 = 9.95 × 10 −4 , ǫ 2 2 = 1.25 × 10 −4 , the probability for the excited state q 1 = 0.0002 is almost negligible. (d) Case D: high temperature, large difference between eigenvalues, medium gap Figure 6. Eigenvalues of the matrix-valued potential V (x) for test cases A to D. In subfigure (a), the mean-field potential λ * (x) (the red curve) is quite close to the ground state λ 0 (x) (the violet curve). Figure 6a presents the eigenvalues λ 0 (x), λ 1 (x) and the mean-field potential function λ * (x) as defined in (1.9) for Case A. With the parameters c = 1 and δ = 1, the system has a large eigenvalue gap. Particularly in this low temperature setting, the mean-field potential λ * (x) is almost identical to the ground-state eigenvalue λ 0 (x). Since the probability for the excited state is very small (q 1 = 0.0002), the three molecular dynamics approximations T mf (τ ), T es (τ ), and T gs (τ ) are similar.
In Figure 7a, the quantum mechanics correlation function curve S qm (τ ) with mass ratio M = 1000 is plotted as a function of correlation time τ , together with the three molecular dynamics approximations S mf (τ ), S es (τ ), and S gs (τ ). The three molecular dynamics correlation function curves are almost on top of each other, as shown in Figure 7b with similarly small errors. This case gives an example where all the three molecular dynamics work analogously, since q 1 ≪ 1, and we note that the error terms ǫ 2 1 , ǫ 2 2 together with 1/M are all very small. Case B: High temperature with small difference between eigenvalues, β = 1, c = 0.1, δ = 1, ǫ 2 1 = 1.9 × 10 −2 , ǫ 2 2 = 3.5 × 10 −3 , the probability for the excited state q 1 = 0.43.
In Figure 6b we observe that with the parameter setting of case B, the mean-field potential λ * (x) lies in between the ground state eigenvalue λ 0 (x) and the excited state eigenvalue   λ 1 (x), indicating that by incorporating the effect of the excited state, the mean-field approximation T mf (τ ) can make a difference from simply using the ground state molecular dynamics.
The improved accuracy of the mean-field molecular dynamics is verified by the Figure 8b, in which we observe a smaller error of mean-field molecular dynamics approximation S qm − S mf L ∞ ([0,τ ]) (the red curve) compared with the molecular dynamics using only the ground state S qm − S gs L ∞ ([0,τ ]) (the violet curve). The excited state molecular dynamics S es (τ ) has the smallest error, manifesting an effective combination of the information from both the ground and the excited eigenstates. Case C: High temperature, small difference between eigenvalues with avoided crossing , β = 1, c = 0.1, δ = 0.01, ǫ 2 1 = 9.1 × 10 −3 , ǫ 2 2 = 9.9 × 10 −3 , the probability for the excited state q 1 = 0.46.
The Case C has a similar parameter setting as the preceding Case B, with the only difference of a smaller parameter δ = 0.01. The small parameter δ leads to a small eigenvalue gap at x = 0, i.e., the two eigenvalues λ 0 (x) and λ 1 (x) almost intersect at this point, as can be seen in the Figure 6c. Compared with Case B, the small eigenvalue gap also makes the probability for the excited state q 1 increase from 0.43 to 0.46 in Case C, with the same inverse temperature β = 1.  The approximate p-auto-correlation function curves with their corresponding maximum errors up to time τ are plotted in the Figures 9a and 9b, respectively. These two figures are quite similar to their corresponding plots in Case B, where the excited state approximation S es has the smallest error, and the mean-field approximation S mf achieves an improved accuracy compared to the ground state approximation S gs .
The similar approximation error of the three molecular dynamics in Case B and Case C can be understood as a result of the relatively small difference between the two eigenvalues λ 0 and λ 1 . For both cases the small parameter c = 0.1 leads to small ǫ 1 and ǫ 2 values, as summarized in Table 1. Case D: High temperature, large difference between eigenvalues with large gap, β = 0.28, c = 1, δ = 1, ǫ 2 1 = 2.01, ǫ 2 2 = 0.42, the probability for the excited state q 1 = 0.30. For this case, we observe from Figure 6d that although the mean-field potential λ * (x) is still in between the ground state λ 0 (x) and the excited state λ 1 (x), the distance between λ * (x) and λ 1 (x) is much larger than that in the previous Case B. Also λ * (x) is much closer to the ground state λ 0 (x) than to the excited state λ 1 (x). The parameter β = 0.28 implies a relatively high temperature, with a considerable contribution from the excited state. Hence we cannot expect the mean-field molecular dynamics to be much better than ground state molecular dynamics. This is also verified by Figure 10b which shows that the error of meanfield molecular dynamics is of the same order as that of ground state molecular dynamics, while the excited state molecular dynamics remains accurate.
The mean-field and ground state molecular dynamics correlations include two approximations: replacing the matrix-valued potential V (x) by a scalar potential λ * (x) or λ 0 (x) respectively, and replacing quantum dynamics with classical dynamics, where S qm,λ * and S qm,λ 0 denote the approximation of auto-correlation function computed with quantum dynamics but using scalar-valued potentials λ * (x) and λ 0 (x), respectively. In the right hand side of (5.11), the first terms S qm (τ ) − S qm,λ * (τ ) and S qm (τ ) − S qm,λ 0 (τ ) correspond to the potential approximations in quantum dynamics, while the second terms S qm,λ * (τ ) − S mf (τ ) and S qm,λ 0 (τ ) − S gs (τ ) are related to classical approximations of quantum dynamics using scalar potentials.
To investigate these two error contributions we compute the correlation function S qm,λ * (τ ) and S qm,λ 0 (τ ) for Case D, using the scalar-valued potential λ * (x) or λ 0 (x) to replace the potential matrix V (x) in the quantum dynamics. The corresponding auto-correlation curves and their maximum error up to time τ are shown in Figures 11a and 11b.
From Figure 11b, we clearly see that the errors S qm − S qm,λ * L ∞ ([0,τ ]) and S qm − S qm,λ 0 L ∞ ([0,τ ]) caused by substituting the potential matrix V (x) with the scalar-valued potentials λ * (x) or λ 0 (x) is of the same order as the total errors S qm − S mf L ∞ ([0,τ ]) and S qm − S gs L ∞ ([0,τ ]) in Figure 10b. Hence we conclude that the main source of error in this case is the simplification by replacing the potential matrix V (x) with a scalar-valued potential, and not the approximation of scalar potential quantum mechanics correlation by classical molecular dynamics.
We also vary the mass ratio M between the heavy particle and the light particle, in order to study the corresponding behaviour of the approximation errors S qm (τ )−S qm,λ * (τ ) L ∞ ([0,τ ]) and S qm (τ ) − S qm,λ 0 (τ ) L ∞ ([0,τ ]) up to time τ = 20 for the mean-field molecular dynamics and ground state molecular dynamics in Case D. As can be seen from the second and fourth column of Table 2, the main error caused by substitution of the potential matrix V with the scalar valued potential λ * or λ 0 varies slightly as the M value changes.  Figure 11. (a) Auto-correlation function P 0 P τ curves S qm , S qm,λ 0 , and S qm,λ * computed by using matrix valued potential V (x), and using scalarvalued potential λ 0 (x) or λ * (x) in the quantum mechanics formula (5.6) in Case D. (b) Maximum error up to time τ in the p-auto-correlation curves computed with quantum mechanics formula using scalar-valued potential λ * (x) or λ 0 (x), by comparing them with the correlation computed from quantum mechanics formula using matrix-valued potential V (x), and with their corresponding molecular dynamics approximations.

M
S qm − S qm,λ * S qm,λ * − S mf S qm − S qm,λ 0 S qm,λ 0 − S gs S qm − S es Case E: High temperature, large difference between eigenvalues with avoided crossing , with β = 1, c = 1, δ = 0.1, ǫ 2 1 = 0.29, ǫ 2 2 = 0.50, the probability for the excited state q 1 = 0. 16. This case has the same parameters as in Section 5.1. In Figure 1 we observe a pattern of the two eigenvalues λ 0 (x) and λ 1 (x) related to avoided crossing of potential surfaces. Our numerical results suggest that for this case all three molecular dynamics are only accurate for short time, as can be seen in Figure 12a. Compared to Case C and Case D, where the excited state dynamics is accurate, the diminished eigenvalue regularity at the avoided crossing may explain the loss of accuracy in Case E.
Apart from the momentum auto-correlation function, we also computed in Case E the correlation function between position observables x 0 and x τ , as plotted in Figure 13. We observe that for short time range (e.g. 0 ≤ τ ≤ 0.1), the error of the ground state molecular dynamics is larger than the error of the mean-field molecular dynamics, which is consistent with the result for equilibrium observables in Section 5.1, in which the ground state molecular dynamics has larger error in approximating the density function µ qm (x) than the mean-field formula. Therefore the mean-field molecular dynamics can improve short-time  Figure 12. (a) Auto-correlation function P 0 P τ curves S qm , S qm,λ 0 , and S qm,λ * computed by using matrix valued potential V (x) and using scalarvalued potential λ 0 (x) or λ * (x) in the quantum mechanics formula (5.6) in Case E. (b) Maximum error up to time τ in the p-auto-correlation curves computed with quantum mechanics formula using scalar-valued potential λ * (x) or λ 0 (x), by comparing them with the correlation computed from quantum mechanics formula using matrix-valued potential V (x), and with their corresponding molecular dynamics approximations.
approximation of position auto-correlation function compared to the ground state molecular dynamics.  Figure 13. Case E: Auto-correlation function X 0 X τ computed by quantum-mechanics formula with M = 100 and by three molecular dynamics formulae.
We also changed the mass ratio M in this case for the momentum auto-correlation function, from M = 100 to M = 50 and to smaller value M = 20. When M becomes smaller, we can expect the error of molecular dynamics approximation becomes larger, since the error includes the O(M −1 ) term. For Case E, since we are only interested in the short time approximation, the time-dependent error term is not much larger than the O(M −1 ) term. Hence the effect of varying the mass ratio M will be considerable. The dependence of the L ∞ -error in momentum auto-correlation approximation on the mass ratio M is summarized in Table 3, from which we observe an improved accuracy in all the three molecular dynamics approximations with an increased M value.  Table 3. Case E: Dependence of the error on the mass ratio M at different correlation times τ .

Conclusion from numerical comparisons.
From the study with equilibrium observables in Section 5.1, we see that by considering the contributions of excited states the classical mean-field approximation of quantum mechanics density at equilibrium achieves a substantial improvement from the approximation which uses only the information from the ground state. The error of the mean-field approximation will decrease as the mass ratio M increases, following the O(M −1 ) relation. For the time-dependent observables, specifically by studying the momentum auto-correlation function, we know from Case A that for a low temperature setting with a large eigenvalue gap, where the probability for an excited state is small, all three molecular dynamics with the mean-field approximation, excited state approximation, or ground state approximation work similarly well.
From Case B and Case D, we observe that the error of the mean-field approximation decreases as the difference between two eigenvalues diminishes (i.e. parameter c becomes small). Furthermore, for Case B with a small difference between two eigenvalues, the meanfield approximation improves the accuracy of molecular dynamics compared to using the ground state only.
With a small eigenvalue difference, even including the avoided crossing in Case C the result is similar to Case B, that is the mean-field approximation is still more accurate than the ground state approximation.
From Case D we know that when the system temperature is high and the difference between two eigenvalues is not small, the excited state approximation outperforms both the mean-field and the ground state molecular dynamics.
From Case E we see that when the difference between the eigenvalues are sufficiently large and when the potential matrix includes avoided crossings, all three molecular dynamics approximations are accurate for a short time range only.
The small error terms ǫ 2 1 and ǫ 2 2 , in Table 1, for the Cases A, B and C is consistent with the actual error being small for the mean-field approximation, while in Cases D and E where the mean-field approximation error is large these error terms are in fact large. Therefore the experiments indicate that the error estimate could be useful to estimate the mean-field approximation error also for realistic problems when the quantum observable is not computable. Figures 11b and 12b together with Table 2 show that in Case D and Case E, where mean-field and ground state approximations are not accurate, the error of the mean-field and ground state molecular dynamics are dominated by the matrix valued potential replaced by a scalar potential on the quantum level, since the classical approximation error of the quantum dynamics for the corresponding scalar potential is clearly smaller. and the eigenfunctions Φ n are approximated with the 2(K + 1)-length vector φ n , as φ n = [φ n,0,1 , φ n,0,2 , φ n,1,1 , φ n,1,2 , · · · , φ n,K,1 , φ n,K,2 ] T .