The classical limit of quantum observables in the conservation laws of fluid dynamics

TL;DR

This paper extends Irving and Zwanzig's classical limit of quantum observables to matrix-valued potentials, enabling molecular dynamics simulations to derive constitutive relations in fluid dynamics, especially at high temperatures.

Contribution

It introduces a matrix formulation for the semi-classical limit of quantum observables in conservation laws, applicable to systems with multiple excited electron states.

Findings

Derived the classical limit of quantum observables for matrix-valued potentials.

Extended the classical limit to high-temperature regimes beyond electron eigenvalue gaps.

Provided a framework for molecular dynamics simulations of stress and heat flux in quantum systems.

Abstract

In the classical work by Irving and Zwanzig [Irving J.H. and Zwanzig R.W., J. Chem. Phys. 19 (1951), 1173-1180 ] it has been shown that quantum observables for macroscopic density, momentum and energy satisfy the conservation laws of fluid dynamics. This work derives the corresponding classical molecular dynamics limit by extending Irving and Zwanzig's result to matrix-valued potentials for a general quantum particle system. The matrix formulation provides the semi-classical limit of the quantum observables in the conservation laws, also in the case where the temperature is large compared to the electron eigenvalue gaps. The classical limit of the quantum observables in the conservation laws is useful in order to determine the constitutive relations for the stress tensor and the heat flux by molecular dynamics simulations. The main new steps to obtain the molecular dynamics limit is to: (i) approximate the dynamics of quantum observables accurately by classical dynamics, by diagonalizing the Hamiltonian using a non linear eigenvalue problem, (ii) define the local energy density by partitioning a general potential, applying perturbation analysis of the electron eigenvalue problem, (iii) determine the molecular dynamics stress tensor and heat flux in the case of several excited electron states, and (iv) construct the initial particle phase-space density as a local grand canonical quantum ensemble determined by the initial conservation variables.

Equations

$$\Phi (x^1,{\sigma }^1,x^2,{\sigma }^2,\dots ,x^N,{\sigma }^N,{\bar{x}}^1,{\bar{\sigma }}^1,{\bar{x}}^2,{\bar{\sigma }}^2,\dots ,{\bar{x}}^J,{\bar{\sigma }}^J,t)\in \mathbb{C},$$

$$\begin{array}{rl}& \Phi (\dots ,{\bar{x}}^i,{\bar{\sigma }}^i,\dots ,{\bar{x}}^j,{\bar{\sigma }}^j,\dots ,t)=-\Phi (\dots ,{\bar{x}}^j,{\bar{\sigma }}^j,\dots ,{\bar{x}}^i,{\bar{\sigma }}^i,\dots ,t)\end{array}$$

$$\frac{|\Phi (x,\sigma ,\bar{x},\bar{\sigma },t){|}^2}{{\sum }_{\sigma \in \Sigma }{\sum }_{\bar{\sigma }\in \bar{\Sigma }}{\int }_{{\mathbb{R}}^{3(N+J)}}|\Phi (x,\sigma ,\bar{x},\bar{\sigma },t){|}^2\mathrm{d}x\mathrm{d}\bar{x}}$$

$$\mathrm{i}{\partial }_t\Phi (x,\bar{x},t)=\hat{\mathcal{H}}\Phi (x,\bar{x},t),$$

$$\hat{\mathcal{H}}=-\sum _{n=1}^N\frac{1}{2M_n}{\Delta }_{x^n}+\hat{W}(x,\bar{x})+V_b(x)+v_b(x,\bar{x}),$$

$$\begin{array}{rl}\hat{W}(x,\bar{x})& =-\frac{1}{2}{\Delta }_{\bar{x}}+\nu (x,\bar{x})\\ \nu (x,\bar{x})& =\sum _{j=1}^J\sum _{k<j}\frac{1}{|{\bar{x}}^k-{\bar{x}}^j|}+\sum _{m=1}^N\sum _{n<m}\frac{Z_n{Z}_m}{|x^n-x^m|}-\sum _{j=1}^J\sum _{n=1}^N\frac{Z_n}{|x^n-{\bar{x}}^j|}.\end{array}$$

$$\begin{split}&v_{b}(x,\bar{x}):=\sum_{j=1}^{J}h(|\bar{x}^{j}-x^{0}|),\\ &h(y)\to\infty\quad\mbox{ as $y\to\infty$},\\ &h(y)=0\quad\mbox{ for $|y|\leq C$,}\\ &V_{b}(x)\to\infty\quad\mbox{ as $|x|\to\infty$},\\ &V_{b}(x)=0\quad\mbox{ for $|x|\leq C$,}\end{split}$$

$$\hat{W}=\sum _{n=1}^N{\hat{W}}^n.$$

$${\hat{W}}^n(x,\bar{x})=-\frac{1}{2}\sum _{k\in {\mathcal{E}}_n}{\Delta }_{{\bar{x}}^k}+\frac{1}{2}\sum _{\ell \in {\mathcal{E}}_n}\sum _{k\ne \ell }\frac{1}{|{\bar{x}}^k-{\bar{x}}^{\ell }|}+\frac{1}{2}\sum _{m\ne n}\frac{Z_n{Z}_m}{|x^n-x^m|}-\sum _{k=1}^J\frac{Z_n}{|x^n-{\bar{x}}^k|},$$

$$(\hat{W}+v_b){\Psi }_k={\lambda }_k{\Psi }_k,$$

$$(x^1,\dots ,x^N,{\bar{x}}^1,\dots ,{\bar{x}}^J)\mapsto (Q{x}^1+\alpha ,\dots ,Q{x}^N+\alpha ,Q{\bar{x}}^1+\alpha ,\dots ,Q{\bar{x}}^J+\alpha ),$$

$$(\hat{W}+v_b)(x,\bar{x})=(\hat{W}+v_b)(Q{x}^1+\alpha ,\dots ,Q{x}^N+\alpha ,Q{\bar{x}}^1+\alpha ,\dots ,Q{\bar{x}}^J+\alpha ).$$

$$\lambda_{k}^{n}(x)=\left\langle{\Psi_{k}(x)},{\Big{(}\hat{W}^{n}+\frac{v_{b}}{N}\Big{)}(x)\Psi_{k}(x)}\right\rangle_{e}\,.$$

$$\left\langle{\Psi_{k}},{\big{(}\hat{W}+v_{b}-(\hat{W}^{n}+v_{b}/N)\big{)}\Psi_{k}\Big{)}}\right\rangle_{e}\,,$$

$${\lambda }_k=\sum _{n=1}^N{\lambda }_k^n.$$

$$\begin{array}{rl}V\psi & =\sum _{k=1}^d{\lambda }_k{\Psi }_k{⟨{\Psi }_k,\psi ⟩}_e,\\ V^n\psi & =\sum _{k=1}^d{\lambda }_k^n{\Psi }_k{⟨{\Psi }_k,\psi ⟩}_e,\end{array}$$

$$V(x)\psi_{k}(x)=\lambda_{k}(x)\psi_{k}(x)\,,\;\;\;\mbox{for all $x\in\mathbb{R}^{3N}$,}$$

$${\lambda }_1(x)<{\lambda }_2(x)<\dots <{\lambda }_d(x),$$

$$\displaystyle\lambda_{1}(x)+{V_{b}}(x)\rightarrow\infty\;\;\;\mbox{ as $|x|\rightarrow\infty$,}$$

$$\displaystyle\lambda_{d}(x)\to\infty\;\;\;\mbox{ as $d\to\infty$},$$

$$\displaystyle\{\lambda_{k}\}_{k=1}^{d}\mbox{ depend continuously on $\delta_{a}$ as $\delta_{a}\to 0+$}\,.$$

$$\begin{split}\partial_{t}\rho(y,t)+\sum_{\ell=1}^{3}\partial_{y_{\ell}}\big{(}\rho(y,t)u_{\ell}(y,t)\big{)}&=0\,,\\ \partial_{t}\big{(}\rho(y,t)u_{j}(y,t)\big{)}+\sum_{\ell=1}^{3}\partial_{y_{\ell}}\big{(}\rho(y,t)u_{j}(y,t)u_{\ell}(y,t)-\sigma_{\ell j}(y,t)\big{)}&=F_{j}(y,t)\,,\\ \partial_{t}E(y,t)+\sum_{\ell=1}^{3}\partial_{y_{\ell}}\big{(}E(y,t)u_{\ell}(y,t)+q_{\ell}(y,t)-\sum_{j=1}^{3}\sigma_{\ell j}(y,t)u_{j}(y,t)\big{)}&=P(y,t)\,,\\ \end{split}$$

$$H(x,p)=\sum _{n=1}^N\frac{1}{2M_n}|p^n{|}^2+\lambda (x)+V_b(x),$$

$$\begin{split}\dot{x}^{n}_{t}&=\frac{1}{M_{n}}p^{n}_{t}\,,\;\;\mbox{$n=1,\dots,N$}\\ \dot{p}_{t}&=-\nabla\lambda(x_{t})-\nabla V_{b}(x_{t})\,,\end{split}$$

$$\lambda (x)=:\tilde{\lambda }(r(x)).$$

$$\begin{split}\int_{\mathbb{R}^{3}}\eta(y)\mathrm{d}y&=1\,,\\ \eta(y)&\geq 0\,,\;\;\;\eta(y)=\eta(-y)\,,\;\;\quad\mbox{for all $y\in\mathbb{R}^{3}$,}\\ \eta(y)&=0\,,\quad\mbox{ for }|y|>\epsilon\,.\end{split}$$

$$\begin{array}{rl}\rho (y,t)& ={\int }_{{\mathbb{R}}^{6N}}\sum _n{M}_n\eta (y-x_t^n)f(x_0,p_0)\mathrm{d}{x}_0\mathrm{d}{p}_0,\end{array}$$

