The chaotic four-body problem in Newtonian gravity - II. An ansatz-based approach to analytic solutions

In this paper, we continue our analysis of the chaotic four-body problem by presenting a general ansatz-based analytic treatment using statistical mechanics, where each outcome of the four-body problem is regarded as some variation of the three-body problem (e.g., when two single stars are produced, called the 2+1+1 outcome, each ejection event is modeled as its own three-body interaction by assuming that the ejections are well separated in time). This is a generalization of the statistical mechanics treatment of the three-body problem based on the density-of-states formalism. In our case, we focus on the interaction of two binary systems, after which we divide our results into three possible outcome scenarios (2+2, 2+1+1, and 3+1). For each outcome, we apply an ansatz-based approach to deriving analytic distribution functions that describe the properties of the products of chaotic four-body interactions involving point particles. To test our theoretical distributions, we perform a set of scattering simulations in the equal-mass point particle limit using FEWBODY . We compare our final theoretical distributions to the simulations for each particular scenario, finding consistently good agreement between the two. The highlights of our results include that binary-binary scatterings act to systematically destroy binaries producing instead a single binary and two ejected stars or a stable triple, the 2+2 outcome produces the widest binaries and the 2+1+1 outcome produces the most


INTRODUCTION
The four-body problem has scarcely been studied analytically.This can be understood, at least in part, upon considering the history of the three-body problem, and its notoriety for being a strong example of chaos in nature.With the twobody problem solved, the temptation to find an analytic solution to the three-body problem attracted many researchers.Generally, the goal was to predict the positions of the particles at any future time, for any set of initial conditions.This eventually led to the understanding that the addition of even one extra particle (relative to the two-body problem) renders the number of variables in the equations of motion greater than the number of equations.The problem is unsolvable, and more particles will only make it worse.Consequently, the four-body problem received little attention for centuries (e.g.Nash & Monaghan 1980).
More recently, computational advances have allowed for numerical studies of the four-body problem (e.g.Harrington 1974;Saslaw, Valtonen & Aarseth 1974;Mikkola 1983Mikkola , 1984a,b;,b;Rasio, McMillan & Hut 1995;Fregeau et al. 2004;Leigh & Geller 2012;Leigh et al. 2016;Ryu, Leigh & Perna 2017a,b,c).Ignoring planetary dynamics, most of these studied scattering interactions between two binary star systems.For example, Mikkola (1983) confirmed that stable triple systems form during encounters between identical binaries.Mikkola (1984a) extended this result to include binaries with different initial orbital energies, finding in the process that significantly more triples form as the ratio of binding energies increases from unity.
The primary astrophysical motivation for this paper is binary-binary scatterings in dense stellar systems, such as open, globular, or nuclear star clusters.In such systems, the rate of binary-binary scatterings can dominate the rate of binary-single scatterings provided the binary fraction satisfies f b ≳ 10% (Sigurdsson & Phinney 1993;Leigh & Sills 2011).In this case, binary-binary scatterings are the dominant cluster heating source and the 4-body problem's scattering outcomes become critical for understanding the thermodynamic evolution of the host star cluster.Moreover, one possible outcome of a binary-binary scattering (which does not occur for binary-single scatterings in the point particle limit) is the dynamical formation of a stable triple star system (e.g.Leigh et al. 2016).Such triple systems are of great interest for their susceptibility to the Kozai-Lidov mechanism, which can create accreting inner binaries or exotic astrophysical transients (e.g.Perets & Fabrycky 2009).The decay products of binarybinary scatterings are thought to have been observed directly, both in the form of stable triples (e.g.Leigh & Sills 2011;Leigh & Geller 2013) and runaway O/B stars (e.g.Hoogerwerf et al. 2001;Oh et al. 2015).
The strongly chaotic nature of the generic four-body problem makes an analytic solution to a set of specific initial conditions impossible (except for some fine-tuned sets of measure zero).However, chaos becomes a useful tool if we are interested in probabilistic distributions of outcomes corresponding to distributions of initial conditions.Following the pioneering work of Monaghan (1976a,b, hereafter the "Monaghan formalism"), we employ statistical mechanics to compute distributions of outcomes for the generic problem of binary-binary scatterings.This paper is a more analytic continuation of our previous work, which examined binary-binary scattering using a large suite of numerical integrations (Leigh et al. 2016, hereafter "Paper I").
In §2, we provide an overview of the geometry and phase space of the problem.In § 3, 4, and 5, we use the Monaghan formalism to compute phase space volumes and parameter distributions for the three possible outcomes of a binarybinary scattering event involving point particles.In § 6, we compute branching ratios for these three outcomes.In § 7 we present the toolkit FEWBODY to describe and explain the numerical simulations we use to test our analytic distributions.
In § 8 we show our results comparing the simulated data with the analytic derivations found in the previous sections.In § 9 we consider the astrophysical implications of our results, and in §10 we summarize our work.

