What good are hard-sphere event-driven molecular dynamics simulations in the face of chaos?

Question

Simple hard-sphere dynamical systems can exhibit chaotic dynamics. Due to finite-precision arithmetic when implemented on a computer, the presence of chaos implies that for a given set of initial data, simulations of such systems on different machines or run with different arithmetic precisions can yield wildly different dynamics for sufficiently long times.

Nonetheless, for e.g. systems of many hard spheres in a box with reflecting boundaries, can it be mathematically proven that qualitative and/or statistical aspects of the system's finite-precision-simulated behavior will approximate its infinite-precision behavior?

I'm particularly curious about simple event-driven simulations of such systems such as that devised by Alder and Wainwright in 1959¹. I would imagine that the answer to the question is yes, because otherwise what's the use of these methods which, as far as I can tell, are ubiquitous? I can't seem to find a good, reasonably specific discussion of this anywhere. References appreciated.

---

1. Alder, B.J. and Wainwright, T.E., 1959. Studies in molecular dynamics. I. General method. The Journal of Chemical Physics, 31(2), pp.459–466.

Answer 1

Yes, it is possible to show that the statistical behavior of the approximate system will reach that of the "exact" system. (This is true even though hard-sphere dynamics do not accurately describe molecular systems!)

The basic premise underlying molecular dynamics is the ergodic theorem, which states that, in the limit of long times, the time average of a statistical measure of the system is equivalent to its ensemble average: $$\lim_{t \to \infty} \bar{A} = \left< A \right>,$$ where $\bar{A}$ is the time average and $\left<A\right>$ is the ensemble average of the quantity $A$.

Now, in principle, the time average must be measured for an infinite period of time to get exact equality, but in practice, the results will begin to converge within experimental error as the time and the system size increases. (The exact amount of time this will take can be measured, as is the decorrelation time required to achieve statistical independence between time samples, which can itself be used as a proxy to determine how long a simulation is required.) It can be shown using backward error analysis that the error of the numerical integration of a Hamiltonian can be made exponentially small; this applies to the integrators in molecular dynamics.

While numerical errors will take you on eventually exponentially divergent trajectories, in principle, so long you're generating the states in a manner consistent with the ensemble one is trying to observe, this should not ultimately impact the results, so long as a sufficient number of samples are collected. This means, for example, if we use a thermostat to ensure constant temperature, we can only equate the results with the results of the canonical ensemble.

An excellent discussion of the ergodic theorem can be found in Tolman's Principles of Statistical Mechanics.

Answer 2

Something that is at least very close to what your are asking for is the Shadowing Lemma, which roughly states that near any numerical solution of a dynamical system, you will find a real solution.

I am no expert on this, so I can only provide this a starting point for further research. In particular I do not know whether this translates to non-smooth dynamical systems (such as hard-sphere molecular dynamics) or whether any statement can be made about statistical properties. Still, it should provide a starting point for a literature search.