$$\begin{array}{rl}{\partial }_t\rho (y,t)& =-{\int }_{{\mathbb{R}}^{6N}}\sum _{n=1}^N{M}_n{\dot{x}}_t^n\cdot \nabla \eta (y-x_t^n)f(z_0)\mathrm{d}{z}_0\\ & =-{\int }_{{\mathbb{R}}^{6N}}\sum _{n=1}^N{p}_t^n\cdot \nabla \eta (y-x_t^n)f(z_0)\mathrm{d}{z}_0\end{array}$$

$$\rho (y,t)u(y,t):={\int }_{{\mathbb{R}}^{6N}}\sum _{n=1}^N\eta (y-x_t^n)p_t^{n}f(z_0)\mathrm{d}{z}_0$$

$$\partial_{t}\rho(y,t)+\sum_{k=1}^{3}\partial_{y_{k}}\big{(}\rho(y,t)u_{k}(y,t)\big{)}=0\,.$$

$$\begin{split}\partial_{t}\big{(}\rho(y,t)u(y,t)\big{)}&=-\int_{\mathbb{R}^{6N}}\sum_{n=1}^{N}M_{n}^{-1}p^{n}_{t}\cdot\nabla\eta(y-x^{n}_{t})p_{t}^{n}f(z_{0})\mathrm{d}z_{0}\\ &\quad-\int_{\mathbb{R}^{6N}}\sum_{n=1}^{N}\eta(y-x^{n}_{t})\nabla_{x^{n}}\lambda(x_{t})f(z_{0})\mathrm{d}z_{0}\\ &\quad-\int_{\mathbb{R}^{6N}}\sum_{n=1}^{N}\eta(y-x^{n}_{t})\nabla_{x^{n}}V_{b}(x_{t})f(z_{0})\mathrm{d}z_{0}\,.\\ \end{split}$$

$$\begin{split}\sum_{n=1}^{N}\eta(y-x^{n})\nabla_{x^{n}}\lambda&=\sum_{n}\sum_{j<k}\eta(y-x^{n})\partial_{r^{jk}}\tilde{\lambda}(r)\,\nabla_{x^{n}}(|x^{j}-x^{k}|)\\ &=\sum_{n<k}\big{(}\eta(y-x^{n})-\eta(y-x^{k})\big{)}\partial_{r^{nk}}\tilde{\lambda}(r)\,\nabla_{x^{n}}(|x^{n}-x^{k}|)\\ \end{split}$$

Full text

{AMS}

35L65, 35Q70, 82C10, 81Q20

1 The purpose of the work and the results

The macroscopic conservation laws for mass, momentum and energy form the basis of continuum fluid mechanics. These conservation laws are formulated in terms of the stress tensor and the heat flux. In order to form a closed system constitutive relations for the stress tensor and the heat flux are used. Such constitutive relations can be determined approximately from measurements or from molecular dynamics simulations. In both cases one seeks approximations of the stress tensor and the heat flux as functions of the density, momentum and energy and their derivatives. The molecular dynamics formulation requires derivation of the stress tensor and the heat flux as functions of the particle dynamics. The derivation of such functional relations is the focus of this work.

The stress tensor and the heat flux were first derived by Irving and Kirkwood, [5], from molecular dynamics systems based on interaction with scalar pair potentials and has later been modified by Noll, [14], and Hardy, [4]. These formulations have been used frequently to numerically determine the constitutive relations, cf. [3]. For instance, the works [18] and [19] include comparisons of different methods to numerically determine the stress tensor in molecular dynamics simulations.

Already in 1951 Irving and Zwanzig, [6], showed that quantum observables for the density, momentum and energy satisfy the conservation laws and derived observables for the stress tensor and the heat flux. Since it is only at the quantum level the particle interaction is determined from fundamental principles their result provides a solid foundation for the basic conservation laws in continuum mechanics. The property that the observables for the density, momentum and energy satisfy the conservation laws does not mean that a closed system of conservation laws is derived, since the derived stress tensor and the heat flux are not determined as constitutive functions of the macroscopic conservation variables. To form a closed system would include the additional step to determine constitutive functions of the conservation variables that approximate the data from molecular dynamics or measurements, which is not studied here.

Irving and Zwanzig used a quantum model with the Hamiltonian given by a sum of kinetic energy and scalar pair potential energy including all particles, i.e. both the nuclei and the electrons. The aim of this work is to extend the derivation by Irving and Zwanzig to a setting with a matrix-valued Hamiltonian consisting of a sum of the kinetic energy of the nuclei (times the identity matrix) and a matrix representing the electron kinetic energy, the electron-electron, electron-nuclei, and nuclei-nuclei interaction. The purpose of having a matrix for the electron part in the Hamiltonian is to replace the time evolution for the electrons by the Schrödinger electron eigenvalue problem. An advantage of including the electron part as a matrix-valued operator is that the classical limit, as the nuclei-electron mass ratio tends to infinity, has been derived rigorously, [16] and [7], and by knowing the classical limit the system can be simulated by ab initio molecular dynamics for nuclei with the potential generated by the electron eigenvalue problem. For instance, one may ask how the observables of the density, momentum, energy, stress tensor and heat flux are effected by the possibility of excited electron states and how these observables should be computed in molecular dynamics simulations. This question is answered in Theorem 5.5 by applying the classical molecular dynamics limit of quantum observables in [7]. The work [7] is for the setting of constant temperature in the canonical ensemble and shows, for example, how the potential is modified also when the difference of the excited and ground state electron eigenvalue is not large compared to the temperature.

The time evolution of the conserved quantum observables uses the ingenious observation by Irving and Zwanzig that, for an observable that is a polynomial of the degree at most two in the momentum coordinate, the commutator of the Hamiltonian operator and that quantum observable becomes equal to the Weyl quantization of the Poisson bracket. Combined with the observation that the observables for density, momentum and energy are polynomials of degree at most two in the momentum coordinate the quantum observables therefore satisfy the same conservation laws as in the derivation based on classical particle dynamics by Irving and Kirkwood. However, in the case of matrix-valued potentials the commutator of the Hamiltonian and the quantum observables for mass, momentum and energy does not reduce to a Poisson bracket since the matrix-valued symbols do not commute in general. In this work we show that for a certain diagonalization, based on a nonlinear eigenvalue problem, these commutators are reduced to a quantization of corresponding Poisson brackets.

To define the energy observable the works [5, 14, 6, 4] use that the potential energy can be split into a sum of potential energies related to each particle as defined by pair interactions. In the matrix-valued case considered here the splitting is required for the eigenvalues of the matrix potential which is not a sum of pair potential interactions. Our splitting is instead obtained by using perturbation theory for eigenvalues.

The pair potential property is also used in the works [5, 14, 6, 4] to reduce forcing terms to divergence of a stress term. Such reduction has been obtained in [1] for general potentials that are invariant with respect to translation and orthogonal transformations by changing to the coordinates depending on all pair distances. This change to the pair distance coordinates is also used here.

The compressible Euler equations have been derived from classical perturbed Newtonian particle dynamics using the relative entropy method in [15]. The classical Newtonian particle dynamics based on short range pair potential interactions is then weakly perturbed in two ways: to avoid unbounded velocities the kinetic energy is modified, for instance as relativistic, and to prove ergodicity with respect to Gibbs distributions the Hamiltonian dynamics is perturbed by a weak noise term that vanishes in the macroscopic hydrodynamic limit. A main accomplishment in [15] is to show that the density solving the Liouville equation that is initially close to a grand canonical Gibbs measure remains close to a grand canonical Gibbs measure at later time, so that the Gibbs measure determines the compressible Euler equations for all times, as long as the solution to the Euler equations remains smooth. The work [15] achieves the mathematically ambitious goal to derive a closed system of conservation laws from microscopic dynamics, which also requires additional assumptions and restricts to a setting with smooth classical solutions to the Euler equations. The Euler equations includes a pressure term that originates from microscopic particle forces. The relative entropy method has also been used to derive the compressible Euler equations, with a certain pressure term, in a scaling limit from a quantum system of fermions under an assumption of ergodicity of the quantum dynamics with respect to the Gibbs measure, see [13]. The objective in our work here is different from [15] and [13], in particular, we derive microscopic expressions for the stress tensor and the heat flux from a general quantum mechanical setting but we do not address the question of deriving a closed system of conservation laws from quantum mechanics.

We formulate the quantum mechanical model and the conservation laws in Section 2. In Sections 3 and 4 we review the derivations of the conservation laws from classical and quantum dynamics, respectively, following the works [5, 14, 4, 6], although reformulated in order to prepare for the new results in Section 5. These derivations are then used to obtain the matrix-valued extension of the quantum dynamics in Section 5. We derive the main result in Theorem 5.5, namely a classical molecular dynamics limit of the quantum observables in the conservation laws, under some assumptions on regularizations: the semiclassical analysis result requires $L^{2}$ bounded symbols, which is not satisfied in the canonical quantum formulation. The theorem therefore assumes that certain regularized $L^{2}$-bounded symbols depend continuously on the regularization parameters. In Section 6 we present an approach for determining an initial phase-space density that matches the initial conservation variables locally. Following [1] we discuss in Section 7 the non-uniqueness question for the stress tensor.

In conclusion, the main result in this work is to formulate the quantum conservation laws using matrix valued symbols and apply recent techniques from semiclassical analysis to determine the molecular dynamics stress tensor and heat flux in the case of several excited electron states. A new ingredient in the formulation is also to define the local energy density by partitioning a general potential, applying perturbation analysis of the canonical electron eigenvalue problem; previous work on molecular dynamics formulations of the conservation laws used empirical pair potentials, [5, 14, 4, 15]. As in [15] we construct the initial particle phase-space density as a local grand canonical quantum ensemble determined by the initial conservation variables.

We think that the presented work is the first ab initio result that determines the stress and the heat flux in a molecular dynamics setting with several excited states and at any temperature. Therefore we believe that this is a valuable first step for further study. For instance, it would be interesting to extend the result by determining conditions that imply the assumed continuous dependence on the regularization parameters.

