---
abstract: 'Numerical solutions to Newton’s equations of motion for chaotic self gravitating systems of more than 2 bodies are often regarded to be irreversible. This is due to the exponential growth of errors introduced by the integration scheme and the numerical round-off in the least significant figure. This secular growth of error is sometimes attributed to the increase in entropy of the system even though Newton’s equations of motion are strictly time reversible. We demonstrate that when numerical errors are reduced to below the physical perturbation and its exponential growth during integration the microscopic reversibility is retrieved. Time reversibility itself is not a guarantee for a definitive solution to the chaotic N-body problem. However, time reversible algorithms may be used to find initial conditions for which perturbed trajectories converge rather than diverge. The ability to calculate such a converging pair of solutions is a striking illustration which shows that it is possible to compute a definitive solution to a highly unstable problem. This works as follows: If you ([i]{}) use a code which is capable of producing a definitive solution (and which will therefore handle converging pairs of solutions correctly), ([ii]{}) use it to study the statistical result of some other problem, and then ([iii]{}) find that some other code produces a solution [S]{} with statistical properties which are indistinguishable from those of the definitive solution, then solution [S]{} may be deemed veracious.'
address: |
    $^{1}$Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands.\
    $^{2}$CIDMA, Departamento de F' isica, Universidade de Aveiro, Campus de Santiago, 3810-193 Aveiro, Portugal.
author:
- 'Simon F. Portegies Zwart$^{1}$ Tjarda C. N. Boekholt$^{2}$'
title: ' Numerical verification of the microscopic time reversibility of Newton’s equations of motion: Fighting exponential divergence'
---

Introduction
============

General analytic solutions to problems in Newtonian dynamics [@Newton:1687] can only be achieved for a single particle, $N=1$, or for $N=2$ [@Euler1760; @Lagrange1772 see [@doi:10.1119/1.4867608] for a historical overview]. Families of periodic solutions exist for $N>2$ [@0951-7715-11-2-011], and in particular the parameter-space search of [@2017arXiv170904775L] has succesully identified more than 1000 new periodic solutions to the restricted 3-body problem, suggesting that the number of such solutions is interminable. The latter is consistent with the existence of periodic solutions within the KAM theorem [@1954Kolmogorov; @1962Moser; @Arnold1963]. For all other solutions approximate methods have to be employed. These approximate solutions are unreliable due to the intrinsic chaotic nature of the problem, leading to exponential growth of small perturbations [@1892mnmc.book.....P; @1899lmnd.book.....P; @1908BSAFR..22..389P]. This notion is not new, as Poincare [@1903mnm..book.....P p.138] already pointed out:\

-   [*Une cause très petite, qui nous échappe, détermine un effet considérable que nous ne pouvons pas ne pas vois, et alors nous disons que cet effect est dû au hasard. Si nous connaissions exactement les lois de la nature et la situation de l’Universe à l’instant initial, nous pourrions prédire exactement la situation de ce même Universe a un instant ulterieur. Mais, lors même que les lois naturelles n’auriaient plus de secret pour nous, nous ne pourrions connâitre la situation initiale qu’approximativement. Si cela nous perment de prévoir la situation ulterieure avec la même approximation, c’est tout ce qu’il nous faut, nous disons que le phénomènon a été prévu, qu’il est régi par des lois; mais il n’en est pas toujours ainsi, il peut, arriver que de petites différences dans les conditions initiales en engendrent de très grandes dans les phénomènes finaux; une petitt erreur sur les premières produirait une erreur énorme sur les derniers. La prédiction devient impossible et nous avons le phénomène fortuit.* ]{}[^1]

As a result, small errors in the temporal or spatial discretizations, or in the numerical integration scheme cause any solution to eventually become invalid [@1964ApJ...140..250M], maybe even within a few dynamical time scales [@1986LNP...267..212D].

The relation between instability and chaoticity manifests itself by the high sensitivity to small changes in the initial conditions [@Hadamard1898; @1908BSAFR..22..389P; @1963JAtS...20..130L]. Divergent behavior is often demonstrated by performing two calculations with a slight offset $\delta$ in phase space. The exponential growth with time of this phase-space distance $\delta(t)$ is expressed in the largest Lyapunov exponent [@lyapunov1892]. For a chaotic system, the associated e-folding time is positive and finite. The rate of divergence of two trajectories in phase space is characterised by this $e$-folding time, which is specific for the $N$-body realization. Since Newton’s equations of motion are time reversible a finite parameter space must exist for which $\delta(t)$ diminishes [@1930Perron]; namely, when an initial divergent pair of trajectories is time reversed, it should return to the moment the perturbation was introduced. Such behavior can only be established when solving Newton’s equations of motion with sufficient accuracy and enough precision to guarantee that the accumulated error and its exponential growth remains below the introduced deviation $\delta(t=0)$. In Appendix A we present a short glossary of terms used in the manuscript to help the reader appreciate our disucssion.