STATISTICAL MECHANICS OF THE FOUR-BODY PROBLEM
In this paper we consider the outcome of the generic, nonhierarchical four-body problem.We assume all stars are point particles and neglect tidal forces, physical collisions, the effects of general relativity, and non-gravitational forces.The four interacting stars have masses ma, m b , mc, and m d , which may differ from each other.The interacting system has a conserved total energy, E0, and a conserved total angular momentum, ⃗ L0.Our approach, which was first developed in the Monaghan formalism to treat the chaotic three-body problem, is to consider the statistical phase space of different outcomes of the four-body problem.By calculating the phase space volume of a single outcome as a function of parameters of interest (e.g.ejection velocities), we can construct distributions of these outcome parameters.By calculating the relative volumes of different outcomes, we can compute branching ratios between them.
There are four possible outcomes of the generic four-body problem.Following a phase of chaotic, non-hierarchical gravitational interactions, the system eventually forms some combination of hierarchical, bound particles and unbound, escaping particles.The four possible combinations are (i) One escaping, unbound star and a hierarchically stable triple (the "3+1" outcome).
(iii) Two surviving binaries which are mutually unbound from each other (the "2+2" outcome).
In isothermal star clusters with Maxwellian velocity distributions, the 1+1+1+1 outcome is almost always energetically forbidden (Leigh et al. 2016), so we ignore it for the remainder of this work.
The primary assumption we make in analyzing the 3+1, 2+1+1, and 2+2 cases is that the strongly chaotic interactions of a non-hierarchical four-body system will uniformly populate the phase space available in each of these outcomes.
The same assumption motivated the statistical mechanical treatment of the three-body problem in Monaghan (1976a,b), and was validated by post-hoc checks using numerical orbit integrations.Because the parameter space of the four-body problem is much larger than that of the three-body problem, we cannot hope to fully cover it with numerical integrations, but we will check our results against a suite of four-body scattering experiments.
Although our primary motivation is to make statistical predictions for the outcomes of binary-binary scatterings, our results should apply equally well to other non-hierarchical four-body systems, such as a strong triple-single scattering, or to other, less probable events (such as a simultaneous encounter between a binary and two single stars).
In the following three sections, we consider the 2+2, 3+1, and 2+1+1 outcomes.It is important to note that variables in these sections may sometimes have the same names but different definitions.We define all variables locally at the beginning of their respective sections.

THE 2+2 OUTCOME
We begin with the 2+2 case, as it has the greatest similarity to the classic 2+1 three-body problem.The total energy of this final state is E0 = Es + EB1 + EB2, where In the above equations, ⃗ r ab is the separation vector between mass a and mass b, ⃗ r cd is the separation vector between mass c and mass d, and ⃗ rs is the separation vector between the two binary centers of mass, as is shown in Fig. 1.
There is likewise a well-defined final state angular momentum, The two final state binaries in this outcome have characteristic semi-major axes a1 and a2.We choose a1 > a2.If a1 ≫ a2, then the 2+2 problem reduces to the standard Monaghan 2+1 formalism with extra degrees of freedom.More specifically, we can apply the loss cone formalism while treating the second binary as a point particle (although we must account for its reservoir of energy and angular momentum).
Operating first in this limit, we simplify further by working also in the low angular momentum limit.In this case, the (3) We eliminate three variables of integration with the following simplification: Under the assumption that a1 ≫ a2, we use the areal loss cone factor Here α is a dimensionless fudge factor of order unity that ultimately must be calibrated from numerical scattering experiments.In the classic 2+1 problem, α ≈ 7. We next evaluate In the last approximate equality, we have taken R ≲ 3a1.The density of states for this outcome is now Further simplification (making use of isotropy) gives It is important to mention that in this derivation we have called me the mass of the object that we consider is the one that escapes (analogous to the three-body case in which the escaper is a single), which in this case is a binary (so mB1 = mB2 = me).For simplicity we will define A := e M, however it is important to remember that the value of me and M will be different for each of the three different outcomes.Eq. 9 encodes the final probability distribution of binary energies (Eq.8 does not because its limits of integration depend on EB1 and EB2).Specifically, these probability distributions are However, in this simple derivation, we have followed (Valtonen & Karttunen 2006), and have neglected angular momentum conservation, rendering our results inaccurate.We now consider the power-law ansatz of (Valtonen & Karttunen 2006), In the zero angular momentum limit of the chaotic threebody problem, n = 3, leading us to speculate that here n = 3 and ν = 1.This, however, would imply a UV divergence in P (|EB2|) as EB2 → ∞.We fix this by truncating at a maximum energy Emax motivated by physics beyond the Newtonian point particle limit of this paper (e.g.physical collisions, tidal interactions, relativistically unstable orbits, etc.).Thus, These approximations will begin to break down if a1 ≈ a2.