2 Problem formulation

2.1 The quantum-mechanical model

We consider derivation of conservation laws from ab initio dynamics for which the starting point is the quantum mechanical model consisting of $N$ nuclei (heavy particles or slow degrees of freedom) and $J$ electrons (light particles or fast degrees of freedom). Each particle has a related position coordinate in $\mathbb{R}^{3}$ and a discrete spin coordinate. The spin coordinate $\bar{\sigma}_{i}$ for each electron takes the value in the set $\{-1/2,1/2\}$, and similarly the spin coordinate $\sigma_{i}$ for a nucleus can take values in a discrete set $\{-s,-s+1,\dots,s\}$, see [10, 2]. The quantum system at time $t$ is then described by a wave function

[TABLE]

with nuclei position coordinates $x=(x^{1},x^{2},\ldots,x^{N})\in\mathbb{R}^{3N}$ and electron position coordinate $\bar{x}=(\bar{x}^{1},\bar{x}^{2},\ldots,\bar{x}^{J})\in\mathbb{R}^{3J}$, and spin coordinates $\sigma=(\sigma^{1},\ldots,\sigma^{N})\in\Sigma\equiv\{-s,-s+1,\ldots,s\}^{N}$ and $\bar{\sigma}=(\bar{\sigma}^{1},\ldots,\bar{\sigma}^{J})\in\bar{\Sigma}\equiv\{-1/2,1/2\}^{J}$. The wave function is required to satisfy the Pauli exclusion principle which implies that it is anti-symmetric with respect to interchanging electron coordinates, namely

[TABLE]

and similarly identical fermion nuclei are also anti-symmetric while identical boson nuclei are symmetric with respect to its nucleon coordinates, see [10, 2]. We note that

[TABLE]

is the probability to find the quantum system in $(x,\sigma,\bar{x},\bar{\sigma})$ at time $t$. In the absence of magnetic fields, the wave function depends on the spin coordinates only parametrically since the Hamiltonian does not depend on the spin coordinates. To simplify the notation we therefore suppress the spin coordinates in the sequel although we include the dependence on the spin in Pauli exclusion principle for the electron wave functions. We assume the atomic units (a.u.) in which the mass and charge of the electron are equal to one, and the Planck constant is $\hbar=1$. We denote by $M_{n}$ the mass of individual nuclei.

The quantum mechanical evolution is described by the Schrödinger equation for a wave function $\Phi:\mathbb{R}^{3N}\times\mathbb{R}^{3J}\times[0,\infty)\rightarrow\mathbb{C}$ satisfying

[TABLE]

with the Hamiltonian operator

[TABLE]

where $V_{b}:\mathbb{R}^{3N}\to\mathbb{R}$ and $v_{b}:\mathbb{R}^{3N}\times\mathbb{R}^{3J}$ are given ad hoc external scalar smooth potentials and $\hat{W}(x,\bar{x})$ is the electronic operator formed by the electron kinetic energy, electron-electron repulsion, nuclei-nuclei repulsion and electron-nuclei attraction, see [2, 9, 10],

[TABLE]

Here $Z_{n}$ denotes the charge of the $n$th nucleus and $\bar{x}^{j}\in\mathbb{R}^{3}$ the coordinate of the electron $j$. The external potentials $V_{b}:\mathbb{R}^{3N}\to\mathbb{R}$ and $v_{b}:\mathbb{R}^{3N}\times\mathbb{R}^{3J}\to\mathbb{R}$ make the system confined and do not perturb the system far away from the boundary by the assumptions

[TABLE]

with a large constant $C$, where $x^{0}:=\sum_{n=1}^{N}x^{n}/N$ and $h:\mathbb{R}\to\mathbb{R}$ is a given smooth function.

We now show how the summation over all particles in the definition of the potential operator (2.3) can be rearranged into a sum over contributions from each nucleus

[TABLE]

We partition the electron index set $\mathcal{E}\equiv\{1,\dots,J\}$ into disjoints subsets $\mathcal{E}_{n}$, $n=1,\dots,N$ such that $\bigcup_{n=1}^{N}\mathcal{E}_{n}=\mathcal{E}$ and the number of elements in $\mathcal{E}_{n}$ is equal to $Z_{n}$, i.e., $\#\mathcal{E}_{n}=Z_{n}$. For the ease of exposition we assume charge neutrality of the electron-nuclei system. Using this partitioning we define

[TABLE]

where the terms including the sums in $\mathcal{E}_{n}$ correspond to the electron kinetic energy and electron-electron repulsion for electrons associated with the nucleus $n$.

2.2 The electronic operator

The electron eigenvalue problem takes the form

[TABLE]

with each function $\Psi_{k}$ in the anti-symmetric wave function subspace of $L^{2}$, based on $\bar{x}$ and the electron spin coordinates. Since $\hat{W}+v_{b}$ depends on the nuclei coordinate $x$ only parametrically, the eigenfunctions and eigenvalues will depend on $x$. We will see that the eigenvalues $\lambda_{k},\ k=1,\ldots$ will be the potentials for the nuclei dynamics which determine the molecular dynamics. The work [5, 6, 14, 4] use the explicit pair interactions of their potentials to define the energy related to each particle and to treat the conservation of total momentum. Our ab initio potentials $\lambda_{k}$ are not given by pair interactions. To derive the conservation laws from the quantum mechanical formulation (2.1), we will in particular use two properties of the electron eigenvalue problem (2.11), namely its invariance with respect to translations and rotations in $\mathbb{R}^{3}$ and its construction by pair interactions.

2.2.1 Rotational invariance

We note that $\hat{W}+v_{b}$ is invariant with respect to the affine transformations defined by

[TABLE]

where $Q\in O(3)$ is an orthogonal transformation of $\mathbb{R}^{3}$ and $\alpha\in\mathbb{R}^{3}$ a translation, i.e.,

[TABLE]

Therefore also the eigenvalues $\lambda_{k}$ and eigenfunctions are invariant with respect to such translations and orthogonal transformations.

2.2.2 Partition of the energy

We define the potential energy related to the nucleus $n$ as

[TABLE]

Here we denote by $\left\langle{u},{w}\right\rangle_{e}=\sum_{\bar{\sigma}\in\bar{\Sigma}}\int_{\mathbb{R}^{3J}}u^{*}(\bar{x},\bar{\sigma})w(\bar{x},\bar{\sigma})\,\mathrm{d}\bar{x}$ the scalar product on the Hilbert space $L^{2}(\mathbb{R}^{3J}\times\bar{\Sigma})$ of electron states. This definition can be motivated by standard perturbation analysis, see [8]: A small perturbation $w$ of the potential $\hat{W}+v_{b}$ yields to leading order that the $k$th eigenvalue $\lambda_{k}=\left\langle{\Psi_{k}},{(\hat{W}+v_{b})\Psi_{k}}\right\rangle_{e}$ is perturbed to \langle{\Psi_{k}},{\big{(}\hat{W}+v_{b}+w\big{)}\Psi_{k}}\rangle. The energy related to the particle $n$, for eigenvalue $k$, can therefore be viewed as the difference of $\lambda_{k}$ and the $k$th eigenvalue for the perturbed potential $\hat{W}+v_{b}-(\hat{W}^{n}+v_{b}/N)$, which corresponds to removing particle $n$ and its electrons and the fraction $1/N$ of the external potential $v_{b}$. This perturbed eigenvalue is, to the leading order,

[TABLE]

which leads to the definition in (2.8). Property (2.6) and (2.8) imply that for all $k$

[TABLE]

In the introductory Section 3 we will consider only one energy surface, often the ground state $\lambda_{1}(x)$, and therefore we drop the index $k$ there, and use the notation $\lambda:=\lambda_{k}$ for a given fixed $k$. In Section 5 we use all electron eigenvalues $\lambda_{k},\ k=1,\ldots d$, to establish molecular dynamics representations of the stress tensor and the heat flux.

2.2.3 Reduction to finite dimensional electron space and regularization of the electron operator

To obtain the classical limits in Section 5 we first regularize the electron operator $\hat{W}+v_{b}$ by replacing all Coulomb terms $1/|z|$ in $\nu$ by $1/\sqrt{|z|^{2}+\delta_{a}}$, with a positive constant $\delta_{a}>0$. We note that this regularization preserves the translation and rotation invariance. Then we introduce a finite dimensional approximation $V$ of $\hat{W}+v_{b}$ in the anti-symmetric electronic wave function space (including also the spin coordinates) by applying a projection onto the finite-dimensional subspace of $L^{2}(\mathbb{R}^{3J}\times\bar{\Sigma})$ spanned by the eigenfunctions $\{\Psi_{1},\dots,\Psi_{d}\}$ of the regularized Hamiltonian $\hat{W}+v_{b}$, with the corresponding eigenvalues $\{\lambda_{1},\ldots,\lambda_{d}\}$, as

[TABLE]

for any wave function $\psi\in L^{2}(\mathbb{R}^{3J}\times\bar{\Sigma})$. Since we have $v_{b}(x,\bar{x})\to\infty$ as $|\bar{x}|\to\infty$ for given $x$ and the potential part in $\hat{W}+v_{b}$ is locally integrable with respect to $\bar{x}$, the spectrum of $\hat{W}+v_{b}$ is discrete, see [12].

The motivation to introduce the reduction (2.10) to a smooth $d\times d$ matrix with distinct eigenvalues is that we then have a precise result of the approximation of quantum observables by molecular dynamics, in particular Theorem 5.4 proved in [7]. Thus $\hat{W}+v_{b}$ is approximated by a Hermitian matrix-valued operator $V:\mathbb{R}^{3N}\to\mathbb{C}^{d\times d}$ by (2.10). We assume that the eigenvalues $\lambda_{1}(x),\lambda_{2}(x),\ldots,\lambda_{d}(x)$ of this matrix-valued potential $V(x)$, i.e., solutions