The combination of hypersensitivity to small perturbations and non-integrability prevents us from manifesting Newton’s reversibility numerically, because the underlying method should be accurate as well as precise in order to arrive at a converged solution. High-order and symplectic numerical solvers tend to be insufficiently accurate, in the sense that reducing the time step tends to interfere with the lack of precision due to the growth of the error introduced by round off [@2015ComAC...2....2B]. The combination of the exponential growth of small perturbations and the inevitability of numerical errors forms the fundamental argument why solving Newton’s equations of motion still is one of the hardest problems in computational physics. Individual numerical trajectories quickly forget their initial conditions[^2]. Instead of integrating a single realization until a converged solution is obtained one often considers ensembles of trajectories in phase space, starting with a random sample of initial realizations close to, and possibly including, the objected realization [@Szebehely1991]. Each of these realizations is subsequently calculated and the objected phase space is anticipated to provide a probability density distribution around the true solution. In principle this leads to a veracious solution, but this is not guaranteed, in which case the ensamble average may well be distinct from the true solution. [@2015ComAC...2....2B] demonstrated that an ensamble of reprehansible solutions was statistically indistringhuisable from the an ensable of converged solutions with idential realizations, a quality we call [*nagh Hoch*]{}.

It is not clear whether the strict time-reversibility of Newton’s equations of motion can be feasibly maintained from a practical numerical point of view. The authors in [@2008MNRAS.388..965L; @2016ttp..book.....V], claim that the chaotic nature of the underlying dynamical processes then prevents us, via the second law of thermodynamics, to reverse time and calculate backward. [@2008MNRAS.388..965L] associate such irreversible dynamical process to the increase of entropy and the arrow of time in systems that, from a theoretical perspective should be strictly deterministic. We instead call this process numerical confusion.

Chaos in $N$-body systems is often confused with a number of side effects, including uncertainties in the initial conditions, round-off, integration errors, spatial and temporal discretization, as was already pointed out by [@1991pscn.proc...47H]. In a chaotic self-gravitating system errors grow exponentially but with sufficient accuracy and precision this does not prevent the system from resulting in a definitive solution. Such calculations are time reversible [@1994CeMDA..60..365A], but time reversibility itself is insufficient to guarantee a definitive result because the method may still fail to resolve close encounters either by lack of accuracy or by lack of precision. Two neighboring solutions on the other hand recover their initial offset when reversed after some finite time. In such a time-reversed evolution trajectories approach each other until the initially introduced phases-space distance $\delta(t=0)$ is reached, after which they diverge again. Such a phase-space offset deminishing simulation leads to exponential convergence within a finite time interval. Most studies in $N$-body dynamics focus on the largest positive Lyapunov exponent to measure the rate of divergence between neighbouring trajectories rather than on converging solutions [@LichtenbergAndLieberman1992]. In the next § we demonstrate that such converging trajectories exist, and we speculate that they are important when the duration of the convergence phase becomes comparable to the simulation time-scale of interest.

Results for the Pythagorean 3-body problem
==========================================