THE 3+1 OUTCOME
In Paper I, we demonstrated that to a good approximation, the distribution of escaper velocities in the 3+1 outcome can be computed with a straightforward application of the 2+1 Monaghan formalism, considering only the binding energy of the inner binary EB.This approach works because the total binding energy, E1, of a hierarchically stable triple (at least in the equal mass case considered in Paper I) is dominated

Strongly Hierarchical Triples
When masses are comparable, the stable triple produced in the 3+1 outcome can be thought of in the following way: the inner binary contains the bulk of the energy E1, while the outer binary contains the bulk of the angular momentum ⃗ L1.In paper I, we showed that the standard 2+1 Monaghan formalism applies reasonably well to the binding energy distributions, so now applying factor A (with their respective values) we found where n = 9/2 if angular momentum conservation is neglected and n = 3 in zero-angular momentum ensembles.We can determine the distribution of outer triple binding energies at a similar level of approximation by assuming a thermal distribution of outer eccentricities eT, i.e. dN/deT = 2et.Then, under the assumption that the outer orbit angular momentum LT ≫ LB, we use the relation e or, equivalently, For fixed LT, Eq. 17 specifies the distribution of |ET|, which varies from 0 to a maximum value ≲ |E1|/2.
We can proceed further under the assumption of small angular momentum in the inner binary; in this limit, the angular distribution of ⃗ LB will be approximately isotropic with respect to ⃗ LT.If we define a reference axis ẑ ∥ ⃗ LB, then the distribution of misalignment angles where cos θ ≡ L1 • LB = L z 1 /L1.In the remainder of this section, we denote the z component of a vector with a superscript z, and the components orthogonal to ẑ with a superscript ⊥.
We now complete this perturbative calculation: having assumed a distribution of θ which is isotropic, we wish to know the distribution of a different misalignment angle, cos ψ ≡ LB • LT = L z T /LT.In general, ψ ≈ θ, but we aim here to quantify the leading order deviation from isotropy in ψ.Since and the quadratic formula then provides The final distribution of interest is dN/dψ = (dN/dθ)(dθ/dψ), which evaluates to where and both cos θ and sin θ can be computed from Eq. 20.

THE 2+1+1 OUTCOME
The 2+1+1 case is in some ways the most distinct from the classic 2+1 problem.It differs not only in additional degrees of freedom, but more fundamentally in its complicated causality.The 2+2, 3+1, and 2+1 scenarios terminate a single impulsive escape, but the 2+1+1 does not, and two different escape events must be considered.In Fig. 3, we show a cartoon of the 2+1+1 outcome.Here masses ma and m b form the survivor binary, while masses mc and m d are escaping on unbound orbits.The order of escape matters: the distribution of parameters for the first escaper (mc) will differ from the parameter distribution for the second escaper (m d ).
We begin by making an approximation of sequential escape: we assume that in general, a metastable triple is formed after the escape of particle C, and that particle D is only ejected after the gravitational influence of C becomes negligible.This approximation lets us apply the standard 2+1 Monaghan formalism in an iterated way.Based on the numerical scattering experiments of Paper I, we believe it to be well justified for low virial ratios (k ≪ 1) but a poor approximation for high virial ratios (k ≈ 1), when both ejected stars are ejected almost simultaneously.We first estimate the distribution of binding energies ET of the metastable triple: As before, E0 is the conserved total energy of the four-body encounter.For a given ET value, we can take the standard 2+1 distribution of binding energies for EB (the binding energy of the final surviving binary), but if we want a distribution of ET values, we need to integrate over Eq. 23: More generally, if we substitute a power law index n for the triple binding energy distribution (9/2 above) we find Likewise, we can apply the results of the standard 2+1 formalism.