[TABLE]

satisfy

[TABLE]

The first assumption is in order to have differentiable eigenvectors, the second condition implies that the system is confined, with respect to the nuclei, and the third condition is used to have a discrete spectrum of $\hat{W}+v_{b}$ and ensure that $V$ is a consistent approximation of $\hat{W}+v_{b}$.

In general, nuclei positions $x$ where electron eigenvalues coincide, e.g. points $x$ where $\lambda_{1}(x)=\lambda_{2}(x)$, form a co-dimension two set in $\mathbb{R}^{3N}$. Such points include so-called conical intersections which are difficult to handle in the classical limit and they are not included in Theorem 5.4. Section 5.4 describes how an extension of Theorem 5.4 to include conical intersection could be possible.

2.3 The conservation laws

The conservation laws or balance laws for mass, momentum and energy, based on the density $\rho:\mathbb{R}^{3}\times[0,\infty)\rightarrow[0,\infty)$, velocity $u:\mathbb{R}^{3}\times[0,\infty)\rightarrow\mathbb{R}^{3}$ and energy density $E:\mathbb{R}^{3}\times[0,\infty)\rightarrow\mathbb{R}$, take the form

[TABLE]

where $\sigma_{\ell j}:\mathbb{R}^{3}\times[0,\infty)\rightarrow\mathbb{R}$ is the $\ell j$-component of the $3\times 3$ stress tensor, $q_{\ell}:\mathbb{R}^{3}\times[0,\infty)\rightarrow\mathbb{R}$ is the $\ell$-th component of the heat flux, $F:\mathbb{R}^{3}\times[0,\infty)\to\mathbb{R}^{3}$ is an external force and $P:\mathbb{R}^{3}\times[0,\infty)\to\mathbb{R}$ is an external energy source. The purpose of this work is to derive these conservation/balance laws from microscopic dynamical systems. First we consider classical systems and then quantum systems.

3 The conservation laws derived from classical particle dynamics

In this section we consider a system of $N$ classical particles, where each particle has mass $M_{n}$, position coordinate $x^{n}$ and momentum coordinate $p^{n}$, $n=1,\ldots,N$. We use the notation $x=(x^{1},x^{2},x^{3},\ldots,x^{N})$ and $p=(p^{1},p^{2},p^{3},\ldots,p^{N})$. The position coordinates $x:[0,\infty)\rightarrow\mathbb{R}^{3N}$ and momentum coordinates $p:[0,\infty)\rightarrow\mathbb{R}^{3N}$ satisfy the classical equations of motion given by the Hamiltonian

[TABLE]

where $\lambda:\mathbb{R}^{3N}\to\mathbb{R}$ is a given interaction potential which we relate to the potential energy surface, i.e. an eigenvalue $\lambda_{k}(x)$ in (2.11), in the next section and $V_{b}:\mathbb{R}^{3}\to\mathbb{R}$ is the external potential. Thus the evolution of the system is given by the solution of Newtonian dynamics

[TABLE]

where $x^{n}_{t}\in\mathbb{R}^{3}$ and $p^{n}_{t}\in\mathbb{R}^{3}$ is the position and momentum, respectively, of the particle $n$ at time $t$. The given potential $\lambda:\mathbb{R}^{3N}\rightarrow\mathbb{R}$ is assumed to be invariant under the Euclidean group of transformations of $\mathbb{R}^{3}$, i.e. $\lambda(x^{1},\ldots,x^{N})=\lambda(Qx^{1}+\alpha,\ldots,Qx^{N}+\alpha)$ for any orthogonal $3\times 3$ matrix $Q$ and any translation $\alpha\in\mathbb{R}^{3}$. The following lemma will be used to represent the potential $\lambda$ as a function of pairwise distances between particles rather than the positions of each particle. It shows that there is a transformation in the Euclidean group, which consists of isometries, that maps all elements of one point set to another point set, provided the distances between points in both sets coincide.

Lemma 3.1**.**

If the two sets of points $\{x^{i}\}_{i=1}^{N}$ and $\{y^{i}\}_{i=1}^{N}$, where $x^{i},y^{i}\in\mathbb{R}^{3}$, satisfy $r^{ij}=|x^{i}-x^{j}|=|y^{i}-y^{j}|$ for $1\leq i,j\leq N$, then there exist an orthogonal matrix $Q\in\mathbb{R}^{3\times 3}$ and a translation vector $\alpha\in\mathbb{R}^{3}$ such that $x^{i}=Qy^{i}+\alpha$ for $1\leq i\leq N$.

Proof 3.2**.**

Let $\bar{x}^{i}:=x^{i}-x^{1},\bar{y}^{i}=y^{i}-y^{1}$, for $1\leq i\leq N$. If $\bar{x}^{i}=\bar{y}^{i}=0$, for all $1\leq i\leq N$, then clearly the claim in the theorem is true. If not, let $i_{1}$ be an index such that $\bar{x}^{i_{1}}\neq 0$ (which also implies that $\bar{y}^{i_{1}}\neq 0$). Let $Q_{1},Q_{2}\in\mathbb{R}^{3\times 3}$ be two orthogonal matrices such that $Q_{1}\bar{x}^{i_{1}}$ and $Q_{2}\bar{y}^{i_{1}}$ both lie on the first positive coordinate axis. Then clearly $Q_{1}\bar{x}^{i_{1}}=Q_{2}\bar{y}^{i_{1}}$.

Define $\bar{\bar{x}}^{i}:=Q_{1}\bar{x}^{i}$ and $\bar{\bar{y}}^{i}:=Q_{2}\bar{y}^{i}$ for all $1\leq i\leq N$. If all $\bar{\bar{x}}^{i}$ and $\bar{\bar{y}}^{i}$ lie on the first coordinate axis then $\bar{\bar{x}}^{i}=\bar{\bar{y}}^{i}$, for $1\leq i\leq N$, since every $\bar{\bar{x}}^{i}$ and $\bar{\bar{y}}^{i}$ have the same distance to $\bar{\bar{x}}^{1}$ in the origin, and $\bar{\bar{x}}^{i_{1}}$. Assume now that there exists an index $i_{2}$ such that $\bar{\bar{x}}^{i_{2}}$ does not lie on the first coordinate axis. Since $\bar{\bar{x}}^{i_{2}}$ and $\bar{\bar{y}}^{i_{2}}$ have the same distance to $\bar{\bar{x}}^{1}$ and $\bar{\bar{x}}^{i_{1}}$, also $\bar{\bar{y}}^{i_{2}}$ does not lie on the first coordinate axis. Let $Q_{3},Q_{4}\in\mathbb{R}^{3\times 3}$ be two orthogonal matrices that are rotations around the first coordinate axis such that $Q_{3}\bar{\bar{x}}^{i_{2}}$ and $Q_{4}\bar{\bar{y}}^{i_{2}}$ are both in the “positive $xy$-plane”, i.e. given as $(a,b,0)$ for $b>0$. This makes $Q_{3}\bar{\bar{x}}^{i_{2}}=Q_{4}\bar{\bar{y}}^{i_{2}}$ since the points are on the same distance to $\bar{\bar{x}}^{1}$ and $\bar{\bar{x}}^{i_{1}}$.

Define $\bar{\bar{\bar{x}}}^{i}:=Q_{3}\bar{\bar{x}}^{i}$ and $\bar{\bar{\bar{y}}}^{i}:=Q_{4}\bar{\bar{y}}^{i}$ for all $1\leq i\leq N$. Since the points $\bar{\bar{\bar{x}}}^{i}$ and $\bar{\bar{\bar{y}}}^{i}$ have the same distance to the points $\bar{\bar{\bar{x}}}^{1}$, $\bar{\bar{\bar{x}}}^{i_{1}}$, and $\bar{\bar{\bar{x}}}^{i_{2}}$, that all lie in the plane spanned by the first two coordinate directions, but not all of them on a straight line, we must either have that $\bar{\bar{\bar{x}}}^{i}=\bar{\bar{\bar{y}}}^{i}$ or $\bar{\bar{\bar{x}}}^{i}=Q\bar{\bar{\bar{y}}}^{i}$, for the reflection in the $xy$-plane Q=\bigl{(}\begin{smallmatrix}1&0&0\\ 0&1&0\\ 0&0&-1\end{smallmatrix}\bigr{)}. There cannot be two points $\bar{\bar{\bar{x}}}^{i}$ and $\bar{\bar{\bar{x}}}^{j}$ that do not lie in the $xy$-plane and satisfy $\bar{\bar{\bar{x}}}^{i}=\bar{\bar{\bar{y}}}^{i}$ and $\bar{\bar{\bar{x}}}^{j}=Q\bar{\bar{\bar{y}}}^{j}$, since then $\bar{\bar{\bar{x}}}^{j}$ and $\bar{\bar{\bar{y}}}^{j}$ would be on different distance from $\bar{\bar{\bar{x}}}^{i}=\bar{\bar{\bar{y}}}^{i}$.

Hence either $\bar{\bar{\bar{x}}}^{i}=\bar{\bar{\bar{y}}}^{i}$ for all $1\leq i\leq N$, or $\bar{\bar{\bar{x}}}^{i}=Q\bar{\bar{\bar{y}}}^{i}$ for all $1\leq i\leq N$. Since $\bar{\bar{\bar{x}}}^{i}$ are obtained from $x^{i}$ by the same set of translations and multiplications by orthogonal matrices for all $1\leq i\leq N$, and likewise for $\bar{\bar{\bar{y}}}^{i}$, the proof is complete.

