On the reliability of N-body simulations

Abstract

The general consensus in the N-body community is that statistical results of an ensemble of collisional N-body simulations are accurate, even though individual simulations are not. A way to test this hypothesis is to make a direct comparison of an ensemble of solutions obtained by conventional methods with an ensemble of true solutions. In order to make this possible, we wrote an N-body code called Brutus, that uses arbitrary-precision arithmetic. In combination with the Bulirsch-Stoer method, Brutus is able to obtain converged solutions, which are true up to a specified number of digits.We perform simulations of democratic 3-body systems, where after a sequence of resonances and ejections, a final configuration is reached consisting of a permanent binary and an escaping star. We do this with conventional double-precision methods, and with Brutus; both have the same set of initial conditions and initial realisations. The ensemble of solutions from the conventional simulations is compared directly to that of the converged simulations, both as an ensemble and on an individual basis to determine the distribution of the errors.We find that on average at least half of the conventional simulations diverge from the converged solution, such that the two solutions are microscopically incomparable. For the solutions which have not diverged significantly, we observe that if the integrator has a bias in energy and angular momentum, this propagates to a bias in the statistical properties of the binaries. In the case when the conventional solution has diverged onto an entirely different trajectory in phase-space, we find that the errors are centred around zero and symmetric; the error due to divergence is unbiased, as long as the time-step parameter, η ≤ 2 -5 and when simulations which violate energy conservation by more than 10% are excluded. For resonant 3-body interactions, we conclude that the statistical results of an ensemble of conventional solutions are indeed accurate.

Dynamical formation and scattering of hierarchical triples: cross-sections, Kozai–Lidov oscillations, and collisions

“…In an individual scattering event, the semi-major axis and eccentricity of all component binaries (both binaries of the triple and, in the case of triple-binary scattering, that of the interloping binary) in the system are fixed, but the argument of pericenter 2 Portegies Zwart & Boekholt (2014) have argued that typical energy conservation standards in N -body studies are extremely conservative and that maintaining energy conservation to better than one part in 10 is sufficient to preserve the statistical properties of dynamical experiments. These results have been confirmed by Boekholt & Portegies Zwart (2015). Since very few of our experiments violate our energy conservation standard, this finding does not change our results, but does imply that we may have slightly overestimated our systematic uncertainties.…”

A Quantification of the Butterfly Effect in Cosmological Simulations and Implications for Galaxy Scaling Relations

“…Since any specific order of summation is arbitrary, no emergent sequence of round-off errors (and hence evolution of the system) is more correct than any other (whether and in what sense the ensemble of solutions to the system represents the true physical solution is a different question, see e.g. Boekholt & Portegies Zwart 2015;Portegies Zwart & Boekholt 2018). This is true even with the simplest set of physics, namely in pure N-body simulations, as well as in pure hydrodynamical simulations, let alone in a combined gravity and hydrodynamics case.…”

The Approximation Tower in Computational Science: Why Testing Scientific Software Is Difficult

“…This is typically done by comparing the orbits obtained for different parameter values and aiming for convergence within the desired accuracy. 5 testing At this point, you're entitled to believe that we've arrived at the top of the approximation tower, but there are still two more floors to go: program transformations. The first level contains transformations performed manually by a software developer, usually for performance.…”