BRANCHING RATIOS
The "branching ratio" defines the probability of obtaining a given outcome, for a given total encounter energy and angular momentum.Hence, for the chaotic four-body problem in the point-particle limit with total energy E < 0, there are three branching ratios to consider.In general, the relative fractions for these different outcomes must be determined using numerical scattering simulations.However, as we are about to show, these branching ratios can also be computed analytically, if all particles are identical.Consider performing N0 simulations of a chaotic four-body interaction involving identical point-particles, with nearly identical initial conditions.Given N0 simulations, we must obtain N0 outcomes.Then, the total number of simulations performed can be written: where N2+1+1, N3+1 and N2+2 correspond to the number of simulations resulting in, respectively, the 2+1+1, 3+1 and 2+2 outcomes.Now, by conservation of energy, we must also find that the total amount of energy put in across all simulations is equal to the total energy we get back out.In other words: Similarly, the total angular momentum must be conserved, ensuring that the following must be true: Each term on the right-hand-sides of Equations 27 and 28 must be broken up in to the individual contributions from each decay product (i.e., single, binary and/or triple star(s)).
The limits of the resulting integrals must then be chosen appropriately.For example, for the 2+1+1 outcome, we have: where the indices 1 and 2 correspond to, respectively, the first and second ejected single stars, and all distributions correspond to those presented for the 2+1+1 outcome.Note as well that f (ES,i) = f (vS,i)/(mS,ivS,i).
If we divide both sides of all three equations by N0, then Equations 26, 27 and 28 constitute three equations, each with the same three unknowns.Hence, this system of equations is solvable.The factor in front of each term corresponds to the branching ratio for that outcome.
We caution that if the particles are not identical, then the formalism presented here for calculating branching ratios for the different outcomes is no longer valid, strictly speaking.We defer this issue to a future paper, along with a more thorough comparison between the predicted branching ratios and the results of numerical scattering simulations.

METHODS
In this section, we describe and present the numerical scattering simulations used to test directly the analytic distribution functions derived in the preceding sections.

Numerical scattering simulations
The numerical scattering simulations used throughout this paper are the same as presented in Leigh et al. (2016).For completeness, we repeat our description of the code and initial set-up here.
We calculate the outcomes of a series of binary-binary (2+2) encounters using the FEWBODY numerical scattering code1 .The code integrates the usual N -body equations in configuration-(i.e.position-) space in order to advance the system forward in time, using the eighth-order Runge-Kutta Prince-Dormand integration method with ninth-order error estimate and adaptive time-step.For more details about the FEWBODY code, we refer the reader to Fregeau et al. (2004).
The outcomes of these 2+2 encounters are studied for the initial virial ratio k = 0, where k is defined as: here the indexes 1 and 2 correspond to the two initial binaries.
The initial kinetic energy corresponding to the centre of mass motion of binary i is: where mi = mi,a + m i,b is the total binary mass and v inf,i is the initial centre of mass velocity for binary i.The initial orbital energy of binary i is: where mi,a and m i,b are the masses of the binary components and ai is the initial orbital separation.Given this definition for the virial ratio, k = 0 corresponds to the binaries starting from rest, and maximizes the fraction of longer-lived chaotic interactions (which is a necessary prerequisite to apply the Monaghan formalism).
All objects are point particles with masses of 1 M⊙.All binaries have ai = 1 AU initially, and eccentricities ei = 0. We fix the impact parameter at b = 0 for all simulations.The angles defining the initial relative configurations of the binary orbital planes and phases are chosen at random.
We use the same criteria as Fregeau et al. (2004) to decide when a given encounter is complete.To first order, this is defined as the point at which the separately bound hierarchies that make up the system are no longer interacting with each other or evolving internally.More specifically, the integration is terminated when the top-level hierarchies have positive relative velocity and the corresponding top-level Nbody system has positive total energy.Each hierarchy must also be dynamically stable and experience a tidal perturbation from other nodes within the same hierarchy that is less than the critical value adopted by FEWBODY, called the tidal tolerance parameter.For this study, we adopt the tidal tolerance parameter δ = 10 −7 for all simulations.2This choice for δ, while computationally expensive, is needed to maximize the accuracy of our simulations, and ensure that we have converged on the correct encounter outcome (see Geller & Leigh 2015 for more details).
Because of the isotropy of our initial conditions, the typical four-body encounter we simulate has L0 > 0. If one considers a binary-binary scattering event where the two initial binaries have isotropically oriented angular momentum vectors of magnitude L1 and L2 (L1 ≥ L2 by assumption), then the total angular momentum ⃗ L0 = ⃗ L1 + ⃗ L2 spans a range of magnitudes from L1 − L2 to L1 + L2 with a distribution The first moment of this distribution is If we specialize now to our initial conditions (equal masses m, equal initial semi-major axes, L ≡ L1 = L2), we find that where we have followed Valtonen & Karttunen (2006) in defining the maximum system angular momentum Lmax ≡ 5 2 G m 5 /|E0|.Their numerical fitting formula for the classic 2+1 problem predicts n = 3 + 18 L2 for ensembles of resonant three-body encounters with angular momentum L = L/Lmax.This gives us a naive expectation of n ≈ 5.6 for our numerical simulations.

