SciPost Phys. Core 6, 016 (2023)

Chaos in the vicinity of a singularity in the three-body problem:
The equilateral triangle experimentin the zero angular momentum limit
Hugo D. Parischewsky1⋆ , Alessandro A. Trani2,3 and Nathan W. C. Leigh1,4
1 Departamento de Astronomía, Facultad de Ciencias Físicas y Matemáticas,
Universidad de Concepción, Chile, Avenida Esteban Iturra s/n Casilla 160-C
2 Department of Earth Science and Astronomy, College of Arts and Sciences,
The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan
3 Okinawa Institute of Science and Technology, 1919-1 Tancha,
Onna-son, Okinawa 904-0495, Japan
4 Department of Astrophysics, American Museum of Natural History,
New York, NY 10024, USA⋆ hugo.parischewsky.z@gmail.com

Abstract
We present numerical simulations of the gravitational three-body problem, in which
three particles lie at rest close to the vertices of an equilateral triangle. In the unperturbed problem, the three particles fall towards the center of mass of the system to forma three-body collision, or singularity, where the particles overlap in space and time. By
perturbing the initial positions of the particles, we are able to study chaos in the vicinity of the singularity. Here we cover both the singular region close to the unperturbed
configuration and the binary-single scattering regime where one side of the triangle isvery short compared to the other two. We make phase space plots to study the regularand ergodic subsets of our simulations and compare them with the outcomes expectedfrom the statistical escape theory of the three-body problem. We further provide fits tothe ergodic subset to characterize the properties of the left-over binaries. We identify
the discrepancy between the statistical theory and the simulations in the regular subset
of interactions, which only exhibits weak chaos. As we decrease the scale of the perturbations in the initial positions, the phase space becomes entirely dominated by regularinteractions, according to our metric for chaos. Finally, we show the effect of general relativity corrections by simulating the same scenario with the inclusion of post-Newtonian
corrections to the equations of motion.
Copyright H. D. Parischewsky et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 26-10-2021
Accepted 17-10-2022
Check forupdatesPublished 21-03-2023
doi:10.21468/SciPostPhysCore.6.1.016

Contents
1 Introduction

2

2 Quantifying Chaos

41

SciPost Phys. Core 6, 016 (2023)

2.1
2.2
2.3

Ergodicity
Lévy flights
Scramble cut

45
6

3 Method
3.1 N-body Integrator
3.2 Initial conditions
3.2.1 Unperturbed initial configuration
3.2.2 Perturbed initial configurations
3.3 Ejection criteria
3.4 Simulations with post-Newtonian corrections

6
6
6
6
7
8
9

4 Theoretical Expectations
4.1 Dividing the phase space into ergodic and regular subsets
4.2 Orbital energy distribution
4.3 Eccentricity distribution
4.4 Escape velocity distribution

9
9
10
10
10

5 Results
5.1 Phase space divided by ejected particle
5.1.1 Full phase space
5.1.2 Zoom-in on the singularity
5.1.3 Zoom-in on the binary-single regime
5.1.4 Ergodic and regular phase space
5.2 Including post-Newtonian corrections and general relativistic effects
5.3 Comparison with theoretical expectations
5.3.1 Distributions of outcome properties
5.3.2 Phase space distribution of outcome properties
5.3.3 Weak chaos and the disagreements between theory and simulations

10
11
11
12
12
13
13
15
15
16
17

6 Discussion & Summary
6.1 Astrophysical implications
6.1.1 Three body binary formation in isotropic clusters
6.1.2 Binary properties

18
18
19
19

References

20

1

Introduction

The gravitational two-body problem was solved analytically by Newton in his PhilosophiæNaturalis Principia Mathematica. His work made it possible to predict exactly the positions andvelocities of two self-gravitating point masses at any point in the future, given their initial
masses, positions and velocities.
After Newton, generations of mathematicians and physicists tried to obtain an equivalently
elegant solution for the three-body problem, without success. It was eventually Poincaré [1]who showed that such a solution would be impossible since the general three-body problem
is an example of chaos in nature.

2

SciPost Phys. Core 6, 016 (2023)

More than a century passed with little further progress, but the introduction of computer
simulations in the late 1900’s by Aarseth and collaborators [2–5] re-vitalised popularity in thethree-body problem, allowing researchers to confront their analytic models directly with computer simulations. This allowed for the development of analytic approximations for regions ofphase space that are regular (i.e., the interactions are prompt and never enter a long-lived resonant state; see [6] and the text below for a more detailed definition). However, only recently
the probabilistic theories of the three-body problem have begun to be explored [7–13].
A prerequisite to having a chaotic three-body interaction is that the system enters a longlived resonant state, where all three particles undergo numerous close approaches before theinteraction ends. If the interaction ends promptly, with the incoming single star interactingwith the initial binary components only once before the interaction ends, the interaction issaid to be regular and analytic solutions can typically be found e.g. [14–16].
There are two main states for the time evolution of a chaotic three-body system in a resonant state. The first state consists of a hierarchy, which is composed of a temporary close
binary system and a temporary single star going on a prolonged excursion with a total energy approaching zero but remaining negative. The time evolution is such that the system
continually breaks apart into such a hierarchy, with the temporary single sometimes going on
long-lived and sometimes short-lived excursions. The second state is such that all three particles are in approximate energy equipartition, leading to a fast and chaotic exchange of energy
and angular momentum e.g. [17, 18]. These "scramble states" are needed for the system tobecome chaotic, such that the particles lose all memory of their initial conditions, recalling
just the total energy, total angular momentum and the particle masses. The time evolution
continues in this way, inter-changing between these two states chaotically. Eventually, if all
bodies are point particles, the interaction ends with one of the particles being ejected with apositive total energy and a finite velocity at spatial infinity [9].
Chaos can be defined as occurring when small perturbations to the initial conditions result in different macroscopic outcomes (e.g., which of the three particles is the one to beejected). The extreme sensitivity to the initial conditions manifests itself through the exponential growth of small perturbations that can exhibit unpredictable and divergent behaviour,yielding completely different outcomes for nearly identical sets of initial conditions [19]. Thisimplies numerical and physical consequences [20] and, as already mentioned, points to theneed to develop a probabilistic theory for chaotic regions of phase space [10, 11, 21, 22].
This extreme sensitivity to the initial conditions can also result in diverging trajectories
through phase space purely due to the accumulation of numerical errors, which act as small
perturbations to the system at each time-step, effectively mimicking and amplifying the effects
of chaos. We depict this behaviour in Figure 1, by showing the trajectories of the particles inposition space for similar sets of initial conditions. Consequently, a high degree of accuracy andprecision is needed over long-time-scale integrations, to ensure that the simulated solutions
are correct, and to remove the concern of numerical errors mimicking the effects of chaos.
The development of sophisticated gravity integrators has proved essential by reducing theaccumulation of numerical errors, and guaranteeing the reliability of orbital integrations overlong periods of times (i.e., many orbital periods) [20,23–25]. Most dynamical systems display
some aspect of chaos, including solar system small bodies e.g. [26–28], small stellar systems
e.g. [6, 10, 11, 29–33], star clusters e.g. [34, 35], galaxies e.g. [36], and so on.
In this paper, our focus is to study chaos in the vicinity of a singularity in the three-bodyproblem. The singularity we consider is a three-body collision at the system centre of mass,
which occurs when the three particles are released from rest, each initially at the vertex of aperfect equilateral triangle. This creates a singularity in the gravitational acceleration, since
the particles overlap in both space and time upon reaching the system centre of mass for thisidealized initial configuration. By perturbing the planar equilateral triangle and releasing the3

