# The paradox of infinitesimal granularity:
Chaos and the reversibility of time in Newton’s theory of gravity

Simon Portegies Zwart<EMAIL_ADDRESS>Leiden Observatory, Leiden
University, PO Box 9513, 2300 RA, Leiden, The Netherlands  Tjarda Boekholt
<EMAIL_ADDRESS>Rudolf Peierls Centre for Theoretical
Physics, Clarendon Laboratory, Parks Road, Oxford, OX1 3PU, UK

###### Abstract

The fundamental laws of physics are time-symmetric, but our macroscopic
experience contradicts this. The time reversibility paradox is partly a
consequence of the unpredictability of Newton’s equations of motion. We
measure the dependence of the fraction of irreversible, gravitational N-body
systems on numerical precision and find that it scales as a power law. The
stochastic wave packet reduction postulate then introduces fundamental
uncertainties in the Cartesian phase space coordinates that propagate through
classical three-body dynamics to macroscopic scales within the triple’s
lifetime. The spontaneous collapse of the wave function then drives the global
chaotic behavior of the Universe through the superposition of triple systems
(and probably multi-body systems). The paradox of infinitesimal granularity
then arises from the superposition principle, which states that any multi-body
system is composed of an ensemble of three-body problems.

During the Great Plague of London in 1665, Isaac Newton went in quarantine for
two years at his home in Woolsthorpe. There he laid down the foundation for
his theory on gravitation. During last year’s COVID-19 pandemic, we have been
working on one of the consequences of Newton’s theory of gravitation and it’s
relation to chaos and time reversibility.

The chaotic problem of three bodies started with Henry Poincaré who argued in
1891 [1] that the problem was unsolvable, but he was wrong. This led to the
Kolmogorov-Arnold-Moser theorem, and a solution based on an infinite time-
series [2, 3], which therefore has little practical use [4]. Montgomery [5]
discovered a family of periodic solutions, One of which is presented in fig.
1.

Newton could, in principle, have realized that stable periodic solutions
existed. Such solutions, called braids, are found by laboriously searching
parameter space using action minimizers under topological constraints [6, 7].
Braids are stable, regular, and once found, they can be integrated relatively
easily [8]. Due to a lack of computer resources, Newton would not have been
able to calculate these solutions, even during his two years of social
distancing.

Figure 1: Two orbits of three-body problems, a periodic braid to the left, and
right a partial solution of a democratic triple system. The orbit to the left
is bounded and time reversible, it was found by [9]. The orbit to the right is
unbounded and not time reversible.

Non-periodic democratic solutions, as depicted in fig. 1, on the other hand,
are impossible to solve by hand [10], and even with modern digital computers
they can only be approximated [11].

## The paradox of infinitesimal granularity

The intrinsic chaotic nature of the democratic N-body problem causes small
errors introduced by computer round-off and integration time discretizations
to grow exponentially [12]. Still, solutions can be acquired by systematically
reducing the time step and increasing the computer’s mantissa until the
solution converges. Finding such a solution is an elaborate iterative process
in which the time step is reduced and the mantissa extended until the result
becomes independent of these.

For this purpose, we developed Brutus, a direct $N$-body solver for self-
gravitating systems of $N$ point-masses to arbitrary precision [13]. In Brutus
round-off is controlled by specifying the word-length $L_{w}$ in which
floating point arithmetic is performed using the MPFR arbitrary-precision
library [14]. To assure that the method of numerical integration does not
introduce a bias in our simulations we adopted the Gragg-Bulirsch-Stoer
algorithm [15], which combines the modified midpoint method [16] for solving
Newton’s ordinary differential equation [17] with Richardson extrapolation
[18] to improve the rate of convergence. We can now control the algorithmic
error by changing an energy tolerance parameter $\epsilon$ and the length of
the mantissa $L_{w}$: Both parameters can be varied simultaneously as in [19]

$L_{w}=4|\log_{10}(\epsilon)|+32\,\textrm{bits}.$ (1)

By repeating calculations with identical realizations of the initial
conditions while systematically changing $\epsilon$ and $L_{w}$ we can achieve
a solution that becomes independent of $\epsilon$ and $L_{w}$ to $p$ decimal
places.

Integrating converged solutions scales $\propto N^{2}$, but requires many
iterations, which progressively become more expensive. Acquiring a converged
solutions to an $N$-body problem is easily $10^{4}$ times more expensive than
regular integration techniques using IEEE 754 double precision arithmetic
[20]. Any non-converged solution to a self-gravitating system under Newton’s
equations of motion should be subject to considerable concern [21], although
statistically, the phase-space coverage of the solutions may be
indistinguishable from the true outcome space [22].