RESULTS
In this section, we compare the results of our numerical scattering simulations to the fitting formulae presented in the previous sections.

Comparing to the simulations
In Figure 4, we present the final outcome distributions after the interaction between our two initial binaries.We separate our results for different semi-major axis ratios (a1/a2, indicated on the x-axis) and show the fraction of simulations resulting in each of our three possible outcomes (i.e., 3+1, 2+1+1 and 2+2).For each combination of a1/a2 we perform 10,000 scattering simulations.The results for the 2+2, 3+1 and 2+1+1 cases are shown in red, black and blue colour bars, respectively.Note that this colour relationship will be used for all figures in this paper.We note that these branching ratios should be computable explicitly using the equations in Section 6, but we defer this to a future paper, since we first need to re-derive the analytic functions in this paper from first principles as in Stone & Leigh (2019) or Ginat & Perets (2021) to obtain the needed angular momentum dependences.We see that for the case in which our initial semi-major axes are relatively similar (1 ≤ a1/a2 ≤ 4) the largest outcome fraction corresponds to the 2+1+1 scenario.This is in agreement with what was shown by Mikkola (1983) and Leigh et al. (2016) for identical initial conditions.The trend changes as the ratio of the semi-major axis increases.For a1/a2 > 8 the scenario that occurs the most corresponds to the 3+1 case, that is a particle of the system is transformed into an escaper leaving behind a dynamically stable triple system.Finally, the formation of the 2+2 case is always the least probable, tending to decrease as the semi-major axis ratio increases.This is because for this to occur, both binaries must form simultaneously, which effectively requires that two stars be ejected in similar directions, with similar ejection times and escape velocities, rendering this outcome improbable.Alternatively, this can be viewed as the wider binary struggling more and more to eject the more compact binary (in analogy to the the three-body case) as the ratio of semi-major axes increases at fixed particle mass.
In Figure 5, we present normalized histograms of the final distributions of binary orbital energies, parameterized using the total encounter energy E0 as z = E0/EB.This is the result of our binary-binary scattering experiments, for different values of the semi-major axis ratio (a1/a2, as in Figure 4).For the 2+2 case, we show the orbital energies of both final binaries, divided according to their final orbital energies using the solid and dashed lines for the compact and wide binary, respectively.For the 3+1 case, we show only the orbital energy of the inner binary of the resulting triple.
For the energy range between 0 ≤ |E0|/|EB| ≤ 1, the analytical distributions reproduce very well what was obtained through the simulations.Especially for the 2+1+1 case, the curve is wholly reproduced, which shows that our assumption that this outcome can be modeled by applying the three-body disintegration scheme twice, by assuming additionally that each escape is well separated in time, work quite well.The same happens for the 3+1 case, assuming that all the interaction energy is held in the inner binary of the triple system, while all the angular momentum is kept in the outer orbit of the same system.In the 2+2 case, we see that the distribution fits well for the compact binary provided |E0|/|EB| ≤ 1, while for the case of the wider binary the distribution fits poorly.This is likely because our assumptions here begin to break down such that very few of these interactions are truly chaotic, and hence very few of our simulations should actually produce results that agree with theory.This is supported by the fact that the good agreement between our results and theory begins to decrease as the ratio of the semi-major axis increases, particularly for a1/a2 ≥ 8.This last point for the 2+2 configuration, and the distributions for which E0/EB > 1, will be taken up again in Section 9.
The binding energy distribution can be used to derive the escape velocity distribution f (ve)dve for the escaping star(s), as given by Equation 37 with me = m b /3 = M/4 for the 3+1 case for all identical particles (Leigh et al. 2016).This gives the following functional form (Equations 7.19 and 7.26 in Valtonen & Karttunen 2006): In Figure 6 we show, for different initial a1/a2 ratios, the distribution of escape velocities for the single star for a 3+1 outcome (black lines), as well as for both single stars for a 2+1+1 outcome (blue lines), where the solid lines show the simulated data and the dotted lines show the analytic fits.Note that to calculate the distribution analytically using Equation 37, we use the values n = 4.5 to account for the angular momentum dependence and me = m b /3 = M/4 for the mass of the escaping particle(s).
Figure 7 shows the distributions of escape velocities for binaries for the 2+2 outcome (red lines).The solid lines show the results of our numerical scattering simulations, whereas the dotted lines show our analytic fits.Note that by conser-  [1,2,4,8,16,32].
vation of linear momentum, the escape velocity distributions are equivalent for both binaries, since we are dealing with the equal particle mass case.
We clearly see that the analytic distribution of the escape velocities shows a clear agreement with what was seen in the simulations.On the other hand, there is a clear tendency for poorer agreement between the theory and the simulations as the semi-major axis ratio increases (as in Figure 5).The velocity range as well as the corresponding outcome fractions for the 3+1 and 2+1+1 cases show very similar distributions as those found in (Leigh et al. 2016), but in our case for n = 4.5.The same is true for the 2+2 case where our distribution behaves as expected for n = 4.5.
In Figure 8, we show the final distributions of orbital eccentricities for every encounter outcome, including the inner orbits of stable triples.The solid black line shows a thermal eccentricity distribution f(e)de = 2e, which matches the simulated data quite well for all orbits and all encounter outcomes.