Numerical verification of the time reversibility of Newton’s equations of motion can be obtained by demonstrating the existence of converging trajectories in phase space, and the ability to calculate any non-colliding initial condition to an arbitrary point in time and backward to recover the initial realization with an accumulated integration error much smaller than the initial offset. In colliding orbits, such as the case in homothetic hyperbolic orbits [@10.2307/24893288; @1974InMat..27..191M], the singularity can only be reached asymptotically which restricts the time-scale over which the solution can be obtained numerically. In Fig.\[Fig:errorproagation\] we demonstrate numerical time-reversibility by presenting the result for the Pythagorean $3$-body problem [@1913AN....195..113B], which was first solved in 1967 [@1967AJ.....72..876S] and which is chaotic [@1994CeMDA..58....1A] but relatively short lived. In Tab.\[ICs\] we present an initial realization for the Pythagorean problem (and the perturbed conditions). In the Appendix B we present a brief summary of the results of the calculations. Snapshots of the converged solution are available in the online supplementary material. The possibility to time-revers the calculations apply to any self-gravitating multi-body system and the choice of the Pythagorean problem is motivated by limited resources and its historical context. We integrate the equations of motion for the $3$-bodies using the [ Brutus]{} arbitrary precise $N$-body code [@2015ComAC...2....2B; @2016MNRAS.461.3576B], which is part of the AMUSE software environment . In [Brutus]{} the integration accuracy is controlled by the tolerance $\epsilon$ in the Bulirsch-Stoer [@springerlink:10.1007/BF01386092] integrator, and the precision is controlled using length of the mantissa $\eta$. By changing these, [Brutus]{} can be tuned until the solution converges. Finding a converged solution is an iterative procedure in which the values of the Bulirsch-Stoer tolerance $\epsilon$ and the length of the mantissa $\eta$ are improved until to a pre-determined precision the phase-space coordinates of the solution becomes independent of these parameters.

![Evolution of the phase-space separation $\delta$ between two solutions of the Pythagorean $3$-body problem [@1913AN....195..113B]. The perturbed trajectories are calculated using an initial offset of $10^{-10}$ in the x-coordinate of the least massive body (see Tab.\[ICs\]). After $t=100$ we reverse all the velocities and continue the integration to $t=200$. The accuracy of the Bulirsch-Stoer integrator is tuned using two parameters: the tolerance parameter $\epsilon$ and the word-length $\eta$. The identical evolution of both systems back to their initial conditions is achieved using a Bulirsch-Stoer tolerance of $10^{-24}$ and a word-length of 136bits (see Tab.\[Final\] for the interim and final results). Less accurate and less precise calculations are presented with the colored curves (see top left inset for the parameters). \[Fig:errorproagation\] ](./delta_vs_t_precise_varprecision.pdf){width="100.00000%"}

  ------ -------------- ---- --- ---
  mass                           
  3      1              3    0   0
  4      -2             -1   0   0
  5      1              -1   0   0
                                 
  3      1.0000000001   3    0   0
  4      -2             -1   0   0
  5      1              -1   0   0
  ------ -------------- ---- --- ---

  : Initial realization of the Pythagorean problem as adopted in our numerical experiment. Each line contains mass, $x$- and $y$-position and $x$- and $y$-velocity in units in which we adopt $G=1$. The $z$-coorindates are zero and omitted from the table. The mit-time and final simulation results are presented in the Appendix B. []{data-label="ICs"}

We integrate the Pythagorean problem for $100$ time-units ($G=1$, the units of mass and length are implicitly defined in Table\[ICs\]): any time and any initial realization would have worked but may require a different tolerance and word-length in the integration scheme [@2015ComAC...2....2B]. Two trajectories are calculated, for one we offset the $x$-position of the least massive body by $\delta =
10^{-10}$ (see Tab.\[ICs\]) but a similarly small offset in any of the Cartesian coordinates would have worked. The three bodies engage in a resonant interaction, which lasts until one particle escapes[^3] at $t \simeq
63.4^{\ref{Footnote:Dissolution}}$, consistent with earlier results [@1967AJ.....72..876S; @1986LNP...267..212D; @1994CeMDA..58....1A; @2008MNRAS.386..425F]. During this period two distinct chaotic behaviors manifest themselves. The system starts with a relatively long $e$-folding time scale of $\sim 28.9$ of the exponential growth of the initial phase-space separation (see Fig.\[Fig:errorproagation\]). After 28.6 time units a transition occurs, possibly initiated by a close three-body encounter (see also [@2008MNRAS.386..425F]), leading to a faster growth of the phase-space distance. In this lap the system exhibits a much shorter $e$-folding time scale of $\sim 2.1$. The unperturbed and the perturbed systems eventually dissolve into a bound pair and a single body[^4]. After this the phase-space separation grows approximately linearly. This evolution is consistent with results of previous studies [@1986LNP...267..212D; @1994CeMDA..58....1A].