The measure of chaos, often quantified with the largest positive Lyapunov
exponent expressed in a time scale, $t_{\rm Ly}$, cannot be determined without
performing extensive numerical analyses. We do not know how chaotic $N$-body
systems are, but estimate values for $t_{\rm Ly}$ by running the same initial
realization twice until a converged solution is achieved: In the second
calculation, a small relative displacement of a random particle by $10^{-10}$
is introduced in one of the principal Cartesian coordinates. We calculate the
systems until a single star escapes. We then turn time backward in the
computer, and calculate the final conditions back to their initial
realization, recovering the $10^{-10}$ initial displacement of the perturbed
solution with respect to the unperturbed solution [19]. When performing this
operation for a large number of randomly generated triples, we can calculate
how many of those converged solutions were successfully run backward as a
function of $\epsilon$. In figure 2 we present this fraction of irregular (or
irreturnable) solutions $f_{\rm irr}$ as a function of $\epsilon$.

Figure 2: The fraction of irreversible systems of converged solutions for
zero-angular momentum-triples ($J=0.0$ in red) and triples with the maximum
amount of angular momentim giving virial initial conditions ($J=0.19$ in
blue). Both curves can be fitted by a power law:
$\log_{10}\,f_{\rm{irr}}\simeq 0.036\log_{10}\,\epsilon+0.25$ for $J=0$, and
$\log_{10}\,f_{\rm{irr}}\simeq 0.016\log_{10}\,\epsilon+0.11$.

We were surprised to find that triples with finite angular momentum turn out
to be more chaotic than their zero-angular momentum counterparts. The
underlying reason hides in the triple’s lifetime; high-angular momentum
triples have a lower rate of close encounters, and those encounters are
typically at a larger mutual distance. As a consequence, low angular momentum
systems tend to have a higher probability of violent, stellar ejecting
encounters causing these systems to live shorter on average compared to their
high angular-momentum equivalents. The shorter lifetime of low-angular
momentum systems cause them to have less time for chaos to propagate, even at
the same Lyapunov time scale.

Another remarkable aspect of figure 2 is that about $\
{\raise-2.15277pt\hbox{$\buildrel>\over{\sim}$}}\ 10$% of virialized triples
are so chaotic that these trajectories are affected by perturbations on a
relative scale of $10^{-70}$ [11]. When considering the size of the universe,
a perturbation of $1.3\cdot 10^{-34}$ Å, or $\sim 8.1\cdot 10^{-10}$ lp would
then still drive exponential divergence to a level of making the system
unpredictable at any scale.

## Global chaos due to Newton’s equations of motion

Newton operates on all scales and independently of the intrinsic nature of the
objects, either being black holes or the balls of a juggler. The interference
of other physical processes is irrelevant so long as the dynamics is driven in
part by Newton’s equations of motion. The fundamental irreversibility of time
in Newton’s equations of motion propagates to all scales and affects every
object [21]. Time irreversibility on macroscopic scales then is a consequence
of the chaotic nature of Newton’s equation of motion, as the spontaneous
collapse of the wave function [23] introduces stochastic processes on a
microscopic scale that propagate to macroscopic scales.

After exploring the three-body problem, we continue with virialized equal-mass
large $N$-systems in a homogeneous unit cube (we adopted dimension-less
coordinates in units with $G=1$ [24, 25]). We explore systems with up to
$N=131072$ ($\equiv 128$ k). The results are presented in fig 3, where we show
the Lyapunov time scale for of ensembles of multi-body systems. These
calculations are continued for about 3 or 4 dynamical times, but the
measurements are carried out on the sub-sample of the data, from $t=1/\pi$
until growth of the phase-space distance between two converged solutions of a
perturbed and a non-perturbed initial realization exceeds 0.1. In some cases,
the exponential growth of the phase-space distance eventually flattens to a
$\log_{10}(\delta)\propto t$ scaling relation (consistent with the expectation
for a relaxation-driven process, see [26]). We acquired converged solutions
for $N=4$, 8, 16, 32, 64, 128, and $N=1024$, and reprehensive [19] solutions
for $N$ up to 128k particles using a regular fourth-order Hermite integration
algorithm using IEEE 754 double precision arithmetic.

The Lyapunov time-scale $\propto\ln(\ln(N))$ [27] with a discontinuity near
$N=32$ due to the transition from crossing-time dominated to relaxation time
dominated dynamics. The distribution width of Lyapunov timescales for larger
$N$ is narrower than for small $N$. For up to $N\sim 128$ k the majority of
systems tend to be more regular than for $N=3$. For larger $N$, however, all
systems tend to be even more chaotic than the most chaotic three-body systems.