SciPost Phys. Core 6, 016 (2023)

particles from rest, they will arrive at the system centre of mass at nearly the same time, butslightly offset. In this way, we can probe the interaction outcomes directly in the vicinity of thesingularity, resolving it down to very small spatial scales. Our experiments are designed to be
entirely planar, with each system composed of three point-particles with equal masses, eachlocated initially at a randomly sampled position for the corresponding vertex for that particle.
This experiment is in the zero angular momentum limit, since there is no angular momentum
for any configuration initially at rest.
In Section 2, we introduce a number of key definitions used throughout this paper to describe chaotic three-body interactions. In Section 3, we present our methods and experimental
set-up, including justification for our choice of gravity integrator TSUNAMI [37] and the issueof regularization, the initial conditions, the number of simulations performed, and so on. In
Section 4 we present the theoretical expressions used to compare analytic theory with our numerical scattering experiments. In Section 5, we present the results of our experiments, witha focus on comparing the simulated distributions for the final binary semi-major axis (i.e.,
the inverse binding energy distribution), eccentricity, single star escaper velocity, and the totalinteraction lifetimes to the theoretical expectations. In Section 6, we summarise our results
and discuss their implications for the three-body problem, chaos in the vicinity of a singularity.
Finally, we discuss the properties of binaries formed via three-body interactions of isolated
single stars in dense isotropic stellar systems.

2

Quantifying Chaos

In this section we introduce several key definitions used throughout the paper to describe
chaotic interactions.
We introduce the concept of ergodicity, its relation to Lévy flights and their combined implications to defining a chaotic interaction in the three body problem. As already discussed,
separating three-body interactions into ergodic and regular subsets is critical for a proper comparison between the simulated data and theoretical predictions. The regular subset tends tocorrespond to prompt interactions, for which analytic methods have already been developed
to quantify the outcome properties as a function of the initial conditions. The ergodic subset
tends to correspond to longer-lived resonant interactions, and require a probabilistic theory toquantify the outcome properties as a function of the initial conditions.

2.1

Ergodicity

When comparing simulations of three-body interactions to theoretical predictions, it is crucial
to first define a criterion for when the interaction formally becomes a chaotic one. This is
typically done by quantifying the amount of time all three particles have comparable energies and are in an approximate state of energy equipartition. If this criterion can be achieved
multiple times, the system gradually loses memory of its initial conditions (i.e., the particles
no longer remember from where they originated in the initial phase space), leaving the system only knowing the total interaction energy, the total angular momentum and the particle
masses. Stone and Leigh [11] found that if the system enters two or more such states, which
the authors term "scrambles", then this removes most of the regular subset, leaving beyond
only the ergodic subset of the interactions. More recently, other authors have begun to explore more sophisticated methods for defining when an interaction formally becomes ergodic
e.g. [13, 38].
We also point out that longer integration times for the simulations tend to correlate withincreased error accumulation. As already discussed, this can mimic the effects of chaos. Hence,
one must be careful to manage the errors properly in order to identify the subset of interactions
4

SciPost Phys. Core 6, 016 (2023)

Perturbed triangle number 195908

0.0

0.0

0.2
0.4
0.6
0.8
1.0

0.5

0.0

X[AU]

0.5

0.4
0.2

0.2

0.6

0.6

0.8

0.8
2.0

1.5

1.0

0.5 0.0

X[AU]

0.5

1.0

1.0

0.5
1.0

1.5

0.5

0.0

0.5

X[AU]

1.0

Perturbed triangle number 5757532

1.0

Y[AU]

0.0

0.2
0.4

P1 Initial position
P2 Initial position
P3 Initial position
P1 Trajectory
P2 Trajectory
P3 Trajectory

0.5

0.0

0.4

Perturbed triangle number 252488

1.0

Y[AU]

1.0

Perturbed triangle number 6320049

0.6

Y[AU]

0.2

Y[AU]

Y[AU]

0.2

0.5
0.0
0.5

1.5
2.0

0.8

0.4

0.4

1.5

Perturbed triangle number 265238

0.6

0.6

1.0

0.5

0.0

X[AU]

1.0

0.5

0.8

0.6

0.4

0.2

X[AU]

0.0

0.2

0.4

Figure 1: Each panel shows the trajectories in position-space for all three particles fordifferent sets of initial conditions. Red, blue, and orange stars represent the initial
positions for each of the particles. The black arrows indicate the directions of ejection
for all three particles, where one of the particles is ejected as a single particle andthe other two form a binary system.

that should correspond to the true end state of the system (i.e., what nature would produce).
Our methods for managing error accumulation in the simulations are described in detail in
Section 3.

2.2

Lévy flights

Lévy flights correspond to very long excursions of one of the particles, where the particle has
nearly zero but still negative total energy and remains barely bound to the system. These
correlate with very long integration times, since the duration of these excursions can be very
long e.g [39, 40].
The nature of Lévy flights remains an open research question. Lévy flights are a class ofrandom walks with a heavy-tailed distribution. In our context, a Lévy flight is an extremely
long excursion following a fast chaotic interaction (see [39–41]), whereas interactions without
exceedingly long Lévy flights tend to follow an exponential drop-off in the cumulative interaction lifetimes, which can be described using a half-life formalism (see [13, 38–40] for a morein-depth discussion of the relevant physics).
The importance of Lévy flights for defining chaotic interactions is still not completely understood. For example, consider a three-body interaction that undergoes a very long excursion or Lévy flight with the energy of the single being very close to zero. The interaction ends
promptly at the end of the Lévy flight, with the single interacting once with the temporary
binary followed by a prompt ejection and the end of the interaction. If this experiment isrepeated with minute perturbations to the initial conditions, one can imagine that these simulations, grouped very closely together in phase space, will all end with a prompt ejection of
the returning single, only altering the properties of the ejected single and the left-over binary

5

SciPost Phys. Core 6, 016 (2023)

very slightly. This would, at least in principle, introduce a small regular subset into three-bodyinteractions that were already defined as being chaotic, having achieved our chosen criterion
for ergodicity. These subsets should occupy sufficiently small phase space volumes, however,that they should constitute ≪1% of our chaotic subset.

2.3

Scramble cut

A scramble is defined as a temporary state in which all three particles have comparable energies
and are not in a hierarchical configuration. We define NS as the number of scrambles, whichis the number of times the system has entered into such a temporary state of approximate
equipartition. If the system achieves approximate energy equipartition, it can lose all memoryof the initial conditions and become ergodic, retaining only memory of the particle masses, thetotal interaction energy and the total angular momentum. Stone and Leigh [11] showed thattwo or more such scrambles are a sufficient criterion to remove most of the regular regions ofphase space, leaving only the ergodic subset. In this paper, if a three-body interaction satisfies
the criterion NS ≥ 4, where NS is the number of scrambles, it is defined to be a chaotic or
ergodic interaction. After applying this cut, we find that 82.1% of our simulations correspondto the chaotic region and 17.9% to the regular region.

3