We stop the calculation at $t = 100$ by which time the phase-space separation between the two solutions has grown by a factor of $10^{9}$ to $\delta \sim 0.1$. We subsequently reverse the velocities of all particles and continue the calculation to $t=200$. At this moment, we compare the resulting positions and velocities with the original initial conditions (see Tab.\[Final\]). In order to realize a reversible calculation in which the initially introduced perturbation grows by 9 orders of magnitude the numerical round-off error and its exponential growth has to remain well within a factor of $10^{18}$ of the initially introduced offset. With a Bulirsch-Stoer tolerance of at most $10^{-24}$ and a word-length of at least 136bits, the final and initial conditions are identical to a sufficiently large mantissa to recover the initial offset between the two initial realizations. In Tab.\[Final\] we present the final conditions for the unperturbed and the perturbed initial realizations.

Discussion
==========

With the calculations of the Pythagorean 3-body problem, we demonstrate that irrespective of the exponential growth of perturbations the equations of motion can be solved accurately although not exactly using floating-point arithmetic. Newton’s laws of motion are time reversible and chaos does not prevent this. Microscopic irreversibility (microscopic in the sense of detailed positions and velocities), as found in previous studies [@hollinger2012nature; @2008MNRAS.388..965L] is then the result of the secular growth of numerical errors. From a theoretical perspective, [@MR103254; @MR103256] already argued that irreversibility is not an intrinsic quality of chaos. Earlier claims of reversibility in $N$-body simulations [@2017arXiv170407715R] using high-order direct integrators with double floating point precision results in an evolution of the phase-space distance comparable to the blue curve in Fig.\[Fig:errorproagation\], but much higher precision is needed in order to recover the initial conditions (black curve).

In $N$-body integrations the numerical error accumulates due to round-off and discretization errors. Irreversibility of these numerical errors will render the entire calculation irreversible. Integer arithmetic, as alternative to floating points, could solve this problem because round-off in that case is deterministic and time symmetric [@1992PhyD...56....1E]. This makes the calculations strict time symmetric, but they remain reprehansible in the sense that they are not necessarily accurate and not precise. We test this hypothesis by repeating the calculations (see Tab.\[ICs\]) using the [Janus]{} integrator in the [ Rebound]{} package [@2017arXiv170407715R], which employs integer arithmetic. We verify that integer arithmetic in the $N$-body problem leads to a time reversible result, but it is insufficiently accurate to resolve close encounters, irrespective of the time symmetry in the round off. This is evidenced by a growth of the phase-space distance that deviates from the curve presented in Fig.\[Fig:errorproagation\]. In our experiment we explored a time step size of 0.1 to $10^{-26}$ but the solution did not converge. We tested the effect of floating-point behavior on the growth of the phase-space distance using different integrators, 4th order predictor-corrector Hermite [@1992PASJ...44..141M], using the implementation of [ph4]{}, which is part of the [AMUSE]{} package [@AMUSE] and the Verlet leap-frog integrator [@PhysRev.159.98] which is part of the [Starlab]{} package and found identical behavior up to the point the curves start to deviate some time after the time-reversed calculation (see the blue curve in fig.\[Fig:errorproagation\]).

The parameter space for converging trajectories is small. In our experiment, where we follow the growth of the perturbation over $9$ orders of magnitude in the normalized $6$-dimensional phase space, convergence is not reached before we extend the mantissa of the calculation to exceed 18 decimal places (for calculating forwards in time and backwards). In order to continue hiding the round-off, the entire calculation had to be performed with 40 significant figures (see Tab.\[Final\]). A naive estimate of the available parameter space for finding such converging solutions is then $\propto
1/10^{18}$. The arrow of time then manifests itself by the larger parameter space for which irreversible perturbations lead to a divergent growth rather than convergence.

Once we are able to time-reverse a relatively simple three-body simulation, it remains unclear to what extend this can be attributed to a complicated multi-body simulation. In a gedanken experiment one could evolve a star cluster to the moment just before the first hard binary forms in the collapsing core. When starting from this point we would calculate backwards in time, the system should return to its initial realization. When at this point we introduce a small perturbation to one of the single objects the system will not return to its initial realization. But we wonder how large a perturbation are we allowed to introduce for the system to still climb out of core collapse, rather than continue to collapse until the actual formation of the first binary?