DISCUSSION
In this section, we discuss possible caveats of our work.
First, we note that the agreement between our simulations and the analytic fits is consistently good for small initial semimajor axis ratios but becomes poorer as this ratio increases.This is because, at least in part, a smaller fraction of the interactions become chaotic, which is a prerequisite to performing this comparison.Hence, we are left with fewer simulations that should be in agreement with theory and must omit those simulations that ended deterministically.Most of the interactions that cause this turn out to be simple exchanges.In other words, in the limit of a large semi-major axis ratio, the probability that the compact binary will simply be exchanged into the wide binary becomes high, liberating one single star in the process and forming a dynamically stable hierarchical triple system.Hence, fitting a three-body interaction, treating the more compact binary as a heavy single particle, is a more appropriate model in this limit.[1,2,4,8].All histograms have been normalized by the total number of simulations that resulted in the corresponding outcome.Note that for the 2+2 case (red color) there is a solid line and another dashed line, which represent the simulated values for the compact and wide binary, respectively.
In Figure 5, we show the simulated distributions of leftover binary orbital energies for the most compact orbit in the final outcome and compare to our analytic fits.We see good agreement between the two for E0/EB < 1, since beyond this limit significant angular momentum is contained in the final orbit and our analytic fits do not account for any angular momentum dependence.For example, in the 3+1 case, the inner binary also contains angular momentum but we assume that it does not.For this case, the angular momentum in the inner and outer orbits ultimately limits the minimum binary energy in the interior because the outer orbit cannot contain enough angular momentum to accommodate a dynamically stable hierarchical triple.This itself explains why the 2+2 distributions also tend to zero, although at lower minimum binary energies, since wider binaries can be accommodated in this case, since there is no requirement mentioned above for dynamic stability in the triplet case.for normalized distributions of escape velocities from the single star (in km/s) for the 3+1 (black color) and 2+1+1 (blue color) outcomes.The different insets show the same semi-major axis ratios, number of simulations, line types and colours as in Figure 5.The dotted black line shows the distribution of escape velocities calculated using Equation 37 for a 3+1 outcome and assuming n = 4.5, corresponding to approximately isotropic scattering.The dotted blue line shows the same thing but for the 2+1+1 case assuming n = 4.5.For both analytical curves we assume me = m b /3 = M/4.Note that for the 2+1+1 case there are 2 solid blue lines, which correspond to the 2 escapers in this scenario.
In the same figure, we can see that in the case of the wide binary in the 2+2 scenario, our approximation does not fit the simulated values.This occurs because the observed values are given in the domain E0/EB > 1, which means that in this binary system there is a large amount of angular momentum.Therefore, since our analytic formalism does not explicitly take into account the angular momentum dependence, it is perhaps not surprising that we do not see good agreement between the simulated data and theory in this domain.
In Section 7, we show that we expect the analytic distributions to match the simulated ones when n ≈ 5.6.When performing our comparisons, we adopt a value of n = 4.5 which shows the best agreement.However, we note that our expected value of n = 5.6 also does a good job of describing the data.
This difference is probably due to the fact that we are not considering the total angular momentum dependence in Figure 7.The same as in Figure 6 but for both binaries in the 2+2 scenario.The different insets show the same semi-major axis ratios, number of simulations, line types and colours as in Figure 5.The dashed red line shows the distribution of escape velocities calculated using Equation 37 and assuming n = 4.5, with me = m b = M/2.Note that there is only one solid red line, since both binaries (wide and compact) have the same escape velocity distribution due to conservation of linear momentum.our derivation.Moreover, for our simulations we incorporated different initial semi-major axes in the four-body interaction, while for the initial derivation we assumed equal semi-major axes.
In Figure 8, we see that the eccentricity distribution tends to be quite similar to the thermal distribution.While this is theoretically expected for the three-body case, assuming a detailed balance between binary creation and destruction (Heggie 1975), we are not aware of any expectation for the four-body case.Nevertheless, the reason we see this agreement is probably the same as argued in (Heggie 1975) for the three-body case.This is because a thermal distribution is expected for ergodized outcomes in three-body interactions, and since in this paper we treat each decomposition of the four-body case as a variation of the three-body decomposition, this distribution makes sense (e.g., the 2+1+1 case is modeled as two sequential disruptions of three-body systems).
On the other hand, we see that the distribution of eccentricities corresponding to the external component of the triple system (green dashed line) shows a distribution that does The solid black line shows the distribution of eccentricities for the inner binary of the triple system for the 3+1 case while the dashed green line shows the same for the outer orbit of the triple system for the 3+1 case.The red solid line shows the eccentricity distribution for the compact binary in the 2+2 case, while the red dashed line shows the distribution of eccentricities for the wide binary in the same case.For comparison, we plot a black dashed line showing a thermal eccentricity distribution f(e)de = 2e.The different insets show the same semi-major axis ratios and number of simulations as in Figure 5.
not quite match the thermal distribution for values close to 1 (highest eccentricities).Here we see a paucity of triples with high outer eccentricities relative to a thermal distribution because stable triple systems cannot exist if the external component has very large eccentricity.Otherwise, the triple system will tend to break up, ending up in the 2+1+1 configuration.As expected, this tendency is repeated in each inset (a1/a2 = 1, 2, 4, 8), but we see a tendency for the eccentricity distribution of the outer orbits of stable triples to flatten as the ratio of semi-major axes increases.This is likely because we begin with all binaries being initially circular, and when there is a large ratio between their semi-major axes a simple exchange interaction is the most likely outcome, and here the outer orbit is more often left unaffected, remaining approximately circular.