MethodIn this section, we describe the numerical simulations used to perform all numerical experiments conducted in this paper. We first present TSUNAMI, the gravity integrator we use tocompute the time evolution of our three body systems, before moving on to describing how
we define our initial conditions and set up our experiments.

3.1

N-body Integrator

We use TSUNAMI [37], a few-body code ideally suited to follow the time evolution of stronglyinteracting systems. TSUNAMI manages the computations using sophisticated regularization
techniques to achieve the required accuracy and precision, which are needed to properly model
very close approaches between particles.
TSUNAMI implements a leapfrog scheme derived from a time-transformed Hamiltonian,using a combination of the logH and TTL schemes [42,43]. Here we use only the logH formulation. Because the leapfrog scheme is accurate only to 2nd order, we increase the accuracyof the integration by employing a Bulirch-Stoer extrapolation scheme [44] with a tolerance of
10−13 . TSUNAMI also includes the post-Newtonian terms of order 1, 2 and 2.5 [45]. For moredetails on the implementation and performance of TSUNAMI, we refer the reader to [46].

3.2

Initial conditions

In the following sub-sections, we describe how the initial conditions for our suites of numericalsimulations are defined.
3.2.1

Unperturbed initial configuration

The initial conditions are generated by perturbing the positions of the particles from the unperturbed triangular configuration. The unperturbed triangular configuration is described here.
We label each particle with the numbers 1,2,3, which also indicate the side of the triangle
that each particle opposes. The initial, unperturbed positions for the particles 1 and 2 areP1 = (x 1 , y1 , z1 ) = (1.5, 1.5, 0) and P2 = (x 2 , y2 , z2 ) = (2.5, 1.5, 0). Note that the triangle
6

SciPost Phys. Core 6, 016 (2023)

Y [AU]

Triangle perturbation
3.0

Equilateral triangle

2.5

Generated triangle

2.0
1.5
1.0
1.0

1.5

2.0

X [AU]

2.5

3.0

Figure 2: Scheme of the initial setup. The equilateral triangle vertices are labeled
as P1 , P2 , and P3 . The generated particles corresponding to the perturbed triangle
are labeled as P1gen , P2gen , and P3gen . The green circles correspond to the boundaries
of each generated particle (rgen = 0.5 au). The insets in the upper right show theequilateral triangle and the perturbed version.

side is 1 au. This represents a perfect equilateral triangle, where each side will have a lengths△ = ((x 1 − x 2 )2 + ( y1 − y2 )2 )0.5 in Cartesian coordinates. The initial position of particle P3
was then calculated as:
x3 =
y3 =

x1 +
y2 +

s△
2


,

2
s△
−

(1)
 s 2 0.5
△

2

.

(2)

We set z = 0 , and focus on purely 2D coplanar configurations.
3.2.2

Perturbed initial configurations

We perturb the initial positions of the particles in the following way. The position of every
particle is chosen within a radius of size rgen from the vertex, as shown in Figure 2. First,
we randomly draw a perturbation radius r within a circle of rgen . Since r represents a radial
2
perturbation, we randomly sample it from a distribution uniform in rgen
. The position of eachparticle is perturbed by the same amount r, but with a radial direction that is independently
drawn for each particle. Specifically, the direction of the perturbation is set by the angles θ1 ,
θ2 and θ3 for particle 1, 2 and 3. The angles are sampled from a uniform distribution in therange (0, 2π).
We generate three different sets of simulations with different perturbation magnitudes, inorder to increase sampling in the vicinity of the singularity and to re-perform our calculations
with post-Newtonian corrections turned on. To explore the entire phase space, we first sample
each generated particle within a rgen = 0.5 au circle positioned at the vertex of the opposing side of the triangle. To explore the vicinity of the singularity we increase our sampling
resolution and use instead rgen = 10−7 au.

7SciPost Phys. Core 6, 016 (2023)

Particle 1
Particle 2
Particle 3
Random Position

Generated triangle 888
Generated triangle 777
Equilateral triangle

Figure 3: Left panel: three perturbed triangles configurations. The green circles
indicate the boundaries where each new particle can be generated. Top-right panel:
initial state of the system corresponding to an equilateral triangle 1 au side. Bottomright panel: the distribution of perturbation radius, defined as the distance from theperturbed position to the vertex of the triangle.

We assume identical point particles and purely Newtonian forces for Set A and Set B. In
Set A we assumes a vertex radius of 0.5 au, while for Set B, the sampling is performed usinga vertex radius of 10−7 au.
We re-run the same simulations in Set C but with post-Newtonian corrections turned on.
These are required very close to the singular regions in the three-body problem, in order toavoid nonphysical solutions (e.g., escape velocities for the single star that exceed the velocity
of light). This requires the introduction of a critical radius to define the conditions for a merger
to occur, which we will come back to in more detail below. The initial particle velocities arealways zero and the particle masses are all set to 1 M⊙ .
Table 1 shows the number of simulations performed for each set (Column 2), the particle
radii (Column 3), and our choices for the sampling radii rgen (Column 4).
As shown in Figure 2, for every new triangle that we generate, we re-compute the lengthsof each side, where S1 S2 and S3 represent the initial lengths of the sides opposing particles 1,
2, and 3, respectively. These choices are made to set a good balance between drop-in time andthe total number of simulations we can perform (i.e., shorter drop-in times result in simulations
that finish quicker, hence we can perform more simulations for a constant total real-time runtime), and to connect each one of the possible initial configurations to the properties of theejected particle and the left-over binary (see Figure 4 and Figure 9).

3.3

Ejection criteria

To determine when a particle is ejected from the system, we check the state of the systemat every time step (each simulation has a maximum integration time of 105 yrs). The mostbound pair is defined as the binary, and at every time-step, we check if the ejected single
remains bound to the binary centre of mass. The integration stops if the escaping single star
has positive total energy and the final hyperbolic binary-single orbital separation is 100 times
that of the most compact binary semi-major axis.
Roughly 0.12% of our simulations in Set A did not finish, i.e the integration time needed

8

SciPost Phys. Core 6, 016 (2023)

Table 1: For Set A and Set B, the initial particle positions are sampled following theprocedure described in Section (3.2.2). For Set C, the initial configurations are shownin Section (3.4). The variables rs , rgen corresponds to the Schwarzschild radius for a
1M⊙ particle, and to the sampling radii respectively.

Set A
Set B
Set C

Number of
Simulations107
5 × 106
5 × 105

ParticleradiusPoint particle
Point particle
500 rs

rgen
[au]
0.5
10−7
0.5

for some systems to fulfill our ejection criteria is longer than our total integration time. Thisis typically due to long-lived Lévy flights, where the return time of the escaper exceeds ourchosen maximum integration time. This represents a negligible fraction of our simulations,and our results remain unchanged without these simulations.

3.4

Simulations with post-Newtonian corrections

We use Set A to define our initial particle positions, but set the radius of each particle to beequal to 500 times the Schwarzschild radius. Particles are assumed to merge using the stickystar approximation: if the particle radii overlap in both time and space, we assume a merger
occurs. We choose 500 Schwarzschild radii for the following reasons. First, when two particles
are closer than this distance, they are practically decoupled from the rest of the system, and thebinary will merge very quickly. Second, the post-Newtonian approximation used to correct the
Newtonian acceleration begins to break down at around this distance [45]. Third, by avoiding
integrating the final part of the inspiral we save a great amount of computational run-time.
This is because, as the period of the merging binary decreases, the post-Newtonian corrections
would start to dominate and the integration time would increase. This scenario is avoided
by choosing 500rs for the particle radius and by adopting the sticky-star approximation todetermine when the particles merge (i.e., when their radii overlap in both time and space).