Acknowledgments {#acknowledgments .unnumbered}
===============

We thank Fouvry Jean-Baptiste, Vincent Icke, Junichiro Makino, Cole Miller, Inti Pelupessy, Dagmar Portegies Zwart, Hanno Rein, Mauri Valtonen and the accurate and precise non-anonymous referee for discussions. This work was supported by the NWO (grant \#621.016.701 \[LGM-II\]) the European Union’s Horizon 2020 research and innovation programme under grant agreement No 671564 (COMPAT project), and Fundaç\~ ao para a Ci\^ encia e a Tecnologia (grant SFRH/BPD/122325/2016), and the Center for Research & Development in Mathematics and Applications (CIDMA) (strategic project UID/MAT/04106/2013) and from ENGAGE SKA, POCI-01-0145-FEDER-022217, funded by COMPETE 2020 and FCT, Portugal.

Appendix A {#AppendixA .unnumbered}
==========

Here we present a small glossary of terms used in this manuscript.

-   [**Accuracy, accurate:**]{} degree to which a calculation is numerically accurate in terms of the discretization of the numerical method. Accuracy is controlled by the time step size of the integration. In the Bulirsch-Stoer method the accuracy is of the order the tolerance $\epsilon$.

-   [**Converged:**]{} A numerical solution that to a certain pre-determined length of the mantissa will not change irrespective of the further increase in the mantissa during the calculation (precision) or a further decrease of the time step (accuracy). The term should be used in conjuction with the length of the matissa of the converged simulation. The iterative process that leads to a converged solution is considered [*converging*]{}, which should not be confused with the term [*convergence*]{}. Converged solutions are considered [*definitive*]{}.

-   [**Convergence:**]{} The gradual decrease in the phase-space distance between two different solutions over a course of time.

-   [**Definitive solution:**]{} The solution for some initial realization that is expected to be numerically indistinghuishable from the true solution. From a mathematical perspective a definitive solution represents a computable non-periodic pseudo-orbit [@1987PhLA..122..399P].

-   [**Initial conditions:**]{} statistical description of boundary values from which, using some random number sequence with a seed value, an [*initial realization*]{} can be generated.

-   [**Initial realization:**]{} actual specific phase-space coordinates that uniquely describes a system of particles. Any initial realization results after some time in one unique final realization.

-   [**nagh Hoch:**]{} The concept that an ensemble of random initial realizations in a wide range of parameters gives statistically the same result as the converged solutions of the same ensemble realizations. This concept is a quality of the numerical method, but we speculate that this quality also applies to real systems. nagh Hoch means as much as “similar appearance” [@KlingonDictionary].

-   [**neighboring solutions:**]{} Two initial realizations that started with an infinitesimal offset integrated over a time scale sufficiently short that the exponential growth of the offset remains below a predetermined limit (typically small compared to the initial size of the system but possibly large compared to the initial offset) [@1992PhLA..166..340F].

-   [**Precision, precise:**]{} degree to which a calculation is numerically reproducible in terms of the length of the mantissa. Precision in Brutus is controlled using length of the mantissa $\eta$.

-   [**Reprehansible reprehensive:**]{} Solution to Newton’s equations of motion for which the accumulation of numerical errors and the system’s response exceeds the exponential growth of the initial offset $\delta$.

-   [**Time reversibility:**]{} The ability of a numerical integrator to recover the initial realization from reversing the final realization. Time reversibility does not guarantee that the solution is converged as it can already be obtained by insisting time symmetry in numerical errors without resolving close encounters.

-   [**True solution:**]{} The unique solution Nature would provide in the limited physical domain (such as Newtonion dynamics in isolation) for a specific initial realization.

-   [**Veracity, veracious:**]{} the concept that a limited and bounded variation of the initial conditions near a pre-selected initial realization gives statistically an indistinguishable ensemble average as a single [*converged*]{} solution. A veracious solution can be perceived as the experimental observation of the true solution [@2003PhRvL..90e4104H]. The presumption that Newtonian simulations are [*veracious*]{} is used for example in stability studies of the Solar System [@2008DDA....39.1514V; @2015ApJ...798....8Z see e.g.].

Appendix B {#AppendixB .unnumbered}
==========

The Pythagorean $3$-body problem [@1913AN....195..113B] starts with three point masses of 3, 4, and 5 in the $x$-$y$–plane and with zero velocities. The three masses are located in the corners of a Pythagorean triangle with (1, 3) for the first mass, (-2, 1) for the second and (1, -1) for the third and most massive object. The values are dimensionless modified $N$-body units [@1986LNP...267..233H]. In Tab.\[ICs\] we present these coordinates, and the initial realization of our perturbed solution.

In Fig.\[Fig:errorproagation\] we present the phase-space separation between the converged solution and a perturbed solution, and in addition a number of less accurate calculations. The less accurate solutions do not recover the initial phase-space distance of $10^{-10}$ but continue to diverge at some moment during the integration. The moment they start to diverge depends on the adopted accuracy.

In Tab.\[Final\] we present the mid-point (at $t=100$) and final conditions (at $t=200$) for the converged unperturbed configuration and the perturbed initial realization. Here we omitted the third coordinate ($z$) because it is $0$ and remains $0$ throughout the calculation. We did not round the numbers but present them as they are represented in the computer in 40 significant figures.

[llll]{} mass & position ($x$, $y$) & velocity ($v_x$, $v_y$)\
\
\
3& 23.177734271287932984376719864262469663010 & 0.5319181795481124436207125296002406419911\
& 68.527351147839055139293277996603512381330 & 1.5806796434841945398880662764620232011290\
4& -7.259955925491511142826576066493924610372 & -1.3069988754994211535542088320251159535050\
& -22.988301724483254153623924485962763879080 & -0.3946057714182357399709966348235988763510\
5& -8.098675822379550876364771064194141771170 & 0.7264481926706694566709395478765597018697\
& -22.725769309116829760676827209503214745270 & -0.6327231689559281319560424580219478896762\
\
3& 1.000000000000029147791398925636660776567& -2.270447738062808857494824732048534406358e-14\
& 3.000000000000036383772882843757984240521& 4.774666017242739047666501771157480010792e-14\
4& -1.999999999999974915829187984060851093363& 4.774665242050519395810021558737175983173e-14\
& -0.999999999999810155587671818578010637774& 2.270447738062669364553550881039426305535e-14\
5& 0.999999999999962443988511031800475014671& -2.457463550802730796268351737568898744656e-14\
& -1.000000000000173705793592251440952990686& -4.681157800795778807649659958570352078546e-14\
\
\
3& 23.190521482793319989765681634995570288100 & 0.5322503567960587903769290142249037912735\
& 68.556547837283584752026354886258680715190& 1.5814556205092701449052881349866961910410\
4& -7.287343535103771287105271550570713456589 & -1.4321791210002474579596089005290112193940\
& -22.995007467016918242398592448181673127520& -0.3552190002153252430432594844029797100956\
5& -8.084438061532974964175191740749603394790 & 0.8263930827225626921415297118904104229181\
& -22.737922728756616257296938973228645901360& -0.6646981721333018925085652934724455130194\
\
3& 1.000000000100049154879867130170780709564 & -3.828886532300641220319197711163974610107e-14\
& 3.000000000000061357649755490016408020941& 8.052003949335830235349211152837151414787e-14\
4& -1.999999999999957698016272046329507322690 & 8.052002826247233760595438150140762918720e-14\
& -0.999999999999679846094415534203355381441& 3.828886532501206909565529941811723471334e-14\
5& 0.999999999999936665485097358523258172108 & -4.144270341617403501659375346577704899564e-14\
& -1.000000000000292937714320866161883344304& -7.894311595602463066508120495426148312190e-14\

chaos.bbl

[^1]: A tiny difference, which we didn’t notice, has considerable repercussions that we can’t ignore, and then we say that this effect is pure chance. If we knew exactly the laws of nature and the conditions of the universe at the initial moment, we could exactly predict the situation of that same universe at any succeeding moment. But even if it were the case that the natural laws had no longer any secret for us, we could still only know the initial situation approximately. If that enabled us to predict the succeeding situation with the same approximation, that is all we require, and we should say that the phenomenon had been predicted, that it is governed by laws. But it is not always so; it may happen that small differences in the initial conditions produce very great ones in the final phenomena. A small error in the former will produce an enormous error in the latter. Prediction becomes impossible, and we have the fortuitous phenomenon.

[^2]: With the term “forget” we mean that the behaviour of a solution (in some chaotic region) becomes statistically independent of its initial conditions (see [@2014ApJ...785L...3P]).

[^3]: A particle is considered to escape when it is receding and its binding energy with respect to the system is negative [see @1983ApJ...268..319H].

[^4]: The dissolution time was determined by fitting an exponential plus a linear model to $\delta$ as a function of time. At $t\simeq 63.4$ the linear behavior overtakes the exponential behavior of $\delta$. The corresponding time scale then indicates when the system’s diverging characteristics are no longer exponential. For other escape criteria the dissolution time may vary.\[Footnote:Dissolution\]