CONCLUSIONS
In this paper, we have derived analytic distribution functions using the density of states formalism and an ansatz-based approach for the outcomes of four-body (i.e., binary-binary) scatterings in the equal-mass point-particle limit.We have further confronted our analytic fits with the results of numerical scattering simulations, and find good agreement.The highlights of our results can be summarized as follows: • We have derived analytic distribution functions (DFs) to describe the properties of the products of chaotic four-body interactions in the equal-mass point particle limit.These DFs include, for the most compact orbit in the left-over binaries and/or triples in the final outcome state, the distributions of orbital energies for the left-over binaries and triples, the distributions of ejection velocities and the orbital parameters of any left-over binaries or triples.We find good agreement between our analytic theory and the simulations for low semi-major axis ratios, since for larger ratios the angular momentum dependence would need to be integrated into our analytic formalism to expect good agreement in this limit.
• For most of the relevant parameter space, binary-binary scatterings act to systematically destroy binaries by either forming two ejected singles or a stable hierarchical triple instead.
• The 2+1+1 outcome (i.e., one binary and two singles are formed) tends to form the most compact binaries.
• The 2+2 outcome (i.e., two binaries are produced) is consistently the least likely outcome for all ratios of the initial binary semi-major axes, and tends to produce the widest binaries.This is because, in order to form two binaries in the end, effectively the more compact final binary must eject the other two stars at about the same time, in similar directions and with comparable ejection velocities.Alternatively, this can be viewed as the wider binary having more and more difficulty in ejecting the more compact binary (in analogy to the three-body case) as the ratio of semi-major axes increases (at fixed particle mass).
• All outcomes of binary-binary scatterings produce binaries with a distribution of eccentricities consistent with being thermal.This is the case except for very large initial semimajor axis ratios, for which we find a flat eccentricity distribution for the inner binaries of dynamically-formed triples (while that for the 2+1+1 outcome remains consistent with thermal).Since it is those binary-binary scatterings with the largest semi-major axis ratios that produce the most triples, we naively expect that the inner binaries of dynamicallyformed triples should show an approximately flat distribution of orbital eccentricities (in part since we assume initially circular orbits).Finally, for the outer orbits of stable triples, we see a slight deviation from a thermal distribution at high eccentricities, since here very high values for the eccentricity are forbidden if the formed triple is to be dynamically stable.
• We have derived a prediction for the distribution of inclination angles between the inner and outer orbital planes of dynamically-formed stable hierarchical triples.We find that it deviates from an isotropic distribution more and more with increasing angular momentum, potentially allowing for an observational signature to test if triples are primarily formed dynamically (including during the star formation phase).
Fig. 1 shows a cartoon of the 2+2 outcome.The final state here is two binaries; the first, composed of masses ma and m b , has total mass mB1 = ma+m b .The second, composed of masses mc and m d , has total mass mB2 = mc + m d .The binaries have individual reduced masses M1 = mam b /mB1 and M2 = mcm d /mB2, and a joint reduced mass of m = mB1mB2/M , where M = mB1 + mB2.