4

Theoretical Expectations

In this section, we describe the theoretical expressions (see Chapter 7 of [47]) for the expecteddistributions of the final binary orbital energy or, equivalently, the final binary semi-major
axis (i.e.,the inverse binding energy distribution), the final binary eccentricity, the single starejection velocity and total system lifetime, for comparison to our numerical scattering results.

4.1

Dividing the phase space into ergodic and regular subsets

In order to properly compare our simulated data to theoretical expectations, the simulations
must first be divided into chaotic and regular subsets. This is done using the criterion definedin Section 2.3. In Section 5, we will perform our comparisons with and without applying thesecuts, to better understand how the chaotic and regular subsets are each contributing to thetotal distribution (i.e., without applying any cuts for ergodicity).

9

SciPost Phys. Core 6, 016 (2023)

4.2

Orbital energy distribution

We use the inverse orbital energy distribution we obtain from our simulations to compare to
the theoretical binary orbital energy distributions provided in [47] and shown in Figure. 8. As
we are working with particles released from rest, the initial kinetic energy is equal to zero. We
can thus write the initial energy in terms of the potential energy alone, and the binary binding
energy in terms of the binary semi-major axis:
m1 m2
m1 m3
m2 m3
+
+
|r − r2 | |r1 − r3 | |r2 − r3 | 1
ma m b.
= −G
2a binE0 = −G
EB





,

(3)
(4)

The variables ma and m b correspond to the particle masses in the final left-over binary
post-interaction.
The distribution used in [47] for the inverse orbital energy is given by Equation 5, wherez = |E0 |/|EB | corresponds to the inverse binary orbital energy. For the planar case, this isf (z) = 3.5z 2.5 .

4.3

(5)

Eccentricity distribution

We compare our simulated results to several different theoretical eccentricity distributions inFigure. 8. These are described below.
The thermal distribution for the binary eccentricity, which assumes a detailed steady-statebalance between binary creation and destruction [48], is given by Equation 6 :f (e) de = 2e de .

(6)

In the planar case, the eccentricity should be distributed as Equation 7
f (e) de = e(1 − e2 )−1/2 de .

(7)

The power-law index on (1 − e2 ) is taken to be −1/2, which is appropriate for the low angularmomentum limit [8, 49].

4.4

Escape velocity distribution

We adopt Equation 7.19 in [47] for the escape velocity distribution shown in Figure. 8. Thisis appropriate to the planar case and is given byf (vs ) dvs

5=(3.5|E0 |7/2 (ms M /mB ))vs dvs(|E0 | + 12 (ms M /mB )vs2 )9/2.

(8)

Results

In this section, we present the results of our numerical scattering experiments. We begin bypresenting the phase space plots for our zero angular momentum equilateral triangle set-up,
before moving on to the final calculated distributions for the post-interaction binary orbitalenergy, orbital eccentricity, escaper velocity and total system lifetime. We then compare theseto the expected theoretical predictions, using our phase space plots to identify the origins of
the disagreements between our theoretical calculations and the simulated data.
10

SciPost Phys. Core 6, 016 (2023)

Region 4
Region 6

Particle 1
Particle 2
Particle 3

Region 5
Region 2

Region 1
Region 3

Figure 4: Phase space plot, color-coded by ejected particle. X- and y-axes indicate theratio between the sides S1 /S2 and S1 /S3 of the perturbed equilateral triangles. The
singularity is at S1 /S2 = S1 /S3 = 1, corresponding to a perfect equilateral triangle.
Each subsequent inset corresponds to a zoom-in on the regions in the phase space.

5.1

Phase space divided by ejected particle

We represent the initial phase space in 2 dimensions, so that each point in this space corresponds to a unique initial configuration (i.e. a single simulation). For our purposes, theparameters chosen for the X and Y axes are the ratios between the lengths of the sides ofeach of the generated triangles, as shown in Figure. 3. Here, S1 , S2 and S3 correspond to thelengths of the sides opposing, respectively, particles P1 , P2 , and P3 . We then plot the ratios
S1 /S2 and S1 /S3 on each axis to construct our phase space plots for the outcomes of the perturbed equilateral triangle experiment. In this coordinate system, the singularity is locatedat the coordinates (1, 1), and corresponds to the initial conditions for a perfect equilateraltriangle.
We colour-code each data point according to which particle is ejected, using the colour
corresponding to that particle (see Figure. 4). In this case, uniform swaths of colour correspond to prompt regular outcomes (i.e., these initial configurations evolve passing through 4
or less scrambles) and multi-coloured patches correspond to chaotic regions of phase space(i.e., systems that experience more than 4 scrambles in their lifetime).
5.1.1

Full phase space

The central panel in Figure 4 shows the phase space limited to 0 < S1 /S2 < 3 and 0 < S1/S3 < 3
for the interactions of Set A.

11

SciPost Phys. Core 6, 016 (2023)

Near the singularity, we see regular swaths of uniform colour, separated by chaotic regionsof phase space. In the regular regions, small perturbations to the initial state (i.e., equilateraltriangles where the particles have similar initial positions) do not lead to a different particle
being ejected. In the multi-coloured regions of our phase space corresponding to chaotic or
ergodic regions, the identity of the ejected particle is very sensitive to the initial conditions,causing different particles to be the ones ejected in spite of very minor changes to the initialconditions.
Moving away from the central singularity, we see emerging more and more chaotic regions of phase space. Moving far from the singularity to initial conditions corresponding to
extremely deformed triangles, the phase space has three arms. In each arm, the initial configuration is approaching the binary-single scattering regime, since the length of one of thesides of the triangle is very small. The arm that goes to the top represents the deformationof the triangle where particle P3 (coloured by orange) is initial far from particles P1 and P2 .
Similarly, the arms that goes to the right and the left bottom represent the deformation of thetriangle where the particles P2 (coloured by blue) and P1 (coloured by red) are far from theothers.

5.1.2

Zoom-in on the singularity

In Set B, the perturbations to the equilateral triangle are very small. Therefore, we can better
appreciate the behaviour of the system very close to the three-body collision.
Each top panel in Figure. 5 is a subsequent zoom-in on the central singularity, with thezoom-ins becoming more extreme going from left to right. We see regular coloured regionswith clear boundaries separated from the chaotic ones. As we zoom-in more and more on
the central singularity, we find that the chaotic regions disappear to eventually become fully
regular in the vicinity of the singularity (i.e., we see only regular patches corresponding touniform colours).
Panel a is a zoom-in on the singularity for Set A. We see regular uniformly coloured regions with clear boundaries separated from the chaotic ones. In panel b (a subsequent zoom-inon the singularity corresponding to panel a) we see a plume of uniform colour emerging in
between the red and yellow swaths of regular regions of uniform colour. This blue plume becomes prominent as we zoom-in more, until eventually we are left with only uniform colouringin the vicinity of the singularity, as shown in panel d. Here we find that the chaotic regions
disappear and the phase space becomes fully regular in the vicinity of the singularity.
5.1.3

Zoom-in on the binary-single regime