To handle conservation of total momentum we will use Newtons third law for pair interactions and we follow the construction in [1] to determine pair interactions in a general potential that is invariant with respect to translations and orthogonal transformations in $\mathbb{R}^{3}$: knowing all $N(N-1)/2$ pair distances $r:=(r^{12},r^{13},\ldots,r^{N-1N}):=(|x^{1}-x^{2}|,|x^{1}-x^{3}|,\ldots,|x^{N-1}-x^{N}|)$ determines $x$ up to a translation and orthogonal transformation in $\mathbb{R}^{3}$ and since $\lambda(x)$ remains the same for such translations and orthogonal transformations the potential is determined by all pair distances, i.e.

[TABLE]

We will use the partial derivatives $\partial_{r^{jk}}\tilde{\lambda}(r^{12},r^{13},\ldots,r^{N-1N})$  . Not all $r\in\mathbb{R}^{N(N-1)/2}$ correspond to particle positions $x\in\mathbb{R}^{3N}$ and there are $N(N-1)/2$ partial derivatives $\partial_{r^{jk}}\tilde{\lambda}(r)$ while the gradient $\nabla\lambda(x)$ only has $3N$ components. Therefore the partial derivatives $\partial_{r^{jk}}\tilde{\lambda}(r)$ are not uniquely defined by $\nabla\lambda(x)$. Section 7 shows how to determine $\partial_{r^{jk}}\tilde{\lambda}(r)$.

To define the observables for density, momentum and energy and their dependence on the space coordinate $y\in\mathbb{R}^{3}$ we use a non-negative smooth mollifier $\eta:\mathbb{R}^{3}\rightarrow\mathbb{R}$, $\eta\in C^{\infty}(\mathbb{R}^{3})$, with a compact support, satisfying

[TABLE]

The macroscopic density $\rho:\mathbb{R}^{3}\times[0,\infty)\rightarrow\mathbb{R}$ is defined by the particle system as

[TABLE]

where $x_{t}^{n}$ is a function of the initial condition $(x_{0},p_{0})$, and $f:\mathbb{R}^{6N}\rightarrow[0,\infty)$ is a given initial particle distribution function normalized so that $\int_{\mathbb{R}^{6N}}f(x_{0},p_{0})\mathrm{d}x_{0}\mathrm{d}p_{0}=1$. Irving and Kirkwood, [5], use this definition with $\eta$ equal to a point mass and a general initial distribution $f$. Noll, [14], formulates the integration with respect to point masses in terms of the one-point and two-point density correlations functions instead and provides precise conditions for the validity of the derivation. Hardy, [4], uses the mollifier $\eta$ but not the integration over the initial particle distribution.

3.1 The conservation of mass

Let $z_{0}=(x_{0},p_{0})$ denote the phase-space coordinate in $\mathbb{R}^{6N}$ and $x\cdot y=\sum_{i=1}^{3}x^{i}y^{i}$ the Euclidean scalar product in $\mathbb{R}^{3}$. Differentiation of the density implies

[TABLE]

and by defining the velocity $u:\mathbb{R}^{3}\times[0,\infty)\rightarrow\mathbb{R}^{3}$ as

[TABLE]

we obtain the conservation law for the mass

[TABLE]

3.2 The conservation of momentum

Differentiation of the momentum yields

[TABLE]

In order to write the second term as a divergence term we follow Noll’s, [14], and Hardy’s method, [4], based on identifying gradients with respect to pair distances and converting the difference in $\eta$ at the the corresponding point to a gradient term: the combination of the pair distance derivative, using the definition of $\tilde{\lambda}$ in (3.18),

[TABLE]

and the difference at the corresponding points

[TABLE]

shows that

[TABLE]

We conclude that the following conservation law for the momentum holds

[TABLE]

To write the conservation law for momentum in the form (2.16) we follow the steps in [14]. Let

[TABLE]

where by (3.20)

[TABLE]

Definition (3.20) implies

[TABLE]

so that the first integral in (3.23) satisfies

[TABLE]

The conservation law for the momentum can therefore be formulated as

[TABLE]

where

[TABLE]

defines the stress tensor using $v^{n}=M_{n}^{-1}p^{n}-u(y,t)$, and

[TABLE]

defines the macroscopic external force.

3.3 The conservation of energy

In order to define the energy density we need to define the potential energy related to each particle. In the case of pair potentials this is straight forward by summing the pair potentials including the particle, as in [5, 14, 4]. Since we have a more general potential, which does not have to be a sum of pair potentials, this step requires a new construction: here we use the potential energy $\lambda^{n}$, related to the particle $n$, introduced in (2.8) and (2.9).

We define the energy density $E:\mathbb{R}^{3}\times[0,\infty)\rightarrow\mathbb{R}$ by

[TABLE]

and differentiate it to obtain

[TABLE]

We have

[TABLE]

so the right hand side in (3.26) becomes

[TABLE]

By using (3.22) we obtain the conservation law for the energy

[TABLE]

In order to write the energy conservation in the standard form (2.16) we use

[TABLE]

The first term in the right hand side becomes $E(y,t)u_{\ell}(y,t)$ under the integration in (3.27), and the last term is part of $\sum_{j}\sigma_{\ell j}(y,t)u_{j}(y,t)$. The term including the factor $v^{n}_{\ell}|u|^{2}$ vanishes upon integration over the initial distribution due to the definition (3.20) and the second term is included in the heat flux. The conservation of energy can therefore be written

[TABLE]

where the heat flux is defined as

[TABLE]

and the external energy source as

[TABLE]

4 The conservation laws derived from quantum mechanics

Irving and Zwanzig [6] derived the conservation laws when the particle system is modeled by the Schrödinger equation for a wave function $\Phi:\mathbb{R}^{3(N+J)}\times[0,\infty)\rightarrow\mathbb{C}$ satisfying

[TABLE]

with the Hamiltonian

[TABLE]

based on the nuclei and electron coordinates now written together as $\tilde{x}=(x,\bar{x})\in\mathbb{R}^{3(N+J)}$ with particle masses of nuclei and electrons denoted by $\tilde{M}_{n}=M_{n},\,n=1,\ldots N$, and $\tilde{M}_{n}=1,\,n=N+1,\ldots N+J$, as defined in (2.3). The wave function $\Phi$ is in an appropriate subset of $L^{2}(\mathbb{R}^{3(N+J)})$ taking anti-symmetry of electron coordinates into account. Irving and Zwanzig used the Wigner function to establish correspondence between classical and quantum observables.

We use instead the related Weyl quantization, which associates to a (Weyl) symbol $A:\mathbb{R}^{6(N+J)}\rightarrow\mathbb{C}$, i.e., a classical observable $A(x,\bar{x},p,\bar{p})$, an operator on $L^{2}(\mathbb{R}^{3(J+N)})$, in fact the Weyl quantization represents an isomorphism between $L^{2}(\mathbb{R}^{3(J+N)})$ and the space of Hilbert-Schmidt operators on $L^{2}(\mathbb{R}^{3(J+N)})$. The Weyl quantized operator $\hat{A}$ associated with the symbol $A$ in the Schwartz space acting on a function $\phi\in L^{2}(\mathbb{R}^{3(J+N)})$ is defined by

[TABLE]

and the definition is extended to more general symbols $A$ by the standard density arguments. For instance, quantization of the Hamiltonian symbol

[TABLE]

yields the operator $\hat{H}$.

We recall that the dependence on spin variables is only parametric and it does not enter the quantization procedure since we consider observables and Hamiltonians that do not depend on spin operators, for example modelling systems in absence of magnetic fields, or spin orbit coupling interactions etc. Therefore we omit spin variables $\sigma$, $\bar{\sigma}$ in the notation. However, we note that the scalar product of wave functions $\phi$, $\psi\in L^{2}(\mathbb{R}^{3(J+N)}\times\Sigma\times\bar{\Sigma})$ is $\langle\phi,\psi\rangle=\sum_{\sigma}\sum_{\bar{\sigma}}\int_{L^{2}(\mathbb{R}^{3(J+N)}\times\Sigma\times\bar{\Sigma})}\phi^{*}(x,\sigma,\bar{x},\bar{\sigma})\psi(x,\sigma,\bar{x},\bar{\sigma})\,\mathrm{d}x\,\mathrm{d}\bar{x}$.

The Schrödinger equation implies the evolution of the wave function $\Phi(\cdot,t)=e^{-\mathrm{i}t\hat{H}}\Phi(\cdot,0)$ and consequently an observable at time $t$ defined by $\langle\Phi(\cdot,t),\hat{A}\Phi(\cdot,t)\rangle$ satisfies

[TABLE]

By defining the evolution of observables as

[TABLE]

differentiation implies the Heisenberg-von Neumann equation

[TABLE]

where $[\hat{B},\hat{C}]=\hat{B}\hat{C}-\hat{C}\hat{B}$ is the commutator. We also obtain

[TABLE]

and $\hat{A}_{0}=\hat{A}$. Let $\hat{f}$ be the Weyl quantization of any initial classical density distribution $f:\mathbb{R}^{6(N+J)}\to\mathbb{R}$. Section 6 presents a precise definition of a density symbol related to the given initial data of the macroscopic density, momentum and energy. The Irving and Zwanzig quantum density observable is then defined by the $L^{2}(\mathbb{R}^{3(J+N)}\times\Sigma\times\bar{\Sigma})$ trace

[TABLE]

where $\{\Phi_{j}\}_{j=1}^{\infty}$ is a basis of the subspace of $L^{2}(\mathbb{R}^{3(N+J)}\times\Sigma\times\bar{\Sigma})$ based on the symmetry conditions of fermions and bosons and the density observable is the quantization of the density symbol

[TABLE]

The quantum momentum and energy observables are analogously defined as

[TABLE]

where $\nu=\sum_{n=1}^{N+J}\nu^{n}$ is a partition with the potential energy related to each particle, analogous to (2.6), now defined as

[TABLE]

Irving and Zwanzig observed the following crucial property.

Lemma 4.1**.**

For any twice differentiable scalar symbol that is a polynomial of degree two in the momentum variable

[TABLE]