Figure 3: Estimate of the Lyapunov timescale as a function of the number of
particles. Here the horizontal axis is not linear, but in $\ln(\ln(N))$ to
illustrate the scaling proposed in [27]. The different symbols and colors
represent different calculations (see legend). The vertical bars, plotted for
Newton’s Hermite only, show the root-mean-square of the dispersion in the
series of solutions. The error bars in the results obtained with Brutus are
statistically indistinguishable from the presented bars. The blue horizontal
dotted line accentuates the boundary where large $N$ systems become on average
as chaotic as the most chaotic 3-body systems in our adopted initial
conditions.

We argue that if wave functions collapse spontaneously, as is the case in
Ghirardi-Rimini-Weber theory [23], microscopic stochastic effects propagate
non-linearly to macroscopic scales through Newtonian dynamics. The motion of a
considerable fraction of triples of self-gravitating bodies, and possibly the
vast majority of systems with $\
{\raise-2.15277pt\hbox{$\buildrel>\over{\sim}$}}\ 128$ k are affected by the
spontaneous collapse of the wave function. As a consequence, the evolution of
self-gravitating multi-body system is eventually stochastic, due to the
proliferation of chaotic processes that propagate from microscopic scales to
macroscopic scales. This connection between the microscopic quantum scale and
macroscopic scales cause the dynamics to be fundamentally time irreversible;
Newton’s equations of motion are incalculable (they cannot be calculated
precisely in a deterministic nature) and inconsequential (there is no point in
achieving such a result).

Changes in the macroscopic distribution of mass in the Galaxy are driven by
chaos on microscopic scales. Due to the spontaneous collapse of the wave
function, perturbations at Planck length scales are magnified exponentially by
the chaotic dynamics of planetary systems, stellar multiples, star clusters
and giant molecular clouds. This growth cannot continue indefinitely, and
converges once the perturbation has approximately reached the size of the
Galaxy [28]. The divergence in orbital energy then starts exponentially but
eventually transforms into a much slower random walk, consistent with the
global relaxation process [26, 29]. Hence, one may think of a Galactic
butterfly effect, where fluctuations lead to unpredictable local events, but
which only gradually affect the global Galactic climate.

## Acknowledgments

It is a pleasure to thank David Albert and Douglas Heggie for discussions.