The bottom left panel of Figure 5 shows the same phase space as Figure 4 of our Set A , but
for smaller values of S1 /S2 and S1 /S3 . The right middle pair of panels shows the binary-single
scattering regime. The bottom right figures show two zoom-in on the regular regions of thebinary-single scattering regime.
The three corners of the phase space in Figure 5 correspond to the binary-single scattering
regime, where two of the bodies are very far from the third. In panel e, we see regular regions
repeating at approximately regular intervals, which are separated by chaotic patches of phasespace. The first three regular regions have roughly the same shape (see panel f). However,the second zoom-in panel reveals two bulges arising from the centre of each regular region(see panel g). The bulges increase in size as we move farther into the extreme binary-single
scattering regime (see panel h). Importantly, the regular structures do not disappear andappear to continue indefinitely down to increasingly small scales.

12

SciPost Phys. Core 6, 016 (2023)

Particle 1
Particle 2
Particle 3

Figure 5: Same as Figure 4, but zoomed on specific regions of the phase space. Top
panels, left to right: each subsequent inset on the top is a zoom-in on the singularity(see panels a,b,c, and d). Bottom panels e and g zoom on the regime of binary-single
scattering, where one side is very small compared to the others. Bottom right panels
f and h: the zoomed-in version of the two different shapes we see emerging in panelse and g.

5.1.4

Ergodic and regular phase space

We can better appreciate the clear distinction between regular regions and the chaotic regions
by separating them in the initial phase space plots of Figure. 4.
Figure 6 shows the region of our phase space close to the singularity for which
S1 /S2 = S1 /S3 (left column) and the binary-single regime S1 ≫ S2 , S3 (middle and right panels). The last two panels show a zoom-in of the middle panels. The top row shows the ergodic
subset, and the bottom row shows the regular subset. As expected, our scramble criterion
neatly separates the regular and chaotic regions, including the triangle-shaped regions aroundthe singularity and the self-repeating structures in the binary-single regime.

5.2

Including post-Newtonian corrections and general relativistic effects

Very near to the singularity, general relativity should begin to become important as the distance between the particles approach the Planck scale. We re-perform our suite of scatteringexperiments, adopting the same sets of initial conditions, but now including post-Newtonian
corrections to account for the effects of general relativity. The results of this exercise are shownin Figure 7, which shows how Figure 5 changes when general relativity is accounted for.
We find that turning on general relativity (i.e turning on post-Newtonian corrections)causes us to lose information in the binary-single regime. This is made evident by the darkgrey data points, which correspond to prompt mergers of two of the particles. This loss of
information (i.e., chaotic and regular zones of the parameter space are replaced by the dark13

SciPost Phys. Core 6, 016 (2023)

Phase space plot
Near singularity

Phase space plot
Binary-single regime

Phase space plot
Binary-single regime zoom in

P2

Ejected particle

P3

P1

Figure 6: Top panels: phase space for the ergodic subset near the singularity (left)
and in the binary-single scattering regime (middle and right). The bottom panel
shows the same regions of the top panels, but for the regular subset instead.

fParticle 1
Particle 2
Particle 3
Collision

h

g
e

a

g

e

Figure 7: Same as Figure 5, but employing post-Newtonian corrections. The darkgrey dots correspond to simulations ending in a two-body collision, or merger.

14

SciPost Phys. Core 6, 016 (2023)

grey area) occurs because the two black holes merge quickly when the distance separating
their centres of mass approaches the sum of their Schwarzschild radii due to the loss of energy
from gravitational waves. In the binary-single regime of the phase space, two particles always
begin very close together at the beginning of the simulations. This is not a problem for the
Newtonian regime, but leads to quick mergers in the post-Newtonian regime. We caution that
newtonian gravity is scale invariant whereas post-Newtonian gravity is not. Hence we do notextend our results to systems with different particle masses.

5.3

Comparison with theoretical expectations

In this section, we present our final simulated distributions of binary orbital energies, orbital
eccentricities, escaper velocities and total system lifetime, and compare the results to the theoretical expectations provided in Section 4. Because the statistical theory is based on purely
Newtonian physics, we compare it only with our Set A, which does not include post-Newtonian
corrections.
5.3.1

Distributions of outcome properties

In the top panels of Figure 8, the histograms show the final distributions of binary eccentricities, ejected particle velocities, and system lifetimes. The distributions are separated into
regular (all simulations with NS ≤ 4) and ergodic (all simulations with NS > 4) subsets, alongwith the entire suite of simulations in Set A. The criterion for separating the two subsets is
NS ≥ 4, as discussed in Section 2.3.
In the binary eccentricity distribution panel (top left panel), both subsets (i.e, ergodic andregular) show a strongly suprathermal distribution. Concerning the escape velocity distribution, the ergodic subset is similar to a log-normal distribution while the regular subset showsan irregular distribution. Both of them peak near 104.7 m/s. Further differences between theergodic and regular subsets emerge from the system lifetime distribution (top right panel).
The regular interactions tend to be shorter than the ergodic ones. For the ergodic subset, thedistribution looks like a log-normal with the peak near 4 yr (log10 (lifetime/yr) ∼ 0.6), while
the regular one peaks near 10 yr (log10 (lifetime/yr) ∼ 1).
In the middle panels of Figure 8, the histograms show the distributions of initial systemenergies, binary orbital energies, and the inverse binding energies. Relative to the theoretical
expectation, we find too many compact binaries (i.e., binary stars with energies E0 /EB ≲ 0.6)
and obtain not enough wide binaries (i.e., binary stars with energies E0 /EB ≳ 0.6). Because
the system has no angular momentum, the total energy depends only on the initial triangle
shape (see Equation 4). The vertical dashed line represents the median value for the initialenergy. This value matches the peak of the regular subset. Each numerical integration ends
in a binary who takes most of the energy of the system, which we calculate using Equation 4.
The dashed line corresponds to the median value for the binary energy and we see this does
not match the peak of any of the subset curves. The peak we see in the regular subset is aconsequence of the mismatching values between the median of the initial energy distribution
and the binary energy distribution.
In the bottom panels of Figure 8, the histograms show the distributions of binary orbital
eccentricities, ejected particle velocities and inverse binding energies, respectively. We compare them with our theoretical expectations (see the equations in Section 4) [9]. The panels
only show the ergodic subset, because we do not expect the regular subset to follow the outcome distributions expected from the ergodic theory. We see an excess of eccentric binaries
relative to a thermal distribution (Equation 6), but the agreement is overall good with the zeroangular momentum distribution (Equation 7). Relative to the theoretical expectation for the15

SciPost Phys. Core 6, 016 (2023)

Ergodic subset Ns > 4
Regular subset Ns 4
All setMed(E0)Med(EB)

Ergodic subset
eq. 6
eq. 7

Ergodic subset
eq. 8

Med(E0)Med(EB)

Ergodic subset
eq. 5

Figure 8: Distributions of outcome parameters. The green lines show the distributions for our entire suite of experiments, blue lines correspond to the regular subset
where NS ≤ 4 , the red line corresponds to the ergodic subset where NS > 4. Bottom
panels show the ergodic subset along with the theoretical curves shown in 4.

ejection velocities (Equation 8), we obtain an excess of high-velocity escapers and a lack of
low-velocity escapers.
5.3.2

Phase space distribution of outcome properties

