Issue 
A&A
Volume 659, March 2022



Article Number  A86  
Number of page(s)  28  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202141789  
Published online  10 March 2022 
Chaos in selfgravitating manybody systems
Lyapunov time dependence of N and the influence of general relativity
^{1}
Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
email: spz@strw.leidenuniv.nl
^{2}
Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, Parks Road, Oxford OX1 3PU, UK
email: tjarda.boekholt@physics.ox.ac.uk
^{3}
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
email: epor@stsci.edu
^{4}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{5}
Department of Physics, Drexel University, Philadelphia, PA 19104, USA
Received:
14
July
2021
Accepted:
24
September
2021
In selfgravitating Nbody systems, small perturbations introduced at the start, or infinitesimal errors that are produced by the numerical integrator or are due to limited precision in the computer, grow exponentially with time. For Newton’s gravity, we confirm earlier results that for relatively homogeneous systems, this rate of growth per crossing time increases with N up to N ∼ 30, but that for larger systems, the growth rate has a weaker scaling with N. For concentrated systems, however, the rate of exponential growth continues to scale with N. In relativistic selfgravitating systems, the rate of growth is almost independent of N. This effect, however, is only noticeable when the system’s mean velocity approaches the speed of light to within three orders of magnitude. The chaotic behavior of systems with more than a dozen bodies for the usually adopted approximation of only solving the pairwise interactions in the EinsteinInfeldHoffmann equation of motion is qualitatively different than when the interaction terms (or cross terms) are taken into account. This result provides a strong motivation for followup studies on the microscopic effect of general relativity on orbital chaos, and on the influence of higherorder crossterms in the Taylorseries expansion of the EinsteinInfeldHoffmann equations of motion.
Key words: gravitation / relativistic processes / celestial mechanics / chaos / stars: kinematics and dynamics / methods: numerical
© ESO 2022
1. Introduction
Soon after the first gravitational Nbody problems were computed (von Hoerner 1960; Aarseth & Hoyle 1964; van Albada 1968), Miller (1964) questioned the validity of such simulations. The nature of his concern was based on the intrinsically chaotic behavior of Newton’s law of gravity. Errors during integration are introduced by the limited precision of the computer, together with the limited accuracy of the numerical integration scheme. The exponential growth of both sources of errors then contributes to the lack of reproducibility in Nbody simulations (Boekholt et al. 2020).
In an attempt to acquire a converged solution to the 25body problem, Hayli (1970) integrated Newton’s equations of motion using various precisions on a 60 bit CDC 6400, a 64 bit IBM 360/365, and an 80 bit HoneywellBull CII 90–80. He found that although identical initial realizations were used, he acquired a different answer for the final positions and velocities of the 25 objects in his calculations, even with the same algorithm and time step. He argued that the difference in precision of the floatingpoint unit was responsible for the lack of reproducibility of the results. Because the IEEE standard for floatingpoint arithmetic (IEEE 754) was only introduced in 1985, the discrepancy identified by Hayli (1970) could also have originated from differences in roundoff in the least significant digit of the various machines. We still encounter this problem in today’s Graphical Processing Units (GPU), which have limited precision (Fortin et al. 2012). These sources of errors, time step, and roundoff make individual solutions to the Nbody problem notoriously unreliable, although statistically, they may still have the correct phasespace distribution characteristics (Boekholt & Portegies Zwart 2015; Hernandez et al. 2020). The common agreement on roundoff among computer manufacturers, compiler designers, and operating systems hides this problem by generating identical output when the same initial realization is run using different hardware and operating systems, so long as the same source code is run on a single core, and compilers are not too different.
It is a common assumption among Nbody practitioners that the microscopic unpredictability and the consequential loss of reproducibility is irrelevant so long as the global phasespace characteristics are still representative of true physics. It remains an article of faith that such a statistical validity holds for the Newtonian Nbody problem (Heggie 1991). Chaos, however, leads to unpredictability due to temporal discretization, roundoff, and uncertainty in the initial realization (Miller 1964). Unpredictability due to chaotic behavior and the consequential loss of reproducibility may then lead to incorrect physical results.
After the pioneering work of Hayli (1970), the problem of characterizing chaos in selfgravitating Nbody systems received little attention until the late 1980s (except possibly in Zadunaisky 1979, in which the focus was on the computer’s precision), when computers became sufficiently powerful to address the problem for larger N (Heggie 1988; Kandrup & Smith 1991, 1992; Kandrup et al. 1992, 1994; Goodman et al. 1993; Gurzadyan & Kocharyan 1994; Fukushige & Makino 1994). Goodman et al. (1993) and Fukushige & Makino (1994) provided excellent overviews of the underlying arguments for this chaotic behavior in few and manybody systems, respectively.
Even today, the chaotic Nbody problem is hard to address adequately using digital computers, and there is still no analytic solution. We address three aspects of the problem here: (1) the veracity of a solution for N = 4–1024, (2) the scaling of the growth of errors to large N ≳ 10^{5}, and (3) the effect of general relativity on the chaotic behavior of Nbody systems. The terminology used in this text is explained in the glossary in Portegies Zwart & Boekholt (2018). Loosely speaking, a reprehensive solution is a solution in which the errors introduced during integration exceed the system size. For a veracious solution, this is not the case. For N = 3, Boekholt & Portegies Zwart (2015) demonstrated that veracious Nbody solutions give statistically indistinguishable results as an ensemble of converged solutions. They called this behavior Nagh Hoch, to signify the importance of consistent statistical ensemble average behavior of veracious solutions to the chaotic selfgravitating Nbody problem. It is not clear if this concept also holds for larger N.
Hernandez et al. (2020) studied Nagh Hoch by running ensembles of reprehensible numerical solutions for planetary systems. They integrated a single 3 ⋅ 10^{−5} mass planet in a circular orbit at a distance of 0.29 (in dimensionless Nbody units, Hénon 1971) from the star, together with a test particle in the same plane at apocenter and with an eccentricity of 0.19 and semimajor axis 0.21. The calculations were conducted for 200 efoling timescales using four different algorithms for solving the equations of motion. They found that for a sufficiently large sample of initial realizations and a tolerably small time step, the results of the various integrators are statistically indistinguishable. We conclude from their simulation results that their adopted reprehensive Nbody algorithms for a system of two planets complies to Nagh Hoch. They argue that a relative integration error ≲0.05 is sufficient to preserve the quality, consistent with Boekholt & Portegies Zwart (2015), who argued that a time step smaller than 2^{−5} is sufficient to preserve ergodicity in the outcome space.
Ideally, one would like to perform large Nbody simulations to a converged solution, but this is unrealistic, even on modern digital computers. We can achieve converged solutions for N up to about 1k particles for several crossing times, leading to a series of veracious solutions. At the moment, however, a longer evolution or a larger number of particles are too costly. However, we demonstrate that the chaotic behavior of these systems is reprehensive and confirm that they show statistically indistinguishable chaotic behavior. We subsequently perform reprehensible Nbody simulations for up to 128k particles.
The interest in scaling to N ≫ 3 is in part motivated by understanding the chaotic nature of galaxies. Dense stellar systems, such as globular clusters [N = 𝒪(10^{6})], are highly chaotic (Parvulesco 1924; Carpintero et al. 1999). Galaxies (N → ∞) are considered collisionless because their relaxation time exceeds the Hubble time (Binney & Tremaine 2008). For sufficiently large N the background potential becomes smooth, and the collisionless assumption becomes increasingly applicable (Muzzio & Mosquera 2004; Muzzio et al. 2009). However, due to the pointparticle granularity of the potential, the microscopic exponential instability remains present in the system (Valluri & Merritt 2000). Even large Nsystems are therefore affected by the chaos in smallN subsystems. It is nonetheless unclear how the microscopic exponential instability propagates to the macroscopic structure of the stellar system as a whole. At which value of N does the system exhibit the transition from chaotic smallN to smooth largeN systems?
Heggie (1988) argued that the dynamics Nbody systems under Newton’s equations of motion are dominated by encounters at an impact parameter of about r/N^{1/2}. Since chaos is driven by encounters, the Lyapunov timescale then has a similar scaling^{1}. In a pioneering study, Goodman et al. (1993) discussed this scaling and found a transition in the chaotic behavior around N ≃ 32. They argued that the Lyapunov timescale is proportional to the dynamical crossing time t_{λ} ∝ γt_{cr}/ln(ln(N)) over all values of N. For large N (≳32), the constant γ is smaller, leading to a weaker dependence on the Lyapunov timescale (see also Kandrup & Smith 1992; Kandrup 1998). Hemsendorf & Merritt (2002) found a similar transition, but argued in favor of scaling t_{λ} ∝ 1/ln(N) over the entire range of N. In this latter study, however, the perturbed particle is evolved in a static background (Hemsendorf & Merritt 2002), and it is not clear if their different scaling resulted from this particular assumption or from the slightly different potential. We therefore extend the range for which the Lyapunov time was determined to N = 128k using reprehensive Nbody solutions, allowing us to test the scaling of the Lyapunov timescale to large N in an actual selfgravitating system. In addition, we chose various initial density profiles: a smooth profile (as in Goodman et al. 1993), a Plummer profile (1911, as in Hemsendorf & Merritt 2002), and a King profile (W_{o} = 12, King 1966, just for fun).
The change in slope near N ≃ 32 pinpoints the transition from chaotic behavior driven by local fewbody relaxation (for N ≲ 32) to farfield manybody relaxation. This transition can globally be understood by comparing the dynamical crossing time t_{cross} with the twobody relaxation time t_{rlx}, which can be approximated by (Spitzer 1971, 1987; Spitzer & Hart 1971b,a)
$$\begin{array}{c}\hfill {t}_{\mathrm{rlx}}\simeq \left(\frac{N}{8ln(N)}\right){t}_{\mathrm{cross}}.\end{array}$$(1)
This relation indicates that for N ≳ 32, the relaxation timescale exceeds the crossing time. It remains unclear, however, which underlying process determines the slope.
Hemsendorf & Merritt (2002) argued in favor of the same scaling for all N and suggested that a galaxy with 10^{12} stars would have a Lyapunov timescale of t_{λ} ≈ t_{cr}/30 or t_{λ} ∼ 8 Myr, whereas Goodman et al. (1993) argued for a Lyapunov timescale that is roughly an order of magnitude larger. Regardless of either of the two scaling relations, the Galaxy would be subject to chaotic motion on a fraction of the crossing timescale, and would therefore not be representable by the collisionless Boltzmann equation (Boltzmann 1872).
For the Solar System, there seems to be a difference as well, this time, in terms of chaotic behavior. In the Newtonian case, resonances between Jupiter and Mercury greatly enhance the eccentricity of the orbit of the latter planet. This may eventually lead to a collision between Mercury and the Sun (Laskar et al. 1992; Milani & Nobili 1992). However, the inclusion of relativistic effect tends to stabilize the system, resulting in the muchreduced probability that Mercury collides with the Sun (Laskar & Gastineau 2009). Based on our results, we tend to agree with the increased stability of the Solar System when considering relativistic mechanics. Interestingly, the global Lyapunov timescale for the Solar System is not dissimilar from the galactic result, being ∼5 Myr (Laskar et al. 1992) to ≳6.8 Myr (Applegate et al. 1986) or even slightly longer (Duncan & Quinn 1993; Batygin et al. 2015).
We argue that ensembles of reprehensive Nbody solutions, from N = 4 to 1024 (1k), give statistically the same chaotic behavior as the converged solutions. This trend holds at least up to N = 1k, for which we acquire converged solutions for 10 Hénon time units (equivalent to four crossing times). It becomes rather unpractical to continue with converged solutions for N ≥ 1k.
Regardless of the largeN behavior of the chaotic selfgravitating systems under Newton’s forces, we also study the relativistic case. Rather than solving Einstein’s field equations directly, we address relativistic dynamics through an expansion to the gravitational field in terms of v/c, which expresses the speed of light (c) in dimensionless Nbody units in terms of the velocity of the particles (v). The zeroth order in this expansion represents Newton’s equations of motion, which does not depend on the velocity. The firstorder postNewtonian term (the 1PN) is proportional to v^{2}/c^{2}, describes the motion of N Schwarzschild black holes, and is known as the EinsteinInfeldHoffmann (EIH) equation (Einstein et al. 1938). Lorentz & Droste (1917) worked out a first generalization for these postNewtonian Nbody equations of motion, but the final formulation was realized by Einstein et al. (1938).
Our study is motivated by the generally adopted view that the radiationinduced dissipation in general relativity quenches the chaotic behavior of the Nbody system. Spyrou (1975, but see also, Wanex 2002; Cornish & Levin 2003; Galaviz 2011; Neilsen et al. 2014) demonstrated that the pairwise postNewtonian expansion to 2.5th order is chaotic in democratic threebody systems. They also studied the gravitationalwave signal for such systems (Gültekin et al. 2006). Chaotic behavior was not demonstrated for relativistic systems with more than three particles.
We adopt the EIH equation and compare the degree of chaos measured in systems from N = 3 to N = 1k with the Newtonian solutions (v/c → 0). Our approach, however, is limited to firstorder postNewtonian correction terms for the EIH equations of motion. The 1PN terms are not dissipative and therefore are not expected to result in less chaotic behavior when compared to Newton’s equations of motion.
The motivation for performing our a study is somewhat academic because it is currently unclear to which degree largeN relativistic systems appear in nature or how frequently they emit gravitational waves in an observable wavelength. Galactic nuclei are probably the most promising places to find multiple supermassive black holes orbited by intermediatemass black holes or other compact objects. Portegies Zwart et al. (2006) argued that the Galactic center would be populated by a steady population of a dozen intermediatemass black holes within a few milliparsec of the supermassive black hole. In addition, there could be many stellarmass black holes and a rich population of neutron stars in the Galactic center (Muno et al. 2004; Rimoldi et al. 2015). These black holes may merge, producing observable gravitationalwave signals (Pretorius 2005; Abbott et al. 2017) and highvelocity recoiling black holes (Campanelli et al. 2007). Similar, but less extreme, situations may be present in the cores of some globular clusters (Banerjee et al. 2010; Banerjee 2021). Banerjee & Kroupa (2011) even speculated that clusters composed only of dark objects might exist, which could be copious sources of gravitational radiation. Recently, Gieles et al. (2021) argued that the globular cluster Palomar 5 may host such a central collection of compact objects and that in the coming 100 Myr, its entire core may be composed of black holes. With the current rather large virial radius, the dynamics in this cluster is not expected to be subject to strong relativistic effects, however.
It would be interesting to determine the gravitationalwave signature of such chaotic systems and maybe even search for them in the data collected by gravitationalwave observatories. An analysis like this was done for hierarchical triples (Galaviz & Brügmann 2011; Meiron et al. 2017; Robson et al. 2018; Lim & Rodriguez 2020; Will 2021), but not for N > 3.
2. Methods
In the gravitational manybody problem, N objects move under the attractive influence and spacetime distortions of each other. We use three independently developed implementations of the forcelaw and integration method. One of these is the code ph4 (see Appendix A, McMillan 2014), which is designed for the problem of N point masses under a Newtonian force law with regular hardware and compilersupported precision (see Sect. 2.1) For chaotic systems, one may desire more control over the precision and accuracy of the integrator because reprehensible Nbody calculations provide insufficient trust due to roundoff and integration errors that grow exponentially with time. On the other hand, veracious calculations are also reprehensible, but are statistically indistinguishable from the converged solutions (Boekholt & Portegies Zwart 2015). We therefore also perform calculations using Brutus (Portegies Zwart & Boekholt 2014), which allows us to control these parameters with a tolerance ϵ and wordlength L_{w} (see Sect. 2.2).
The relativistic calculations are performed up to the first postNewton order using the pairwise approximation and the full EIH equations of motion. The relativistic Nbody code is called Hermite_GRX, and it is described in Sect. 2.3. The equations of motion in ph4 and Hermite_GRX are integrated using the Hermite algorithm (Makino 1991) or an adaptation of it to accommodate the velocity dependence of the acceleration. In Brutus, the equations of motion are solved with a secondorder Verlet (1967) scheme.
We finally use these three methods to study chaos in the largeN limit. Large here means ∼1k for the postNewton case and converged Newton solutions, and up to N = 128k for reprehensible Newtonian solutions. Each of the codes is interfaced as a community code to the Astronomical Multipurpose Software Environment (Portegies Zwart & McMillan 2018). We describe these implementations in Appendix A.3.3.
2.1. Regular Nbody calculations
According to Newton’s laws of motion, the acceleration a_{i} on particle i is given by the sum over all other particles (Newton 1687):
$$\begin{array}{c}\hfill {\mathit{a}}_{i}={\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{n}\frac{G{m}_{j}}{{r}_{\mathit{ij}}^{3}}{\mathit{r}}_{\mathit{ij}}.}\end{array}$$(2)
Here m_{i} is the mass of particle i, and r_{ij} is the relative position vector from particle i to j: r_{ij} ≡ r_{i} − r_{j}. Newton’s constant is G = 𝒪(10^{−10} N kg^{−2} m^{−2}), but in our calculations, G ≡ 1 (Hénon 1971; Heggie et al. 1986).
Integration was performed in 64bits using IEEE 754 Standard for floatingpoint arithmetic under the Linux operating system with kernel version 5.8.048generic. CPU calculations were performed on an 192 core Intel Xeon E78890 v4 workstation running at 2.20 GHz. Parallelization was realized using the Message Passing Interface (version Open MPI 4.0.3) (Portegies Zwart et al. 2008). For simulations with N > 1k we adopted the GPU version of ph4, which uses the Sapporo GPU library (Portegies Zwart et al. 2007; Gaburov et al. 2009) running on Xeon E2176M CPU with Quadro P2000 MaxQ GPU. The GPU has 4 GB GDDR5 RAM, but is not equipped with error correction, potentially leading to nonIEEE compliant errors in the calculations.
2.2. Arbitrarily precise Nbody calculations
The Hermite algorithm is a fourthorder scheme that reaches a relative energy error dE/E ≃ 10^{−15} with a timestep parameter η ≃ 10^{−4} for a distribution of equalmass objects in a virialized homogeneous distribution in space. All our calculations were performed with η = 0.01. Roundoff and timestep errors become important for smaller time steps when integrating for more than ∼10 crossing times, resulting in a systematic growth of the energy error. This growth scales ∝η^{2}. Although small, typically about 1/10^{16} (in relative coordinates), these errors drive the eventual irreproducibility of the simulations through exponential growth. This irreproducibility is undesirable when one is interested in studying chaotic motion. We therefore also performed calculations with Brutus (Boekholt & Portegies Zwart 2015), an Nbody code that allows us to integrate any Nbody system to arbitrary precision. In Brutus we control the different sources of error by adopting the GraggBulirsch–Stoer algorithm (Bulirsch 1964; Gragg 1965). In this algorithm, one performs a single step using a Verlet integrator (Verlet 1967), then repeats that same step using half the step size. Now the relative error between the two solutions can be determined by taking the absolute value of the differences in each coordinate. If this error is smaller than some predetermined tolerance, the result is accepted, and the next step is calculated. Otherwise, the same step is repeated with a quarter time step. This procedure is repeated until the relative error between two subsequent solutions is smaller than the tolerance (see Press et al. 1992).
We controlled the discretization error with arbitraryprecision arithmetic, using the GMP (Granlund & The GMP development team 2012, 2015) and MPFR libraries (Fousse et al. 2007) instead of conventional double precision. This allowed us to control the roundoff error by changing the number of digits, which we express in a wordlength L_{w}. A wordlength L_{w} = 64 bits then corresponds approximately to the usual 16decimal place precision in standard IEEE 754 floatingpoint operations on current regular microprocessors.
In practice, we only specified the tolerance and calculated the word length L_{w} ∈ ℤ with (Boekholt & Portegies Zwart 2015),
$$\begin{array}{c}\hfill {L}_{w}=\phantom{\rule{0.333333em}{0ex}}\text{int}(324{log}_{10}(\u03f5)).\end{array}$$(3)
A converged solution to n decimal places is achieved by iteratively repeating a calculation that started with one selected realization of the initial conditions with lower tolerance for each subsequent calculation. This process is repeated until the first n decimal places of the final phasespace coordinates of two subsequent iterations lead to identical values for the first n digits in the positions and velocities of all particles. When the numerical solution has achieved this state of convergence, it is deemed to be definitive (Portegies Zwart & Boekholt 2018).
Calculation typically started with tolerance ϵ = 10^{−5} (L_{w} = 52) for N ≤ 64, reducing the tolerance by a factor 10^{5} upon subsequent calculations. For high values of N, we started with ϵ = 10^{−20} (L_{w} = 112), reducing the tolerance by a factor 10^{10} upon subsequent calculations. With ϵ = 10^{−40} (L_{w} = 192), all solutions up to N = 1k have converged to n = 3 decimal places (the adopted convergence limit).
Calculation time with Brutus scales ∝N^{2}, but due to the expense of the large mantissa calculations, the offset in computer time is long. In addition, several recalculations may be needed before a converged solution is achieved. In Fig. 1 we present the scaling of the calculations with Brutus (ochre symbols show the mean of the computing time for that particular N). The lines to higher values indicate the upper limit for a single most expensive calculation for each value of N. These still scale roughly proportional to N^{2}, but are often ∼100 times more costly than the average.
Fig. 1. Scaling of the various integration methods as a function of the number of particles (N). The classic Newton integrations scale as ∝N^{2}, as indicated with the solid and dashed dark blue lines. The pairwise first expansion scales similarly, but tends to be slightly slower than the pure Newton expansion. The EinsteinInfeld Hoffmann equations to first order scale ∝N^{3}, making large calculations that include the crossterms unpractical. The scaling presented for Brutus is based on the calculations presented here. The bullet points indicate the mean timescale for acquiring a converged solution, and the line pointed upward ends at the single most expensive calculation in our sample of simulations for that particular N. For ph4, we included the regular implementation as well as the GPUenabled version (to the right), running on an Intel Xeon CPU E5620 operating at 2.40 GHz and NVIDIA G96 (Quadro FX580), running on a generic 64bit Ubuntu Linux kernel 2.6.3532. 
2.3. EinsteinInfeldHoffmann solver
Between 1907 and 1915, Einstein developed general relativity (see Weinstein 2012; Poisson & Will 2014, for an interesting read on the history of this development) and viewed gravity as the result of the curvature of spacetime (Einstein & Grossmann 1914; Einstein 1915a,b; Hilbert 1915). The Einstein field equations dictate the gravitational interaction between particles, but these equations are nonlinear and notoriously hard to solve. Schwarzschild (1916) found the first nontrivial solution to the Einstein field equations: the Schwarzschild metric describes a pointlike particle. A rotating black hole was first described analytically as a solution to the field equations by Kerr (1963). Today, there are rather standard software implementations to solve for general relativistic dynamical systems (Mewes et al. 2018; BabiucHamilton et al. 2019), even including magnetic fields (Mewes et al. 2020). However, simulating multiple black holes in a relativistic context is somewhat expensive in terms of computer time.
A numerically cheaper solution is the EinsteinInfeldHoffmann equations (Will 2014, 2021; Poisson & Will 2014), in which the acceleration of body i, a_{i} is given by
$$\begin{array}{cc}& {\mathit{a}}_{i}={\displaystyle \sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}^{3}}{\mathit{r}}_{\mathit{ij}}}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}+\frac{1}{{c}^{2}}{\displaystyle \sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}^{3}}{\mathit{r}}_{\mathit{ij}}[4\frac{G{m}_{j}}{{r}_{\mathit{ij}}}+5\frac{G{m}_{i}}{{r}_{\mathit{ij}}}+\sum _{k\ne i,j}\frac{G{m}_{k}}{{r}_{\mathit{jk}}}}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}+4{\displaystyle \sum _{k\ne i,j}\frac{G{m}_{k}}{{r}_{\mathit{ik}}}\frac{1}{2}\sum _{k\ne i,j}\frac{G{m}_{k}}{{r}_{\mathit{jk}}^{3}}({\mathit{r}}_{\mathit{ij}}\xb7{\mathit{r}}_{\mathit{jk}})}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{v}_{i}^{2}+4{\mathit{v}}_{i}\xb7{\mathit{v}}_{j}2{v}_{j}^{2}+\frac{3}{2}{({\mathit{v}}_{j}\xb7{\mathit{n}}_{\mathit{ij}})}^{2}]\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\frac{7}{2{c}^{2}}{\displaystyle \sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}}\sum _{k\ne i,j}\frac{G{m}_{k}}{{r}_{\mathit{jk}}^{3}}{\mathit{r}}_{\mathit{jk}}}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}+\frac{1}{{c}^{2}}{\displaystyle \sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}^{3}}{\mathit{r}}_{\mathit{ij}}\xb7(4{\mathit{v}}_{i}3{\mathit{v}}_{j}){\mathit{v}}_{\mathit{ij}}.}\hfill \end{array}$$(4)
Here n_{ij} = r_{ij}/r_{ij}, v_{ij} ≡ v_{i} − v_{j}, and ${\widehat{\mathit{r}}}_{\mathit{ij}}={\mathit{r}}_{\mathit{ij}}/{\mathit{r}}_{\mathit{ij}}$ is the unit vector along r_{ij}. To conserve energy, an addition term has to be introduced that depends on the 1 PN approximation,
$$\begin{array}{cc}& E=\frac{1}{2}{\displaystyle \sum _{i}{m}_{i}({v}_{i}^{2}\sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}})}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}+\frac{1}{{c}^{2}}{\displaystyle \sum _{i}{m}_{i}[\frac{3}{8}{v}_{i}^{4}+\frac{3}{2}{v}_{i}^{2}\sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}}}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}+\frac{1}{2}{\displaystyle \sum _{j\ne i}\sum _{k\ne i}\frac{{G}^{2}{m}_{j}{m}_{k}}{{r}_{\mathit{ij}}{r}_{\mathit{ik}}}}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\frac{1}{4}{\displaystyle \sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}}(7{\mathit{v}}_{i}\xb7{\mathit{v}}_{j}+({\mathit{v}}_{i}\xb7{\widehat{\mathit{r}}}_{\mathit{ij}})({\mathit{v}}_{j}\xb7{\widehat{\mathit{r}}}_{\mathit{ij}}))],}\hfill \end{array}$$(5)
$$\begin{array}{cc}& \mathit{P}={\displaystyle \sum _{i}{m}_{i}{\mathit{v}}_{i}+\frac{1}{2{c}^{2}}\sum _{i}{m}_{i}{\mathit{v}}_{i}({v}_{i}^{2}\sum _{j\ne i}\frac{G{m}_{j}}{{r}_{\mathit{ij}}})}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\frac{G}{2{c}^{2}}{\displaystyle \sum _{i}\sum _{j\ne i}\frac{{m}_{i}{m}_{j}}{{r}_{\mathit{ij}}}({\mathit{v}}_{i}\xb7{\widehat{\mathit{r}}}_{\mathit{ij}}){\widehat{\mathit{r}}}_{\mathit{ij}}.}\hfill \end{array}$$(6)
Equation (4) gives the full EIH equations of motion. The first term in Eq. (4) (zerothorder term in the Taylor expansion) is identical to Eq. (2) and represents Newton’s acceleration. The other terms reflect postNewtonian corrections. Several of them depend on velocity, and the penultimate term contains the accelerations of the other particles, making it expensive to compute.
Equation (4) is first order, but higherorder postNewtonian corrections exist, although only under the assumption of pairwise interactions. The threebody Hamiltonian, and therefore the corresponding equations of motion, are known in closed form to second postNewtonian order 𝒪(v^{4}/c^{4}) (Schäfer 1987; Lousto & Nakano 2008), and the twobody equations of motion up to 3.5 PN order, or 𝒪(v^{7}/c^{7}) (Futamase & Itoh 2007; Itoh 2009).
Due to the summations over pairs of particles in Eq. (4), the motion of one particle due to a second particle depends on the other particles in the system. As a consequence, the EIH equations of motion scale as 𝒪(N^{3}), rather than the usual scaling to 𝒪(N^{2}) for Newton’s case. This scaling is confirmed in Fig. 1.
We implemented the pairwise and the full EIH equations of motion to 1PN order using a fourthorder Hermite predictorcorrector scheme (see Sect. 2.1). We refer to Hermite_GR1P as the pairwise equations of motion, and to HermiteGRX for the EIH solution to 1PN order. To illustrate the working of the various implementations, we present Figs. 2–4 for orbits of the N = 2, N = 4, and N = 16 Newtonian case, 1PN pairwise equations of motion, and for the full EIH equations of motion to first order. For these simulations, we adopted black hole masses of 10^{6} M_{⊙}, in a 1 pc cube. Equivalent to specifying the mass of the system, we can also use the relative speed of light v/c. In Hénon units, in which G = M = 1, Newton’s kinetic energy of a system of N bodies with total mass M is E_{kin} = 0.5Mσ^{2} = 1/4, with a velocity dispersion σ^{2} = 1/2 (Heggie et al. 1986).
Fig. 2. Orbital evolution of two black holes of masses 10^{6} M_{⊙} with an initial separation of 0.819 au integrated for half a day for three integrations, pure Newton expansion, Newton with expansions to first order, and the full firstorder EinsteinInfeldHoffmann equations. The first two are precisely on top of each other, and we plotted the firstorder pairwise solution last. The last two solutions are identical because the crossterms do not lead to deviations from the firstorder expansions. The postNewton orbits are not closed, as in the Newtonian case (green). No separate scaling of v/c is applied here because the system is initialized in physical units. 
Fig. 3. Orbital evolution of four of 10^{6} M_{⊙} black holes with the same integrators as in Fig. 2. 
All the initial conditions in this study are virialized according to Newton’s equations of motion (see Sect. 2.4), and therefore in Hénon units, the scaled velocity $\mathit{v}=1/\sqrt{2}$, which sets the scaling of our Nbody simulations. We specify the relative scaling with respect to mass or size by changing the speed of light, or v/c. In the numerical implementation, this parameter is specified through the parameter ζ, which is the reciprocal of v/c (see Appendix A.3.8). For clarity, in the main paper, we only use v/c as free parameter.
2.4. Initial conditions
Initial conditions were generated in standard IEEE double precision (L_{w} = 64). This introduces a discrepancy with the lowN ( ≤ 64) experiments, for which the initial iteration was performed with L_{w} = 52, but if convergence is already achieved for ϵ = 10^{−10} (L_{w} = 72), roundoff at the 13th mantissa (the limiting precision for L_{w} = 52) cannot have propagated to the first 3 decimal places.
We adopted the initial conditions from Goodman et al. (1993). All objects then have the same mass and are distributed in a unit cube in phase space (position and velocity) using dimensionless Nbody units. After generating random positions and velocities, the system was moved to the centerofmass frame and scaled to virial equilibrium for the Newtonian solution. The simulation runs with a finite speed of light also used the identical initial realizations as for the Newtonian case. The slight deviations from virial equilibrium in the relativistic initial conditions have no effect on our results because we started measuring the phasespace distance at t = 1, and by that time, the system was well virialized. We validated and confirmed this statement by recalculating all simulations for N = 16 with v/c = 0.01, for which we scaled the initial conditions to virial equilibrium for v/c = 0.01 by adapting masses and velocities according to Buchdahl (1964), however. The difference between the Newtonian and the relativistically virialized initial conditions was negligible.
The number of runs performed varied for each code and the number of particles, as listed in Table 1. In the second column, we list the number of runs performed with ph4. For the other runs, with Hermite_GRX and Brutus, the same initial realizations were adopted, but sometimes this was a subset. These calculations were performed up to N = 1k using the same number of runs with the same initial conditions as in Goodman et al. (1993).
Number of simulations performed per implementation and number of particles.
We performed an additional series of simulations using ph4 and Hermite_GRX, but with realizations generated using a Plummer (1911) sphere and a King model (W_{o} = 12, King 1966).
The main reason not to perform larger simulations including the EIH equations is their unfavorite scaling of the computer time with N, which we depict in Fig. 1. We estimate approximately one year of integration for N ∼ 10^{4} with the current CPU implementation.
2.5. Measuring the Lyapunov timescale
Measuring the Lyapunov timescale for the gravitational Nbody system is not trivial. Several methods for deriving this quality have been proposed. One method uses the geodesicdeviation vectortechnique (Weinberg 1972; Nieto et al. 2003) for two nearby orbits with projection operations and with time as an independent variable (Wu & Huang 2003), and the twonearby realizations without projection operations and with time as an independent variable. We adopted the last, which may be more expensive to calculate, but is considerably simpler for large N, and it is least affected by underlying assumptions. This same technique was adopted in Goodman et al. (1993), which means that our analysis at least starts from the same assumptions. Strictly speaking, the Lyapunov timescale is defined properly from some starting point until the system dissolves (Urminsky & Heggie 2009; Mel’nikov et al. 2013). Because this definition is rather unpractical, particularly for large N, we stopped the calculations at 10 Nbody time units (equivalent to Goodman et al. 1993).
The degree of chaos in the simulation was measured using the evolution of the phasespace distance between two almost identical initial realizations (see Sect. 2.4). The second realization was constructed by increasing the Cartesian x coordinate of a randomly selected particle with a value of 10^{−7} (in dimensionless Nbody units). Just to emphasize, this initial displacement is 10 million times shorter than the size of the initial extent of the Nbody system. The perturbed realization is therefore not in strict equilibrium, but deviates from Newton’s equilibrium potential energy by 𝒪(10^{−7}/N^{2}).
We integrated both initial realizations to 10 Nbody time units while saving a snapshot every 0.1 Nbody time units, resulting in 100 snapshots per run. The phasespace distance was determined by taking the difference in position and velocity between the same particles in each snapshot, and summing them,
$$\begin{array}{c}\hfill ln(\delta )=\frac{1}{2}ln[{\displaystyle \sum}{({\mathit{r}}_{b}{\mathit{r}}_{a})}^{2}+{({\mathit{v}}_{b}{\mathit{v}}_{a})}^{2}].\end{array}$$(7)
This leads to a phasespace distance as a function of time. To calculate the Lyapunov exponent, we only used the data from t = 1 to a maximum of either t = 10, or the first moment in which the phasespace distance exceeded 0.1. The choice of starting the Lyapunov timescale measurements at t = 1 guarantees that the system is in virial equilibrium even in the most relativistic cases. We subsequently performed a leastsquares fit to the phasespace distance evolution. The fitted slope (in log space) to this phasespace distance evolution gives the Lyapunov exponent. The Lyapunov timescale t_{λ} is the reciprocal of the Lyapunov exponent.
This procedure is slightly different than what was used in Goodman et al. (1993), who adopted t_{λ} = 9/(ln(δ_{t = 10})−ln(δ_{t = 1})), but results in a better estimate of the global Lyapunov timescale. We stopped our measurement when δ ≥ 0.1 because due to conservation of the phasespace characteristics, the system then grows on a relaxation timescale, rather than on a Lyapunov timescale, and δ saturates when it becomes on the order of unity (Hut & Heggie 2002). For very chaotic systems, the procedure adopted by Goodman et al. (1993) leads to an underestimate of the Lyapunov exponent and therefore to an overestimate of t_{λ}, as is the case for N ≳ 1k King models (see, e.g., in Fig. 13).
3. Results
3.1. Chaos in largeN Newtonian systems
In Fig. 5 we show an example for an 1kbody system, starting with the initial conditions of Goodman et al. (1993). The gray square in the middle represents these initial conditions; particles, according to Goodman et al. (1993), are initialized in a unit cube. One calculation (bullet points) gives the result of the unperturbed solution. The perturbed solution is not shown, but the colors of the particles give the phasespace distance between the final perturbed and unperturbed solutions. The black bullet point toward the top right corner of the gray area identifies the (randomly selected) particle for which the initial xcoordinate was increased by 10^{−7}. The least (red) and most (blue) chaotic particles are represented as lines. The overplotted thin black curves show the orbit of the perturbed solution.
Fig. 5. Distribution of phasespace distances in a cluster with 1024 particles. Units are dimensionless Nbody units (Hénon 1971). The gray shaded region indicates the initial conditions in a virialized unit cube. One particle, indicated with the black bullet point, is displaced by 10^{−7} along the Cartesian x coordinate. The final conditions (at t = 10) of the unperturbed particles are represented with the bullet points. The color and size of the points represents the phasespace distance measured over the duration of the simulation (10 Nbody time units), and ranges over about four orders of magnitude. The majority of objects experience considerable change in their orbits, but some are hardly perturbed. Calculations were performed using Brutus until convergence to 3 decimal places, which requires a tolerance of τ = 10^{−40}. 
Figure 5 illustrates Miller’s (1964) point that a small perturbation in a single object leads to large variations in the final phasespace distribution. In fact, most objects experience a strong variation, whereas only a minority of objects are hardly affected. In Fig. 6 we plot the distribution of phasespace distances (log_{10}(δ)) for the calculation of Fig. 5.
Fig. 6. Distribution of phasespace distances for individual particles δ_{i} in the simulations with N = 4 (red) and those with N = 1024 (blue) after integrating for t = 10 Nbody time units. The data for N = 4 are the result of 200 runs. For N = 1024, we adopted the run used in Fig. 5. Calculations were performed using Brutus until the solution was converged. 
The degree to which particles are affected by a small initial perturbation depends on the number of particles in the system. This is illustrated in Fig. 6, where we show the distribution of phasespace distance between a perturbed and an unperturbed solution for the same simulation as in Fig. 5 using 1k particles, and compare this distribution with 200 simulations of N = 4. The smallN systems (red histogram) exhibits a much weaker response to a perturbing particle than the largeN systems (blue); fewbody systems are less chaotic than largeN systems (at least for this selection of initial conditions, and under Newton’s forces).
The different behavior for smallN systems compared to largeN systems motivated Goodman et al. (1993) and Hemsendorf & Merritt (2002) to conduct their analysis and study the source of chaos in small versus large Nbody systems. In Fig. 7 we show the results of Goodman et al. (1993) and compare them with converged solutions using Brutus up to N = 1k and reprehensible solutions using ph4 for up to N = 128k. The consistency between the results obtained by Goodman et al. (1993) (red), Brutus (blue), and ph4 (ochre) gives us confidence in the validity of the nonconverged (reprehensible) Nbody solutions by Goodman et al. (1993) and using the regular Hermite algorithm implemented in ph4 without going through the elaborate process of reaching a converged solution for N > 1k.
Fig. 7. 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 Goodman et al. (1993). The different symbols and colors represent different calculations (see legend). The vertical bars, plotted for Newton’s Hermite only, show the rootmeansquare 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 scaling we observe in Fig. 7 is consistent with that found by Goodman et al. (1993) over the entire range they explored, from N = 4 to N = 512. We therefore conclude that 1) reprehensible simulations are adequate for studying shorttimescale Lyapunov exponent measurements for relatively homogeneous systems, and 2) the Lyapunov timescale t_{λ} ∝ γt_{cr}/ln(ln(N)), with γ = −1.39 for N ≲ 32 and shallower with γ = −0.498 for N ≳ 32.
3.2. Chaos in largeN relativistic systems
To study the degree of chaos in the relativistic regime, we used the EIH equations of motion with initial realizations (masses, positions, and velocities) identical to those used in the Newtonian simulations. Therefore the latter initial conditions are in virial equilibrium for the Newtonian case, but not for the highly relativistic cases. Whereas the Newtonian Nbody initial realizations and calculations were scale free, we have to relax this assumption when introducing the speed of light. Since our calculations scale with mass M, size R, or velocity v, we have the option to qualify the scaling by just changing v/c. Here v/c → 0 corresponds to the Newtonian case (because in that case, c → ∞); higher values of v/c indicate a more relativistic regime.
We started by confirming that for v/c → 0, we reproduce the results from Fig. 7. As long as v/c ≲ 10^{−4}, the results of the relativistic integration can hardly be distinguished from the Newtonian case (see also Fig. 9).
This is further illustrated in Fig. 8, where we present two histograms for N = 4 (red) and N = 64 (blue) for v/c = 0.010 (top panel) and v/c = 0.002 (bottom). In the top panel, the distribution in phasespace distance for both N = 4 and N = 64 are comparable, although the dispersion for N = 64 is somewhat larger. When we reduce the speed of light (expressed as the parameter v/c), both distributions move toward higher values of δ. For v/c = 0.002 (bottom panel), the systems, though still somewhat relativistic, have mean and median values that already approach the Newtonian values. We recall that with this adopted scaling, the system would correspond to a cluster with a total mass of ∼1 M_{⊙} at a size scale of ∼436 au. Such a cluster of stars appears insufficiently relativistic for the degree of chaos in the equations of motion to be affected noticeably.
Fig. 8. Distribution of phasespace distances for individual particles in the simulations with N = 4 (red) and for N = 64 (blue) after integrating for t = 10 Nbody time units. The data for 200 runs were used for each histogram. Calculations were performed using Hermite_GRX using v/c = 0.01 (top panel) and for v/c = 0.002 (bottom panel). 
However, for v/c = 0.0005, which corresponds to a size scale four times larger, or equivalently, to a system four times more massive, the measurements in the Lyapunov timescale start to deviate from the Newtonian case. When we further decrease v/c ≲ 10^{−4}, the mean and median values of both distributions are statistically indistinguishable from the Newtonian case. The dispersion, particularly noticeable for N = 64, remains skewed to low values of δ.
In Fig. 9 we present the median Lyapunov timescale for N = 4 and N = 64 as functions of v/c. For the asymptotic Newtonian case, v/c → 0, the Lyapunov timescale converges to the median for the Newtonian case. The dispersion in the distribution in the relativistic case remains somewhat larger, however, even for v/c → 0, for which ⟨t_{Ly}⟩ = 1.61 ± 2.36, compared to the Newtonian case, for which ⟨t_{Ly}⟩ = 1.71 ± 1.68.
Fig. 9. Lyapunov timescale as a function of v/c for N = 4 (green) and N = 64 (blue). The Newtonian case (run with ph4) is presented as arrows in orange. The vertical bars, only for the green points, indicate the dispersion in the simulation results. The short horizontal dotted green line indicates the lowest value for the Lyapunov timescale reached for v/c = 10^{−3} for N = 4. 
We suspect that these small systematic effects (which are not statistically significant) could result from a few encounters sufficiently close to be affected by general relativity. For the extreme relativistic case, v/c ≳ 0.001, the Lyapunov timescale rises quickly to ⟨t_{Ly}⟩ = 3.89 ± 0.64. We present the results for v/c ≳ 0.02 (two green points to the right), even though these are beyond the regime where the 1PN Taylorseries expansion to the EIH equations of motion is valid (see Sect. 4.2). We therefore limit further analysis to v/c ≲ 0.010.
For N = 4, we observed a minimum in the Lyapunov timescale for v/c ∼ 10^{−3} (signified by the horizontal dotted green line in Fig. 9). The change in behavior for less and more relativistic systems might be interpreted as a signature that the adopted Taylor expansion starts to break down, but for v/c ≳ 10^{−3}, the postNewtonian terms should still be valid. We expect the lowN configurations to break down earlier when they are evolved with time because they are more relaxation dominated than the largeN systems.
With a typical distribution in velocities matching a truncated Maxwellian, a small fraction (∼2.9%) of the systems has a velocity that exceeds the mean dispersion by factor of 3. Even for a value of v/c ≳ 0.05, the fraction of stars with a velocity v ≳ 0.3c is smaller than 1/10^{7}, and it is unlikely that when integrating over only 10 Hénon units, the postNewtonian Taylorseries expansion breaks down.
In Fig. 10 we present measurements for the Lyapunov timescale for the postNewtonian equations of motion. For N ≲ 4, the EIH equations of motion as well as the pairwise 1PN terms show similar chaotic behavior in the sense that the relativistic system is less sensitive to initial perturbations than the Newtonian case. For N ≳ 20, the pairwise 1PN terms result in smaller Lyapunov timescales compared to Newton’s equations of motion, whereas the EIH equations of motion continue to result in a rather large Lyapunov timescale compared to the Newtonian case. The Lyapunov timescales for the EIH equations of motion are roughly twice as large as for the Newtonian case for v/c = 0.005 and roughly four times as large for v/c = 0.010.
Fig. 10. Lyapunov timescale as a function of N for Goodman et al. (1993) (red bullets) compared to the various relativistic solutions. In blue we present the solutions using 1PN pairwise terms with a scaling to the speed of light of v/c = 0.010 and twice this value. For the EIH equations of motion we show the case for v/c = 0.010 in green, twice and four times this value. The top blue curve fits t_{λ} ≃ 0.255 + 13.43e^{−0.371N} for v/c = 0.005, and t_{λ} ≃ 0.051 + 5.50e^{−0.115N} for v/c = 0.010. 
The effect of N on the degree of chaos in the equations of motion is further illustrated in Fig. 11. Here we show for N = 4, N = 16, and for N = 64 the ratio of the pairwise 1PN solution as a function of the full EIH equations of motion. The relativistic calculations were performed with v/c = 0.010.
Fig. 11. Estimate for the Lyapunov timescale for the pairwise 1PN terms as a function of the full EIH equations of motion for three choices of N, 4 (dark blue), 16 (aquamarine), and 64 (green). Both, the pairwise 1PN calculation and the one including the cross terms, are represented as fraction of the Newtonian solution. All relativistic simulations adopt v/c = 0.010. 
For N = 4, the Newtonian case shows a smaller average Lyapunov timescale than the two relativistic cases. This was also confirmed in the threebody simulations by Boekholt et al. (2021), who found no significant difference in the chaotic behavior of Newtonian versus relativistic systems. When N increases, the distribution of the Lyapunov timescale for the EIH equations of motion continues to be large compared to the Newtonian case, but the value for the pairwise 1PN terms tends to drop to less than 1/10th of the Newtonian solution. We conclude that if one is interested in the dynamical behavior of N ≳ 4 black holes, the pairwise 1PN terms do not reliably represent the chaotic behavior expected for such a relativistic system. The pairwise 1PN terms address pairs of compact objects, whereas the EIH equations of motion should give a more reliable representation for relativistic systems of N ≥ 3. We expect this difference to persist for the higherorder Taylor expansion terms.
4. Discussion
4.1. Consequences of the 1PN terms
Phasespace volume is preserved in a solution to a conservative Hamiltonian system, but shrinks in a dissipative system (Shivamoggi 2014), as is the case in general relativity. The contraction of the phasespace volume gives rise to an attractor, and as a consequence, can have bounded trajectories (Bergé et al. 1987). It is not a priori clear, however, how dissipation in an otherwise conservative system affects chaotic motion (Lakshmanan & Rajasekar 2003). We find that in the conservative 1PN regime, the chaotic behavior of Nbody systems is already affected. In particular, for N ≳ 10, simulations that only include the 1PN pairwise terms behave differently than when the 1PN crossterms for the EIH equations of motion are incorporated into the simulations. The behavior of relativistic systems with v/c ≳ 0.005 and for N ≳ 10 is considerably less chaotic than their lessrelativistic v/c ≲ 10^{−3} and Newtonian counterparts.
4.2. Validity of the postNewtonian terms
In this study, we rely on the postNewtonian expansion of the EIH equations of motion. Ideally, we would have adopted full general relativity in our Nbody calculations, but this is somewhat beyond the scope of our study and is numerically challenging.
In an attempt to quantify the validity of the postNewtonian expansion adopted here, we compared the apsidal motion of the orbitaveraged evolution for a twobody system with total mass M, semimajor axis a, and eccentricity e. We write the relative velocity in a circular orbit in terms of the gravitational radius, r_{g} = GM/c^{2}, and the speed of light as
$$\begin{array}{c}\hfill c\simeq \sqrt{GM/(10{r}_{g})}.\end{array}$$(8)
The Taylorseries expansion then starts to break down for $\mathit{v}\equiv c/\sqrt{10}\sim 0.3c$ (Will 2011). During our Nbody calculations, we kept track of this velocity to ensure that the Taylor series expansion in our calculations remained reliable. However, this safety check does not guarantee that our results are not affected, particularly for high values of v/c.
In the regime in which the Taylorseries expansion of the EIH equations of motion breaks down, the 1PN terms adopted here are insufficient to capture the correct physical behavior. In this case, the 2PN terms become essential for the correct physical interpretation of the numerical results. By definition, the 2PN terms are smaller than the 1PN terms because the former scale as v/c and the latter as v^{2}/c^{2}. On the other hand, both terms approach each other for more relativistic systems, with v → c. It is somewhat tricky to give an absolute measure when the 2PN terms should be used in addition to the 1PN terms. In a general Nbody problem, stars may approach each other at a short distance with relatively high velocities with respect to c. When such an encounter is a onetime event in the nondissipative limit, the lack of precision in the PN terms is not expected to make a great difference in the eventual results (Will 2011). Reprehensive simulations are therefore sufficient to derive the largest global positive Lyapunov exponent for the system.
In an attempt to quantify the relative importance of 1PN with respect to the 2PN terms, we compared the apsidal motion of the orbitaveraged evolution for a twobody system with a total mass M, semimajor axis a, and eccentricity e. We write (Iorio 2020)
$$\begin{array}{c}\hfill {\dot{\omega}}_{2PN}/{\dot{\omega}}_{1PN}\simeq ({r}_{g}/a)(1/12)(28{e}^{2})/(1{e}^{2}).\end{array}$$(9)
In Fig. 12 we show as a function of v/c the relative drift in the apsidal motion for the 1PN and 2PN terms for two bodies in a circular orbit at 100r_{g},
$$\begin{array}{c}\hfill v/c=(1/c)\sqrt{GM/a}\sqrt{(1+e)/(1e)}.\end{array}$$(10)
Fig. 12. Relative importance in the apsidal motion of the 1PN terms in comparison to the 2PN terms. 
The boundary at which the postNewtonian expansion is no longer reliable is indicated by the dashed vertical line, near v/c ≃ 0.4, which happens when the two objects approach within 10r_{g}. If we compare this boundary to the range in v/c in Fig. 9, we find that all lie below the boundary, and we therefore argue that the increase in the Lyapunov timescale toward the relativistic regime is physical and not a numerical artifact.
4.3. Other initial density profiles
The initial conditions adopted in Goodman et al. (1993) have a homogeneous phasespace distribution within the adopted limits, and they do not represent any observed stellar systems (Portegies Zwart et al. 2010). Clusters of stars are better represented with a Plummer (1911) distribution or a King (1966) model. For this reason, we also performed a series of calculations with these distributions. One series of Newtonian calculations used Plummer models and King models with dimensionless depth of W_{0} = 12 (which is rather concentrated), and one set of calculations used the EIH equations of motion for the King model case (with identical initial realizations).
In Fig. 13 we present the Lyapunov timescales for these simulations as functions of N. The Newtonian Plummer case shows a slightly smaller Lyapunov timescale than the homogeneous distribution used in Goodman et al. (1993). The Plummer distribution is consistent with the initial conditions adopted by Hemsendorf & Merritt (2002), and their results are consistent with our results for the Plummer sphere, see Fig. 13 (black points).
Fig. 13. Estimated Lyapunov timescale for the Newtonian case with particles distributed in a Plummer sphere (blue) and a King model (yellow). In addition, we show results for the King model, but for the EIH equations of motion with v/c = 0.010 (green). Here the xaxis is in ln(ln(N)). 
The Newtonian King model with W_{0} = 12 is considerably more chaotic than the homogeneous initial realization adopted in Goodman et al. (1993), at least given that for N ≳ 10^{3} we were unable to measure a reliable Lyapunov timescale because the phasespace distance grew beyond δ = 0.1 within 1 Nbody time unit. King models with a central potential depth of W_{0} = 12 turn out to be considerably more chaotic than Plummer models (for N ≳ 40), while Plummer models are expected to behave more chaotically than a homogeneous distribution of particles. This is not a complete surprise because the choice of W_{0} = 12 places the model near core collapse (at least in its density profile). The growth of an initial phasespace distance between two subsequent calculations with almost identical initial relations is then dominated by fewbody interactions in the core. One could argue that the entire chaotic behavior of the star cluster is driven by fewbody interactions in the cluster center. Since some complex threebody interactions are fundamentally unpredictable (Boekholt et al. 2020), the dynamical evolution of the entire cluster will be unpredictable.
The extreme relativistic case, with v/c = 0.01, shows a similar characteristic again as the homogeneous initial realizations, but a rather different scaling when adopting the King model. In the Newtonian case, King models tend to have considerably smaller Lyapunov timescales, but when extremely relativistic, they tend to be more regular than the Newtonian case Fig. 10.
5. Conclusions
We have numerically analyzed the rate at which neighboring solutions of the equations of motion for N selfgravitating bodies diverge in the Newtonian regime, but also with the 1postNewtonian expansion terms for the pairwise approximation and the EinsteinInfeldHoffmann equations of motion. Our results can be interpreted as the rate of growth of the error in an Nbody solution, caused by uncertainties in the initial conditions or errors produced numerically during integration.
Our Newtonian simulations were repeated with higher precision and accuracy until a converged solution was achieved. Due to computer limitations, this was performed for N up to 1024 particles. For large particle numbers (and for the relativistic simulations), we used reprehensible Nbody solutions, and we confirmed them to be veracious for N up to 1k particles.
The motivation to study the growth of errors stems from our desire to understand the role of chaos in these systems. The macroscopic distribution of material in the Galaxy may not be affected by microscopic chaos. But the Galaxy is built up of small subsystems of stars, each of which exhibits chaotic behavior, and the range of the Newtonian force law causes chaos in these microscopic systems to propagate to the Galaxy at large. Chaos in the Galaxy is then governed by the chaos in small N subsystems and not by the global slow (on timescales longer than the dynamical timescale) variations of the orbits of stars in a smooth potential.
We confirm the earlier result of Goodman et al. (1993) and Hemsendorf & Merritt (2002) that the divergence in Nbody systems grows exponentially, with an efolding timescale on the order of the crossing time. Our results agree with the t_{λ} ∝ t_{cross}/ln(ln(N)) scaling of Goodman et al. (1993) and are inconsistent with a t_{cross}/ln(N) scaling.
Our conclusions are listed below.