there holds

[TABLE]

where $a_{0}(\tilde{x})\in\mathbb{R},a_{n}(\tilde{x})\in\mathbb{R}^{3},a_{n,m}(\tilde{x})\in\mathbb{R}^{{3}\times 3}$.

The result is also known in the mathematics literature, cf. [11, Remarks 2.6.9 and 2.7.6], and we include a proof of the lemma in Section 4.1, since it is important for this work. The lemma shows that for a symbol that is a polynomial of degree at most two in the momentum variable the quantum evolution of the observable reduces to the classical evolution based on the Poisson bracket. The quantum observables for density, momentum and energy are based precisely on symbols which are degree zero, one and two, respectively, in the momentum variable. Irving and Zwanzig therefore conclude that the quantum observables satisfy analogous conservation laws as those for classical particle dynamics, namely: differentiation of the density and using (4.33), Lemma 4.1, and the definition of $\hat{p}$ imply the conservation of mass

[TABLE]

where the first equality follows by the definition of the macroscopic density (4.34), the second by the Heisenberg-von Neumann dynamics of quantum observables (4.33), the third by Lemma 4.1, the forth by $\{H,\rho_{0}\}=-\sum_{n}\nabla\eta(y-x^{n})p^{n}=-{\rm div}_{y}p_{0}$ using (4.35) and the fifth by the definition of quantum time evolution (4.32).

Similarly differentiation of the momentum and energy establish the conservation laws for the momentum

[TABLE]

and the energy

[TABLE]

4.1 Proof of Lemma 4.1

We have the composition rule $\hat{B}\hat{C}=\hat{D}$ where

[TABLE]

see [20, Theorem 4.11]. Therefore

[TABLE]

and, letting $\tilde{M}$ denote the diagonal matrix with $\tilde{M}_{n}$ in the diagonal, we obtain

[TABLE]

Evaluation at the point $(\tilde{x}^{\prime},\tilde{p}^{\prime})=(\tilde{x},\tilde{p})$ yields

[TABLE]

The differentiation to the second order becomes

[TABLE]

so that $\mathrm{II}(\tilde{x},\tilde{x},\tilde{p},\tilde{p})=0$.

Since the symbol $A$ is a polynomial of degree two in the momentum variable we have

[TABLE]

which together with the Taylor expansion of the exponential proves the lemma.

4.2 Regularization of the observables

The use of the semiclassical analysis in the next section requires conditions on the observables, which are not satisfied for $\eta(y-x^{n})$. Therefore we replace in the proof of Theorem 5.5 the function $\eta(y-x^{n})$ by

[TABLE]

where $\zeta:\mathbb{R}\to[0,\infty)$ is a smooth cut-off function satisfying

[TABLE]

for a small positive constant $\delta_{p}>0$  .

5 The classical limit of the quantum conservation laws