Figure 9 shows the initial phase space plots colour-coded by the properties of the final particles.
The use of the initial phase space plots is particularly useful to identify the regular regions,which are likely the cause of the disagreement between the simulations and the statisticaltheory.
In the zero angular momentum limit, our naive expectation is to preferentially produce
compact binaries with high eccentricities (i.e., binaries with low orbital angular momentum).
But this is not always observed in our experiments. This occurs because the escaper tends tocarry away some significant angular momentum, requiring the binary to retain the remaining
reservoir of angular momentum (but going in opposite directions in order to cancel out tozero). This is clearly observed in Figure 9 upon comparing the ejection velocities to the orbital
properties of the binary. The final binary semimajor axes in the top right panel of Figure 9 are
expressed in units ofGm2
anorm = −,(9)
2E0
where anorm is the semimajor axis corresponding to the initial energy E0 of the system.
Focusing on the exterior boundaries of the regular regions surrounding the central singularity, we observe higher ejection velocities corresponding to higher eccentricities and wider
binaries. Conversely, focusing on the interior boundaries of the same regular regions, we see
higher ejection velocities correlating with low eccentricities and more compact binaries. Focusing on the singularity, we observe compact binaries correlated with high eccentricities and highejection velocities. In the limit of the escaper velocity becoming very low, we tend to observe
16

SciPost Phys. Core 6, 016 (2023)

1.0
0.8

0.6

0.6

0.4

0.4

0.2

0.2

ab

0.8

0.0

3

5.0

2
4.5

1
0

4.0
1

Figure 9: Color-coded phase space plot showing the final binary semi-major axes inunits of anorm (top left panel), orbital eccentricities (top right panel), escaper ejection
velocities (lower right panel) and total system lifetimes (lower left panel).

either high eccentricities and wide binaries or low eccentricities and compact binaries, since
the escaper is no longer able to carry away much angular momentum, or LS ∼ 0. Conversely,as the escaper velocity reaches its maximum near the singularity, we observe more compact
binaries and higher eccentricities, since the binary is left with little to no remaining angularmomentum. It is for these reasons that we observe a supra-thermal eccentricity distribution(e.g., [11]) and a binary orbital energy distribution that prefers compact binaries.
The shortest-lived systems (i.e., those that satisfy the ejection criteria during the early integration times of the system) are closest to the singularity, with the external borders of themain regular regions reaching the longest lifetimes, and correlating with low ejection velocities, high eccentricities, and compact binaries. Very near to the singularity, the very shortinteraction lifetimes likely indicate regular interactions, which is consistent with the progressive zoom-in on the singularity shown in Figure 5. The low lifetime regions are surrounded
by a narrow strip of very long-lived Lévy flight interactions (i.e., one of the particles goes on
a long excursion with near zero but still negative total energy before the ejection criteria issatisfied), consistent with the results of [38].
5.3.3

Weak chaos and the disagreements between theory and simulations

In this section, we summarize our findings regarding the disagreements between our theoretical expectations and simulated results, as presented in the previous sections. This is donein an effort to identify their origins using the post-interaction distributions of binary orbitalenergies, orbital eccentricities and escaper velocities. In this section, we divide our total suiteof simulations into ergodic and regular subsets. To do this, we select our data following thethe scramble number criterion discussed previously in Section 2.3 (NS > 4), and isolate mostof the regular regions from the chaotic ones. Approximately 17.9% of the simulations belong
to the regular regions and 82.1% belong to the ergodic regions, revealing that most of theparameter space is chaotic.
17

SciPost Phys. Core 6, 016 (2023)

One of the main inconsistencies of the data with our theoretical expectations is the bumpat z = 0.71 in the inverse binding energy distribution (see the right-hand middle panel of
Figure. 8). In order to identify the origin of this bump, we plot separately the distributions forthe ergodic and regular subsets for both E0 and EB . We find a misalignment between the peaksof the regular E0 and EB distributions relative to the medians (see middle insets of Figure. 8).
The origin of this misalignment can be traced back to the slightly asymmetric distribution of
initial binary orbital energies from which we sampled our initial conditions, which slightly
over-populate the hard end of the distribution relative to the soft end.
As a consequence of this misalignment, we expect a bump in the regular subset of theinverse binding energy distribution, as seen in the right-most middle panel of Figure. 8. In
fact, the ratio between the medians of the E0 and EB distributions matches the peak in the
inverse binding energy distribution for the entire suite of simulations, therefore confirmingthat the origin of the misalignment is in the regular, and not in the ergodic, subset.

6

Discussion & Summary

In this paper, we consider an idealized initial configuration of the general three-body problem,in order to study chaos in the vicinity of a singularity in the general N-body problem. The
experiment we perform is to construct an equilateral triangle, with the initial positions of
the three identical particles located at each vertex. For a perfect equilateral triangle, all threeparticles will eventually collide at the centre of the triangle in a three-body collision, if released
from rest at the same time and all particles have equal masses. This would give rise to a
singularity in the Newtonian acceleration because the particles would occupy the same positionin space and time at the system centre of mass. By perturbing the triangle and re-performing
the experiment, the particles arrive at the origin slightly offset in space and time. The larger
the perturbations, the farther is the parameter space from this singularity.
We perform in total 107 simulations using the TSUNAMI code [37]. The large suite ofsimulations proves necessary in order to generate phase space plots with sufficient resolution
to be interpretable (i.e., the chaotic and regular zones of the parameter space are clearly
identified). With our final suite of numerical simulations, we generate the final distributionsof binary orbital energy, orbital eccentricity, escaper velocity, and total interaction lifetime.
We then compare these distributions to the standard theoretical expectations provided in theliterature. We find significant differences in the final orbital properties of our binaries, relative
to the theoretical expectations, including more eccentric orbits and a higher fraction of morecompact binaries. Using our phase space plots, we explain the origins of these differences,which primarily come down to the total angular momentum carried away by the escaping
single star. For example, in the limit of v∞ → 0, the left-over binary tends to be more circular,
since it contains the entire angular momentum reservoir of the interaction.
Finally, using our results for this idealized experiment in the zero-angular momentum limit,
we make predictions for the expected properties of binaries formed via three-body interactionsof initially isolated single stars in very dense environments. We expect these predictions to
be most applicable to star clusters with isotropic stellar velocity distributions and little to norotation.

6.1

Astrophysical implications

In this section, we discuss the properties of binaries formed via three-body interactions in dense
environments, and compare our results to theoretical expectations for the ergodic subset ofour simulations.

18

SciPost Phys. Core 6, 016 (2023)

6.1.1

Three body binary formation in isotropic clusters

In isotropic star clusters, binaries formed from three-body interactions of initially all single
stars should have a low total angular momentum. This is because there is little to no angularmomentum in the cluster, and the stellar orbits are such that three single stars that meetnear the cluster centre of mass should be on nearly radial orbits, causing the three stars toapproach their common centre of mass with small impact parameters. In the opposite limit ofhigh angular momentum and three-body binary formation in clusters with significant rotation,
we would expect the opposite. That is, the three stars should all approach their common centreof mass with large impact parameters, causing the total angular momentum to be large, since
the orbits tend to be less radial and more tangential due to rotation. In this paper, we focus
on the zero angular momentum limit, and defer the rotating case to a future paper.
We expect binary formation via 1+1+1 scatterings to be the dominant mode of binary
formation in very dense environments [50]. The rate is given by [51]: dn 
bdt= 1.91 × 10−133  m 5  10 kms−1 9npc−3 yr−1 ,104 pc−3
M⊙
vrms