For a homogeneous distribution of equalmass particles in virial equilibrium, the efolding timescale for the growth of an initial perturbation in an Nbody system, the socalled Lyapunov timescale, scales for small systems of N ≲ 32 as t_{λ}/t_{cross} = (0.88 ± 0.12)−(1.39 ± 0.13)ln(ln(N)). For larger systems of N ≳ 32, the Lyapunov timescale scales as t_{λ}/t_{cross} = ( − 0.094 ± 0.129)−(0.498 ± 0.066)ln(ln(N)).

For an initial Plummer distribution, the efolding timescale is smaller than the homogeneous initial realizations by about a factor of 5 but preserves the same trend, or for N ≳ 32, it fits t_{λ}/t_{cross} = ( − 0.475 ± 0.018)−(0.528 ± 0.010)ln(ln(N)).

For more concentrated models, such as a King model with W_{0} = 12, the scaling is comparable to the slope observed in the homogeneous unitcube or the Plummer distribution, with t_{λ}/t_{cross} = (1.346 ± 0.110)−(2.212 ± 0.107)ln(ln(N)), but extending somewhat farther, to N ≃ 64. For larger N, the slope is much steeper, t_{λ}/t_{cross} = (4.970 ± 2.03)−(4.813 ± 1.147)ln(ln(N)), indicating that these systems are chaotic for large N on a timescale smaller than a crossing time.