This section first extends the formulation of quantum conservation laws to the case with matrix valued potentials in Subsection 5.1, and then in Subsection 5.2 this formulation is used to derive molecular dynamics limits of the stress tensor and the heat flux consistent with the quantum conservation laws as formulated in Theorem 5.5. Therefore we will redefine the notation for the Hamiltonian $H$ and rescale the Weyl quantization $A\mapsto\hat{A}$ and the composition operator ${\#}$.

5.1 Quantum conservation laws with matrix valued potentials

The aim here is to consider the quantum evolution where the electron part is considered matrix valued and the Weyl quantization is only in the nuclei part. For simplicity all nuclei masses are assumed to be equal, $M_{n}=M$, and we assume that the nuclei electron mass ratio $M\gg 1$ is large. The case with individual masses is treated by a change of variables in Section 5.3. To obtain the classical molecular dynamics limit, as $M\to\infty$, we restrict the electron operator to the finite dimensional $d\times d$ matrix $V$ in (2.11). We also change the time scale so that the nuclei dynamics has a limit, when $M\to\infty$, as follows: the wave function related to (4.30) can be written as $\Phi:\mathbb{R}^{3N}\times[0,\infty)\rightarrow\mathbb{C}^{d}$ and it solves the Schrödinger equation

[TABLE]

with the change of the time scale $\tau=M^{-1/2}t$. The classical limit obtained as $M\rightarrow\infty$ is well behaved in this time scale, in the sense that the nuclei move a distance of order one in time one, as we shall see in Theorem 5.4. The Hamiltonian is now defined as

[TABLE]

with the Hermitian $d\times d$ matrix $V$ defined in (2.10), the external potential $V_{b}$ satisfying (2.4) and ${\rm I}$ the $d\times d$ identity matrix. In this time scale the corresponding Heisenberg-von Neumann equation takes the form

[TABLE]

for the $d\times d$ matrix valued symbol $A_{\tau}(x,p)$. This form of matrix valued symbols are suitable for studying the $M\rightarrow\infty$ limit of the observables we have in the conservation laws, since the corresponding Weyl quantizations are based on highly oscillatory Fourier integral operators that only use the nuclei coordinates. We also assume that the eigenvalues $\lambda_{1}(x),\lambda_{2}(x),\ldots,\lambda_{d}(x)$ of $V(x)$, defined by (2.11), satisfy (2.12)-(2.15). The new Weyl quantization takes the form

[TABLE]

which differs from (4.31) by the scaling $M^{1/2}$. Although this form of matrix valued symbols $H$ and $A$ is useful to obtain the classical limit as $M\to\infty$, matrix valued symbols introduce a complication: the important property in Lemma 4.1 that the commutator with respect to the conservation observables reduces to the quantization of the Poisson bracket does not hold for matrix valued symbols, since these matrices do not commute in general, unless the symbols are diagonal. A main tool to determine the classical limit is therefore to diagonalize $H$ and the observables $A_{\tau}$ in the conservation laws, based on the composition operator $\#$, as follows.

The symbol $C$ for the product of two Weyl operators $\hat{A}\hat{B}=\hat{C}$ is determined by

[TABLE]

see [20, Theorem 4.11] which now includes the scaling $M^{1/2}$ as compared to (4.40). Let $\Psi:\mathbb{R}^{3N}\rightarrow\mathbb{C}^{d^{2}}$ satisfy that $\Psi(x)$ is a unitary matrix for every $x$ with the Hermitian transpose $\Psi^{*}(x)$, and define ${\bar{A}}:\mathbb{R}^{3N}\times[0,\infty)\rightarrow\mathbb{C}^{d^{2}}$ by

[TABLE]

so that

[TABLE]

Then

[TABLE]

and consequently

[TABLE]

The composition rule (5.42) yields $\Psi^{*}\hat{H}\Psi=(\Psi^{*}\#H\#\Psi)^{\widehat{}}$. The next step is to determine $\Psi$ so that

[TABLE]

is diagonal or approximate diagonal. Having $\bar{\bar{H}}$ diagonal implies that $\hat{\bar{\bar{H}}}$ is diagonal and then $\hat{\bar{A}}$ remains diagonal if it is initially diagonal, since then

[TABLE]

The composition rule (5.42) with

[TABLE]

implies that

[TABLE]

as verified in [[7], Lemma 3.1]. Therefore the aim is to choose the unitary matrix $\Psi$ so that it is an approximate solution to the nonlinear eigenvalue problem

[TABLE]

where $\bar{\bar{\Lambda}}$ is diagonal. A solution, $\Psi$, to this nonlinear eigenvalue problem is an $\mathcal{O}(M^{-1})$ perturbation of the eigenvectors to $V(x)$ provided the eigenvalues do not cross and $M$ is sufficiently large. The work [7, (3.18)] shows that (5.43) has an approximate solution $\Psi$ that satisfies

[TABLE]

with $\bar{\Lambda}$ diagonal, based on the following iteration. Let $\Psi_{0}=[\psi_{1}\ \psi_{2}\ \ldots\ \psi_{d}]$ be the matrix of eigenvectors to $V$. The approximate eigenvectors $\Psi$ are the normalized eigenvectors of

[TABLE]

and $\bar{\Lambda}$ are the corresponding eigenvalues. Since regular perturbation theory shows $\|\Psi-\Psi_{0}\|_{C^{1}(\mathbb{R}^{3N})}=\mathcal{O}(M^{-1})$ we obtain

[TABLE]

where the remainder $d\times d$ matrix $r_{0}(x)$ is small

[TABLE]

with the diagonal matrix

[TABLE]

We also need a partition of the eigenvalues $\bar{\Lambda}=\sum_{n=1}^{N}\bar{\Lambda}^{n}$ related to the potential energy for each particle similar to (2.8), now including also the small nonlinear part. In fact, also the nonlinear part has a natural composition into particle contributions, now based on the sensitivity of the eigenvectors with respect to position $x^{n}$. With $V^{n}$ defined by (2.10) we have

[TABLE]

and define for $n=1,\ldots,N$

[TABLE]

which implies

[TABLE]

We will use the observables defining density, momentum and energy as follows. Let

[TABLE]

and as before its time evolution is determined by the Heisenberg-von Neumann equation

[TABLE]

with the solution

[TABLE]

which shows that the time evolution also can be written as

[TABLE]

The momentum and energy density symbols are defined as

[TABLE]

The next step is to derive the quantum conservation/balance laws by studying the evalution of the observables for density, momentum and energy. Let $\hat{f}={\Psi}\hat{\bar{f}}{\Psi}^{*}$ be the Weyl quantization of a given density symbol $\bar{f}(x,p)\in\mathbb{R}^{d\times d}$ as described precisely in Section 6.

Lemma 5.1**.**

Assume that the eigenvalues $\lambda_{k},\ k=1,\ldots,d$ of $V$ are distinct and there are positive constants $C$ and $c$ such that

[TABLE]

then

[TABLE]

as $M\to\infty$.

As compared to (4.37), the quantum continuity equation (5.49) includes a $\mathcal{O}(M^{-3/2})$ remainder term, due to non perfect diagonalization in (5.44).

Proof 5.2**.**

We have

[TABLE]

where $\|r_{0}\|_{L^{\infty}(\mathbb{R}^{3N})}=\mathcal{O}(M^{-2})$ by (5.45).

The diagonal form of $\bar{H}$ combined with the property that the symbols $\bar{\rho}_{0}$, $\bar{p}_{0}$, $\bar{E}_{0}$ and $\bar{H}$ are polynomials of degree at most two as functions of the momentum coordinate imply by Lemma 4.1, rewritten in the new scaling, the reduction of the corresponding quantum commutators to classical Poisson brackets. In the new scaling Lemma 4.1 takes the form: assume that $\bar{H}$ and $\bar{A}$ are diagonal matrices where each component of $\bar{A}$ is a polynomial of degree at most two in $p$, as in (4.36), then ${\rm i}M^{1/2}[\hat{\bar{H}},\hat{\bar{A}}]=\{\bar{H},\bar{A}\}^{\widehat{}}$  . Therefore, we have as in (4.37)

[TABLE]

and by (5.48)

[TABLE]

By taking the trace in (5.50) the conservation law for the mass becomes

[TABLE]

The next step is to estimate the remainder term including $r_{0}$. Cauchy’s inequality in the Hilbert-Schmidt inner product $\mathrm{Tr}\,(\hat{A}^{*}\hat{B})$ implies that the remainder term has the estimate

[TABLE]

The Weyl quantization satisfies

[TABLE]

where $A(z)B(z)$ is the matrix product of the two $d\times d$ matrices $A(z)$ and $B(z)$, with the second trace acting on matrices, see [16] and [7, Lemma 3.1]. This isometry between Hilbert-Schmidt operators and $L^{2}(\mathbb{R}^{6N},\mathbb{C}^{d\times d})$ functions also extends the Weyl quantization from symbols in the Schwartz class to $L^{2}(\mathbb{R}^{6N},\mathbb{C}^{d\times d})$, see [16]. Lemma 5.3 implies

[TABLE]

and we obtain by (5.45) and (5.52)

[TABLE]

*which by (5.51) and (5.52) proves the lemma. *

The work [7, Lemma 3.11] proves

Lemma 5.3**.**

Assume $D:\mathbb{R}^{6N}\rightarrow\mathbb{C}^{d\times d}$ belong to $L^{2}(\mathbb{R}^{6N})$ and $A:\mathbb{R}^{3N}\rightarrow\mathbb{C}^{d\times d}$ depends only on the $x$-coordinate (or only on the $p$-coordinate) and is bounded in $L^{\infty}(\mathbb{R}^{3N})$ then

[TABLE]

The conservation of momentum and energy are also based on the reduction from commutators to Poisson brackets, in Lemma 4.1, as follows

[TABLE]

As for the conservation of mass in Lemma 5.1 we obtain the conservation/balance law for the momentum

[TABLE]

and the conservation/balance law for the energy

[TABLE]

provided $\bar{p}_{0},\bar{S}_{0},\bar{F}_{0},\bar{E}_{0},\bar{Q}_{0},\bar{P}_{0}$ are all bounded in $L^{2}(\mathbb{R}^{6N})$.

In Section 6 we will motivate an initial density $\hat{f}=\hat{\Psi}\hat{\bar{f}}\hat{\Psi}^{*}$ as a local grand canonical Gibbs density, where $\bar{f}$ is diagonal, with the local temperature and chemical potential determined by the macroscopic density and energy. That is, $\hat{f}$ is diagonalized by the same transformation as $\hat{H}$. The traces in the quantum conservation laws (5.49) (5.55) and (5.56) can then be written as

[TABLE]

where $\bar{A}$ is diagonal and equal to $\bar{\rho},\bar{p},\bar{E},\bar{S}$ and $\bar{Q}$, respectively. The next section presents the classical limit of these traces.

5.2 The molecular dynamics limit of the quantum conservation laws

The work [7] proves in Theorem 3.7 the following which provides the classical limit of the quantum observables in the conservation laws. A related result, with different assumptions, is in [16]. The proof is based on Weyl’s law, see [16, 7] namely that quantum observables have the classical representation

[TABLE]

for any $\bar{A}_{jj}(0,\cdot)\in L^{2}(\mathbb{R}^{6N})$ and $\bar{f}_{jj}\in L^{2}(\mathbb{R}^{6N})\cap L^{1}(\mathbb{R}^{6N})$.

Theorem 5.4**.**

Assume that $V$ satisfies the coercivity condition (2.13), the $d\times d$ matrices $\bar{A}_{0}$ and $\bar{f}$ are diagonal, the $d\times d$ matrix-valued Hamiltonian $H$ has distinct eigenvalues, and that there is a constant $C$ such that

[TABLE]

hold, then there is a constant $c$, depending on $C$, such that the canonical ensemble average satisfies the error estimate

[TABLE]

as the mass ratio $M\rightarrow\infty$, where $z^{j}_{\tau}=(x_{\tau},p_{\tau})$ is the solution to the Hamiltonian system

[TABLE]

based on the Hamiltonian $\bar{H}_{jj}(z)=|p|^{2}/2+\bar{\lambda}_{j}(x)+V_{b}(x)$, with initial data $(x_{0},p_{0})=z_{0}$.

We note that the classical limit can be written

[TABLE]

where the probability, $q_{j}^{*}$, to be in electron state $j$ is

[TABLE]

To apply Theorem 5.4 to the quantum observables (5.49), (5.55) and (5.56) for macroscopic density, momentum and energy the momentum variable $\eta$ needs to be regularized, since e.g. the momentum symbol $\bar{p}_{0}=\sum_{n}\eta(y-x^{n})p^{n}{\rm I}$ is not in $L^{2}(\mathbb{R}^{6N})$. Therefore we regularize all symbols by replacing $\eta(y-x^{n})$ by $\eta_{\delta_{p}}:=\eta(y-x^{n})\zeta(|p^{n}|^{2})$, given in (4.41), and denote the quantum observables using $\eta_{\delta_{p}}$ instead of $\eta$ in (5.49), (5.55) and (5.56) as $\mathrm{Tr}\,(\hat{\bar{\rho}}_{\tau}\hat{\bar{f}})_{\delta_{p}}$ (replacing $\mathrm{Tr}\,(\hat{\bar{\rho}}_{\tau}\hat{\bar{f}})$ and similarly for the other observables). {assumption} Assume that, for any regularization $\delta_{a}$ and dimension $d$, the quantum observables

[TABLE]

and their derivatives with respect to $\tau$ and $y$ have limits as $\delta_{p}\to 0+$, with the limits based on $\eta=\eta_{0}$.

Theorem 5.4 and the assumed continuous dependence on the regularization parameters can be used to show a consistency result, namely that as the nuclei electron mass ratio $M$ tends to infinity in the quantum conservation laws (5.49), (5.55) and (5.56), using the splitting $p_{\tau}^{n}=v_{\tau}^{n}+u(y,\tau)$ for the fluxes $\bar{S}$ and $\bar{Q}$, we obtain the following limit in the form of a macroscopic conservation/balance law based on a certain stress tensor and heat flux defined by molecular dynamics including several electron eigenvalues. We note that the diagonal terms in the flux terms $\bar{S}$ and $\bar{Q}$ are the same as in the classical dynamics (3.23) and (3.27).

Theorem 5.5**.**

Assume that the approximate electron operator $V$ satisfies(2.12)-(2.15) and Assumption 5.2 and the assumptions in Theorem 5.4 hold, with $\bar{A}$ diagonal and equal to $\bar{\rho}_{0},\bar{p}_{0},\bar{E}_{0},\bar{S}_{0}$ and $\bar{Q}_{0}$, then as the nuclei-electron mass ratio $M\to\infty$ the quantum conservation laws (5.49), (5.55) and (5.56) have a classical molecular dynamics limit that satisfies

[TABLE]

where

[TABLE]

and $\bar{\sigma}(z^{j};y,{\tau})$ and $\bar{q}(z^{j};y,{\tau})$ are defined in (3.24) and (3.28), respectively, now using $\lambda=\bar{\lambda}_{j}$, and the matrix valued symbols

[TABLE]

Proof 5.6**.**

The quantum observables (5.60), based on the approximate electron operator $V$ with finite $d$ and positive parameter $\delta_{a}$, satisfy the assumptions in Theorem 5.4 with $\eta$ regularized as $\eta_{\delta_{p}}$. These regularized quantum observables therefore have a classical limit as $M\to\infty$: given small positive $(\delta_{a},\delta_{p})$ and large finite $d,M$, the leading order terms in the classical approximation is by Theorem 5.4 arbitrary close to the terms in (5.62) while the error term $c/M$ can be made sufficiently small using sufficiently large $M$.

To verify that the observables satisfy the conservation laws, we need to take the limit $\delta_{p}\to 0+$ in both the quantum observables (5.60) and the leading order classical term (5.62), since the quantum conservation laws (5.49), (5.55) and (5.56) are given with $\delta_{p}=0$. As $\delta_{p}\to 0+$, the regularized observables in (5.60), based on the positive $\delta_{a}$, converge by Assumption 5.2 to the observables that satisfy the conservation laws (5.49), (5.55) and (5.56). Combined with the continuous dependence on $\delta_{p}\to 0+$ in the leading order terms of the molecular dynamics approximation (5.62) we obtain the consistency results (5.61)-(5.63).

The result of Theorem 5.5 still depends on the regularization parameters $\delta_{a}$ and $d$, while the ab initio model corresponds to $\delta_{a}=0$ and $d=\infty$. The molecular dynamics observables in Theorem 5.5 are consistent with the ab initio quantum model in the sense that if the molecular dynamics observables in Theorem 5.5 and the $M\to\infty$ limits of the quantum observables in (5.60) have limits as $\delta_{a}\to 0+$ and $d\to\infty$, with the limits based on $\delta_{a}=0$ and $d=\infty$, then the corresponding molecular dynamics and quantum limits are equal.

Theorem 5.5 proves in particular that the observables are determined by a weighted average with the probability $q^{*}_{j}$ to be in state $j$. Using this probability model, the conservation laws (5.61) are consistent with the derivation of the conservation laws from the classical dynamics in Section 3.

5.3 General nuclei masses

The general case of individual nuclei masses and a diagonal mass matrix $M$ can be treated by rescaling the nuclei position coordinates as $M_{1}^{1/2}x^{\prime}=M^{1/2}x$, which transforms the Hamiltonian into

[TABLE]

In these transformed coordinates the classical limit is by Theorem 5.4 based on Lemma 4.1 applied to \{\bar{H}(x^{\prime},p^{\prime}),\bar{A}\big{(}x(x^{\prime}),p(p^{\prime})\big{)}\} with $\bar{A}$ given by the conservation variables (5.63), which take the form \{\bar{H}_{jj}(x^{\prime},p^{\prime}),\bar{A}_{jj}\big{(}x(x^{\prime}),p(p^{\prime})\big{)}\}=\{|p|^{2}/2+M_{1}M^{-1}\bar{\lambda}_{j}(x),\bar{A}_{jj}(x,p)\} in the original variables. Therefore the transformed Hamiltonian system

[TABLE]

provides the untransformed system

[TABLE]

which is equivalent to the standard form

[TABLE]

with individual masses in the diagonal mass matrix $M$.

5.4 Coinciding eigenvalues

The assumption on distinct eigenvalues in Theorem 5.4 is used in the proof to obtain differentiable eigenvalues and eigenvectors by regular perturbation theory in (5.44). Here we indicate how to relax the assumption on distinct eigenvalues for a special example.

The first step is to perturb $V$ to obtain distinct eigenvalues. We have $V=\Psi\Lambda\Psi^{*}$ where $\Psi$ is the matrix with the eigenvectors as columns and $\Lambda$ is the diagonal matrix of eigenvectors $(\lambda_{1},\ldots,\lambda_{d})$. Assume, for example, there are points $x^{\prime}$ such that eigenvalues $\lambda_{1}(x^{\prime})=\lambda_{2}(x^{\prime})$ coincide and the other pairs of eigenvalues do not coincide. Since we have $\lambda_{1}\leq\lambda_{2}<\lambda_{3}$ we can regularize $\breve{\lambda}_{i}(x):=\int_{\mathbb{R}^{3N}}\lambda_{i}(y)\frac{e^{-|x-y|^{2}/2\delta_{e}}}{(2\pi\delta_{e})^{d/2}}{\rm d}y$ for $i=1,2$ so that for small $\delta_{e}>0$ we obtain distinct eigenvalues $\breve{\lambda}_{1}<\breve{\lambda}_{2}<\lambda_{3}$, which are smooth since $\lambda_{i}$ is Lipschitz continuous. Replace $V$ by the new matrix $\breve{V}:=\Psi\breve{\Lambda}\Psi^{*}$ where $\breve{\Lambda}$ is the diagonal matrix with $\lambda_{i}$ replaced by $\breve{\lambda}_{i}$ for $i=1,2$. Since the eigenvalues of $\breve{V}$ are distinct and smooth.

It remains to study the eigenvectors. The eigenvectors corresponding to $\lambda_{1}$ and $\lambda_{2}$ are typically not continuous at conical intersections $x$, where $\lambda_{1}(x)=\lambda_{2}(x)$, see e.g. [17, Section 12.2.3]. The proof of Theorem 5.4 uses the regularity of the eigenvectors of $V$ only to obtain (5.44). The approximate nonlinear diagonalization becomes surprisingly simple near a conical intersection for the case of a real valued $2\times 2$ matrix V=\left[\begin{array}[]{cc}v_{11}&v_{12}\\ v_{12}&v_{22}\end{array}\right]. Its eigenvalues are $\lambda_{\pm}=\frac{v_{11}+v_{22}}{2}\pm r$ where $r:=\sqrt{(\frac{v_{11}-v_{22}}{2})^{2}+v_{12}^{2}}$ with the corresponding eigenvectors

[TABLE]

Here $r(\cos\alpha,\sin\alpha)=(\frac{v_{11}-v_{22}}{2},v_{12})$, so

[TABLE]

and when $r(x)=0$ there is a conical intersection with $\lambda_{+}(x)=\lambda_{-}(x)$. Let $\tilde{v}:=\frac{v_{11}-v_{22}}{2}$, then $\tan\alpha=v_{12}/\tilde{v}$ so $\alpha=\tan^{-1}v_{12}/\tilde{v}+n\pi$ and $\nabla\alpha(x)=\frac{\tilde{v}\nabla v_{12}-v_{12}\nabla\tilde{v}}{\tilde{v}^{2}+v_{12}^{2}}$, which implies

[TABLE]

and consequently

[TABLE]

in fact solves the nonlinear eigenvalue problem (5.43) exactly. We note that the nonlinear eigenvalue $\bar{\lambda}_{k}$ has the additional term $|\nabla\alpha|^{2}/(16M)$ which is large if $r$ is small, i.e. in the potential landscape of $\bar{\lambda}_{k}$ there is a mountain around the conical intersection that becomes lower as M gets larger. These observations could be the first step to extend Theorem 5.4 to include conical intersections.

6 The initial particle density

Assume we know the initial data \big{(}\rho(\cdot,0),\rho u(\cdot,0),E(\cdot,0)\big{)} for the macroscopic conservation laws, although in practise this data can be hard to determine, e.g. for the flow in a river. Current molecular dynamics simulations can only use a small fraction of the number of particles in a real system. Therefore we need an initial particle density that is related to a larger ensemble. We seek a density that has the property that the marginal distribution of a subsystem weakly coupled to a larger heat bath system is the same as the whole system. Under certain assumptions stated in [7] the classical Gibbs density is the only density with this property. Given the local values of the macroscopic conserved variables the goal here is therefore to determine a local grand canonical Gibbs density

[TABLE]

where

[TABLE]

based on local values of the temperature $T(y)$ and the chemical potential $\mu(y)$, such that

[TABLE]

Weyl’s law given by the quantum-classical representation (5.57) combined with (6.65) show that the equation

[TABLE]

defines the initial velocity $u(y,0)$. It remains to verify if varying $T(y)$ and $\mu(y)$ yield large enough sets to match the initial data for $\rho$ and $E$.

If $|\mu(y)|\gg 1$ we roughly get $\mathcal{H}\simeq-\sum_{n=1}^{N}\eta(y-x^{n})\mu(y)$. Laplace principle implies that as $\mu(y)\to-\infty$ the grand canonical density will sample the minimum of the observable $\sum_{n=1}^{N}\eta(y-x^{n})$ and as $\mu(y)\to\infty$ the microscopic particle density will sample the maximum of $\sum_{n=1}^{N}\eta(y-x^{n})$. Therefore varying $\mu(y)$ from $-\infty$ to $\infty$ will change the local density

[TABLE]

from nearly vacuum to arbitrary high macroscopic density.

The temperature is related to the microscopic kinetic energy and we have

[TABLE]

We see that the factor $\eta(y-x^{n})\frac{|p^{n}|^{2}}{2M_{n}}e^{-\eta(y-x^{n})|p^{n}|^{2}/(T(y)2M_{n})}$ upon integration with respect to $p^{n}$ will be proportional to $T(y)$. The other terms in $E(y,0)$ will have upper and lower bounds uniform in $T$. Therefore by varying the temperature a large open set of macroscopic energies can be attained.

Finally, we note that minimizing the entropy $\sum_{j=1}^{d}\int_{\mathbb{R}^{6N}}\bar{f}_{jj}(z)\log\bar{f}_{jj}(z)\mathrm{d}z$ under the constraints that the value of the macroscopic density is $\rho(y,0)$ and the macroscopic energy is $E(y,0)$ yields the probability density

[TABLE]

for $c=1/\sum_{j=1}^{d}\int_{\mathbb{R}^{6N}}e^{-\mathcal{H}(z,j;y)/T(y)}\mathrm{d}z$ with the Lagrange multipliers $\mu_{1}(y)=1/T(y)$ and $\mu_{0}(y)=\mu(y)/T(y)$. The Gibbs density (6.64) is therefore consistent with this constrained minimization.

7 The partial derivatives $\partial_{r^{jk}}\tilde{\lambda}(\tilde{x})$

The function $\tilde{\lambda}$ is defined by (3.18) on the set $M\subset\mathbb{R}^{N(N-1)/2}$ consisting of all vectors $r:=(r^{12},r^{13},\ldots,r^{N-1N})$ such that there exist particle positions $x^{i}\in\mathbb{R}^{3}$, $i=1,\ldots,N$ and $r^{ij}=\lvert x^{i}-x^{j}\rvert$ for $1\leq i<j\leq N$. The goal here is to define the gradient $\nabla_{r}\tilde{\lambda}$ such that the chain rule \nabla\lambda(x)=\big{(}\partial r/\partial x\big{)}^{T}\nabla_{r}\tilde{\lambda}(r) is valid on $M$. We therefore need to solve the linear equation $Av=b$ for $v=\nabla_{r}\tilde{\lambda}$, where we have used the notation A=\big{(}\partial r/\partial x\big{)}^{T} and $b=\nabla\lambda$. Since this is an underdetermined linear system we choose the solution $v$ that minimizes the $\ell^{2}$ norm. The stationary point to the Lagrangian $L(v,y)=\lvert v\rvert^{2}/2+y\cdot(Av-b)$, where $y$ is the Lagrange multiplier, has the solution

[TABLE]

Since generically the matrix $A$ has full rank, equation (7.66) contains a computable expression for $v=\nabla_{r}\tilde{\lambda}$, which can be used in the expression for the stress tensor in (3.24).