(10)

where n is the number density, m is the average mass and vrms is the 3D root-mean-squarevelocity. Assuming m = 0.5 M⊙ , a mass-to-light ratio of 2 and multiplying by the core volume
of a typical globular cluster (i.e., assuming a core radius of 1 pc), we find using values takenfrom the Harris globular cluster catalog [52] for NGC 362, NGC 6397, NGC 6256 and NGC
6522 that the time between binary formation events due to 1+1+1 scatterings of all single
stars are, respectively, ∼ 609 Gyr, 22 Myrs, 2.2 Gyrs and 5.5 Gyrs. Note that these are all postcore collapse clusters. Hence, although often rare, binary formation due to 1+1+1 scatteringsof all single stars do occur in the densest Milky Way globular clusters well within the lifetimeof the cluster. Post-core collapse clusters are the most likely to harbour binaries formed via
1+1+1 interactions.
6.1.2

Binary properties

For the ergodic subset, stars are ejected at higher velocities relative to the theoretical expectation. Hence, the single star can take away a significant fraction of the total angular momentum
reservoir. In this limit, we produce more wide binaries relative to the theoretical expectation,
along with very high eccentricities. In the limit that the single star takes away no angular momentum (i.e., its final velocity at infinity is zero), we find more compact binaries and smaller
eccentricities relative to the theoretical expectations.
Overall, our results suggest that binaries formed from three-body interactions in isotropic
star clusters should end up with orbital parameters that correspond to rapid coalescence forbinaries composed of two black holes due to the emission of gravitational waves (GWs). We
tend to find either compact circular binaries or wide very eccentric binaries. In both cases, themerger times will be short, since the timescale for a merger to due gravitational wave emission
is proportional to the binary orbital separation ∝ a4 (i.e., more compact binaries merge faster)and also proportional to (1 − e2 )7/2 for the orbital eccentricity [53]. Nonetheless, dependingon the scale of the system, even the high eccentricity might not be enough to merge a binary.
This is because post-Newtonian corrections are not scale-free, and we can always scale up the
distance so that the pericenter passage is not close enough to ensure a rapid merger, evenat very high eccentricities. This will be investigated in the future using dedicated numericalexperiments with post-Newtonian corrections.

19

SciPost Phys. Core 6, 016 (2023)

Acknowledgements
NWCL gratefully acknowledges the generous support of a Fondecyt Iniciación grant 11180005,
as well as support from Millenium Nucleus NCN19-058 (TITANs) and funding via the
BASAL Centro de Excelencia en Astrofisica y Tecnologias Afines (CATA) grant PFB-06/2007.
NWCL also thanks support from ANID BASAL project ACE210002 and ANID BASAL projects
ACE210002 and FB210003. A.A.T. acknowledges support from JSPS KAKENHI Grant Numbers
17H06360, 19K03907 and 21K13914.

Data Availability
The TSUNAMI code, the initial conditions and the simulation data underlying this article will
be shared on reasonable request to the corresponding author.

References[1] H. Poincaré, Les méthodes nouvelles de la mécanique céleste (1856).
[2] V. Szebehely and C. F. Peters, Complete solution of a general problem of three bodies, Astron.
J. 72, 876 (1967), doi:10.1086/110355.
[3] S. J. Aarseth and F. Hoyle, Dynamical evolution of clusters of galaxies, I, Mon. Not. R.
Astron. Soc. 126, 223 (1963), doi:10.1093/mnras/126.3.223.
[4] S. J. Aarseth, Dynamical evolution of clusters of galaxies, II, Mon. Not. R. Astron. Soc. 132,
35 (1966), doi:10.1093/mnras/132.1.35.
[5] S. J. Aarseth, Dynamical Evolution of clusters of galaxies–III, Mon. Not. R. Astron. Soc.
144, 537 (1969), doi:10.1093/mnras/144.4.537.
[6] P. Hut and J. N. Bahcall, Binary-single star scattering. I - Numerical experiments for equalmasses, Astrophys. J. 268, 319 (1983), doi:10.1086/160956.
[7] J. J. Monaghan, A statistical theory of the disruption of three-body systems – I: Low angularmomentum, Mon. Not. R. Astron. Soc. 176, 63 (1976), doi:10.1093/mnras/176.1.63.
[8] J. J. Monaghan, A statistical theory of the disruption of three-body systems – II: High angularmomentum, Mon. Not. R. Astron. Soc. 177, 583 (1976), doi:10.1093/mnras/177.3.583.
[9] M. J. Valtonen, The general three-body problem in astrophysics, Vistas Astron. 32, 23(1988), doi:10.1016/0083-6656(88)90395-9.
[10] N. W. C. Leigh and S. Wegsman, Illustrating chaos: A schematic discretization of the generalthree-body problem in Newtonian gravity, Mon. Not. R. Astron. Soc. 476, 336 (2018),
doi:10.1093/mnras/sty192.
[11] N. C. Stone and N. W. C. Leigh, A statistical solution to the chaotic, non-hierarchical threebody problem, Nature 576, 406 (2019), doi:10.1038/s41586-019-1833-8.
[12] B. Kol, Flux-based statistical prediction of three-body outcomes, Celest. Mech. Dyn. Astr.
133, 17 (2021), doi:10.1007/s10569-021-10015-x.

20

SciPost Phys. Core 6, 016 (2023)[13] V. Manwadkar, B. Kol, A. A. Trani and N. W. C. Leigh, Testing the flux-based statistical prediction of the three-body problem, Mon. Not. R. Astron. Soc. 506, 692 (2021),
doi:10.1093/mnras/stab1689.
[14] P. Hut and J. N. Bahcall, Binary-single star scattering. I - Numerical experiments for equalmasses, Astrophys. J. 268, 319 (1983), doi:10.1086/160956.
[15] P. Hut and B. Paczynski, Effects of encounters with field stars on the evolution of low-mass
semidetached binaries, Astrophys. J. 284, 675 (1984), doi:10.1086/162450.
[16] D. C. Heggie, Chaos in the N -body problem of stellar dynamics, in Predictability, stability, and chaos in N -body dynamical systems, Springer, Boston, Massachusetts, US, ISBN
9781468459999 (1991), doi:10.1007/978-1-4684-5997-5_4.
[17] Z. P. Anosova and V. V. Orlov, Dynamical evolution of triple systems, Astrophys. Space Sci.
142, 233 (1988) doi:10.1007/BF00656215.
[18] M. Valtonen and S. Mikkola, The Few-Body problem in astrophysics, Annu. Rev. Astron.
Astrophys. 29, 9 (1991), doi:10.1146/annurev.aa.29.090191.000301.
[19] T. C. N. Boekholt, S. F. Portegies Zwart and M. Valtonen, Gargantuan chaotic gravitational
three-body systems and their irreversibility to the Planck length, Mon. Not. R. Astron. Soc.
493, 3932 (2020), doi:10.1093/mnras/staa452.
[20] J. Laskar, A numerical experiment on the chaotic behaviour of the Solar System, Nature
338, 237 (1989), doi:10.1038/338237a0.
[21] G. Boeing, Visual analysis of nonlinear dynamical systems: Chaos, fractals, self-similarity
and the limits of prediction, Systems 4, 37 (2016), doi:10.3390/systems4040037.
[22] N. W. C. Leigh, et. al., The chaotic four-body problem in Newtonian gravityI. Identical point-particles, Mon. Not. R. Astron. Soc. 463, 3311 (2016),
doi:10.1093/mnras/stw2178.
[23] G. J. Sussman and J. Wisdom, Chaotic evolution of the solar system, Science 257, 56(1992), doi:10.1126/science.257.5066.56.
[24] T. Ito and K. Tanikawa, Long-term integrations and stability of planetary orbits in
our Solar system, Mon. Not. R. Astron. Soc. 336, 483 (2002), doi:10.1046/j.13658711.2002.05765.x.
[25] W. B. Hayes, Is the outer Solar System chaotic?, Nat. Phys. 3, 689 (2007),
doi:10.1038/nphys728.
[26] J. Wisdom, S. J. Peale and F. Mignard, The chaotic rotation of Hyperion, Icarus 58, 137(1984), doi:10.1016/0019-1035(84)90032-0.
[27] T. C. N. Boekholt, F. I. Pelupessy, D. C. Heggie and S. F. Portegies Zwart, The originof chaos in the orbit of comet 1P/Halley, Mon. Not. R. Astron. Soc. 461, 3576 (2016),
doi:10.1093/mnras/stw1504.
[28] A. C. M. Correia, Chaotic dynamics in the (47171) Lempo triple system, Icarus 305, 250(2018), doi:10.1016/j.icarus.2018.01.008.
[29] S. Portegies Zwart and T. Boekholt, On the minimal accuracy required for simulating selfgravitating systems by means of direct N -body methods, Astrophys. J. Lett. 785, L3 (2014),
doi:10.1088/2041-8205/785/1/L3.
21

