{'fraction_non_alphanumeric': 0.051949751373985865, 'fraction_numerical': 0.018188955770740645, 'mean_word_length': 3.9412905728695202, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Classical Langevin dynamics derived from quantum mechanics

The classical work by Zwanzig [J. Stat. Phys. 9 (1973) 215-220] derived Langevin dynamics from a Hamiltonian system of a heavy particle coupled to a heat bath. This work extends Zwanzig's model to a quantum system and formulates a more general coupling between a particle system and a heat bath. The main result proves that ab initio Langevin molecular dynamics, with a certain rank one friction matrix determined by the coupling, approximates for any temperature canonical quantum observables, based on the system coordinates, more accurately than any Hamiltonian system in these coordinates, for large mass ratio between the system and the heat bath nuclei.

Langevin dynamics for (unit mass) particle systems with position coordinates X L : [0, ∞) → R N and momentum coordinates P L : [0, ∞) → R N , defined by dX L (t) = P L (t)dt dP L (t) = −∇λ X L (t) dt −κP L (t)dt + (2κT ) 1/2 dW (t) , (1.1) is used for instance to simulate molecular dynamics in the canonical ensemble of constant temperature T , volume and number of particles, where W denotes the standard Wiener process with N independent components. The purpose of this work is to precisely determine both the potential λ : R N → R and the friction matrixκ ∈ R N ×N in this equation, from a quantum mechanical model of a molecular system including weak coupling to a heat bath. Molecular systems are described by the Schrödinger equation with a potential based on Coulomb interaction of all nuclei and electrons in the system. This quantum mechanical model is complete in the sense that no unknown parameters enter -the observables in the canonical ensemble are determined from the Hamiltonian and the temperature. The classical limit of the quantum formulation, yields an accurate approximation of the observables based on the nuclei only, for large nuclei-electron mass ratio M. Ab initio molecular dynamics based on the electron ground state eigenvalue can be used when the temperature is low compared to the first electron eigenvalue gap. A certain weighted average of different ab initio dynamics, corresponding to each electron eigenvalue, approximates quantum observables for any temperature, see [10], also in the case of observables including time correlation and many particles. The elimination of the electrons provides a substantial computational reduction, making it possible to simulate large molecular systems, cf. [18].
In molecular dynamics simulations one often wants to determine properties of a large macroscopic system with many particles, say N ∼ 10 23 . Such large particle systems cannot yet be simulated on a computer and one may then ask for a setting where a smaller system has similar properties as the large. Therefore, we seek an equilibrium density that has the property that the marginal distribution for a subsystem has the same density as the whole system. In [10] it is motivated how this assumption leads to the Gibbs measure, i.e. the canonical ensemble; this is also the motivation to use the canonical ensemble for the composite system in this work, although some studies on heat bath models use the microcanonical ensemble for the composite whole system.
Langevin dynamics is often introduced to sample initial configurations from the Gibbs distribution and to avoid to simulate the dynamics of all heat bath particles. The friction/damping parameter in the Langevin equation is then typically set small enough to not perturb the dynamics too much and large enough to avoid long sampling times.
The purpose of this work is to show that Langevin molecular dynamics, for the non heat bath nuclei, with a certain friction/damping parameter determined from the Hamiltonian, approximates the quantum system, in the canonical ensemble for any temperature, more accurately than any Hamiltonian dynamics (for the non heat bath particles) in the case the system is weakly coupled to a heat bath of many fast particles. Our heat bath model is based on the assumption of weak coupling -in the sense that the perturbation in the system from the heat bath is small and vice versa -which we show leads to Zwanzig's model for nonlinear generalized Langevin equations in [27], with a harmonic oscillator heat bath. Zwanzig also derives a pure Langevin equation: he assumes first a continuous Debye distribution of the eigenvalues of the heat bath potential energy quadratic form; in the next step he lets the coupling to the heat bath have a special form, so that the integral kernel for the friction term in the generalized Langevin equation becomes a Dirac-delta measure. Our derivation uses a heat bath based on nearest neighbour interaction on an infinite cubic lattice, so that the continuum distribution of eigenvalues is rigorously obtained by considering a difference operator with an infinite number of nodes. Our convergence towards a pure Langevin equation is not based on a Dirac-delta measure for the integral kernel but obtained from the time scale separation of the fast heat bath particles and the slower system particles, which allows a general coupling and a bounded covariance matrix for the fluctuations. The fast heat bath dynamics is provided either from light heat bath particles or a stiff heat bath potential. By a stiff heat bath we mean that the smallest eigenvalue of the Hessian of the heat bath potential energy is of the order χ −1 , where χ ≪ 1. We show that the friction/damping coefficient in the pure Langevin equation is determined by the derivative of the forces on the heat bath particles with respect to the system particle positions. We also prove that the observables of the system coordinates in the system-heat bath quantum model can be approximated using this Langevin dynamics with accuracy O(m log m −1 + (mM) −1 ), where M ≫ 1 is the system nuclei-electron mass ratio and m ≪ 1 is the heat bath nuclei -system nuclei mass ratio; the approximation by a Hamiltonian system yields the corresponding larger error estimate O(m 1/2 + (mM) −1 ), in the case of light heat bath particles. In this sense, our Langevin equation is a better approximation. The case with a stiff heat bath system has the analogous error estimate O(χ 2δ log χ −1 + M −1 ) with the correct friction/damping parameter, while a Hamiltonian system gives the larger error O(χ 2δ−1/2 + M −1 ), where χ δ ≪ χ 1/4 measures the coupling between the heat bath and the system. Our main assumptions are: • the coupling between the system and the heat bath is weak and localized, • either the heat bath particles are much lighter than the system particles or the heat bath is stiff, • the harmonic oscillator heat bath is constructed from nearest neighbour interaction on an infinite cubic lattice in dimension three, • the heat bath particles are initially randomly Gibbs distributed (conditioned by the system particle coordinates), and • the system potential energy and the observables are sufficiently regular. The system and the heat bath is modelled by a Hamiltonian where the system potential energy is perturbed byV (x, X) with system particle positions X ∈ R N and heat bath particle positions x ∈ R n . Our assumption of weak coupling is formulated as the requirement that the potentialV satisfies min x∈R nV (x, X) =V a(X), X = 0 . Theorem 3.7 proves that the obtained minimizer a(X) ∈ R n determines the friction matrixκ in (1.1) by the rank one N × N matrix whereV and a have limits as n → ∞; hereV ′′ denotes the Hessian ofV with respect to x, the brackets ·, · denote the scalar product in R n andc is a constant (related to the density of states).
In a setting when the temperature is small compared to the difference of the two smallest eigenvalues of the Hamiltonian symbol, with given nuclei coordinates, it is well known that ab initio molecular dynamics is based on the ground state electron eigenvalue as the potential λ in (1.1), cf. [18]. When the temperature is larger, excited electron states influence the nuclei dynamics. The work [10] derives a molecular dynamics approximation of quantum observables, including time correlations, as a certain weighted average of ab initio observables in the excited states, and these ideas are put into the context of this work in Section 5.
The analysis of particles in a heat bath has a long history, starting with the work by Einstein and Smoluchowski. The Langevin equation was introduced in [14] to study Brownian motion mathematically, before the Wiener process was available. Early results on the elimination of the heat bath degrees of freedom to obtain a Langevin equation are [7,6] and [27], [26,Section 9.3]. The work [7,6] include in addition a derivation of a quantum Langevin equation, which is based on an operator version of the classical Langevin equation. Our study of the classical Langevin equation from quantum mechanics is not related to this quantum Langevin equation. We start with an ab initio quantum model of the system coupled to a heat bath in the canonical ensemble and use its classical limit to eliminate the electron degrees of freedom. Then the classical system coupled to the heat bath is analyzed by separation of time scales. The separation of time scales of light bath particles and heavier system particles was first used in [15], see also [20,23] to determine a Fokker-Planck equation for the system particle, from the Liouville equation of the coupled system using a formal expansion in the small mass ratio. Section 8 in [22] presents a proof, and several references to related work, where the Langevin equation is derived from the generalized Langevin equation, using exponential decay of the kernel in the memory term of the generalized Langevin equation.
Our contribution employs the separation of time scales approach previously used in [15,20,23], but a novelty is that we here present mathematical proofs for the weak convergence rate of Langevin dynamics towards quantum mechanics. Our work also differs from [22,Section 8], wherein the kernel of the generalized Langevin equation is assumed to be such that by adding a finite dimensional variable the system becomes Markovian and the kernel tends to a point mass. For instance, our kernel vanishes as m → 0 + and we use the density of heat bath states to determine the kernel. Using the precise information from the density of states in the case of heat bath nearest neighbor interaction in a cubic lattice we obtain a positive definite friction matrixκ, while if the nearest neighbor interaction would be related to a lattice in dimension four the friction matrix would vanish, see Remark 3.4.
The main new ideas in our work are the first principles formulation from quantum mechanics, the weak coupling condition as a minimization, the precise use of the density of heat bath states, that the error estimate uses stability of the Kolmogorov backward equation for the Langevin equation evaluated along the dynamics of the coupled system, and formulation of numerical schemes and numerical results related to Langevin dynamics approximation of particles systems. Section 2 formulates a classical model for the system and the heat bath, including the weak coupling, and derives the corresponding generalized Langevin equation, based on a memory term with a specific integral kernel, in a classical molecular dynamics setting. The generalized Langevin equation is analyzed in Section 3, with subsections related to the dissipation, fluctuations and approximation by pure Langevin dynamics in the case of light heat bath particles. The main result of Section 3 is Theorem 3.7, where a certain Langevin dynamics is shown to approximate a classical system weakly coupled to a heat bath. Section 4 extends the analysis to the case of stiff heat baths and derives a corresponding approximation result in Theorem 4.1. Section 5 relates the classical model to a quantum formulation and provides background and error estimates on quantum observables in the canonical ensemble approximated by classical molecular dynamics. The main result of the work is Theorem 5.2, which proves that for any temperature canonical quantum observables based on the system coordinates can be accurately approximates by the Langevin dynamics obtained in Theorems 3.7 and 4.1. Section 6 includes numerical results of the system particle density and autocorrelation of the coupled system and heat bath approximated by the particle density and autocorrelation for the Langevin dynamics.

The model of the system and the heat bath
We consider in this section a classical model of a molecular system, with position coordinates X ∈ R N and momentum coordinates P ∈ R N , coupled to a heat bath, with position and momentum coordinates x ∈ R n and p ∈ R n , respectively, represented by the Hamiltonian where λ : R N → R is the potential energy for the system andV : R n × R N → R is the potential energy for the heat bath including the coupling to the system. The parameter m is the mass ratio between heat bath nuclei and system nuclei. We have set the time scale so that the system nuclei mass is one. In Section 5 we show that this model is the classical limit of a quantum model and study the accuracy of the classical approximation, also in the case with system nuclei that have different masses.
We study small perturbations of the equilibrium bath state x = a(X) ∈ R n where min x∈R nV (x, X) =V a(X), X . Taylor expansion around the equilibrium yields for some ξ ∈ [0, 1], whereV xx is the Hessian matrix in R n×n , with respect to the x coordinate, and the notation ·, · is the standard scalar product in C n . We assume that the coupling between the system and the heat bath is weak, which means that the perturbation in the system from the heat bath must be small. As the perturbation |x − a(X)| → 0 we therefore require that Weak coupling also means that the perturbation in the heat bath from the system is small, so that the influence on the HessianV ′′ xx from X is negligible, and we assume therefore thatV ′′ xx a(X) + ξ(x − a), X = C (2.2) where C is constant, symmetric and positive definite. Hence the vibration frequencies for x are assumed to be constant for all x and X. We use the notation C =V ′′ below.
The assumptions (2.1) and (2.2) of weak coupling lead to the Hamiltonian H : This model (2.3) is of the same form as the model of interaction with a heat bath introduced and analysed by Zwanzig in the seminal work [27], although here the motivation with weak coupling and several system particles is different. The Hamiltonian (2.3) yields the dynamicṡ and the change of variables (mV ′′ ) 1/2 r t := p t implieṡ Define ϕ t := x t − a(X t ) + ir t to obtaiṅ where the third equation uses the notation v · w for the standard scalar product in R N .
Assume that x 0 and p 0 are Gaussian with the distributions provided by the marginals of the Gibbs density e −H(X,P,x,p)/T R 2N+2n e −H(X,P,x,p)/T dXdP dxdp , that is, the momentum p is multivariate normal distributed with mean zero and covariance matrix mT I and independent of x, which is multivariate normal distributed with mean a(X) and covariance matrix T (V ′′ ) −1 . Consequently the initial data can be written using the orthogonal eigenvectors ν ′ m , normalized as ν ′ k , ν ′ k = 1, and eigenvalues µ m of V ′′V and the definition γ m := γ r m + iγ i m (2.7) with γ r m and γ i m , m = 1, . . . , n, independent and normal distributed real scalar random numbers with mean zero and variance T /µ m .
Duhamel's principle shows that which implies a form of Zwanzig's generalized Langevin equation with non Markovian friction term given by the integral and a noise term including the stochastic initial data ϕ 0 . We study two different cases: • either the mass ratio m ≪ 1 is small, or • the smallest eigenvalue ofV ′′ is large of the order χ −1 while the coupling derivative ∇a is small of size χ δ with δ > 1/4. In these cases, both the friction and the noise terms are based on highly oscillatory functions, which will make these contributions small, as explained in the next section. To simplify the analysis, we assume also that ∇a(·) is constant. (2.9) 3. Analysis of the generalized Langevin equation for m ≪ 1 In this section we first study the dissipation term ȧ(X s ), ∇a ds and the fluctuation term Re V ′′ e −iV ′′1/2 m 1/2 t ϕ 0 , ∇a in (2.8), as the number of bath particles tend to infinity and the mass ratio, m, between light heat bath nuclei and heavier system nuclei is small. Then in Section 3.3 we prove an priori estimate of sup s<t E[|X s | 2 + |Ẋ s | 2 + |Ẍ s | 2 ]. Section 3.4 uses the limit terms to construct a Langevin equation, namely the Itô stochastic differential equation dX L (t) = P L (t)dt with a certain symmetric friction matrix κ ∈ R N ×N and a Wiener process W : [0, ∞) × Ω → R N , with N independent components; here Ω is the set of outcomes for the process (X, P ) : [0, ∞) × Ω → R 2N . Finally, we use the solution of the Kolmogorov backward equation for the Langevin dynamics along a solution path (X t , P t ) of (2.8) to derive an error estimate of the approximation, namely E[g(X t , P t )] − E[g(X L (t), P L (t)] = O(m log m −1 ), for any given smooth bounded observable g : R 2N → R and equal initial data (X(0), P (0)) = (X L (0), P L (0)).
3.1. The dissipation term and a precise heat bath. The change of variables We require this dissipation term to be small, so that the coupling to the heat bath yields a small perturbation of the dynamics for X and P . IfV ′′ and ∇a are of order one, the mass m needs to be small, or if ∇a is small we can haveV ′′ large. The case with small m is studied in this section and the case with largeV ′′ is in Section 4. A small mass also requires the integrand to decay as τ → ∞. We will use Fourier analysis to study the decay of the kernel by writing lim n→∞ K n as an integral, which is the next step in the analysis. The kernel K n is based on the equilibrium heat bath position derivative ∂ X ℓ a and this derivative is determined by the derivative of the force on bath particle x j , with respect to X ℓ , by (2.3) as which, for fixed ℓ and j, by assumptions (2.2) and (2.9) is independent of x and X. We assume thatV is constructed so that this force derivative is localized in the sense Let ν k (j) := ν ′ k (j)/ max j |ν ′ k (j)| be the set of orthogonal eigenvectors ofV ′′ that are normalized to one in the maximum norm.The eigenvalue representation ofV ′′ in (2.6) and the definition If n is finite and we make a tiny perturbation ofV ′′ so that all ω k become rational, the function K n will be periodic in τ and consequently it will not decay for large τ . To obtain a decaying kernel we will therefore consider a heat bath with infinite number of particles, n = ∞. The next step is consequently to study the limit of K n , as n → ∞, which requires a more precise formulation of the heat bath.

A precise heat bath. We use the periodic lattice
in dimension three to form the equilibrium positions for the heat bath. The position deviation from the equilibrium, namelyx j := x j − a j (X) ∈ R for each particle j ∈ En, then determines the potential by nearest neighbour interaction in the lattice with periodic boundary conditionsx j+ne i =x j , where j = (j 1 , j 2 , j 3 ), e 1 = (1, 0, 0), e 2 = (0, 1, 0), e 3 = (0, 0, 1),n 3 = n and c is a positive constant independent of n. The small positive constant η n is introduced to make the potential strictly convex for finite n, while it vanishes asymptotically and satisfies n 1/2 η n → ∞ and η n → 0 + as n → ∞ . (3.6) We will see in (3.19) that the zero limit in (3.6) is needed in our model to obtain non zero friction matrices κ. The lower bound in (3.6) implies by (3.3) that ∂ X ℓ a ℓ ∞ remains bounded, provided β ℓ,· ℓ 1 is bounded. We note thatV ′′ consists of the standard finite difference matrix with mesh size one related to the Laplacian in R 3 and a small positive definite perturbation η 2 n I, where I is the corresponding identity matrix. The minimum of the potential is obtained for the position deviationsx j = 0, j ∈ En, which we may view as the n heat bath particles located on the n different lattice points in En. This heat bath model can be extended to positionsx j ∈ R 3 , see Remark 3.1.
we still have the same eigenvalues if c 2 k = c 2 , so the model does not change in principle. If on the other hand the constants c k are different for the different components ofx j the spectrum may change and we obtain a different heat bath model. 3.1.2. The limit friction matrix. We can take the limit asn → ∞ in (3.4), while r = k/n is constant, to obtain an integral (3.10) The change of variables ω i = 2 1/2 c 1 − cos(2πr i ) sgn(r i ) yields, with the spherical coordinate ω = (cos α cos θ, sin α cos θ, sin θ)ω, For the purpose of simplifying later proofs we will assume that β ℓ (ω) is two times differentiable and that it vanishes at the boundary of ω( We also note that the constant c and the density of states are related by We are now ready to formulate the limit as n → ∞ in the friction term based on the constant N × N friction matrix κ. (3.14) Remark 3.3 (Dirchlet boundary condition). If we replace the periodic boundary conditions in the heat bath model (3.7) with homogenous Dirichlet conditions, we have instead . . ,n/2 − 1, and k i = 1, . . . ,n , We can write the eigenvectors as functions of k i /(n + 1) since and split the sum over k i into one sum over odd k i , where the eigenvector is based on cosine functions, and one sum over even k i , where the eigenvector is based on sine functions. With these changes, the derivation of κ follows as in the case with periodic boundary conditions.

Proof of Lemma 3.2.
To study the decay as τ → ∞ of the kernel K ∞ , we use (3.12) and the shorthand f (ω) := f (ω, ·, ·) and integrate by parts Since f ω has its support in ω = |ω| ≤ c, cf. (3.12), and the second derivative For a given t > 0 and bounded function h : R N ×R N → R we next study the expectation of In the case t/ √ m > 1, we split the integral with 1 < τ (3.16) The magnitude of the second integral is bounded in expectation using (3.4), (3.10) and (3.15): The expectation and the limit n → ∞ of the first integral in the right hand side of (3.16) can be written as where, using (3.15), so that with τ * = m − 1 2 the error in (3.14) is bounded by The Fourier transform of f (ω) with respect to ω is integrable and since also f (ω) is continuous, we have by the Fourier inversion property and (3.9) That is, the friction matrix, κ, in the Langevin equation is determined by the ∂ X ℓderivative of the sum of forces on all bath particles and we have proved Lemma 3.2.
where S d is a positive constant related to the dimension d. We see that under the assumption β ℓ (0 + )β ℓ ′ (0 + ) > 0, it is only in dimension d = 3 that this heat bath generates a positive definite friction matrix κ. For d < 3 we obtain κ = ∞ and for d > 3 we have κ = 0.
If we change the heat bath potential energy to be based on any circulant matrix in each dimension, we have the requirement lim ω→0 r 2 (ω) dr dω ω 2 = constant (3.20) to obtain a positive definite friction matrix κ, in dimension d = 3. This limit for the eigenvalues ω 2 implies thatV ′′ becomes a difference quotient approximation of the Laplacian with mesh size one. We conclude that (3.20) leads to a choice ofV ′′ in (3.5) that is another discretization of the Laplacian or discretizations that tends to the Laplacian as n → ∞.

3.2.
The fluctuation term. The initial distribution of ϕ 0 = k (γ r k +iγ i k )ν ′ k , determined by the Gibbs distribution of the Hamiltonian system in (2.5) and (2.7), shows that all γ r k and γ i k are independent and normal distributed with mean zero and variance T /ω 2 k . This initial data provides the fluctuation term in (2.8), namely which has the special property that its covariance satisfies that is, the covariance is T times the friction integral kernel, as observed in [27]. We repeat a proof for our setting here, where we also show that the limit of ζ t as n → ∞ is well defined.
Proof of (3.22). Let e −it is the set of normalized real valued eigenvectors ofV ′′ , defined in (2.5). The random coefficients γ ′ k = γ ′ k,r + iγ ′ k,i are based on γ ′ k,r and γ ′ k,i which are independent normal distributioned with mean zero and variance T /ω 2 k , for k = (k 1 , k 2 , k 3 ) and k i = 1, . . . ,n. Use the orthonormal real valued basis {ν ′ k } to obtain We conclude that ζ is a Gaussian process with mean zero and covariance (3.22). The fluctuation-dissipation property (3.22), the limit (3.10) and the bound (3.15) show that the covariance matrix has a limit as n → ∞: Therefore ζ has a well defined limit, in the L 2 -space with norm ( τ 0 E[|ζ t | 2 ]dt) 1/2 , for any finite time τ as n → ∞.
3.3. The system dynamics. We can write the system dynamics (2.4) and (2.8) as We assume that λ is three times continuously differentiable and sup X D 2 λ(X) + sup X D 3 λ(X) is bounded (3.25) and apply Cauchy's inequality on (3.24) to obtain ≤T trace(K∞(0)) ds .
Here D k λ denotes the set of all partial derivatives of order k and if k = 1 or k = 2 we identify it with the gradient and the Hessian, respectivly. Integration yields the Gronwall inequality which by (3.26) and the mean square oḟ X L (0) = X(0) , To analyse the approximation we use that, for any infinitely differentiable function g : R 2N → R with compact support, (3.27) the expected value u(z, s) = E[g X L (t * ), P L (t * ) | X L (s), P L (s) = z] , s < t * and z ∈ R 2N , is well defined, since the Langevin equation (3.1), for Z t := X L (t), P L (t) , has Lipschitz continuous drift and constant diffusion coefficient. The assumption (3.25) implies that the stochastic flows ∂Zt ∂Zs and ∂ 2 Zt ∂Z 2 s also are well defined, which implies that the function u has bounded and continuous derivatives up to order two in z and to order one in t, see [12] and [8], and solves the corresponding Kolmogorov equation ∂ s u(X, P, s) + P · ∇ X u(X, P, s) − ∇λ(X) + √ mκP · ∇ P u(X, P, s) u(X, P, t * ) = g(X, P ) .

(3.29)
and where the rank one friction matrix is determined by the forceF from (3.3) and (3.8)

(3.31)
We see also that the alternative dynamics, with any friction coefficientκ ≥ 0, dX t =P t dt dP t = −∇λ(X t )dt −κP t dt + (2κT ) 1/2 dW t approximates (2.4) and (2.8) with the larger error Proof of (3.30). We have (3.33) and (3.29) shows that the second term cancels √ mκP t · ∇ P u(X t , P t , t) to leading order in (3.28). It remains to show that the scalar product of ∇ P u(X t , P t , t) and the first term in the right hand side (3.33) cancels the last term in (3.28).
To analyze the dependence between ∂ P u(X t , P t , t) and Re e −itV ′′ /m 1/2 ϕ 0 ,V ′′ ∇a from the initial data ϕ 0 we will study how a small perturbation of ϕ 0 influences ∂ P u(X t , P t , t).

(3.35)
Proof of the claim. The linearized problem corresponding to (2.8) becomes which shows that (3.34) satisfies the linearized equation (3.36), and the claim is proved. The main reason we use assumption (2.9), namely that ∇a(X) is constant, is to obtain this linearized equation for G, which otherwise would include the term Re ϕ,V ′′ D 2 a that would introduce the fast time scale of ϕ in G, which we now avoid.

The difference of the two Green's functions satisfy by (3.34) the perturbation representation
with analogous splittings for and we obtain then as in the proof of (3.22)

Analysis of the generalized Langevin equation for χ ≪ 1
In this section we study a time scale separation for system particles and fast heat bath particles due to a stiff heat bath obtained by changing the scaling tō k are the eigenvalues ofṼ ′′ , The nonlinear generalized Langevin equation then becomes . The analysis in Section 3 can be applied by replacing m by χm;V ′′ byṼ ′′ ; and a byã and we obtain where the friction matrix is determined by the forceF ℓj = lim n→∞ (Ṽ ′′ ∂ X ℓã )(j) as (4.4)

Molecular dynamics approximation of a quantum system
The purpose of this section is to present a molecular dynamics approximation for observables of a quantum particle system consisting of nuclei and electrons coupled to a heat bath. The observables may include correlations in time. The first subsection provides background to quantum observables approximated by molecular dynamics in the canonical ensemble. The next subsection combines these quantum approximation results of with the classical Langevin approximation in Theorems 3.7 and 4.1.

Canonical quantum observables approximated by molecular dynamics.
The quantum formulation is based on wave functions Φ : R N × R n → C d and the Hamiltonian where X ∈ R N and x ∈ R n are the nuclei coordinates of the system and heat bath positions, respectively, and M s and M b are the diagonal matrices of the mass of the system nuclei and heat bath nuclei, respectively, measured in units of the electron mass. The functions V : R N → C d 2 and V b : R n ×R N → C d 2 are finite difference approximations of the electron kinetic energy and nuclei-nuclei, nuclei-electron and electron-electron interactions, related to the system and the heat bath. The matrix I is the identity on C d . This simplification to replace the Laplacians for the electron kinetic energy by difference approximations makes it easier to derive the classical limit. Another simplification is to change coordinates M 1/2X = M 1/2 s X and M 1/2x = M 1/2 b x and letx = (X,x) ∈ R N +n , where M ≫ 1 is a reference nuclei-electron mass ratio. In these coordinates the Hamiltonian takes the formĤ ). We will use the eigenvalues λ k (x) ∈ R and eigenvectors ψ k (x) of the Hermitian matrix v(x) defined by We assume that the eigenvalues satisfy The first assumption is in order to have differentiable eigenvectors and the second condition implies that the spectrum ofĤ is discrete, see [4].
The aim here is to study canonical quantum observables, including correlations in time, namely where {Φ n } ∞ n=1 is a normalized basis of L 2 (dx), e.g. the set if normalized eigenfunctions toĤ and (f, g) := R N+n f * (x)g(x)dx. An operatorB is the Weyl quantization that maps L 2 (R N +n ) to [L 2 (R N +n )] d and is defined, from a d × d matrix valued symbol B : R N +n → C d 2 in the Schwartz class, bŷ For instance, we have |p| 2 2 + v(x) =Ĥ. The time dependent operatorB τ is defined bŷ which implies the von Neumann-Heisenberg equation where [Ĥ,B τ ] :=ĤB τ −B τĤ is the commutator. The example of the observable for the diffusion constant 1 6τ uses the time-correlationx(τ ) ·x(0) whereÂ τ =x τ I andĈ 0 =x 0 I and A main tool to determine the classical limit is to diagonalize (5.3), which is based on the following composition of Weyl quantizations: the symbol C for the product of two Weyl operatorsÂB =Ĉ is determined by see [28]. Assume that Ψ : R N +n → C d 2 and Ψ(x) is any unitary matrix with the Hermitian transpose Ψ * (x) and defineĀ : The composition rule (5.4) impliesΨ * ĤΨ = (Ψ * #H#Ψ) andĀ τ = Ψ * #A τ #Ψ. The next step is to determine Ψ so thatH := Ψ * #H#Ψ is almost diagonal. HavingH diagonal implies thatĤ is diagonal and thenÂ remains diagonal if it initially were diagonal, since then d dτÂ The composition rule (5.4) with as verified in [10, Lemma 3.1]. Therefore, the aim is to choose the unitary matrix Ψ so that it becomes an approximate solution to the nonlinear eigenvalue problem v + 1 4M Ψ∇Ψ * · ∇ΨΨ * Ψ = ΨΛ (5.5) in the sense that v The following result for approximating non equilibrium quantum observables by classical molecular dynamics observables is proved in [10].

5.2.
Langevin dynamics derived from quantum mechanics. Assume that the potential V (X) has the eigenvalues λ j (X) and eigenvectors ψ j (X), j = 1, . . . , d. The eigenvalues of the potential v(x(X, x)) = V (X)+V b (x, X) will to leading order in V b (·, X) be given by λ j (X)+ψ * j V b (x, X)ψ j . We assume now that all coupling potentials ψ * j V b (x, X)ψ j satisfy the weak coupling assumptions (2.1) and (2.2), namely ′′ xx ψ j is a constant matrix as in (2.2) and each coupling ψ * j V b ψ j for j = 1, . . . , d yields one equilibrium position a j (X) and one coupling matrixV ′′ j . The eigenvalues of the Hamiltonian symbol are then Let m be the mass ratio for a reference heat bath nuclei to a reference system nuclei. Theorems 3.7 and 4.1 show the classical dynamics provided by the HamiltoniansH jj are accurately approximated by the Langevin dynamics for m ≪ 1. The next step is to relate this approximation by Langevin dynamics also to quantum observables in the canonical ensemble. In particular we need to obtain the Gibbs distribution of the heat bath particles from the quantum observables in the canonical ensemble.
Define the probability, q j , to be in electron state j as then the molecular dynamics observable becomes a sum of observables in the different electron states with initial Gibbs distribution, namely For instance if only the ground state matters we have q 1 = 1 and q k = 0, k = 2, 3, 4, . . . , d.
We simplify by letting the system nuclei have the same mass M and the heat bath nuclei the same mass mM. Then the diagonalized Hamiltonian can by (2.3) be written asH wherez = (X, P, x, p) ∈ R N × R N × R n × R n . Assume that the observablesĀ jj andB jj only depend on the system coordinates X and P , then the classical molecular dynamics approximation of the canonical quantum observables in Theorem 5.1 satisfies where the expected value is with respect to the Gibbs measure of the heat bath coordinates conditioned on the initial system coordinates. We note that the Gibbs measure (5.11) is the invariant measure used to sample the initial heat bath configurations in theorems 3.7 and 4.1. Therefore the combintion of theorems 3.7, 4.1 and 5.1 show that canonical observables for a quantum system coupled to a heat bath can be accurately approximated by Langevin dynamics, with the friction coefficient determined by (3.31) and (4.4), including several electron surfaces λ j , j = 1, . . . , d.  and (X j t , P j t ), for t > 0, is the solution to the Langevin equation dX j t = P j t dt dP j t = −∇λ j X j t dt −mκ j P j t dt + (2mκ j T ) 1/2 dW t , with initial data (X j 0 , P j 0 ) = (X 0 , P 0 ), and c j set by the density of heat bath states at zero frequency in (3.13). The probability q j to be in electron state j is determined by (5.10) and the expected value is with respect to the Wiener process W , with N independent components.
Since this contradicts that β is twice differentiable, an assumption that was used in the proof of Lemma 3.2, let us demonstrate that said lemma also applies in the current setting with jump-discontinuous β: where the last equality follows from ∞ 0 sin τ τ dτ = π/2 andẊ t being continuous with respect to t. We conclude that To prove the convergence (3.14), observe first that Introducing the mesh τ k = 2kπ anḋ it follows that for any τ ∈ [τ k , τ k+1 ), there exists a random s ∈ [τ k , τ k+1 ) such that By the splitting (3.16), and (3.14) follows by the same reasoning as in the proof of Lemma 3.2.
For the heat bath dynamics we construct a splitting scheme which for a uniform mesh t k = k∆t computes the position at every timestep (X 1 (t 0 ), X 1 (t 1 ), . . .) and the momentum at every half-timestep (P 1 (t 0 ), P 1 (t 1/2 ), P 1 (t 1 ), . . .). The damping term's integral is approximated as follows: We make use of the above approximations of the damping term integral in the following splitting scheme for the heat bath dynamics: P 1 k ∆Si n−k + n−1 k=0 P 1 k+1/2 ∆Si n+1/2−k X 1 n+1 = X 1 n + P 1 n+1/2 ∆t P 1 n+1 = P 1 n+1/2 − ∆t 2 X 1 n+1 − ζ 1 (t n+1 ) + 12π √ m n k=0 P 1 k ∆Si n+1/2−k + P 1 k+1/2 ∆Si n−k . (6.5) The mean-zero Gaussian vector ζ = (ζ 1 (t 0 ), ζ 1 (t 1 ), . . . , ζ 1 (t N )) is sampled by computing the square root of the Toeplitz matrix with first row vector (T K ∞ (0 + / √ m), T K ∞ (t 1 / √ m), . . . , T K ∞ (t N / √ m)) and multiplying the square root matrix, say √ K, with an (N +1)−vector ξ of iid standard normals components: ζ = √ Kξ. See [11] for further details on sampling of Gaussian processes, [5,24,16,19,3,1,17] for numerical methods for Langevin dynamics and [2,9,13] and [16,Chapter 8.7] for an alternative numerical method for generalized Langevin equations based on a truncated prony series approximation of the kernel K ∞ . 6.3. Observations. By setting the temperature to T = 3, the stationary distribution of the exact Langevin dynamics (6.3) becomes N(0, I) for any m > 0, with I denoting the identity matrix in R 2 . In order to reduce the computational challenges of long time numerical integration, we sample the initial data from this stationary, i.e., (ξ X , ξ P ) ∼ N(0, I), as we assume this yields initial data for the numerical dynamics (X 1 L,0 , P 1 L,0 ) and (X 1 0 , P 1 0 ) that are very close to their respective stationary distributions. The numerical computations are performed using a timestep ∆t = 0.005, N = 4000 integration steps relating to the final time t N = 20, and M = 2 × 10 5 solution realizations of the respective dynamics. Figures 1 and 2 show good correspondence between the final time marginal distributions of the heat bath dynamics (X 1 N , P 1 N ) and the Langevin dynamics (X 1 L,N , P 1 L,N ) over a range of m-values (all marginals being approximately N(0, 1)-distributed). Figures 3  and 4 show that the auto-correlation for both the position and the momentum of the heat bath dynamics converges to the corresponding ones for the Langevin dynamics as m → 0 + . The computations of the auto-correlation functions are made under the assumption that both kinds of dynamics are wide-sense stationary. For any m > 0 this property holds for the Langevin dynamics, and in the limit m → 0 + it also holds for the heat bath dynamics.    In addition to our above observations, we believe it would be of great interest to obtain numerical verification for the heat bath dynamics weak convergence rate O(m 1/2 ) in (3.32). But, most likely due to computational constraints, we are unable to achieve this currently since even at the quite computationally demanding level of generating M = 200000 heat bath dynamics sample paths, it seems that the sample error dominates errors pertaining to the parameter m.