## References

  * [1] H. Poincaré,  _Mémoires et observations. Sur le problème des trois corps_ ,  Bulletin Astronomique, Serie I 8, 12 (1891).
  * [2] K. F. Sundman,  _Mémoire sur le problème des trois corps_ ,  Acta Math. 36, 105 (1912).
  * [3] K. F. Sundman,  _Recherches sur le problème des trois corps. Nouvelles recherches su r le problème des trois corps_ ,  Acta Soc. Sc. Fennicae 34, No. 6, 1-43 (1907); 35, No. 9, 1-27 (1909). (1909).
  * [4] M. Henkel,  _Sur la solution de sundman du problème des trois corps_ ,  Philosophia Scientiae 5(2), 161 (2001).
  * [5] R. Montgomery,  _The n -body problem, the braid group, and action-minimizing periodic solutions_ ,  Nonlinearity 11(2), 363 (1998).
  * [6] R. J. Vanderbei,  _New Orbits for the n-Body Problem_ ,  Annals of the New York Academy of Sciences 1017(1), 422 (2004),  doi:10.1196/annals.1311.024,  eprint astro-ph/0303153.
  * [7] D. Yan, R. Liu, X. Hu, W. Mao and T. Ouyang,  _New Phenomena in the Spatial Isosceles Three-Body Problem with Unequal Masses_ ,  International Journal of Bifurcation and Chaos 25(12), 1550169-313 (2015),  doi:10.1142/S0218127415501692.
  * [8] D. C. Heggie,  _A new outcome of binary-binary scattering_ ,  MNRAS 318(4), L61 (2000),  doi:10.1046/j.1365-8711.2000.04027.x,  eprint astro-ph/9604016.
  * [9] X. Li, Y. Jing and S. Liao,  _Over a thousand new periodic orbits of a planar three-body system with unequal masses_ ,  Publ. Astr. Soc. Japan 70(4), 64 (2018),  doi:10.1093/pasj/psy057,  eprint 1709.04775.
  * [10] P. Hut and J. N. Bahcall,  _Binary-single star scattering. I - Numerical experiments for equal masses_ ,  ApJ 268, 319 (1983),  doi:10.1086/160956.
  * [11] T. C. N. Boekholt, S. F. Portegies Zwart and M. Valtonen,  _Gargantuan chaotic gravitational three-body systems and their irreversibility to the Planck length_ ,  MNRAS 493(3), 3932 (2020),  doi:10.1093/mnras/staa452,  eprint 2002.04029.
  * [12] M. Valtonen, J. Anosova, K. Kholshevnikov, A. Mylläri, V. Orlov and K. Tanikawa,  _The Three-body Problem from Pythagoras to Hawking_ ,  doi:10.1007/978-3-319-22726-9 (2016).
  * [13] S. F. Portegies Zwart and T. C. Boekholt,  _Numerical verification of the microscopic time reversibility of newton’s equations of motion: Fighting exponential divergence_ ,  Communications in Nonlinear Science and Numerical Simulation 61, 160 (2018),  doi:https://doi.org/10.1016/j.cnsns.2018.02.002.
  * [14] S. Portegies Zwart and T. Boekholt,  _On the minimal accuracy required for simulating self-gravitating systems by means of direct N-body methods_ ,  ApJL in press, ArXiv e-prints 1402.6713 (2014),  eprint 1402.6713.
  * [15] R. Bulirsch and J. Stoer,  _Fehlerabschätzungen und extrapolation mit rationalen funktionen bei verfahren vom richardson-typus_ ,  Numerische Mathematik 6, 413 (1964),  10.1007/BF01386092.
  * [16] W. Gragg,  _On extrapolation algorithms for ordinary initial value problems_ ,  Journal of the Society for Industrial and Applied Mathematics 2(3), 384 (1965).
  * [17] I. Newton,  _Philosophiae Naturalis Principia Mathematica_ , vol. 1 (1687).
  * [18] L. Richardson,  _The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, with an Application to the Stresses in a Masonry Dam_ ,  Phil. Trans. R. Soc. Lond. A 210(459-470), 307 (1911).
  * [19] S. F. Portegies Zwart and T. C. N. Boekholt,  _Numerical verification of the microscopic time reversibility of Newton’s equations of motion: Fighting exponential divergence_ ,  Communications in Nonlinear Science and Numerical Simulations 61, 160 (2018),  doi:10.1016/j.cnsns.2018.02.002,  eprint 1802.00970.
  * [20] S. F. Portegies Zwart, T. C. N. Boekholt, E. H. Por, A. S. Hamers and S. L. W. McMillan,  _Chaos in self-gravitating many-body systems. Lyapunov time dependence of N and the influence of general relativity_ ,  A&A 659, A86 (2022),  doi:10.1051/0004-6361/202141789,  eprint 2109.11012.
  * [21] R. H. Miller,  _Irreversibility in Small Stellar Dynamical Systems._ ,  ApJ 140, 250 (1964),  doi:10.1086/147911.
  * [22] T. Boekholt and S. Portegies Zwart,  _On the reliability of N-body simulations_ ,  Computational Astrophysics and Cosmology 2, 2 (2015),  doi:10.1186/s40668-014-0005-3,  eprint 1411.6671.
  * [23] G. C. Ghirardi, A. Rimini and T. Weber,  _Unified dynamics for microscopic and macroscopic systems_ ,  Phys. Rev. D 34, 470 (1986),  doi:10.1103/PhysRevD.34.470.
  * [24] M. Hénon,  _Monte Carlo Models of Star Clusters (Part of the Proceedings of the IAU Colloquium No. 10, held in Cambridge, England, August 12-15, 1970.)_ ,  Ap&SS 13, 284 (1971).
  * [25] D. C. Heggie and R. D. Mathieu,  _Standardised Units and Time Scales_ ,  In P. Hut and S. L. W. McMillan, eds., _The Use of Supercomputers in Stellar Dynamics_ , vol. 267 of _Lecture Notes in Physics, Berlin Springer Verlag_ , p. 233,  doi:10.1007/BFb0116419 (1986).
  * [26] P. Hut and D. C. Heggie,  _Orbital Divergence and Relaxation in the Gravitational N-Body Problem_ ,  Journal of Statistical Physics 109(5-6), 1017 (2002),  doi:10.1023/A:1020472526203,  eprint astro-ph/0111015.
  * [27] J. Goodman, D. C. Heggie and P. Hut,  _On the Exponential Instability of N-Body Systems_ ,  ApJ 415, 715 (1993),  doi:10.1086/173196.
  * [28] M. Valluri and D. Merritt,  _Orbital Instability and Relaxation in Stellar Systems_ ,  In V. G. Gurzadyan and R. Ruffini, eds., _The Chaotic Universe_ , pp. 229–246,  doi:10.1142/9789812793621\\_0014 (2000), eprint astro-ph/9909403.
  * [29] M. Hemsendorf and D. Merritt,  _Instability of the Gravitational N-Body Problem in the Large-N Limit_ ,  ApJ 580(1), 606 (2002),  doi:10.1086/343027,  eprint astro-ph/0205538.