SciPost Phys. Core 6, 016 (2023)[30] C. A. Martínez-Barbosa et. al., The evolution of the Sun’s birth cluster and the search
for the solar siblings with Gaia, Mon. Not. R. Astron. Soc. 457, 1062 (2016),
doi:10.1093/mnras/stw006.
[31] N. Leigh and A. M. Geller, Small-N collisional dynamics: Pushing into the realm of
not-so-small N , Mon. Not. R. Astron. Soc. 425, 2369 (2012), doi:10.1111/j.13652966.2012.21689.x.
[32] N. W. C. Leigh, A. M. Geller, M. M. Shara, J. Garland, H. Clees-Baron and A. Ahmed,
Small-N collisional dynamics - III: The battle for the realm of not-so-small-N, Mon. Not. R.
Astron. Soc. 471, 1830 (2017), doi:10.1093/mnras/stx1704.
[33] N. W. C. Leigh, et. al., Small-N collisional dynamics - IV. Order in the realm of not-so-smallN , Mon. Not. R. Astron. Soc. 480, 3062 (2018), doi:10.1093/mnras/sty2046.
[34] R. H. Miller and E. N. Parker, Formation of massive quasi-stellar objects., Astrophys. J.
140, 50 (1964), doi:10.1086/147891.
[35] J. Goodman and P. Hut, Binary-single-star scattering. V - Steady state binary distribution in a homogeneous static background of single stars, Astrophys. J. 403, 271 (1993),
doi:10.1086/172200.
[36] M. Valluri and D. Merritt, Orbital instability and relaxation in stellar systems, in
The chaotic Universe, World Scientific, Singapore, ISBN 9789812793621 (2000),
doi:10.1142/9789812793621_0014.
[37] A. A. Trani, M. S. Fujii and M. Spera, The Keplerian three-body encounter. I. Insights on theorigin of the S-stars and the G-objects in the galactic center, Astrophys. J. 875, 42 (2019),
doi:10.3847/1538-4357/ab0e70.
[38] V. Manwadkar, A. A. Trani and N. W. C. Leigh, Chaos and Lévy flights in the three-bodyproblem, Mon. Not. R. Astron. Soc. 497, 3694 (2020), doi:10.1093/mnras/staa1722.
[39] N. W. C. Leigh, M. M. Shara and A. M. Geller, When does a star cluster become a multiple
star system? - I. Lifetimes of equal-mass small-N systems, Mon. Not. R. Astron. Soc. 459,
1242 (2016), doi:10.1093/mnras/stw735.
[40] T. Ibragimov, N. W. C. Leigh, T. Ryu, T. Panurach and R. Perna, When do star clusters
become multiple star systems? II. Towards a half-life formalism with four bodies, Mon. Not.
R. Astron. Soc. 477, 4213 (2018), doi:10.1093/mnras/sty712.
[41] I. I. Shevchenko, Hamiltonian intermittency and Lévy flights in the three-body problem,
Phys. Rev. E 81, 066216 (2010), doi:10.1103/PhysRevE.81.066216.
[42] S. Mikkola and K. Tanikawa, Algorithmic regularization of the few-body problem, Mon.
Not. R. Astron. Soc. 310, 745 (1999), doi:10.1046/j.1365-8711.1999.02982.x.
[43] S. Mikkola and K. Tanikawa, Explicit symplectic algorithms for time-transformed Hamiltonians, Celest. Mech. Dyn. Astron. 74, 287 (1999), doi:10.1023/A:1008368322547.
[44] J. Stoer and R. Bulirsch, Introduction to numerical analysis, Springer, New York, US, ISBN
9781475722727 (1983), doi:10.1007/978-1-4757-2272-7.
[45] L. Blanchet, Gravitational radiation from post-Newtonian sources and inspiralling compact
binaries, Living Rev. Relativ. 17, 2 (2014), doi:10.12942/lrr-2014-2.

22

SciPost Phys. Core 6, 016 (2023)[46] A. A. Trani and M. Spera, Modeling gravitational few-body problems with TSUNAMI and
OKINAMI, Proc. Int. Astron. Union 16, 404 (2023), doi:10.1017/S1743921322001818.
[47] W. C. Saslaw, The three-body problem by Mauri Valtonen and Hannu Karttunen, Celest.
Mech. Dyn. Astr. 96, 345 (2006), doi:10.1007/s10569-006-9049-2.
[48] D. C. Heggie, Binary evolution in stellar dynamics, Mon. Not. R. Astron. Soc. 173, 729(1975), doi:10.1093/mnras/173.3.729.
[49] E. M. Standish, Jr., The dynamical evolution of triple star systems, Astron. Astrophys. 21,
185 (1972).
[50] B. Reinoso, N. W. C. Leigh, C. M. Barrera-Retamal, D. Schleicher, R. S. Klessen and A. M.
Stutz, The mean free path approximation and stellar collisions in star clusters: Numerical
exploration of the analytic rates and the role of perturbations on binary star mergers, Mon.
Not. R. Astron. Soc. 509, 3724 (2021), doi:10.1093/mnras/stab3254.
[51] L. Spitzer, Dynamical evolution of globular clusters, Princeton University Press, Princeton,
New Jersey, US, ISBN 9780691635361 (1987).
[52] W. E. Harris, A new catalog of globular clusters in the Milky Way, (arXiv preprint)
doi:10.48550/arXiv.1012.3224.
[53] P. C. Peters, Gravitational radiation and the motion of two point masses, Phys. Rev. 136,
B1224 (1964), doi:10.1103/PhysRev.136.B1224.

23