Figure 4 .
Figure 4.The fraction of different outcomes for our binary-binary scattering simulations as a function of the ratio of the initial semimajor axes of the binaries, or a 1 /a 2 .We vary the ratio as a 1 /a 2 =[1, 2, 4, 8, 16, 32].

Figure 5 .
Figure 5.The distributions of final binary binding energies are shown for each encounter outcome, parameterized using the total encounter energy E 0 as z = E 0 /E B .The colours are the same as in Figure 4.The solid lines represent the values of the simulations while the dotted lines show the values obtained analytically.The black vertical dashed line shows the ratio E 0 /E B = 1.Each panel shows the distributions for a different value of the initial semimajor axis ratio where a 1 /a 2 =[1, 2, 4, 8].All histograms have been normalized by the total number of simulations that resulted in the corresponding outcome.Note that for the 2+2 case (red color) there is a solid line and another dashed line, which represent the simulated values for the compact and wide binary, respectively.

Figure 6 .
Figure6.Comparison between simulations and analytic results for normalized distributions of escape velocities from the single star (in km/s) for the 3+1 (black color) and 2+1+1 (blue color) outcomes.The different insets show the same semi-major axis ratios, number of simulations, line types and colours as in Figure5.The dotted black line shows the distribution of escape velocities calculated using Equation 37 for a 3+1 outcome and assuming n = 4.5, corresponding to approximately isotropic scattering.The dotted blue line shows the same thing but for the 2+1+1 case assuming n = 4.5.For both analytical curves we assume me = m b /3 = M/4.Note that for the 2+1+1 case there are 2 solid blue lines, which correspond to the 2 escapers in this scenario.

Figure 8 .
Figure 8.The distributions of final binary orbital eccentricities are shown for each encounter outcome.The solid blue line shows the distribution of eccentricities for the binary for the 2+1+1 case.The solid black line shows the distribution of eccentricities for the inner binary of the triple system for the 3+1 case while the dashed green line shows the same for the outer orbit of the triple system for the 3+1 case.The red solid line shows the eccentricity distribution for the compact binary in the 2+2 case, while the red dashed line shows the distribution of eccentricities for the wide binary in the same case.For comparison, we plot a black dashed line showing a thermal eccentricity distribution f(e)de = 2e.The different insets show the same semi-major axis ratios and number of simulations as in Figure5.