Chaos in the vicinity of a singularity in the Three-Body Problem: The equilateral triangle experiment in the zero angular momentum limit

We present numerical simulations of the gravitational three-body problem, in which three particles lie at rest close to the vertex of an equilateral triangle. In the unperturbed problem, the three particles fall towards the center of mass of the system to form a 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 when one side of the triangle is very short compared to the other two. We make phase space plots to study the regular and ergodic subsets of our simulations and compare them with the outcomes expected from the statistical escape theory of the three-body problem. We further provide fits to the 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 regular interactions, 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.


Introduction
The gravitational two-body problem was solved analytically by Newton in his PhilosophiaeNaturalis Principia Mathematica.His work made it possible to predict exactly the positions and velocities 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.
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][3][4][5] re-vitalised popularity in the three-body problem, allowing researchers to confront their analytic models directly with computer simulations.This allowed for the development of analytic approximations for regions of phase 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][8][9][10][11][12][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 the interaction ends.If the interaction ends promptly, with the incoming single star interacting with the initial binary components only once before the interaction ends, the interaction is said to be regular and analytic solutions can typically be found e.g.[14][15][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 to become 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 a positive 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 be ejected).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].This implies numerical and physical consequences [20] and, as already mentioned, points to the need 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 in position space for similar sets of initial conditions.Consequently, a high degree of accuracy and precision 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 the accumulation of numerical errors, and guaranteeing the reliability of orbital integrations over long periods of times (i.e., many orbital periods) [20,[23][24][25].Most dynamical systems display some aspect of chaos, including solar system small bodies e.g.[26][27][28], small stellar systems e.g.[6,10,11,[29][30][31][32][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-body problem.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 a perfect 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 this idealized initial configuration.By perturbing the planar equilateral triangle and releasing the particles from rest, they will arrive at the system centre of mass at nearly the same time, but slightly offset.In this way, we can probe the interaction outcomes directly in the vicinity of the singularity, 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, each located 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 issue of 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, with a 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 total interaction 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.

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 to correspond 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 to quantify the outcome properties as a function of the initial conditions.

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 with increased 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  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.

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 of random 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][40][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][39][40] for a more in-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 is repeated 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 very slightly.This would, at least in principle, introduce a small regular subset into three-body interactions 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.

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 N S as the number of scrambles, which is 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 memory of the initial conditions and become ergodic, retaining only memory of the particle masses, the total interaction energy and the total angular momentum.Stone and Leigh [11] showed that two or more such scrambles are a sufficient criterion to remove most of the regular regions of phase space, leaving only the ergodic subset.In this paper, if a three-body interaction satisfies the criterion N S ≥ 4, where N S 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 correspond to the chaotic region and 17.9% to the regular region.

Method
In 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 to compute the time evolution of our three body systems, before moving on to describing how we define our initial conditions and set up our experiments.

N-body Integrator
We use TSUNAMI [37], a few-body code ideally suited to follow the time evolution of strongly interacting 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 accuracy of 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 more details on the implementation and performance of TSUNAMI, we refer the reader to [46].

Initial conditions
In the following sub-sections, we describe how the initial conditions for our suites of numerical simulations are defined.

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 are P 1 = (x 1 , y 1 , z 1 ) = (1.5, 1.5, 0) and P 2 = (x 2 , y 2 , z 2 ) = (2.5, 1.5, 0).Note that the triangle side is 1 au.This represents a perfect equilateral triangle, where each side will have a length 5 in Cartesian coordinates.The initial position of particle P 3 was then calculated as: We set z = 0 , and focus on purely 2D coplanar configurations.

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 r gen from the vertex, as shown in Figure 2. First, we randomly draw a perturbation radius r within a circle of r gen .Since r represents a radial perturbation, we randomly sample it from a distribution uniform in r 2 gen .The position of each particle 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 the range (0, 2π).
We generate three different sets of simulations with different perturbation magnitudes, in order 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 r gen = 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 r gen = 10 −7 au.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 using a 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 to avoid 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 are always 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 r gen (Column 4).
As shown in Figure 2, for every new triangle that we generate, we re-compute the lengths of each side, where S 1 S 2 and S 3 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 and the 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 the ejected particle and the left-over binary (see Figure 4 and Figure 9).

Ejection criteria
To determine when a particle is ejected from the system, we check the state of the system at every time step (each simulation has a maximum integration time of 10 5 yrs).The most bound 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 for some systems to fulfill our ejection criteria is longer than our total integration time.This is typically due to long-lived Lévy flights, where the return time of the escaper exceeds our chosen maximum integration time.This represents a negligible fraction of our simulations, and our results remain unchanged without these simulations.

Simulations with post-Newtonian corrections
We use Set A to define our initial particle positions, but set the radius of each particle to be equal 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 the binary 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 500r s for the particle radius and by adopting the sticky-star approximation to determine when the particles merge (i.e., when their radii overlap in both time and space).

Theoretical Expectations
In this section, we describe the theoretical expressions (see Chapter 7 of [47]) for the expected distributions 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 star ejection velocity and total system lifetime, for comparison to our numerical scattering results.

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 defined in Section 2.3.In Section 5, we will perform our comparisons with and without applying these cuts, to better understand how the chaotic and regular subsets are each contributing to the total distribution (i.e., without applying any cuts for ergodicity).

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: The variables m a 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, where z = |E 0 |/|E B | corresponds to the inverse binary orbital energy.For the planar case, this is (5)

Eccentricity distribution
We compare our simulated results to several different theoretical eccentricity distributions in Figure .8.These are described below.
The thermal distribution for the binary eccentricity, which assumes a detailed steady-state balance between binary creation and destruction [48], is given by Equation 6: In the planar case, the eccentricity should be distributed as Equation 7f (e) de = e(1 − e 2 ) −1/2 de .
The power-law index on (1 − e 2 ) is taken to be −1/2, which is appropriate for the low angular momentum limit [8,49].

Escape velocity distribution
We adopt Equation 7.19 in [47] for the escape velocity distribution shown in Figure .8.This is appropriate to the planar case and is given by

Results
In this section, we present the results of our numerical scattering experiments.We begin by presenting 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 orbital energy, orbital eccentricity, escaper velocity and total system lifetime.We then compare these to the expected theoretical predictions, using our phase space plots to identify the origins of the disagreements between our theoretical calculations and the simulated data.Each subsequent inset corresponds to a zoom-in on the regions in the phase space.

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, the parameters chosen for the X and Y axes are the ratios between the lengths of the sides of each of the generated triangles, as shown in Figure .3. Here, S 1 , S 2 and S 3 correspond to the lengths of the sides opposing, respectively, particles P 1 , P 2 , and P 3 .We then plot the ratios S 1 /S 2 and S 1 /S 3 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 located at the coordinates (1, 1), and corresponds to the initial conditions for a perfect equilateral triangle.
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).

Full phase space
The central panel in Figure 4 shows the phase space limited to 0 < S 1 /S 2 < 3 and 0 < S 1 /S 3 < 3 for the interactions of Set A.
Near the singularity, we see regular swaths of uniform colour, separated by chaotic regions of phase space.In the regular regions, small perturbations to the initial state (i.e., equilateral triangles 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 initial conditions.
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 the sides of the triangle is very small.The arm that goes to the top represents the deformation of the triangle where particle P 3 (coloured by orange) is initial far from particles P 1 and P 2 .Similarly, the arms that goes to the right and the left bottom represent the deformation of the triangle where the particles P 2 (coloured by blue) and P 1 (coloured by red) are far from the others.

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 the zoom-ins becoming more extreme going from left to right.We see regular coloured regions with 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 to uniform 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-in on 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 colouring in 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.

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 S 1 /S 2 and S 1 /S 3 .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 the binary-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 phase space.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 and appear to continue indefinitely down to increasingly small scales.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 panels e and g.

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 S 1 /S 2 = S 1 /S 3 (left column) and the binary-single regime S 1 ≫ S 2 , S 3 (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 around the singularity and the self-repeating structures in the binary-single regime.

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 scattering experiments, 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 shown in 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 dark-  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 not extend our results to systems with different particle masses.

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.

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 N S ≤ 4) and ergodic (all simulations with N S > 4) subsets, along with the entire suite of simulations in Set A. The criterion for separating the two subsets is N S ≥ 4, as discussed in Section 2.3.
In the binary eccentricity distribution panel (top left panel), both subsets (i.e, ergodic and regular) show a strongly suprathermal distribution.Concerning the escape velocity distribution, the ergodic subset is similar to a log-normal distribution while the regular subset shows an irregular distribution.Both of them peak near 10 4.7 m/s.Further differences between the ergodic 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, the distribution looks like a log-normal with the peak near 4 yr (log 10 (lifetime/yr) ∼ 0.6), while the regular one peaks near 10 yr (log 10 (lifetime/yr) ∼ 1).
In the middle panels of Figure 8, the histograms show the distributions of initial system energies, 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 E 0 /E B ≲ 0.6) and obtain not enough wide binaries (i.e., binary stars with energies E 0 /E B ≳ 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 initial energy.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 a consequence 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 zero angular momentum distribution (Equation 7).Relative to the theoretical expectation for the  ejection velocities (Equation 8), we obtain an excess of high-velocity escapers and a lack of low-velocity escapers.

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 statistical theory.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 to carry 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 to zero).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 of where a norm is the semimajor axis corresponding to the initial energy E 0 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 high ejection velocities.In the limit of the escaper velocity becoming very low, we tend to observe  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 L S ∼ 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 angular momentum.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 the main regular regions reaching the longest lifetimes, and correlating with low ejection velocities, high eccentricities, and compact binaries.Very near to the singularity, the very short interaction 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 is satisfied), consistent with the results of [38].

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 done in an effort to identify their origins using the post-interaction distributions of binary orbital energies, orbital eccentricities and escaper velocities.In this section, we divide our total suite of simulations into ergodic and regular subsets.To do this, we select our data following the the scramble number criterion discussed previously in Section 2.3 (N S > 4), and isolate most of 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 the parameter space is chaotic.One of the main inconsistencies of the data with our theoretical expectations is the bump at 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 for the ergodic and regular subsets for both E 0 and E B .We find a misalignment between the peaks of the regular E 0 and E B 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 the inverse binding energy distribution, as seen in the right-most middle panel of Figure .8. In fact, the ratio between the medians of the E 0 and E B distributions matches the peak in the inverse binding energy distribution for the entire suite of simulations, therefore confirming that the origin of the misalignment is in the regular, and not in the ergodic, subset.

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 three particles 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 position in 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 10 7 simulations using the TSUNAMI code [37].The large suite of simulations 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 distributions of binary orbital energy, orbital eccentricity, escaper velocity, and total interaction lifetime.We then compare these distributions to the standard theoretical expectations provided in the literature.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 more compact 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 interactions of 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 no rotation.

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 of our simulations.

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 angular momentum in the cluster, and the stellar orbits are such that three single stars that meet near the cluster centre of mass should be on nearly radial orbits, causing the three stars to approach their common centre of mass with small impact parameters.In the opposite limit of high 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 centre of 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]: where n is the number density, m is the average mass and v rms is the 3D root-mean-square velocity.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 taken from 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 scatterings of all single stars do occur in the densest Milky Way globular clusters well within the lifetime of the cluster.Post-core collapse clusters are the most likely to harbour binaries formed via 1+1+1 interactions.

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 for binaries 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, the merger times will be short, since the timescale for a merger to due gravitational wave emission is proportional to the binary orbital separation ∝ a 4 (i.e., more compact binaries merge faster) and also proportional to (1 − e 2 ) 7/2 for the orbital eccentricity [53].Nonetheless, depending on 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, even at very high eccentricities.This will be investigated in the future using dedicated numerical experiments with post-Newtonian corrections.

Figure 1 :
Figure 1: Each panel shows the trajectories in position-space for all three particles for different 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 and the other two form a binary system.

Figure 2 :
Figure 2: Scheme of the initial setup.The equilateral triangle vertices are labeled as P 1 , P 2 , and P 3 .The generated particles corresponding to the perturbed triangle are labeled as P 1 gen , P 2 gen , and P 3 gen .The green circles correspond to the boundaries of each generated particle (r gen = 0.5 au).The insets in the upper right show the equilateral triangle and the perturbed version.

Figure 3 :
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 the perturbed position to the vertex of the triangle.

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

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

Figure 5 :
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 panels e and g.

Figure 6 :Figure 7 :
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.

Figure 8 :
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 N S ≤ 4 , the red line corresponds to the ergodic subset where N S > 4. Bottom panels show the ergodic subset along with the theoretical curves shown in 4.

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