If a small perturbation is introduced into a single particle of a large Nbody system, all particles are affected within a few crossing times.

For selfgravitating systems with v/c ≲ 10^{−3}, the phasespace mixing of relativistic Nbody systems is indistinguishable from the Newtonian case. This limit is already reached for a total of ten black holes of 10 M_{⊙} confined to a spherical volume of radius 10^{−3} pc.

For highly relativistic systems, v/c ≳ 0.002, the EIH equations of motion to 1PN are considerably less chaotic than their Newtonian counterpart over all values of N. The Lyapunov timescale scales with t_{λ}/t_{cross} = 6.63 ± 1.68 − ln(3.72±2.04ln(N)).

For small N (≲10), the pairwise 1PN terms give similar phasespace mixing to the EIH equations of motion.

For N > 4, the pairwise 1PN corrected equations of motion become considerably more sensitive to perturbations in the initial conditions compared to the EIH equations of motion. The former show considerably shorter Lyapunov timescales compared to their Newtonian counterparts, whereas the latter has even longer Lypaunov timescales.

We conclude that the Galaxy is intrinsically chaotic on a very short timescale because of the chaotic behavior of microscopic fewbody interactions in the centers of star clusters. The chaotic behavior of these smallN systems propagates on a local crossing timescale to the entire star cluster, affecting the orbits of neighboring stars and clusters, and eventually, the entire Galaxy.

The pairwise terms for N > 3 give a different dynamical behavior for relativistic Nbody systems compared to the full EIH equations of motion.
Considering the effect of the full EIH equations of motion on a relativistic cluster of compact objects, and the potential consequences for observations with laser interferometric gravitational wave observatories, we look forward to implementing and study the effect of higherorder crossterms in general relativistic Nbody simulations. We do realize, however, that these calculations do not have the most favorable scaling of the computer time with respect to N.
Acknowledgments
It is a pleasure to thank Clifford Will, ln(a) Sellentin, and Arend Moerman for discussions. This project was supported by funds from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement No 638435 (GalNUC), and by the National Science Foundation in the U.S. under grant AST1814772. This work was performed using resources provided by the Academic Leiden Interdisciplinary Cluster Environment (ALICE), and using LGMII (NWO grant # 621.016.701). Energy consumption of this calculation: The calculations using Brutus are elaborate and took about 10^{7} CPU seconds. The other two sets of calculations are comparable in expense, totaling about a year of single CPU usage. Using the tool http://greenalgorithms.org/, we calculated our energy consumption to be about 3.32 MWh. At Dutch electricity rates, this would produce about 1.8 kiloton CO_{2}, but since the computers used are powered by either Dutch wind or Norwegian hydroelectric power (through certificates) the net CO2 emission should be negligible. Public data: The source code, input files, simulation data, and data processing scripts for this manuscript are available at figshare under DOI 10.6084/m9.figshare.xxxx. Software used for this study: This work would have been impossible without the following public opensource packages and libraries: Python (van Rossum 1995), matplotlib (Hunter 2007), numpy (Oliphant 2006), MPI (Gropp et al. 1996; Gropp 2002), and AMUSE (Portegies Zwart et al. 2018, available for download at https://www.amusecode.org/). Sapporo GPU library (Portegies Zwart et al. 2007; Gaburov et al. 2009), MPFR library (Fousse et al. 2007) of the GMP library (Granlund & The GMP development team 2012). A (python notebook) tutorial for students is available at https://github.com/amusecode/Tutorial. All the Nbody codes used in this study are available in the AMUSE repository at amusecode.org.
References
 Aarseth, S. J., & Hoyle, F. 1964, Astrophys. Norvegica, 9, 313 [NASA ADS] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12 [Google Scholar]
 Antonini, F., Murray, N., & Mikkola, S. 2014, ApJ, 781, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Antognini, J. M. O., & Thompson, T. A. 2016, MNRAS, 456, 4219 [NASA ADS] [CrossRef] [Google Scholar]
 Applegate, J. H., Douglas, M. R., Gursel, Y., Sussman, G. J., & Wisdom, J. 1986, AJ, 92, 176 [NASA ADS] [CrossRef] [Google Scholar]
 BabiucHamilton, M., Brandt, S. R., Diener, P., et al. 2019, The Einstein Toolkit, to find out more, visit http://einsteintoolkit.org [Google Scholar]
 Banerjee, S. 2021, MNRAS, 503, 3371 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., & Kroupa, P. 2011, ApJ, 741, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371 [Google Scholar]
 Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120 [CrossRef] [Google Scholar]
 Bédorf, J., Gaburov, E., & Portegies Zwart, S. 2015, Comput. Astrophys. Cosmol., 2, 8 [CrossRef] [Google Scholar]
 Bergé, P., Pomeau, Y., & Vidal, C. 1987, Order Within Chaos (Wiley) [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
 Boekholt, T., & Portegies Zwart, S. 2015, Comput. Astrophys. Cosmol., 2, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Boekholt, T. C. N., Portegies Zwart, S. F., & Valtonen, M. 2020, MNRAS, 493, 3932 [NASA ADS] [CrossRef] [Google Scholar]
 Boekholt, T. C. N., Moerman, A., & Portegies Zwart, S. F. 2021, Phys. Rev. D, 104, 083020 [CrossRef] [Google Scholar]
 Boltzmann, L. 1872, Wiener Berichte, 275 [Google Scholar]
 Buchdahl, H. A. 1964, ApJ, 140, 1512 [NASA ADS] [CrossRef] [Google Scholar]
 Bulirsch, R., & Stoer, J. 1964, Numer. Math., 6, 413 [CrossRef] [Google Scholar]
 Campanelli, M., Lousto, C. O., Zlochower, Y., & Merritt, D. 2007, Phys. Rev. Lett., 98, 231102 [Google Scholar]
 Carpintero, D. D., Muzzio, J. C., & Wachlin, F. C. 1999, Celest. Mech. Dyn. Astron., 73, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Cornish, N. J., & Levin, J. 2003, Phys. Rev. D, 68, 024004 [NASA ADS] [CrossRef] [Google Scholar]
 de Elía, G. C., Zanardi, M., Dugaro, A., & Naoz, S. 2019, A&A, 627, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Lagrange, J. L. 1772, Chapitre II: Essai sur le Problème des Trois Corps [Google Scholar]
 Duncan, M. J., & Quinn, T. 1993, ARA&A, 31, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Efimov, S. S., & Sidorenko, V. V. 2020, Cosmic Res., 58, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Einstein, A. 1915a, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), 778 [Google Scholar]
 Einstein, A. 1915b, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), 799 [Google Scholar]
 Einstein, A., & Grossmann, M. 1914, Z. Math. Phys., 63, 215 [Google Scholar]
 Einstein, A., Infeld, L., & Hoffmann, B. 1938, Ann. Math., 65 [NASA ADS] [CrossRef] [Google Scholar]
 Euler, L. 1760, Nov. Comm. Acad. Imp. Petropolitanae, 10, 207 [Google Scholar]
 Fortin, P., Gouicem, M., & Graillat, S. 2012, in 20th Euromicro International Conference on Parallel, Distributed and Networkbased Processing, 407 [CrossRef] [Google Scholar]
 Fousse, L., Hanrot, G., Lefèvre, V., Pélissier, P., & Zimmermann, P. 2007, ACM Trans. Math. Softw., 33, 13es [CrossRef] [Google Scholar]
 Fukushige, T., & Makino, J. 1994, ApJ, 436, L111 [NASA ADS] [CrossRef] [Google Scholar]
 Funato, Y., Hut, P., McMillan, S., & Makino, J. 1996, ApJ, 112, 1697 [CrossRef] [Google Scholar]
 Futamase, T., & Itoh, Y. 2007, Liv. Rev. Rel., 10, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Gaburov, E., Harfst, S., & Portegies Zwart, S. 2009, New Astron., 14, 630 [Google Scholar]
 Galaviz, P. 2011, Phys. Rev. D, 84, 104038 [NASA ADS] [CrossRef] [Google Scholar]
 Galaviz, P., & Brügmann, B. 2011, Phys. Rev. D, 83, 084013 [NASA ADS] [CrossRef] [Google Scholar]
 Georgakarakos, N. 2008, Celest. Mech. Dyn. Astron., 100, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Gieles, M., Erkal, D., Antonini, F., Balbinot, E., & Peñarrubia, J. 2021, Nat. Astron., 5, 957 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., Heggie, D. C., & Hut, P. 1993, ApJ, 415, 715 [Google Scholar]
 Gragg, W. B. 1965, SIAM J. Numer. Anal., 2, 384 [NASA ADS] [Google Scholar]
 Granlund, T., & The GMP development team 2012, GNU MP: The GNU Multiple Precision Arithmetic Library, 5th edn. http://gmplib.org/ [Google Scholar]
 Granlund, T., & The GMP development team 2015, GNU MP 6.0 Multiple Precision Arithmetic Library (United Kingdom: Samurai Media Limited) [Google Scholar]
 Gropp, W. 2002, in MPICH2: A New Start for MPI Implementations, eds. D. Kranzlmüller, J. Volkert, P. Kacsuk, & J. Dongarra (Berlin, Heidelberg: Springer, Berlin Heidelberg), 7 [Google Scholar]
 Gropp, W., Lusk, E., Doss, N., & Skjellum, A. 1996, Parallel Comput., 22, 789 [Google Scholar]
 Gültekin, K., Miller, M. C., & Hamilton, D. P. 2006, ApJ, 640, 156 [CrossRef] [Google Scholar]
 Gurzadyan, V. G., & Kocharyan, A. A. 1994, J. Phys. Math. Gen., 27, 2879 [NASA ADS] [CrossRef] [Google Scholar]
 Hamers, A. S. 2021, MNRAS, 500, 3481 [Google Scholar]
 Hamers, A. S., Portegies Zwart, S. F., & Merritt, D. 2014, MNRAS, 443, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Hayli, A. 1970, A&A, 7, 249 [NASA ADS] [Google Scholar]
 Heggie, D. C. 1988, in Longterm Dynamical Behaviour of Natural and Artificial Nbody Systems, ed. A. D. Roy, 329 [CrossRef] [Google Scholar]
 Heggie, D. C. 1991, in Predictability, Stability, and Chaos in NBody Dynamical Systems, eds. S. Roeser, & U. Bastian, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C., & Mathieu, R. D. 1986, in The Use of Supercomputers in Stellar Dynamics, eds. P. Hut, & S. L. W. McMillan (Berlin Springer Verlag), Lect. Notes Phys., 267, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Hemsendorf, M., & Merritt, D. 2002, ApJ, 580, 606 [NASA ADS] [CrossRef] [Google Scholar]
 Hénon, M. 1971, Ap&SS, 13, 284 [Google Scholar]
 Hernandez, D. M., Hadden, S., & Makino, J. 2020, MNRAS, 493, 1913 [NASA ADS] [CrossRef] [Google Scholar]
 Hilbert, D. 1915, Gott. Nachr., 27, 395 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Hut, P., & Heggie, D. C. 2002, J. Stat. Phys., 109, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P., Makino, J., & McMillan, S. 1995, ApJ, 443, L93 [NASA ADS] [CrossRef] [Google Scholar]
 Iorio, L. 2020, Universe, 6, 53 [CrossRef] [Google Scholar]
 ISO 1998, ISO/IEC 14882:1998: Programming languages  C++, 732, Available in electronic form for online purchase at http://webstore.ansi.org/ and http://www.cssinfo.com/ [Google Scholar]
 Itoh, Y. 2009, Phys. Rev. D, 80, 124003 [NASA ADS] [CrossRef] [Google Scholar]
 Kandrup, H. E. 1998, Ann. New York Acad. Sci., 867, 320 [NASA ADS] [Google Scholar]
 Kandrup, H. E., & Smith, H., Jr. 1991, ApJ, 374, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Kandrup, H. E., & Smith, H., Jr. 1992, ApJ, 386, 635 [NASA ADS] [CrossRef] [Google Scholar]
 Kandrup, H. E., Smith, H. Jr., & Willmes, D. E. 1992, ApJ, 399, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Kandrup, H. E., Mahon, M. E., & Smith, H. Jr. 1994, ApJ, 428, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Kepler, J. 1609, Astron. Nova, 1 [Google Scholar]
 Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237 [Google Scholar]
 King, I. R. 1966, AJ, 71, 64 [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
 Kustaanheimo, P., & Stiefel, E. 1965, J. Reine. Angew. Math., 218, 204 [Google Scholar]
 Lakshmanan, M., & Rajasekar, S. 2003, Chaos in Conservative Systems (Berlin, Heidelberg: Springer, Berlin Heidelberg), 191 [Google Scholar]
 Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Quinn, T., & Tremaine, S. 1992, Icarus, 95, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. 1962, Planet. Space Sci., 9, 719 [Google Scholar]
 Lim, H., & Rodriguez, C. L. 2020, Phys. Rev. D, 102, 064033 [NASA ADS] [CrossRef] [Google Scholar]
 Lorentz, H., & Droste, J. 1917, in Verslag Koninklijker Akademie van Wetenschchappen, 26, 392 [Google Scholar]
 Lousto, C. O., & Nakano, H. 2008, CQG, 25, 195019 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, J. 1991, ApJ, 369, 200 [Google Scholar]
 Makino, J. 2002, New Astron., 7, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, J., & Taiji, M. 1996, Comput. Phys., 10, 352 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, J., & Taiji, M. 1998, in Scientific simulations with specialpurpose computers : The GRAPE systems, eds. J. Makino, & M. Taiji (Chichester; Toronto: John Wiley& Sons), c1998 [Google Scholar]
 McMillan, S. L. W. 1986, in The Use of Supercomputers in Stellar Dynamics, eds. P. Hut, & S. L. W. McMillan, (Berlin Springer Verlag), Lect. Notes Phys., 267, 156 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, S. L. W. 2014, in AAS/Division of Dynamical Astronomy Meeting, AAS/Div. Dyn. Astron. Meeting, 45, 303.01 [NASA ADS] [Google Scholar]
 Meiron, Y., Kocsis, B., & Loeb, A. 2017, ApJ, 834, 200 [NASA ADS] [CrossRef] [Google Scholar]
 Mel’nikov, A. V., Orlov, V. V., & Shevchenko, I. I. 2013, Astron. Rep., 57, 429 [CrossRef] [Google Scholar]
 Merritt, D. 2013, Dynamics and Evolution of Galactic Nuclei [Google Scholar]
 Mewes, V., Zlochower, Y., Campanelli, M., et al. 2018, Phys. Rev. D, 97, 084059 [Google Scholar]
 Mewes, V., Zlochower, Y., Campanelli, M., et al. 2020, Phys. Rev. D, 101, 104007a [NASA ADS] [CrossRef] [Google Scholar]
 Mikkola, S., & Merritt, D. 2006, MNRAS, 372, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Mikkola, S., & Merritt, D. 2008, AJ, 135, 2398 [Google Scholar]
 Mikkola, S., & Tanikawa, K. 1999, MNRAS, 310, 745 [NASA ADS] [CrossRef] [Google Scholar]
 Milani, A., & Nobili, A. M. 1992, Nature, 357, 569 [Google Scholar]
 Miller, R. H. 1964, ApJ, 140, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Montgomery, R. 1998, Nonlinearity, 11, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Moore, C. 1993, Phys. Rev. Lett., 3675 [NASA ADS] [CrossRef] [Google Scholar]
 Muno, M. P., Baganoff, F. K., Bautz, M. W., et al. 2004, ApJ, 613, 326 [NASA ADS] [CrossRef] [Google Scholar]
 Muzzio, J. C., & Mosquera, M. E. 2004, Celest. Mech. Dyn. Astron., 88, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Muzzio, J. C., Navone, H. D., & Zorzi, A. F. 2009, Celest. Mech. Dyn. Astron., 105, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Naoz, S., Kocsis, B., Loeb, A., & Yunes, N. 2013a, ApJ, 773, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013b, MNRAS, stt302 [Google Scholar]
 Neilsen, D., Jay, J., & Morgan, T. 2014, in APS April Meeting Abstracts, APS Meeting Abstracts, 2014, M15.008 [Google Scholar]
 Newton, I. 1687, Philosophiae Naturalis Principia Mathematica, 1 [Google Scholar]
 Nieto, J. A., Saucedo, J., & Villanueva, V. M. 2003, Phys. Lett. A, 312, 175 [CrossRef] [Google Scholar]
 Nitadori, K., & Makino, J. 2008, New Astron., 13, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Oliphant, T. E. 2006, A guide to NumPy,, 1 (USA: Trelgol Publishing) [Google Scholar]
 Parvulesco, C. 1924, Bull. Astron., 5, 72 [NASA ADS] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
 Poisson, E., & Will, C. M. 2014, Gravity [Google Scholar]
 Portegies Zwart, S., & Boekholt, T. 2014, ApJ, 785, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., & Boekholt, T. C. N. 2018, Commun. Nonlinear Sci. Numer. Simul., 61, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S., & McMillan, S. 2018, Astrophysical Recipes; The art of AMUSE [Google Scholar]
 Portegies Zwart, S. F., Hut, P., Makino, J., & McMillan, S. L. W. 1998, A&A, 337, 363 [NASA ADS] [Google Scholar]
 Portegies Zwart, S. F., Baumgardt, H., McMillan, S. L. W., et al. 2006, ApJ, 641, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., Belleman, R. G., & Geldof, P. M. 2007, New Astron., 12, 641 [Google Scholar]
 Portegies Zwart, S., McMillan, S., Groen, D., et al. 2008, New Astron., 13, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S., McMillan, S., Harfst, S., et al. 2009, New Astron., 14, 369 [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L., van Elteren, A., Pelupessy, F. I., & de Vries, N. 2013, Comput. Phys. Commun., 184, 456 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S., van Elteren, A., Pelupessy, I., et al. 2018, AMUSE: the Astrophysical Multipurpose Software Environment [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992. Numerical recipes in C. The art of scientific computing, 2nd edn. (Cambridge: University Press), c1992 [Google Scholar]
 Pretorius, F. 2005, Phys. Rev. Lett., 95, 121101 [NASA ADS] [CrossRef] [Google Scholar]
 Rimoldi, A., Rossi, E. M., Piran, T., & Portegies Zwart, S. 2015, MNRAS, 447, 3096 [NASA ADS] [CrossRef] [Google Scholar]
 Robson, T., Cornish, N. J., Tamanini, N., & Toonen, S. 2018, Phys. Rev. D, 98, 064012 [NASA ADS] [CrossRef] [Google Scholar]
 Schäfer, G. 1987, Phys. Lett. A, 123, 336 [CrossRef] [Google Scholar]
 Schwarzschild, K. 1916, Abh. Konigl. Preuss. Akad. Wissenschaften Jahre 1906,92, Berlin, 1907, 1916, 189 [Google Scholar]
 Shivamoggi, B. K. 2014, Chaos in Dissipative Systems (Dordrecht, Netherlands: Springer), 189 [Google Scholar]
 Spitzer, L. 1971, in Pontificiae Academiae Scientiarum Scripta Varia, Proceedings of a Study Week on Nuclei of Galaxies, held in Rome, April 13–18, 1970, Amsterdam: North Holland, and New York: American Elsevier, 1971, ed. D. J. K. O’Connell, 443 [Google Scholar]
 Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton, NJ: Princeton University Press), 191 [Google Scholar]
 Spitzer, L. J., & Hart, M. H. 1971a, ApJ, 164, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, L. J., & Hart, M. H. 1971b, ApJ, 166, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Spyrou, N. 1975, ApJ, 197, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Stephan, A. P., Naoz, S., & Gaudi, B. S. 2021, ApJ, 922, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Stiefel, E. L., & Scheifele, G. 1975, Linear and regular celestial mechanics. Perturbed twobody motion. Numerical methods. Canonical theory [Google Scholar]
 Tokovinin, A. 2014, AJ, 147, 87 [CrossRef] [Google Scholar]
 Urminsky, D. J., & Heggie, D. C. 2009, MNRAS, 392, 1051 [NASA ADS] [CrossRef] [Google Scholar]
 Valluri, M., & Merritt, D. 2000, in The Chaotic Universe, eds. V. G. Gurzadyan, & R. Ruffini, 229 [NASA ADS] [CrossRef] [Google Scholar]
 van Albada, T. S. 1968, Bull. Astron. Inst. Neth., 19, 479 [Google Scholar]
 van Rossum, G. 1995, Extending and embedding the Python interpreter, Report CSR9527, pubCWI, pubCWI:adr [Google Scholar]
 Verlet, L. 1967, Phys. Rev., 159, 98 [CrossRef] [Google Scholar]
 von Hoerner, S. 1960, Z. Astrophys., 50, 184 [NASA ADS] [Google Scholar]
 von Zeipel, H. 1910, Astron. Nachr., 183, 345 [Google Scholar]
 Waldvogel, J. 2006, Celest. Mech. Dyn. Astron., 95, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Waldvogel, J. 2008, Celest. Mech. Dyn. Astron., 102, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Wanex, L. F. 2002, PhD Thesis, University of Nevada, Reno, USA [Google Scholar]
 Weinberg, S. 1972, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity [Google Scholar]
 Weinstein, G. 2012, ArXiv eprints [arXiv:1202.2791] [Google Scholar]
 Will, C. M. 2011, Proc. Nat. Acad. Sci., 108, 5938 [NASA ADS] [CrossRef] [Google Scholar]
 Will, C. M. 2014, Phys. Rev. D, 89, 044043 [NASA ADS] [CrossRef] [Google Scholar]
 Will, C. M. 2021, Phys. Rev. D, 103, 063003 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, X., & Huang, T.Y. 2003, Phys. Lett. A, 313, 77 [CrossRef] [Google Scholar]
 Zadunaisky, P. E. 1979, Celest. Mech., 20, 209 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Implementation of the numerical method
The Newtonian solver and the postNewtonian terms are implemented in C and C++ and interfaced with the Astrophysics Multipurpose Software Environment (AMUSE for short, Portegies Zwart & McMillan 2018). In this appendix, we discuss the simple Hermite predictorcorrector Newtonian Nbody solver called ph4 and the postNewtonian solver called Hermite_GRX. The former solves Newton’s equations of motion quite accurately, but is unable to achieve converged solutions. The other code adopts the postNewtonian approach in which we address the pairwise EIH equation as well as the socalled crossterms (Einstein et al. 1938). All equations are implemented to 1PN order.
The Newtonian code, ph4, is optimized for parallel operations using the Message Passing Interface (MPI, Gropp et al. 1996; Gropp 2002), and for GPU using the Sapporo library (Portegies Zwart et al. 2007; Gaburov et al. 2009). The postNewtonian implementation, Hermite_GRX is parallelized using hyperthreading, but not using MPI, and it does not support GPU operations. In this code, however, fewbody interactions are regularized using quaternions.
In the following sections, we discuss the various implementations and optimizations. We also perform some test calculations to demonstrate the efficiency and accuracy of the various implementations.
A.1. Fourthorder Hermite integration scheme
Here we describe the fourthorder Hermite predictor corrector implementation in ph4 and Hermite_GRX briefly. The first solves Newton’s equations of motion; the second also includes various solvers for addressing the postNewtonian expansion terms.
A.1.1. Predict, evaluate, and correct scheme
The Hermite integration scheme is a family of implicit numerical methods for solving ordinary differential equations. Introduced by Makino (1991), the fourthorder integration scheme is written
$$\begin{array}{cc}\hfill y(t+h)& =y(t)+\frac{{y}^{(1)}(t)+{y}^{(1)}(t+h)}{2}h\hfill \\ \hfill & +\frac{{y}^{(2)}(t){y}^{(2)}(t+h)}{12}{h}^{2}+\mathcal{O}\left({h}^{5}\right),\hfill \end{array}$$(A.1)
which has a local truncation error of 𝒪(h^{5}), resulting in a global truncation error of 𝒪(h^{4}). We denote the ith derivative with respect to t using ( ⋅ )^{(i)} or with Einstein’s convention. The sixth and eightorder schemes, derived by Nitadori & Makino (2008), are not implemented here.
Because the scheme is implicit, a fixedpoint iteration to solve eq. A.2 is needed,
$$\begin{array}{cc}\hfill {y}_{[i+1]}(t+h)& =y(t)+\frac{{y}^{(1)}(t)+{y}_{[i]}^{(1)}(t+h)}{2}h\hfill \\ \hfill & +\frac{{y}^{(2)}(t){y}_{[i]}^{(2)}(t+h)}{12}{h}^{2}.\hfill \end{array}$$(A.2)
Here we used a truncated Taylor expansion around t as the initial (boundary) condition
$$\begin{array}{c}\hfill {y}_{[0]}(t+h)=y(t)+h{y}^{(1)}(t)+\frac{{h}^{2}}{2}{y}^{(2)}(t).\end{array}$$(A.3)
This sequence converges to the limit y(t + h), which ends the current time step. In practice, a single iteration suffices when the time step h is small.
Each integration step then consists of

prediction of the positions and velocities at the next time step t + h,
$$\begin{array}{cc}& {\stackrel{\sim}{\mathit{r}}}_{i+1}={\mathit{r}}_{i}+{\mathit{v}}_{i}h+{\textstyle \frac{1}{2}}\stackrel{\sim}{{\mathit{a}}_{i}}{h}^{2}+{\textstyle \frac{1}{6}}{\mathit{j}}_{i}{h}^{3},\hfill \end{array}$$(A.4a)
$$\begin{array}{cc}& {\stackrel{\sim}{\mathit{v}}}_{i+1}={\mathit{v}}_{i}+{\mathit{a}}_{i}h+{\textstyle \frac{1}{2}}{\mathit{j}}_{i}{h}^{2}.\hfill \end{array}$$(A.4b)
Here r, v, a, and $\mathit{j}=\frac{\text{d}\mathit{a}}{\text{d}t}\equiv \dot{a}$ represent vectors for the position, velocity, acceleration, and jerk, respectively. The jerk j is dotted, which in this case does not indicate a time derivative. Predicted values are indicated with the $\stackrel{\sim}{(\xb7)}$.

Acceleration and jerk are calculated using the predicted positions and velocities (eq. A.4).

A subsequent correction is applied to the position and velocity at the next time step using the predicted accelerations and jerks,
$$\begin{array}{cc}& {\mathit{v}}_{i+1}={\mathit{v}}_{i}+{\textstyle \frac{1}{2}}({\mathit{a}}_{i}+{\stackrel{\sim}{\mathit{a}}}_{i+1})h+{\textstyle \frac{1}{12}}({\mathit{j}}_{i}{\stackrel{\sim}{\mathit{j}}}_{i+1}){h}^{2},\hfill \end{array}$$(A.5a)
$$\begin{array}{cc}& {\mathit{r}}_{i+1}={\mathit{r}}_{i}+{\textstyle \frac{1}{2}}({\mathit{v}}_{i}+{\mathit{v}}_{i+1})h+{\textstyle \frac{1}{12}}({\mathit{a}}_{i}{\stackrel{\sim}{\mathit{a}}}_{i+1}){h}^{2}.\hfill \end{array}$$(A.5b)
The corrected velocities increase the order of the method to 𝒪(h^{4}). In such a predict, evaluate, and correct (PEC) scheme, the fixedpoint iteration can be described as P(EC)^{n} for n iterations.
A.1.2. Variable time step
We use variable but shared time steps. After every step, a new timestep size is determined based on the minimum interparticle collision timescale, calculated from unaccelerated linear motion and the freefall time,
$$\begin{array}{c}\hfill h=\eta \underset{i,j\ne i}{min}(\frac{{\mathit{r}}_{\mathit{ij}}}{{\mathit{v}}_{\mathit{ij}}},\frac{{\mathit{r}}_{\mathit{ij}}}{({m}_{i}+{m}_{j}){\mathit{a}}_{\mathit{ij}}}).\end{array}$$(A.6)
Here r_{ij}, v_{ij}, a_{ij} are the relative distance, velocity, and acceleration between particles i and j, and m_{i} is the mass of particle i. The minimum is taken over each pair of particles (i, j), and over the two estimates of the collision time. Here the timestep parameter η is introduced to control the timestep size and therewith the accuracy (and speed) of the integration scheme. Ler values of η generally correspond to smaller errors and a longer integration wallclock time. The default value in AMUSE, η = 0.03, generally leads to acceptable accuracy at a reasonable speed: in many cases, η = 0.1 probably suffices (Portegies Zwart & Boekholt 2014). For safety, we adopted η = 0.01 for our calculations.
The adopted variable time step removes the timesymmetric properties of the integration. The fundamental idea behind timesymmetrization is to prevent systematic drift in any conserved quantity. Time reversibility then introduces the same drift with opposite sign.
A timesymmetric algorithm exhibits the same drift in both directions of time, resulting in identical absolute drifts when integrating forward and backward with time. This is a desirable quality of an integrator because we consider Nature to conserve energy and angular momentum (see also Portegies Zwart & Boekholt 2018, for a discussion on the arrow of time due to the chaotic behavior of selfgravitating systems and uncertainties on the smallest scales).
One can reintroduce timesymmetry by selecting a symmetric time step, for example, by taking the average of some function at either side of the integration step (Hut et al. 1995),
$$\begin{array}{c}\hfill h={\textstyle \frac{1}{2}}(k({t}_{b})+k({t}_{e})).\end{array}$$(A.7)
Here k(t) is a function to determine the step size at the beginning t_{b} and at the end t_{e} = t_{b} + h of the integration step. This implicit expression requires fixedpoint iteration to evaluate
$$\begin{array}{c}\hfill {h}_{[0]}=k({t}_{b}),\end{array}$$(A.8a)
and the eventual time step when the sequence converges becomes
$$\begin{array}{c}\hfill {h}_{[i+1]}={\textstyle \frac{1}{2}}(k({t}_{b})+k({t}_{b}+{h}_{[i]})).\end{array}$$(A.9a)
Generally, the sequence converges in a single iteration (Hut et al. 1995).
A.1.3. Splitting the jerk
Calculating the jerk is expensive in terms of computer time because it requires three passes over all particles. To avoid evaluating the jerk directly, we use a central numerical derivative,
$$\begin{array}{c}\hfill \mathit{j}(t)=\frac{\mathit{a}(t+h)\mathit{a}(th)}{2h}+\mathcal{O}\left({h}^{2}\right).\end{array}$$(A.10)
Here the accelerations are calculated using the Taylor expanded positions and velocities of the particles,
$$\begin{array}{cc}& \mathit{r}(t\pm h)=\mathit{r}(t)\pm h\mathit{v}(t)+{\textstyle \frac{1}{2}}{h}^{2}\mathit{a}(t)+\mathcal{O}\left({h}^{3}\right),\hfill \end{array}$$(A.11a)
$$\begin{array}{cc}& \mathit{v}(t\pm h)=\mathit{v}(t)\pm h\mathit{a}(t)+\mathcal{O}\left({h}^{2}\right).\hfill \end{array}$$(A.11b)
The numerical calculation of the jerk is equally expensive as the analytic calculation because it requires two additional acceleration calculations per jerk. We reduce the computational complexity by the time step h of the previous integration step and the backward derivative. This allows us to reuse the previous steps’ positions and velocities for calculating the jerk at the current time step.
Using a firstorder derivative instead of the analytical expression for the jerk leads to a reduced accuracy, but this is corrected for by splitting the acceleration into two parts: the Newtonian part, and a perturbing part,
$$\begin{array}{c}\hfill \mathit{a}(t)={\mathit{a}}_{\mathrm{Newton}}(t)+{\mathit{a}}_{\mathrm{pert}}(t).\end{array}$$(A.12)
We note here that we already introduced a perturbation, which in the EIH equations of motion will be the postNewtonian terms. The Newtonian jerk can now be calculated analytically and at negligible cost compared to calculating the perturbing acceleration. The perturbing jerk is calculated from the numerical backward derivative,
$$\begin{array}{c}\hfill \mathit{j}(t)={\mathit{j}}_{\mathrm{Newton}}(t)+\frac{{\mathit{a}}_{\mathrm{pert}}(t){\mathit{a}}_{\mathrm{pert}}(th)}{h}+\mathcal{O}\left({h}^{2}\right).\end{array}$$(A.13)
This operation increases the memory requirement by storing two accelerations for each particle.
The algorithm is made to be selfstarting by defining the perturbing acceleration of the previous integration step a_{i − 1, pert}, as it depends on the jerk of the first iteration j_{i}. For this, we chose a_{i − 1, pert} = a_{i, pert}, so that j_{i, pert} = 0. This decreases the local truncation error of the first step to 𝒪(h^{4}), but its impact on the results is small because only one step is taken.
A.2. Regularization
Regularizing the equations of motion for two bodies (or more) in a close encounter improves computational performance and accuracy. The main reason to introduce regularization, however, is to prevent a devision by zero for extremely close approaches between particles (Kustaanheimo & Stiefel 1965; Mikkola & Tanikawa 1999). Regularizing the postNewtonian expressions is harder than the regular Newtonian case because of the velocity dependence on the acceleration (Mikkola & Merritt 2006, 2008). Here we derive the regularized equations of motion in postNewtonian fewbody encounters using quaternions (Waldvogel 2006), but we start with a brief overview on quaternions
A.2.1. Quaternions
Quaternions are an extension of complex numbers to three complex base quaternions i, j, and k (Waldvogel 2006). A quaternion u is constructed from four real numbers u_{ℓ} ∈ ℝ for ℓ = 0, 1, 2, 3,
$$\begin{array}{c}\hfill \mathit{u}={u}_{0}+{u}_{1}\mathit{i}+{u}_{2}\mathit{j}+{u}_{3}\mathit{k},\end{array}$$(A.14)
with the multiplicative identities
$$\begin{array}{c}\hfill \mathit{i}\mathit{i}=\mathit{j}\mathit{j}=\mathit{k}\mathit{k}=\mathit{i}\mathit{j}\mathit{k}=1,\end{array}$$(A.15)
from which we derive the other products,
$$\begin{array}{c}\hfill \{\begin{array}{cc}\mathit{i}\mathit{j}\hfill & =\mathit{j}\mathit{i}=\mathit{k},\hfill \\ \mathit{j}\mathit{k}\hfill & =\mathit{k}\mathit{j}=\mathit{i},\hfill \\ \mathit{k}\mathit{i}\hfill & =\mathit{i}\mathit{k}=\mathit{j}.\hfill \end{array}\end{array}$$(A.16)
From these, the noncummatative property of quaternion multiplication is evident. We define the conjugate of quaternion u,
$$\begin{array}{c}\hfill \overline{\mathit{u}}={u}_{0}{u}_{1}\mathit{i}{u}_{2}\mathit{j}{u}_{3}\mathit{k},\end{array}$$(A.17)
which leads to the definition of the norm
$$\begin{array}{c}\hfill {\mathit{u}}^{2}=\mathit{u}\overline{\mathit{u}}=\overline{\mathit{u}}\mathit{u}={u}_{0}^{2}+{u}_{1}^{2}+{u}_{2}^{2}+{u}_{3}^{2}.\end{array}$$(A.18)
The star conjugate of u then is
$$\begin{array}{c}\hfill {\mathit{u}}^{\star}={u}_{0}+{u}_{1}\mathit{i}+{u}_{2}\mathit{j}{u}_{3}\mathit{k}.\end{array}$$(A.19)
We associate the vector r = (r_{0}, r_{1}, r_{2})∈ℝ^{3} to quaternion u as
$$\begin{array}{c}\hfill \mathit{u}={r}_{0}+{r}_{1}\mathit{i}+{r}_{2}\mathit{j}.\end{array}$$(A.20)
When describing realworld coordinates, quaternions have a vanishing component in k (Waldvogel 2006).
A.2.2. Equations of motion
Regularization is applied to pairs of particles. They, particle i = 1 and i = 2, are located at phasespace coordinate r_{i} with velocities v_{i}, and masses m_{i}. The equations of motion are
$$\begin{array}{c}\hfill {\ddot{\mathit{r}}}_{i}=\frac{G{m}_{j}}{{r}^{3}}({\mathit{r}}_{j}{\mathit{r}}_{i})+{\mathit{a}}_{i}.\end{array}$$(A.21)
The acceleration a_{i} for particle i consists of all perturbing (postNewtonian) and Newtonian accelerations, excluding the Newtonian acceleration between particle i and j. The center of mass
$$\begin{array}{c}\hfill {\mathit{r}}_{\mathrm{cm}}=\frac{{m}_{1}{\mathit{r}}_{1}+{m}_{2}{\mathit{r}}_{2}}{{m}_{1}+{m}_{2}},\end{array}$$(A.22)
and the relative position
$$\begin{array}{c}\hfill {\mathit{r}}_{12}={\mathit{r}}_{2}{\mathit{r}}_{1},\end{array}$$(A.23)
we rewrite these equations of motion as
$$\begin{array}{c}\hfill \ddot{\mathit{r}}=\frac{\mu}{{r}^{3}}\mathit{r}+\mathit{P},\end{array}$$(A.24)
and
$$\begin{array}{c}\hfill {\ddot{\mathit{r}}}_{\mathrm{cm}}=\frac{{m}_{1}{\mathit{a}}_{1}+{m}_{2}{\mathit{a}}_{2}}{{m}_{1}+{m}_{2}}.\end{array}$$(A.25)
Here μ = G(m_{1} + m_{2}) is the gravitational parameter, and P = a_{2} − a_{1} the relative perturbing acceleration.
The above equations of motion have a singular point at r = 0. We integrate the center of mass separately from the relative position, and rewrite the equation of motion in such in a way as to remove this singular point. This is achieved by remapping the world position to a regularized position quaternion u, from which we calculate the world position
$$\begin{array}{c}\hfill \mathit{r}=\mathit{u}{\mathit{u}}^{\star}.\end{array}$$(A.26)
The quaternion r now has vanishing component k and can therefore be transformed into a vector, from which it follows that
$$\begin{array}{c}\hfill r=\mathit{r}={\mathit{u}}^{2}=\mathit{u}\overline{\mathit{u}}.\end{array}$$(A.27)
A (nonunique) solution to the inverse of eq. A.26 is
$$\begin{array}{c}\hfill \widehat{\mathit{u}}=\frac{\mathit{r}+\mathit{r}}{\sqrt{2(\mathit{r}+{r}_{0})}}.\end{array}$$(A.28)
The position vector, r, is almost entirely oriented in the negative rdirection, and the denominator is close to zero. Without loss of generality, we can avoid large numerical errors by swapping indices i = 1 and 2, resulting in the negation of r.
The regularization time τ is
$$\begin{array}{c}\hfill \mathrm{d}t=r\phantom{\rule{3.33333pt}{0ex}}\mathrm{d}\tau .\end{array}$$(A.29)
The equations of motion for the regularized position is written in regularized time:
$$\begin{array}{c}\hfill {\mathit{u}}^{(2)}={\textstyle \frac{1}{2}}b\mathit{u}+{\textstyle \frac{1}{2}}r\mathit{P}{\overline{\mathit{u}}}^{\star}.\end{array}$$(A.30)
Here b is the binding energy of the binary,
$$\begin{array}{c}\hfill b=\frac{\mu}{r}\frac{1}{2}{\dot{\mathit{r}}}^{2}\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}=\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}\frac{\mu 2{\mathit{u}}^{(1)}{}^{2}}{{\mathit{u}}^{2}}.\end{array}$$(A.31)
For an unperturbed twobody system, P = 0, the binding energy b is constant, and the equation of motion describes the harmonic oscilator. For a perturbed twobody system, the expression resembles a perturbed harmonic oscilator. The greatest advantage of this approach is that this equation of motion has no singular points, not even for r = u^{2} = 0. This improves the performance for small perturbations when numerically integrating the equations of motion of a highly eccentric binary, and it improves the accuracy for perturbed binaries.
The binding energy (eq. A.31) is not regularized, and errors continue to increase in close encounters (Funato et al. 1996). This problem is mitigated by also integrating the binding energy numerically. When the binding energy changes slowly with time as a function of the perturbing accelerations, we find
$$\begin{array}{c}\hfill {b}^{(1)}=\u27e8{\mathit{r}}^{\prime},\mathit{P}\u27e9.\end{array}$$(A.32)
Here the ⟨(⋅),(⋅)⟩ is the vectorial scalar product. An initial condition for u^{(1)} is
$$\begin{array}{c}\hfill {\widehat{\mathit{u}}}^{(1)}={\textstyle \frac{1}{2}}\mathit{v}{\overline{\mathit{u}}}^{\star}.\end{array}$$(A.33)
This expression gives a small correction to the original derivation by Waldvogel (2006) and Waldvogel (2008): it can be verified by substitution in eq. A.31, to derive the binding energy in world coordinates. For the reciprocal eq. A.33, the righthand side multiplication with u^{⋆} leads to
$$\begin{array}{c}\hfill \mathit{v}={\textstyle \frac{2}{r}}\dot{\mathit{u}}{\mathit{u}}^{\star}\equiv {\textstyle \frac{2}{r}}{\mathit{u}}^{(1)}{\mathit{u}}^{\star}.\end{array}$$(A.34)
One can also formulate the KustaanheimoStiefel (hereafter KS, Kustaanheimo & Stiefel 1965) regularized equations of motion in terms of the perturbing potential (Stiefel & Scheifele 1975), which can be advantageous in some cases, although it only applies to cases when the potential is independent of velocity (not the case for general relativity). Moreover, Stiefel and Steifele (1975, see pages 30 and 31) argued that it is numerically more efficient to formulate the equations of motion in terms of the total energy, not just the Kepler energy (see their eq. (A.31)). The equations of motion can also be cast into the form of regular elements (Stiefel & Scheifele 1975, their pages 90 and 91) that are advantageous for perturbed twobody systems because the regular elements remain exactly constant for nonperturbed systems (i.e., the twobody system is integrated analytically).
A.2.3. Numerical integration
We solve the equations of motion using the Hermite scheme in the PEC formulation (see sect. A.1.1), with Δτ ≡ κ. The resulting intergration scheme is described below.

Predict the regularized position, regularized velocity, and binding energy at the end of the current integration step,
$$\begin{array}{cc}& {\stackrel{\sim}{\mathit{u}}}_{i+1}={\mathit{u}}_{i}+{\mathit{u}}_{i}^{(1)}\kappa +{\textstyle \frac{1}{2}}{\mathit{u}}_{i}^{(2)}{\kappa}^{2}+{\textstyle \frac{1}{6}}{\mathit{u}}_{i}^{(3)}{\kappa}^{3},\hfill \end{array}$$(A.35a)
$$\begin{array}{cc}& {\stackrel{\sim}{\mathit{u}}}_{i+1}^{(1)}={\mathit{u}}_{i}^{(1)}+{\mathit{u}}_{i}^{(2)}\kappa +{\textstyle \frac{1}{2}}{\mathit{u}}_{i}^{(3)}{\kappa}^{2},\hfill \end{array}$$(A.35b)
$$\begin{array}{cc}& {\stackrel{\sim}{b}}_{i+1}={b}_{i}+{b}_{i}^{(1)}\kappa +{\textstyle \frac{1}{2}}{b}_{i}^{(2)}{\kappa}^{2}.\hfill \end{array}$$(A.35c)
Then predict the centerofmass position and velocity from eq. A.4.

Evaluate the regularized acceleration, regularized jerk, and the derivatives of the binding energy at the end of the current integration step, according to eqs. A.30 and A.32,
$$\begin{array}{cc}& {\mathit{u}}^{(3)}={\textstyle \frac{1}{2}}({b}^{(1)}\mathit{u}b{\mathit{u}}^{(1)}+{r}^{(1)}\mathit{P}{\overline{\mathit{u}}}^{\star}+r{\mathit{P}}^{(1)}{\overline{\mathit{u}}}^{\star}+r\mathit{P}{({\overline{\mathit{u}}}^{(1)})}^{\star}),\hfill \end{array}$$(A.36a)
$$\begin{array}{cc}& {b}^{(2)}=\u27e8{\mathit{r}}^{(2)},\mathit{P}\u27e9\u27e8{\mathit{r}}^{(1)},{\mathit{P}}^{(1)}\u27e9.\hfill \end{array}$$(A.36b)
The testparticle integrator by Hamers et al. (2014) adopts a similar integration scheme for the KS coordinates.

Correct the regularized position, regularized velocity, and binding energy at the end of the current integration step,
$$\begin{array}{cc}& {\mathit{u}}_{i+1}^{(1)}={\mathit{u}}_{i}^{(1)}+{\textstyle \frac{1}{2}}({\mathit{u}}_{i}^{(2)}+{\stackrel{\sim}{\mathit{u}}}_{i+1}^{(2)})\kappa +{\textstyle \frac{1}{12}}({\mathit{u}}_{i}^{(3)}{\stackrel{\sim}{\mathit{u}}}_{i+1}^{(3)}){\kappa}^{2},\hfill \end{array}$$(A.37a)
$$\begin{array}{cc}& {\mathit{u}}_{i+1}={\mathit{u}}_{i}+{\textstyle \frac{1}{2}}({\mathit{u}}_{i}^{(1)}+{\mathit{u}}_{i+1}^{(1)})\kappa +{\textstyle \frac{1}{12}}({\mathit{u}}_{i}^{(2)}{\stackrel{\sim}{\mathit{u}}}_{i+1}^{(2)}){\kappa}^{2},\hfill \end{array}$$(A.37b)
$$\begin{array}{cc}& {b}_{i+1}={b}_{i}+{\textstyle \frac{1}{2}}({b}_{i}^{(1)}+{\stackrel{\sim}{b}}_{i+1}^{(1)})\kappa +{\textstyle \frac{1}{12}}({b}_{i}^{(2)}{\stackrel{\sim}{b}}_{i+1}^{(2)}){\kappa}^{2}.\hfill \end{array}$$(A.37c)
The corrected velocity in the corrector for the position makes the scheme fourth order (Mikkola & Merritt 2006).
A.2.4. Regularized timestep considerations
The regularized time step is determined using
$$\begin{array}{c}\hfill s=\eta min(\sqrt{\frac{{\mathit{u}}^{(2)}\mathit{u}}{{\mathit{u}}^{(3)}{\mathit{u}}^{(1)}}},\frac{{\mathit{u}}^{(1)}}{{\mathit{u}}^{(2)}}),\end{array}$$(A.38)
which is a variation on the time steps suggested by Funato et al. (1996). To advance model time, we need to convert this κ into h, which, according to Funato et al. (1996), is done with
$$\begin{array}{c}\hfill h\equiv T(\kappa )={t}_{{\textstyle \frac{1}{2}}}^{(1)}\kappa +{\textstyle \frac{1}{24}}{t}_{{\textstyle \frac{1}{2}}}^{(3)}{\kappa}^{3}+{\textstyle \frac{1}{1920}}{t}_{{\textstyle \frac{1}{2}}}^{(5)}{\kappa}^{5}.\end{array}$$(A.39)
Here ${t}_{{\textstyle \frac{1}{2}}}^{(i)}$ is the ith derivative of t with respect to τ from the begin time τ_{b} to $\tau ={\tau}_{b}+{\textstyle \frac{1}{2}}\kappa $ , given by
$$\begin{array}{cc}& {t}_{{\textstyle \frac{1}{2}}}^{(1)}={\mathit{u}}^{2},\hfill \end{array}$$(A.40a)
$$\begin{array}{cc}& {t}_{{\textstyle \frac{1}{2}}}^{(3)}={\mathit{u}}^{(2)}\overline{\mathit{u}}+2{{\mathit{u}}^{(1)}}^{2}+\mathit{u}{\overline{\mathit{u}}}^{(2)},\hfill \end{array}$$(A.40b)
$$\begin{array}{cc}& {t}_{{\textstyle \frac{1}{2}}}^{(5)}={\mathit{u}}^{(4)}\overline{\mathit{u}}+4{\mathit{u}}^{(3)}{\overline{\mathit{u}}}^{(1)}+6{{\mathit{u}}^{(2)}}^{2}+4{\mathit{u}}^{(1)}{\overline{\mathit{u}}}^{(3)}+\mathit{u}{\overline{\mathit{u}}}^{(4)}.\hfill \end{array}$$(A.40c)
Here u and all its derivatives are evaluated at half timesteps, using a Taylor expansion at the beginning of the integration step τ_{b}, using the approximated regularized snap and crackle,
$$\begin{array}{cc}& {s}_{b}\equiv {\mathit{u}}_{b}^{(4)}=6\frac{{\mathit{u}}_{b}^{(2)}{\mathit{u}}_{e}^{(2)}}{{\kappa}^{2}}\frac{4{\mathit{u}}_{b}^{(3)}+2{\mathit{u}}_{e}^{(3)}}{\kappa},\hfill \end{array}$$(A.41a)
$$\begin{array}{cc}& {c}_{b}\equiv {\mathit{u}}_{b}^{(5)}=12\frac{{\mathit{u}}_{b}^{(2)}{\mathit{u}}_{e}^{(2)}}{{\kappa}^{3}}+6\frac{{\mathit{u}}_{b}^{(3)}+{\mathit{u}}_{e}^{(3)}}{{\kappa}^{2}}.\hfill \end{array}$$(A.41b)
Equation A.39 has 𝒪(κ^{7}), ensuring an accuracy of 𝒪(κ^{6}) in the final integrated time. The inverse of this transformation cannot be found analytically. We use NewtonRapson iteration,
$$\begin{array}{cc}& {\kappa}_{[0]}=\frac{h}{{\mathit{u}}_{b}{}^{2}},\hfill \end{array}$$(A.42a)
$$\begin{array}{cc}& {\kappa}_{[i+1]}={\kappa}_{[i]}\frac{T({\kappa}_{[i]})b}{{\mathit{u}}_{{\textstyle \frac{1}{2}}}{}^{2}}.\hfill \end{array}$$(A.42b)
Here ${\mathit{u}}_{{\textstyle \frac{1}{2}}}$ is given by a thirdorder Taylor expansion. In practice, convergence to machine precision is reached in four iterations.
A.2.5. Selecting the regularized particle pair
Selection of which particles are to be regularized is done by calculating a regularization criterion for each pair of particles. We then sort the resulting list of pairs, only keeping the N_{reg} highestgraded pairs. Here N_{reg} is a free parameter. We then regularize all particle pairs if they have the smallest norm of the relative acceleration
$$\begin{array}{c}\hfill {Q}_{\mathrm{reg}}=\frac{G({m}_{1}+{m}_{2})}{{r}^{2}}.\end{array}$$(A.43)
The resulting complete integration steps are listed below.

Select pairs of particles that need to be regularized and resolve their step.

Predict all unregularized particles and all regularized pairs of particles.

Convert coordinates of regularized particles into world coordinates.

Evaluate acceleration and jerks for all pairs of particles, excluding the Newtonian interaction between regularized pairs.

For regularized pairs, calculate derivatives for the acceleration, jerk, and binding energies.

Correct all unregularized particles and all regularized pairs of particles.

Convert coordinates of regularized particles into world coordinates to synchronize the whole system of particles.
To symmetrize the global time step, we first need to symmetrize the time steps in their own coordinate system. Regularized time steps need to be symmetrized in regularized time and then need to be transformed into a world time step.
A.3. Implementation
A.3.1. ph4
The Nbody integrator ph4 uses a fourthorder Hermite scheme similar to that described in §A.1.1, with some differences in detail. Its origin lies in the kira integrator, which was part of the Starlab software suite (Portegies Zwart et al. 1998), but most of the complicating elements in kira, such as treatments of arbitrary multiples and close encounters, and stellar and binary evolution, have been removed to be replaced in the AMUSE model by separate community modules communicating at the Python level. The data structures in ph4 are deliberately kept simple, making for a robust module that is easy to manage as a standalone tool, and ph4 facilitates parallelism as well as GPU acceleration, as discussed below.
Although ph4 uses individual block time steps internally, its basic mode of operation, like that of most AMUSE modules (see section A.3.3), is to take an Nbody system, typically synchronized at some initial time t_{0}, and integrate it forward to some new time t_{1}. The ph4 module uses individual block time steps (McMillan 1986), with essentially the same logic as described by Makino (1991) and used in Starlab. Each particle i has its own current time t_{i} and time step δt_{i}. By rounding all steps down to powers of 2, we open the possibility that many particles can be advanced simultaneously. Specifically, it is often the case that at any stage of the calculation, multiple “iparticles” have the same next time t_{next} = t_{i} + δt_{i}, allowing a onetime parallel prediction of all field positions and velocities,
$$\begin{array}{cc}\hfill {\stackrel{\sim}{\mathit{r}}}_{j}& ={\mathit{r}}_{j}+{\mathit{v}}_{j}{h}_{j}+{\textstyle \frac{1}{2}}{\mathit{a}}_{j}{h}_{j}^{2}+{\textstyle \frac{1}{6}}{\mathit{j}}_{j}{h}_{j}^{3}\hfill \\ \hfill {\stackrel{\sim}{\mathit{v}}}_{j}& ={\mathit{v}}_{j}+{\mathit{a}}_{j}{h}_{j}+{\textstyle \frac{1}{2}}{\mathit{j}}_{j}{h}_{j}^{2}\hfill \end{array}$$
(where h_{j} = t_{next} − t_{j}) and parallel computation of the predicted accelerations and jerks of the iparticles,
$$\begin{array}{cc}\hfill {\stackrel{\sim}{\mathit{a}}}_{i}& ={\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{n}\frac{G{m}_{j}}{{r}_{\mathit{ij}}^{3}}{\stackrel{\sim}{\mathit{r}}}_{\mathit{ij}},}\hfill \end{array}$$(A.44)
$$\begin{array}{cc}\hfill {\stackrel{\sim}{\mathit{j}}}_{i}& ={\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{n}\frac{G{m}_{j}}{{r}_{\mathit{ij}}^{3}}[{\stackrel{\sim}{\mathit{v}}}_{\mathit{ij}}+\frac{3({\stackrel{\sim}{\mathit{v}}}_{\mathit{ij}}\xb7{\stackrel{\sim}{\mathit{r}}}_{\mathit{ij}}){\stackrel{\sim}{\mathit{r}}}_{\mathit{ij}}}{{r}_{\mathit{ij}}^{2}}],}\hfill \end{array}$$(A.45)
where, as before, ${\stackrel{\sim}{\mathit{r}}}_{\mathit{ij}}={\stackrel{\sim}{\mathit{r}}}_{i}{\stackrel{\sim}{\mathit{r}}}_{j}$ and ${\stackrel{\sim}{\mathit{v}}}_{\mathit{ij}}={\stackrel{\sim}{\mathit{v}}}_{j}{\stackrel{\sim}{\mathit{v}}}_{j}$.
The corrector step in ph4 differs from the more elegant version described in §A.1.1, with the formulation chosen to allow the use of the traditional Aarseth (1985) time step formula. Using the available derivative information at the beginning and end of the step, we can estimate the next two derivatives in the Taylor series for the position and the velocity (Makino 1991), snap and crackle,
$$\begin{array}{cc}\hfill {\mathit{s}}_{i}\phantom{\rule{3.33333pt}{0ex}}\equiv \phantom{\rule{3.33333pt}{0ex}}{\ddot{\mathit{a}}}_{i}& =\frac{6({\mathit{a}}_{i}{\stackrel{\sim}{\mathit{a}}}_{i})\delta {t}_{i}(4{\mathit{j}}_{i}+2{\stackrel{\sim}{\mathit{j}}}_{i})}{\delta {t}_{i}^{2}}\hfill \end{array}$$(A.46)
$$\begin{array}{cc}\hfill {\mathit{c}}_{i}\phantom{\rule{3.33333pt}{0ex}}\equiv \phantom{\rule{3.33333pt}{0ex}}{\stackrel{\mathbf{\u20db}}{\mathit{a}}}_{i},& =\frac{12({\mathit{a}}_{i}{\stackrel{\sim}{\mathit{a}}}_{i})+6\delta {t}_{i}({\mathit{j}}_{i}+{\stackrel{\sim}{\mathit{j}}}_{i})}{\delta {t}_{i}^{3}},\hfill \end{array}$$(A.47)
(see also Eqs. A.41a and A.41b), leading to the correction
$$\begin{array}{cc}\hfill {\mathit{r}}_{i}& ={\stackrel{\sim}{\mathit{r}}}_{i}+\frac{\delta {t}_{i}^{4}}{24}{\mathit{s}}_{i}+\frac{\delta {t}_{i}^{5}}{120}{\mathit{c}}_{i}\hfill \end{array}$$(A.48)
$$\begin{array}{cc}\hfill {\mathit{v}}_{i},& ={\stackrel{\sim}{\mathit{v}}}_{i}+\frac{\delta {t}_{i}^{4}}{6}{\mathit{s}}_{i}+\frac{\delta {t}_{i}^{5}}{24}{\mathit{c}}_{i}.\hfill \end{array}$$(A.49)
The new time step is (Aarseth 1985)
$$\begin{array}{c}\hfill \delta {t}_{i}=\eta \sqrt{\frac{{\stackrel{\sim}{\mathit{a}}}_{i}{\mathit{s}}_{i}+{\stackrel{\sim}{\mathit{j}}}_{i}{}^{2}}{{\stackrel{\sim}{\mathit{j}}}_{i}{\mathit{c}}_{i}+{\mathit{s}}_{i}{}^{2}}}\end{array}$$(A.50)
(where η is an accuracy parameter), rounded down to a power of 2. We note in passing that ph4 departs from kira in the use of a novel and more efficient blockscheduling algorithm, which reorders the block step (t_{next}) list and hence determines the next time step in 𝒪(1) steps per i particle.
The specialization of ph4 that causes it to perform well, that is, the removal of most of the complicating physics, in principle also limits its range of applicability. The code can be coupled to other physics solvers or other gravity modules on different scales through the AMUSE framework, see section A.3.3.
A.3.2. Hermite_GRX
We implemented the EIH equations of motion in regularized and nonregularized forms in standard C++11 (ISO 1998). Our implementation includes the correction terms for EIH equations of motion, and we include the calculation for the energy and linear momentum for validation.
The code is parallelized using threads, but not with MPI, and it is not GPUenabled. We implemented the postNewtonian cross terms using two particle sets: one set of N particles that are affected by general relativity, and one set of n particles that is purely Newtonian. In principle, all particles can be considered relativistic, or all can be Newtonian. In the first case, the code performance scales with N^{3}, otherwise with the usual N^{2}. In general, the code scales ∝N^{3} + nN^{2} + n^{2}. In the case of a galactic nucleus, N = 1 represents the supermassive black hole and the rest of the particles n for the other stars.
We implemented

1PN EIH for the full EIH equations of motion, resulting in an 𝒪(N^{3}) time complexity.

1PN Pairwise, which neglects the acceleration dependence of the velocity in the EIH equations of motion, resulting in a 𝒪(N^{2}) algorithm.

1PN GC Crossterms for the integration of supermassive black holes in galactic nuclei in which one massive object includes the cross terms with the lowmass objects, but the lowmass objects are not relativistic.
In Hermite_GRX, we implemented several numerical schemes, which include

Hermite, for the standard Hermite predictorcorrector integrator with variable but shared time step (Makino 1991).

SymmetrizedHermite, which is a timesymmetrized version of the Hermite integrator (Hut et al. 1995).

RegularizedHermite, which includes regularized close approaches for pairs of particles using KS regularization. This implementation still uses the standard Hermite scheme for time integrations.

SymmetrizedRegularizedHermite, which adopts RegularizedHermite, but with symmetrized time step (Funato et al. 1996).
A.3.3. Implementation in the Astrophysics Multipurpose Software Environment
The Astrophysics Multipurpose Software Environment is a numerical framework for multiscale and multiphysics simulations (Portegies Zwart & McMillan 2018). AMUSE uses numerical implementations for a wide variety of physical processes, including gas dynamics, star formation, stellar evolution, gravitational dynamics, circumstellar disk evolution, and radiative transport. Other physical processes, such as stellar and binary evolution, the Galactic tidal field or hydrodynamical processes can be accommodated through the AMUSE framework. One complication in multibody dynamics is the formation of substructures, such as binaries and triples.
The treatment of such local condensations, but also of close encounters in unsoftened systems, is a much lowerlevel issue and requires special treatment. AMUSE does contain Nbody modules that include specialized treatment of close encounters, but most do not, and because the guiding principle behind AMUSE is to separate functionality as much as possible, ph4 and ph4 like many other modules, relies on the multiples module to manage close encounters and any longlived binary and multiple systems that arise. In multiples, a binary or stable multiple, once identified, is treated as an unperturbed object, possibly with the inclusion of secular internal evolution terms until it has a close encounter with another object in the system. At that point, the interaction is treated as a fewbody scattering, which eventually results in the creation of new stable objects that are then reinserted into the ph4 or ph4 integration. The multiples module is described in more detail in (Portegies Zwart & McMillan 2018, Sect. 4.5).
Hermite_GRX can also be combined with multiples through the AMUSE framework, but due to the local regularization strategies using quaternions (see section A.2), this is not always necessary.
In AMUSE, at least two independently developed implementations for each of the domainspecific solvers are available. This Noah’s Ark approach (Portegies Zwart et al. 2009) allows the user to swap one simulation code for another without any further changes to the runtime environment (codes, scripts, or underlying hardware) and without the need to recompile. It is nonintrusive in the sense that underlying numerical implementations do not require any modifications or recompilation.
The environment is tuned for running on highperformance architectures. It also includes support for GPU and massive taskbased parallelism (using messagepassing parallelism or open multiprocessing parallelism).
A.3.4. Stopping conditions
In the AMUSE framework, codes are interfaced to allow the generation of a homogeneous and selfconsistent simulation environment for performing multiscale and multiphysics simulations. Many of the codes in AMUSE are not build for this purpose, but for operating within a specific domain and parameter range. Due to the interaction with other codes, the underlying simulation engines (called community codes) may be forced to operate outside their usual domain range. If this goes well, the particular community code crashes with the appropriate memory core dump. In AMUSE, however, such an interrupt is caught by the framework without resulting in a crash.
In a general simulation environment, this would be the moment another code takes over to continue the calculation in a different part of parameter space, another temporal or spatial domain, or including different physics. In AMUSE, this problem is addressed by introducing stopping conditions to interrupt the particular simulation domain that the underlying code is unsuited to handle.
In ph4 and Hermite_GRX, we implemented three different stopping conditions. When any of these interrupts is initiated, the code returns control to the AMUSE framework, where the event can be handeled appropriately. By default, the stopping conditions are turned off. The three stopping conditions are described below.

collision_detection All particles have a property called radius, which is used to check for collisions between pairs of particles. At each integration step, we check whether two stars approach each other within the sum of their radii. For this, we assume that the particles within that internal codetime step move in a straight line,
$$\begin{array}{c}\hfill \mathit{r}(s)={\mathit{r}}_{b}+({\mathit{r}}_{e}{\mathit{r}}_{b}).\end{array}$$(A.51)
Here r_{b} and r_{e} are the positions at the beginning and end of the current time step, respectively. If the relative distance between two particles is smaller than the sum of their radii for some value of s, the interrupt collision_detection is initiated.

Wall_clock_time_out_detection The interrupt is initiated when the code takes too long in terms of the wallclock time to evolve to a specified model time.

maximum_number_of_integration_steps_detection. This is initiated when the evolution to the required model time takes more than a predetermined number of time steps to reach the end time of the simulation.
A.3.5. Parallelization by messagepassing in Hermite_GRX
Solving the EIH equations of motion scales ∝N^{3}, making it a rather slow Nbody code. In addition to regularization, we speed the code up by parallelizing it. We do this by multithreading with C++11. Communication between threads can be done in shared memory, which is implemented through the standard library.
A.3.6. Parallelization by messagepassing in ph4
The most important departure of ph4 from kira is the use of MPI parallelism and GPU acceleration in ph4 to increase performance on parallel architectures and GPUsupported systems. Although kira in Starlab was designed to run on a single processor, it can operate in parallel using MPI (Portegies Zwart et al. 2008), but its normal mode of operation is with specialpurpose GRAPE6 (Makino & Taiji 1996; Makino 1998) or GPU through the sapporo library.
Fig. A.1. Integration of an N = 10 Plummer sphere with a virial radius of 1 mpc and a total mass of 10^{7} equalmass black holes for 200 years. 
A significant difference from the kira formulation is that in ph4, all global (jdata) calculations are implemented as MPI parallel tasks with an arbitrary number of workers, and each worker is optionally GPU accelerated using the sapporo2 library (Gaburov et al. 2009; Bédorf et al. 2015). This allows an arbitrary number of GPUs to be configured per MPI worker (subject to the availability of hardware), with all options settable from the Pythonlevel interface to ph4. The resulting boost in speed makes the GPUaccelerated version of ph4 one of the bestperforming direct (N^{2}) Nbody codes in the AMUSE suite (see Fig 10 in Portegies Zwart et al. (2013)). We illustrate the higher speed of the GPUaccelerated version of ph4 in Fig. 1.
A.3.7. Parallelization by messagepassing in Brutus
In Brutus, the iparallellization scheme is implemented, that is, in the doubleforce loop, the outer forloop is parallelized (Makino 2002). The numbers in arbitraryprecision data type are converted into an array of characters that are subsequently communicated through MPI. When it is received, the data type is converted back into arbitraryprecision variables.
A.3.8. Speed of light
Hermite_GRX, ph4, and Brutus operate internally in dimensionless Nbody units, for which G = 1, the total mass M = 1, and the virial radius R = 1. When we scale to physical units, like in a 10^{6} M_{⊙} star cluster with a 1 pc virial radius, we have to introduce scaling between the physical units and the Nbody units. In AMUSE, this is done with a converter (Portegies Zwart & McMillan 2018).
In Hermite_GRX, we use the relative speed of light parameter ζ in terms of the mean velocity in Nbody units ($\frac{1}{2}\sqrt{2}$), and size
$$\begin{array}{c}\hfill v=\sqrt{\frac{\mathit{Gm}}{r}}.\end{array}$$(A.52)
The speed of light in physical units also has to be converted into Nbody units. In Hermite_GRX, we do this by defining the relative speed of light, or the relativisticality of the conditions, as
$$\begin{array}{c}\hfill \zeta =\frac{1}{2}\sqrt{2}/v.\end{array}$$(A.53)
For m = 10^{6} M_{⊙} and r = 1 au, we obtain ζ = 𝒪(2.4 · 10^{−5}), and for r = 1 pc, we obtain ζ = 𝒪(0.01). Here the system becomes Newtonian for ζ → 1, and lower values of ζ mean a more relativistic system. The assumption of the Tailorseries expansion tends to break down for ζ ≳ 0.3 (see sect. 4.2 or better, read (Will 2011)).
In the main paper we express the speed of light in terms of the velocity dispersion v/c, which is more natural. The conversion from ζ to v/c is ${v}_{\mathrm{nbody}}=1/\sqrt{2}$:
$$\begin{array}{c}\hfill v/c={v}_{\mathrm{nbody}}/{c}_{\mathrm{nbody}}={v}_{\mathrm{nbody}}/(\zeta c).\end{array}$$(A.54)
A.3.9. Example
In Listing 2 and 1, we showed a rudimentary AMUSE script for calculating the evolution of N black holes to 1PN EIH equations of motion (including the cross terms). Listing 2 showed how to start the code, and Listing 1 showed the main event loop. Here we adopted the symmetrized regularized Hermite integrator with the full EIH equations of motion to first order. We adopted a radius of the particles of 10GM/c^{2}.
def gravity (bodies, converter, integrator, dt_param =0.1, t_end = 200 units . yr): gravity = HermitePN (converter) gp = gravity . parameters gp . dt_param = dt_param gp . perturbation = ’1PN_EIH ’ gp . integrator = integrator gp . light_speed = constants . c gp . num_threads = 1 gravity. particles . add _ particles (bodies) sc = gravity . stopping_conditions . \ collision _ detection sc . enable () Etot_init = gravity . kinetic_energy +\ gravity . potential_energy from_gravity = gravity . particles . \ new_channel_to (bodies) dt = 0.001 * t_end while gravity. get_time () < t_end : time = gravity . model_time + dt gravity . evolve_model (time) from_gravity . copy () Ekin = gravity . kinetic_energy Epot = gravity . potential_energy Etot = Ekin + Epot dE = (Etot_init − Etot)/Etot print ("T =", gravity . get_time (), "M =", bodies . mass . sum (), "E = ", Etot, "Q = ", −Ekin / Epot) print ("dE =", dE) if sc . is_set (): print ("Collision detected =", len (sc . particles (0))) resolve_collision (sc, gravity, bodies) gravity . stop ()
mass = 1. e+6  units .MSun size = 0.001 units . pc conv = nbody_system . nbody_to_si (mass, size) bodies = new_plummer_model (10, conv) bodies. radius = Schartzschield (bodies . mass) integrator = ’ SymmetrizedRegularizedHermite ’ gravity (bodies, integrator, converter, dt_param = o . dt_param, t_end = o . t_end)
Appendix B: Validation of the code
B.1. Twobody systems
B.1.1. Integrator performance
The Newtonian solution for the equations of motion for two particles was described by Kepler (1609). In the absence of general relativity, the solution is static, with the exception of the mean anomaly. As a first test, we check the conservation of these theoretically conserved Keplerian elements for one orbit.
We varied the initial eccentricity e_{0} and timestep parameter η and adopted masses of M = 10^{6}M_{⊙}, m = 50M_{⊙}, an initial semimajor axis a_{0} = 1 mpc, and an initial eccentrity e_{0} ∈ {0.1, 0.5, 0.9}. For the timestep parameter η ∈ {0.03, 0.01, 0.003, 0.001}, and we integrated for one orbital period (Kepler 1609),
$$\begin{array}{c}\hfill P=2\pi \sqrt{\frac{{a}_{0}^{3}}{G(M+m)}}.\end{array}$$(B.1)
In figs. B.1 and B.2 we show the relative errors in energy and eccentricity in these integrations for the unregularized and the regularized Hermite integrator, respectively. For each calculation, the error in the energy and eccentricity reaches a maximum near pericenter.
Fig. B.1. Relative error in the energy (top row of panels) and eccentricity (bottom row) for integrating a twobody system using the Hermite integrator without postNewtonian terms for initial eccentricities e_{0} ∈ {0.1, 0.5, 0.9} (from left to right), for timestep parameters η ∈ {0.03, 0.01, 0.003, 0.001} (top to bottom in blue, green, red, and orange, respectively). 
Fig. B.2. Relative integration errors in the energy (top row of panels) and eccentricity (bottom row) for one orbit using the regularized Hermite integrator for Newton’s equations of motion for initial eccentricities e_{0} ∈ {0.1, 0.5, 0.9} (left to right panels, respectively), for various time step parameters η ∈ {0.03, 0.01, 0.003, 0.001} (in colour: blue, green, red and orange, respectively). 
The relative error in the regularized Hermite integration in fig. B.2 is one orders of magnitude smaller than the unregularized integrator error (see Figure B.1). As intended in its design, the regularized Hermite integrator performes equally in terms of conserving energy and angular momentum compared to the nonregularized integrator for loweccentricity orbits, and considerably better for eccentric orbits.
In fig. B.3 we present the secular drift in energy as a function of η for initial eccentricities e_{0} ∈ {0.01, 0.1, 0.5, 0.9, 0.99, 0.999, 0.9999}. The same set of initial conditions were adopted in Hamers et al. (2014). These integrations were performed for t_{end}/P = 10^{3} and t_{end}/P = 3.4 × 10^{5}. As expected for a fourthorder integrator, the integration error scales with 𝒪(η^{4}). The regularized Hermite integrator outperforms the other integrators in terms of energy conservation for high eccentricity. The nonregularized integrators have a secular growth of the energy error. We therefore prefer the regularized integration for the longterm evolution of highly eccentric orbits.
Fig. B.3. Relative energy errors for integration for t_{end}/P = 10^{3} and t_{end}/P = 3.4 × 10^{5} for initial eccentricities e_{0} ∈ {0.01, 0.1, 0.5, 0.9, 0.99, 0.999, 0.9999} as a function of the timestep parameter η. The various integrators are indicated with colors (see legend). 
B.1.2. PostNewtonian corrections
General relativity changes the dynamics of astronomical systems. This results in variations in the evolution of the orbital elements for twobody systems.
The secular evolution in the argument of periastron forms one of the major tests for general relativity. The osculating elements, instantaneous orbital parameters under the influence of a perturbing acceleration, can be derived from the perturbing acceleration using Lagrange’s planetary equations (de Lagrange 1772; Merritt 2013).
We numerically integrated Lagrange’s planetary equations for the firstorder postNewtonian perturbation using Euler’s method (Euler 1760). We decreased the time step until the solution convergenced. For the earlier adopted binary, we used initial eccentrities e_{0} ∈ {0.1, 0.5, 0.9}. The resulting theoretical osculating elements are presented in fig. B.4. The postNewtonian terms become more important for larger eccentricities because near periastron, the relative velocity of the particles becomes large while the distance becomes smaller for higher eccentricities.
Fig. B.4. Osculating elements as a function of time for a binary with a stellar mass 50 M_{⊙} in orbit around a 10^{6} M_{⊙} supermassive black hole in a a = 1 mpc orbit with an initial eccentricity e_{0} ∈ {0.5, 0.7, 0.9}. Only the argument of periastron, ω, shows a secular variation (bottom panel). 
Direct numerical integration of the equations of motion should reproduce these osculating elements, and in particular, the secular change in the argument of periastron. We used the same binary as before and a postNewtonian perturbation to integrate one orbit. For a twobody system, the EIH equations of motion reduce to the pairwise approximation. In fig. B.5 we present the relative error of the osculating elements integrated using the regularized Hermite and compare them to the theoretical prediction. We also show the integration error in total energy, including the postNewtonian energy (eq. 5).
Fig. B.5. Relative errors in semimajor axis, eccentricity, argument of pericenter and postNewtonian energy for a relativistic binary composed of a 50 M_{⊙} star in orbit around a 10^{6} M_{⊙} supermassive black hole in a a = 1 mpc orbit with an initial eccentricity e_{0} ∈ {0.5, 0.7, 0.9}. The initial eccentricity was chosen to be e_{0} ∈ {0.5, 0.7, 0.9}. The integration was done using a regularized Hermite integrator, using a time step parameter of η ∈ {0.03, 0.01, 0.003, 0.001} to show convergence (blue, green, red, and orange, respectively). 
The error in the osculating element remains finite, even for very low values of η. The discrepancy is largest near pericenter and smallest near apocenter. The maximum relative error we observe in fig. B.5 in the first postNewtonian correction is several orders of magnitude smaller than the theoretical predictions, making an implementation error improbable. The discrepancy between the theoretical value and the numerical results may well be caused by roundoff, in particular since the postNewtonian corrections require quite a large number of operations per step and the time step is small. With a timestep parameter η = 10^{−3} and ∼10^{3} operations per postNewtonian evaluation, we expect the roundoff error to grow by some six orders of magnitude over one orbital period. With the adopted 16 mantissa implementation, we then arrive at a mean error of about 𝒪(10^{6}), which is consistent with the observed energy error (bottom row of panels in fig. B.5).
For validation and verification, we recalculated the same initial conditions using ARCHAIN (Mikkola & Merritt 2008). The results are indistinguishable from our implementation. The discrepancy between the numerical Nbody result and the converged semianalytic solution then manifests itself in two independently developed codes. We argue that the conserved energy corresponding to the equations of motions truncated to first postNewtonian order contains some second postNewtonian terms that are ignored. Because these terms have order 𝒪(v^{4}/c^{4}), they tend to be important for higher eccentricity and near pericenter, which is precisely what we observe in our simulations.
We caution about judging the accuracy of an Nbody simulation based on energy conservation alone, in particular when considering the secondorder postNewtonian terms (Portegies Zwart & Boekholt 2018). On the other hand, the secular evolution of the energy does not seem to be affected. For the 1PN terms (including the cross terms) adopted in the main paper, the enery is conserved.
Appendix C: Relativistic von ZeipelLidovKozai effect
There are only a few known solutions to the threebody problem. In addition to several semianalytic solutions to resonant cases, such as we find in periodic braids (Moore 1993; Montgomery 1998), there is also a theory about the general behavior of hierarchical threebody systems. In this section, we focus on the latter, in particular since there is a rich body of literature about the associated phenomena observed in hierarchical triples. We refer to this theory as von ZeipelLidovKozai cycles (von Zeipel 1910; Lidov 1962; Kozai 1962), and there is copious literature about the theory (Efimov & Sidorenko 2020; Hamers 2021), its deeper consequences (de Elía et al. 2019), or the observational aspects (Stephan et al. 2021). Here we adopt the von ZeipelLidovKozai effect of testing the threebody methods for Newtonian and relativistic dynamics.
C.1. Newtonian von ZeipelLidovKozai problem
Numerical and analytical dynamical stability arguments (Georgakarakos 2008) indicate that most triples with comparable masses and mutual distances are dynamically unstable and ultimately decay into a binary and a single star. Counterexamples exist, however, in which the three stars form a stable and periodic braid. Some of these solutions are even stable to secondorder postNewtonian order (Lousto & Nakano 2008). Although dynamically unstable triples are rare, they are of considerable theoretical importance. From an observational perspective, they are also of interest because they lead to relatively highvelocity stars or stellar mergers.
Except for braids, stable triples are always hierarchical in the sense that they can be described as an inner binary and a third body that orbits the center of mass of the inner binary at a distance much larger than the separation of the inner binary, as shown schematically in Figure C.1. Such hierarchical triples are rather common, and Tokovinin (2014) estimated their fraction among solartype stars in the solar neighborhood ∼13 %.
Fig. C.1. Schematic image of a hierarchical binary, consisting of the inner binary with masses m_{1} and m_{2}, orbiting each other with a semimajor axis a_{1} and an eccentricity e_{1}, orbited by a third mass m_{3}, with semimajor axis a_{2} and eccentricity e_{2}. This figure is not to scale, as typically a_{2} ≫ a_{1}. 
An important aspect in the dynamics of hierarchical triples is the periodic exchange of orbital angular momentum between the inner and outer binary. Such coupling is only effective when the relative inclination between the two orbital planes of the inner and outer orbits exceed some critical value
$$\begin{array}{c}\hfill {i}_{\mathrm{rel},\mathrm{crit}}=arccos\left(\sqrt{{\textstyle \frac{3}{5}}}\right)\approx 39.{2}^{\xb0}.\end{array}$$(C.1)
In this case, the inner binary and the relative inclination evolve periodically. To first nonzero order, this periodicity conserves
$$\begin{array}{c}\hfill {L}_{z}\propto \sqrt{1{e}_{1}^{2}}cos({i}_{\mathrm{rel}})=\text{c}onstant.\end{array}$$(C.2)
During such von ZeipelLidovKozai cycles, the orbital energies remain constant, which leads to constant semimajor axes. During such a cycle, the eccentricity of the inner binary can reach values as high as e_{1} ∼ 1 − 10^{−6}, as we show in Figures C.3 and C.4. Such highly eccentric orbits are easily subject to tidal effects or the emission of gravitational waves, and could lead to stellar collisions (Antognini & Thompson 2016). Fortunately, such higheccentricity encounters are not expected to naturally occur in large Nbody systems, except in the presence of hierarchical multiple subsystems. If such high eccentricities are relevant, the entire integration scheme, the Taylor expansion adopted for the postNewtonian terms, and the possibility of tidal effect should all be reconsidered.
One great advantage of von ZeipelLidovKozai cycles is the possibility of deriving the secular evolution analytically by averaging over the inner and outer orbits. Here we assume that the orbital parameters vary slowly compared to the outer orbit: the timescale on which the orbital angular momentum of the inner binary varies is small compared to the inner orbital period (Antonini et al. 2014). Direct numerical integration of the equations of motion of hierarchical triples remains important for validating the underlying assumption on the system’s hierarchy.
In Figure C.2 we present the result of several such numerical integration of some von ZeipelLidovLidov cycles together with the relation between e_{1} and i_{rel}. The lowestorder approximation for von ZeipelLidovKozai cycles is the result of the quadrupole term in the multipole expansion that is used in the derivation. The corresponding timescale is approximately (Naoz et al. 2013b,)
$$\begin{array}{c}\hfill {t}_{\mathrm{quad}}^{\mathrm{Newton}}\sim \frac{2\pi {a}_{2}^{3}{(1{e}_{2}^{2})}^{{\textstyle \frac{3}{2}}}\sqrt{{m}_{1}+{m}_{2}}}{{a}_{1}^{{\textstyle \frac{3}{2}}}{m}_{3}\sqrt{G}},\end{array}$$(C.3)
Fig. C.2. A few KozaiLidov cycles with their distinctive signature of higheccentricity (low 1 − e_{1}) spikes. On the right, the eccentricity and relative inclination are plotted against each other from t = 0 to t = 10000 yr, including the analytical prediction (based on the initial conditions and conservation of linear momentum in Equation C.2). The thickness of this line is due to numerical integration errors. Initial condition is a binary of masses m_{1} = 1 M_{Jup} and m_{2} = 1M_{⊙}, semimajor axis a_{1} = 0.005 AU, and eccentricity e_{1} = 0.001, orbited by a third body of mass m_{3} = 10^{6}M_{⊙} at semimajor axis a_{2} = 51.4 AU, eccentricity e_{2} = 0.7, and relative inclination i_{rel} = 95°. Integration was done using a regularized Hermite integrator with a timestep parameter of η = 0.01. 
in which the eccentricity reaches a maximum value of
$$\begin{array}{c}\hfill {e}_{1,max}=\sqrt{1{\textstyle \frac{5}{3}}{cos}^{2}({i}_{\mathrm{tot}}}).\end{array}$$(C.4)
The latter expression is only valid in the limit in which the octopole terms vanish, the test particle quadrupole order limit (Naoz et al. 2013a), in which a_{2} ≫ a_{1}. When the octupolelevel terms become important, they can be seen as a modulation of von ZeipelLidovKozai cycles. The importance of the octupolelevel variations can be quantified by considering the ratio of the octupole to quadrupolelevel coefficients,
$$\begin{array}{c}\hfill \frac{{C}_{3}}{{C}_{2}}=\frac{15}{4}\frac{{\u03f5}_{M}}{{e}_{2}}.\end{array}$$(C.5)
Here C_{2} and C_{3} are the quadrupole and octupolelevel coefficients given by Naoz et al. (2013a), and ϵ_{M} is the relative importance of the octupolelevel term in the secularized Hamiltonian,
$$\begin{array}{c}\hfill {\u03f5}_{M}=\left(\frac{{m}_{1}{m}_{2}}{{m}_{1}+{m}_{2}}\right)\left(\frac{{a}_{1}}{{a}_{2}}\right)\left(\frac{{e}_{2}}{1{e}_{2}^{2}}\right).\end{array}$$(C.6)
This suggests that octupolelevel variations are important for eccentric innerbinaries with highmass components. We note here that ϵ_{M} is independent of the mass of the third (outer) body m_{3}, but depends on its orbital parameters a_{2} and e_{2}. The timescale of the octupole variation can be defined in a similar fashion,
$$\begin{array}{c}\hfill {t}_{\mathrm{oct}}^{\mathrm{Newton}}\sim \frac{4}{15}{\u03f5}_{M}^{1}{t}_{\mathrm{quad}}^{\mathrm{Newton}}.\end{array}$$(C.7)
Naoz et al. (2013a) demonstrated that this octupole variation can have consequences on the maximum eccentricity reached during the evolution of the system. The induction of variations in the relative inclination i_{tot} over time can lead to a flip in the inner binary’s orbit. The maximum eccentricity is reached at the moment the flip occurs (see eq. C.4).
We demonstrated in sect. B that the unregularized Hermite scheme is prone to introducing integration errors for highly eccentric orbits. This may pose a problem while integrating triples that are subject to von ZeipelLidovKozai cycles, for which the inner binary can become highly eccentric. We illustrated this in fig. C.3, wherethe inner binary reaches an eccentricity in excess of 1 − 10^{−6} at t ∼ 4 Myr, leading to a relative integration error of about eight orders of magnitude larger than when integrating a circular orbit. The result of regularized Hermite scheme for the same triple is presented in fig. C.4, showing a considerably better conservation of energy, and it is also faster.
Fig. C.3. Kozai cycles integrated using the Hermite integrator with a timestep parameter of η = 0.001. Initial condition was a binary with masses m_{1} = 1M_{Jup} and m_{2} = 1M_{⊙} with semimajor axis a_{1} = 6 AU and eccentricity e_{1} = 0.001 orbited by a third body of mass m_{3} = 40M_{Jup} at a distance a_{2} = 100 AU and eccentricity e_{2} = 0.6. The initial relative inclination was i_{tot} = 65°. These initial conditions are identical to Figure 3 of Naoz et al. (2013a). The integration errors in both the energy and the semimajor axis of the inner binary rapidly rise at t ∼ 4 Myr, where the inner binary reaches maximum eccentricity. 
Fig. C.4. Kozai cycles integrated using the regularized Hermite integrator with a timestep parameter of η = 0.003. Initial conditions are identical to those in Figure C.3, using a similar wallclock run time. The large integration errors at high eccentricities are now absent, as regularization is applied to the inner binary. 
C.2. Relativistic von ZeipelLidovKozai problem
Triple star systems that are subject to von ZeipelLidovKozai cycles may be affected by general relativity, in particular, when this leads to highly eccentric orbits. Since these cycles are the result of secular resonances between an inner and an outer orbit, variations in the parameters of the inner binary on a timescale similar to the secular resonance tend to quench the effect and reduce the extremes in the von ZeipelLidovKozai cycles.
The timescale on which the inner binary orbit is affected by general relativity is (Naoz et al. 2013b,)
$$\begin{array}{c}\hfill {t}_{\mathrm{i}}^{1\mathrm{PN}}\sim 2\pi \frac{{a}_{1}^{{\textstyle \frac{5}{2}}}{c}^{2}(1{e}_{1}^{2})}{3G{({m}_{1}+{m}_{2})}^{{\textstyle \frac{3}{2}}}}.\end{array}$$(C.8)
We define the relative (dimensionless) parameter ℛ as the ratio between the 1PN terms and Newtonian quadrupole timescales for a circular inner orbit,
$$\begin{array}{c}\hfill \mathcal{R}={\frac{{t}_{\mathrm{i}}^{1\mathrm{PN}}}{{t}_{\mathrm{quad}}^{\mathrm{Newton}}}}_{{e}_{1}=0}=\frac{1}{3}\frac{{({a}_{1}/{\mathcal{R}}_{1})}^{4}}{{({a}_{2}/{\mathcal{R}}_{3})}^{3}}\frac{1}{{({m}_{3}/{m}_{1})}^{2}{(1{e}_{2})}^{{\textstyle \frac{3}{2}}}}.\end{array}$$(C.9)
Here ℛ_{1} and ℛ_{2} are the gravitational radii of the inner binary and the outer orbiting tertiary body. They are given by
$$\begin{array}{c}\hfill {\mathcal{R}}_{1}=\frac{G({m}_{1}+{m}_{2})}{{c}^{2}},\end{array}$$(C.10a)
and
$$\begin{array}{c}\hfill {\mathcal{R}}_{2}=\frac{G{m}_{3}}{{c}^{2}},\end{array}$$(C.11a)
respectively.
A maximum in the eccentricity of the inner orbit is induced when the timescale for which the inner orbit evolves due to general relativity is on the same order as one von ZeipelLidovKozai cycle due to classical Newtonian resonance. The criteria for this to happen are ℛ ∼ 1 and m_{3} ≫ m_{1} > m_{2} (Naoz et al. 2013b,).
In fig. C.5 we compare Hermite_GRX with the results presented in Naoz et al. (2013b, see their fig.5). For the initial conditions, we used a star with planet that orbit a supermassive black hole. The inner binary, the star and planet, have masses m_{1} = 1M_{⊙} and m_{2} = 0.001M_{⊙}, with semimajor axis a_{1} = 10^{5} ℛ_{1}, eccentricity e_{1} = 0.001, and inclination i_{rel} = 65° with respect to the outer orbit. The third body is a supermassive black hole with mass m_{3} = 10^{6}M_{⊙} and orbits the inner binary with a semimajor axis a_{2} and eccentricity e_{2} = 0.7 for various values of ℛ. The system is evolved to ${t}_{\mathrm{end}}=s{t}_{\mathrm{quad}}^{\mathrm{Newton}}$ for three value of s, where we adopted a timestep parameter η = 0.003 for s = 10 and s = 100, and η = 0.002 for s = 300. The maximum eccentricity reached during this time interval for various values of ℛ is plotted in fig. C.5. Some minor deviation from the theoretical curve, in particular for s = 10 is caused by our incomplete sampling because we only determined the eccentricity of the inner orbit near apocenter. In addition, the distance between the third body to the inner binary also introduces variations in the eccentricity of the inner orbit.
Fig. C.5. Maximum eccentricity reached after ${t}_{\mathrm{end}}=N{t}_{\mathrm{quad}}^{\mathrm{Newton}}$, where N ∈ {10, 100, 300} are denoted by the colors. The initial conditions are m_{1} = 1M_{⊙}, m_{2} = 0.001M_{⊙}, m_{3} = 10^{6}M_{⊙}, a_{1} = 10^{5} ℛ_{1}, i_{tot} = 65°, e_{1} = 0.001, e_{2} = 0.7, and a_{2} was varied to simulate different ℛ. An eccentricity excitation around ℛ = 0.65 is evident. Simulations were performed using a regularized Hermite integrator with a timestep parameter of η = 0.003 for N = 10 and N = 100 and η = 0.002 for N = 300. Each data point is an independent simulation that took a wallclock time of ∼20 min for N = 10, ∼3 hours for N = 100, and ∼12 hours for N = 300. 
In fig. C.5 we show the resonantlike eccentricity excitation discussed in Naoz et al. (2013b,). Our calculations are integrated directly, whereas Naoz et al. (2013b,) conducted an orbitaveraged integration, using test particles with 1PN including terms only up to $\mathcal{O}({a}_{1}^{2})$ and $\mathcal{O}({a}_{2}^{2})$ and the term, whereas we perform a direct Nbody integration of the full EIH equations of motion to 1PN with finite masses.
The first peak in the resonant structure of the eccentricity in fig. C.5 is shifted to ℛ ∼ 0.65 with respect to to ℛ ∼ 0.55 in fig. 5 of Naoz et al. (2013b,). The calculations performed by Naoz et al. (2013b,) adopted orbit averaging, which gives a considerable speedup compared to direct integration. On the other hand, however, this may lead to missing the collision because the maximum eccentricity is reached in a time interval that is shorter than the orbital period of the outer binary. Orbit averaging over the outer orbit then lacks the resolution to resolve the maximum eccentricity in the inner orbit, whereas in the direct Nbody integration presented in fig. C.5, we do resolve the evolution of the eccentricity of the inner orbit.
Naoz et al. (2013b,) further discussed the possibility of orbital flips when including relativistic effects, even in the absence of considerable variations at the Newtonian octupole moments, which form the usual cause of orbital flips. In their fig. 7, they presented a specific case for an inner binary with m_{1} = 10M_{⊙}, m_{2} = 8M_{⊙}, a_{1} = 10 AU, and e_{1} = 0.001 that is orbited by a m_{3} = 30M_{⊙} tertiary body with semimajor axis a_{2} = 502 AU and eccentricity e_{2} = 0.7 and inclined by i_{rel} = 94° to the plane of the inner binary. According to eq. C.6, the high mass ratio of the inner binary suppresses the octupolelevel (Newtonian) effects. By integrating these initial conditions, including the EIH equations of motion using the regularized Hermite integrator with a timestep parameter η = 0.0003, we do not observe such a orbital flip, as we show in fig. C.6.
Fig. C.6. KozaiLidov cycles for identical initial conditions as in Naoz et al. (2013b,), showing no orbital flips, in contrast to the orbital flips shown in that paper. The initial conditions of this direct numerical integration consist of a similarmass inner binary (m_{1} = 10M_{⊙}, m_{2} = 8M_{⊙}, a_{1} = 10 AU, e_{1} = 0.001) orbited by a third body (m_{3} = 30M_{⊙}, a_{2} = 502 AU, e_{2} = 0.7), inclined by i_{rel} = 94° relative to the inner binary. Integration was done using Newtonian (blue line) and the EIH equations of motion (green striped line) using a regularized Hermite integrator with a timestep parameter η = 0.0003. 
The origin of this discrepancy is not so clear. The Newtonian case does not show orbital flips, and we see no direct argument for the presence of orbital flips when adopting the EIH equations of motion. On the other hand, when integrating the EIH equations of motion, we acquire a considerable error in the total energy when the binary reaches its highest eccentricity. This is caused by the large perturbation of the postNewtonian terms when the inner binary reaches pericenter. As a result, this system is hard to integrate numerically. Our integration with a timestep parameter η = 0.0003 integrated for 1 Myr took ∼5 days on a regular workstation and reached a minimum relative inclination of i_{rel, min} = 93°.
Overall, our regularized Hermite integrator performs well and gives results that are consistent with previous calculations. Discrepancies with secular evolution calculations can be understood from the lack of innerorbit resolution in the latter. We therefore see no reason to doubt our implementation of the regularized and postNewtonian terms.
All Tables
All Figures
Fig. 1. Scaling of the various integration methods as a function of the number of particles (N). The classic Newton integrations scale as ∝N^{2}, as indicated with the solid and dashed dark blue lines. The pairwise first expansion scales similarly, but tends to be slightly slower than the pure Newton expansion. The EinsteinInfeld Hoffmann equations to first order scale ∝N^{3}, making large calculations that include the crossterms unpractical. The scaling presented for Brutus is based on the calculations presented here. The bullet points indicate the mean timescale for acquiring a converged solution, and the line pointed upward ends at the single most expensive calculation in our sample of simulations for that particular N. For ph4, we included the regular implementation as well as the GPUenabled version (to the right), running on an Intel Xeon CPU E5620 operating at 2.40 GHz and NVIDIA G96 (Quadro FX580), running on a generic 64bit Ubuntu Linux kernel 2.6.3532. 

In the text 
Fig. 2. Orbital evolution of two black holes of masses 10^{6} M_{⊙} with an initial separation of 0.819 au integrated for half a day for three integrations, pure Newton expansion, Newton with expansions to first order, and the full firstorder EinsteinInfeldHoffmann equations. The first two are precisely on top of each other, and we plotted the firstorder pairwise solution last. The last two solutions are identical because the crossterms do not lead to deviations from the firstorder expansions. The postNewton orbits are not closed, as in the Newtonian case (green). No separate scaling of v/c is applied here because the system is initialized in physical units. 

In the text 
Fig. 3. Orbital evolution of four of 10^{6} M_{⊙} black holes with the same integrators as in Fig. 2. 

In the text 
Fig. 4. Orbital evolution of 16 of 10^{6} M_{⊙} black holes with the same integrators as in Fig. 2. 

In the text 
Fig. 5. Distribution of phasespace distances in a cluster with 1024 particles. Units are dimensionless Nbody units (Hénon 1971). The gray shaded region indicates the initial conditions in a virialized unit cube. One particle, indicated with the black bullet point, is displaced by 10^{−7} along the Cartesian x coordinate. The final conditions (at t = 10) of the unperturbed particles are represented with the bullet points. The color and size of the points represents the phasespace distance measured over the duration of the simulation (10 Nbody time units), and ranges over about four orders of magnitude. The majority of objects experience considerable change in their orbits, but some are hardly perturbed. Calculations were performed using Brutus until convergence to 3 decimal places, which requires a tolerance of τ = 10^{−40}. 

In the text 
Fig. 6. Distribution of phasespace distances for individual particles δ_{i} in the simulations with N = 4 (red) and those with N = 1024 (blue) after integrating for t = 10 Nbody time units. The data for N = 4 are the result of 200 runs. For N = 1024, we adopted the run used in Fig. 5. Calculations were performed using Brutus until the solution was converged. 

In the text 
Fig. 7. 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 Goodman et al. (1993). The different symbols and colors represent different calculations (see legend). The vertical bars, plotted for Newton’s Hermite only, show the rootmeansquare of the dispersion in the series of solutions. The error bars in the results obtained with Brutus are statistically indistinguishable from the presented bars. 

In the text 
Fig. 8. Distribution of phasespace distances for individual particles in the simulations with N = 4 (red) and for N = 64 (blue) after integrating for t = 10 Nbody time units. The data for 200 runs were used for each histogram. Calculations were performed using Hermite_GRX using v/c = 0.01 (top panel) and for v/c = 0.002 (bottom panel). 

In the text 
Fig. 9. Lyapunov timescale as a function of v/c for N = 4 (green) and N = 64 (blue). The Newtonian case (run with ph4) is presented as arrows in orange. The vertical bars, only for the green points, indicate the dispersion in the simulation results. The short horizontal dotted green line indicates the lowest value for the Lyapunov timescale reached for v/c = 10^{−3} for N = 4. 

In the text 
Fig. 10. Lyapunov timescale as a function of N for Goodman et al. (1993) (red bullets) compared to the various relativistic solutions. In blue we present the solutions using 1PN pairwise terms with a scaling to the speed of light of v/c = 0.010 and twice this value. For the EIH equations of motion we show the case for v/c = 0.010 in green, twice and four times this value. The top blue curve fits t_{λ} ≃ 0.255 + 13.43e^{−0.371N} for v/c = 0.005, and t_{λ} ≃ 0.051 + 5.50e^{−0.115N} for v/c = 0.010. 

In the text 
Fig. 11. Estimate for the Lyapunov timescale for the pairwise 1PN terms as a function of the full EIH equations of motion for three choices of N, 4 (dark blue), 16 (aquamarine), and 64 (green). Both, the pairwise 1PN calculation and the one including the cross terms, are represented as fraction of the Newtonian solution. All relativistic simulations adopt v/c = 0.010. 

In the text 
Fig. 12. Relative importance in the apsidal motion of the 1PN terms in comparison to the 2PN terms. 

In the text 
Fig. 13. Estimated Lyapunov timescale for the Newtonian case with particles distributed in a Plummer sphere (blue) and a King model (yellow). In addition, we show results for the King model, but for the EIH equations of motion with v/c = 0.010 (green). Here the xaxis is in ln(ln(N)). 

In the text 
Fig. A.1. Integration of an N = 10 Plummer sphere with a virial radius of 1 mpc and a total mass of 10^{7} equalmass black holes for 200 years. 

In the text 
Fig. B.1. Relative error in the energy (top row of panels) and eccentricity (bottom row) for integrating a twobody system using the Hermite integrator without postNewtonian terms for initial eccentricities e_{0} ∈ {0.1, 0.5, 0.9} (from left to right), for timestep parameters η ∈ {0.03, 0.01, 0.003, 0.001} (top to bottom in blue, green, red, and orange, respectively). 

In the text 
Fig. B.2. Relative integration errors in the energy (top row of panels) and eccentricity (bottom row) for one orbit using the regularized Hermite integrator for Newton’s equations of motion for initial eccentricities e_{0} ∈ {0.1, 0.5, 0.9} (left to right panels, respectively), for various time step parameters η ∈ {0.03, 0.01, 0.003, 0.001} (in colour: blue, green, red and orange, respectively). 

In the text 
Fig. B.3. Relative energy errors for integration for t_{end}/P = 10^{3} and t_{end}/P = 3.4 × 10^{5} for initial eccentricities e_{0} ∈ {0.01, 0.1, 0.5, 0.9, 0.99, 0.999, 0.9999} as a function of the timestep parameter η. The various integrators are indicated with colors (see legend). 

In the text 
Fig. B.4. Osculating elements as a function of time for a binary with a stellar mass 50 M_{⊙} in orbit around a 10^{6} M_{⊙} supermassive black hole in a a = 1 mpc orbit with an initial eccentricity e_{0} ∈ {0.5, 0.7, 0.9}. Only the argument of periastron, ω, shows a secular variation (bottom panel). 

In the text 
Fig. B.5. Relative errors in semimajor axis, eccentricity, argument of pericenter and postNewtonian energy for a relativistic binary composed of a 50 M_{⊙} star in orbit around a 10^{6} M_{⊙} supermassive black hole in a a = 1 mpc orbit with an initial eccentricity e_{0} ∈ {0.5, 0.7, 0.9}. The initial eccentricity was chosen to be e_{0} ∈ {0.5, 0.7, 0.9}. The integration was done using a regularized Hermite integrator, using a time step parameter of η ∈ {0.03, 0.01, 0.003, 0.001} to show convergence (blue, green, red, and orange, respectively). 

In the text 
Fig. C.1. Schematic image of a hierarchical binary, consisting of the inner binary with masses m_{1} and m_{2}, orbiting each other with a semimajor axis a_{1} and an eccentricity e_{1}, orbited by a third mass m_{3}, with semimajor axis a_{2} and eccentricity e_{2}. This figure is not to scale, as typically a_{2} ≫ a_{1}. 

In the text 
Fig. C.2. A few KozaiLidov cycles with their distinctive signature of higheccentricity (low 1 − e_{1}) spikes. On the right, the eccentricity and relative inclination are plotted against each other from t = 0 to t = 10000 yr, including the analytical prediction (based on the initial conditions and conservation of linear momentum in Equation C.2). The thickness of this line is due to numerical integration errors. Initial condition is a binary of masses m_{1} = 1 M_{Jup} and m_{2} = 1M_{⊙}, semimajor axis a_{1} = 0.005 AU, and eccentricity e_{1} = 0.001, orbited by a third body of mass m_{3} = 10^{6}M_{⊙} at semimajor axis a_{2} = 51.4 AU, eccentricity e_{2} = 0.7, and relative inclination i_{rel} = 95°. Integration was done using a regularized Hermite integrator with a timestep parameter of η = 0.01. 

In the text 
Fig. C.3. Kozai cycles integrated using the Hermite integrator with a timestep parameter of η = 0.001. Initial condition was a binary with masses m_{1} = 1M_{Jup} and m_{2} = 1M_{⊙} with semimajor axis a_{1} = 6 AU and eccentricity e_{1} = 0.001 orbited by a third body of mass m_{3} = 40M_{Jup} at a distance a_{2} = 100 AU and eccentricity e_{2} = 0.6. The initial relative inclination was i_{tot} = 65°. These initial conditions are identical to Figure 3 of Naoz et al. (2013a). The integration errors in both the energy and the semimajor axis of the inner binary rapidly rise at t ∼ 4 Myr, where the inner binary reaches maximum eccentricity. 

In the text 
Fig. C.4. Kozai cycles integrated using the regularized Hermite integrator with a timestep parameter of η = 0.003. Initial conditions are identical to those in Figure C.3, using a similar wallclock run time. The large integration errors at high eccentricities are now absent, as regularization is applied to the inner binary. 

In the text 
Fig. C.5. Maximum eccentricity reached after ${t}_{\mathrm{end}}=N{t}_{\mathrm{quad}}^{\mathrm{Newton}}$, where N ∈ {10, 100, 300} are denoted by the colors. The initial conditions are m_{1} = 1M_{⊙}, m_{2} = 0.001M_{⊙}, m_{3} = 10^{6}M_{⊙}, a_{1} = 10^{5} ℛ_{1}, i_{tot} = 65°, e_{1} = 0.001, e_{2} = 0.7, and a_{2} was varied to simulate different ℛ. An eccentricity excitation around ℛ = 0.65 is evident. Simulations were performed using a regularized Hermite integrator with a timestep parameter of η = 0.003 for N = 10 and N = 100 and η = 0.002 for N = 300. Each data point is an independent simulation that took a wallclock time of ∼20 min for N = 10, ∼3 hours for N = 100, and ∼12 hours for N = 300. 

In the text 
Fig. C.6. KozaiLidov cycles for identical initial conditions as in Naoz et al. (2013b,), showing no orbital flips, in contrast to the orbital flips shown in that paper. The initial conditions of this direct numerical integration consist of a similarmass inner binary (m_{1} = 10M_{⊙}, m_{2} = 8M_{⊙}, a_{1} = 10 AU, e_{1} = 0.001) orbited by a third body (m_{3} = 30M_{⊙}, a_{2} = 502 AU, e_{2} = 0.7), inclined by i_{rel} = 94° relative to the inner binary. Integration was done using Newtonian (blue line) and the EIH equations of motion (green striped line) using a regularized Hermite integrator with a timestep parameter η = 0.0003. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.