Three-Body Binary Formation in Clusters: Analytical Theory

Binary formation in clusters through triple encounters between three unbound stars, 'three-body' binary formation, is one of the main dynamical formation processes of binary systems in dense environments. In this paper, we use an analytical probabilistic approach to study the process for the equal mass case and calculate a probability distribution for the orbital parameters of three-body-formed binaries, as well as their formation rate. For the first time, we give closed-form analytical expressions to the full orbital parameter distribution, accounting for both energy and angular momentum conservation. This calculation relies on the sensitive dependence of the outcomes of three-body scatterings on the initial conditions: here we compute the rate of three-body binaries from ergodic interactions, which allow for an analytical derivation of the distribution of orbital parameters of the binaries thus created. We find that soft binaries are highly favoured in this process and that these binaries have a super-thermal eccentricity distribution, while the few hard three-body binaries have an eccentricity distribution much closer to thermal. The analytical results predict and reproduce simulation results of three-body scattering experiments in the literature well.


INTRODUCTION
Binary stars are among the most important astrophysical systems, playing a key part in stellar evolution, and in the dynamics of dense stellar systems.They have a major role in the evolution of globular clusters (e.g.Heggie & Hut 2003;Binney & Tremaine 2008), and effectively serve to stop corecollapse (e.g.Hut 1996;Goodman & Hut 1989).Among the leading channels for the formation of binaries in clusters are primordial binaries formed in the first epoch of star formation in the cluster (e.g.Shu et al. 1987), tidal captures (e.g.Fabian et al. 1975;Press & Teukolsky 1977;Generozov et al. 2018), and dynamical formation of binaries, via three-body interaction -so-called 'three-body' binaries (e.g.Mansbach 1970;Aarseth & Heggie 1976;Stodolkiewicz 1986;Goodman & Hut 1993;Portegies Zwart & McMillan 2000;Heggie & Hut 2003;Atallah et al. 2024), which are the topic of this paper -see, e.g.Pooley et al. (2003) for observational evidence for this process.In such an interaction, three unbound bodies exchange energy, such that at the end two of them remain bound, with the third one serving as a cata-⋆ E-mail: yb.ginat@physics.ox.ac.uk lyst, an energy reservoir where excess energy is deposited.Incidentally, it has been proposed recently that if sufficient amounts of gas are present in a cluster, this could accelerate binary formation by serving as a dissipation mechanisman alternative energy reservoir (Rozner et al. 2023).
In this paper, we wish to study three-body binary formation from an analytical perspective.We seek to offer a complimentary calculation to those of Aarseth & Heggie (1976), who approached the problem in an impulsive approximation, and Goodman & Hut (1993), who used detailed balance to convert an analysis of the inverse problem -that of the ionisation of a binary by a single star, to three-body binary formation.Here, instead, we would like to harness the chaotic nature of the three-body problem (e.g.Valtonen & Karttunen 2006) to derive probabilistic results on threebody binary formation.While previous analytical works on the probabilistic solution of the three-body problem (e.g.Monaghan 1976a,b;Valtonen & Karttunen 2006;Stone & Leigh 2019;Ginat & Perets 2021a,b;Kol 2021) have studied the negative energy case, i.e. a binary-single encounter, here we study the positive-energy case -the formation of a binary from three initially unbound stars.This approach allows us to calculate not only the binary formation rate but also the joint eccentricity and semi-major axis distribution, in a way that accounts for angular momentum conservation as well as energy conservation, during the three-body interaction.
The three-body binary formation can be automatically accounted for in N -body simulations of globular clusters, if close encounters are well resolved, and not smeared by artificial softening (e.g.van Albada 1968;Aarseth 1969;Breen & Heggie 2012a,b;Tanikawa 2013;Wang et al. 2016;Park et al. 2017;Kumamoto et al. 2019;Arca Sedda et al. 2023), but not necessarily in Monte-Carlo simulations (e.g.Giersz 1998;Ivanova et al. 2005;Morscher et al. 2015;Geller et al. 2019;Hong et al. 2020;Rodriguez et al. 2022;Weatherford et al. 2023), which often resort to implementing a prescription for three-body binary formation (Ivanova et al. 2005).The former simulations are often limited, as modelling realistic globular clusters, requires expensive long-term simulations, as such clusters are old, highly populated, and diverse -and thus too computationally expensive -to be accurately modelled by direct N -body simulations.Moreover, few-body simulations do not provide an inherent understanding of the formation process and its dependence on the properties of the encounters.It is thus important to derive an analytic understanding of the processes.Moreover, such analysis can test and replace semi-analytic prescriptions used in Monte-Carlo simulations, which are more versatile than direct Nbody ones.
The paper is structured as follows: we start by calculating the probability distribution of binary parameters, given fixed values of the conserved quantities -the energy and angular momentum (as well as the masses) in §2.This allows us to study the properties of binaries formed this way as a function of these parameters.Then, in §3, we consider the binary formation rate in a system, e.g. a globular cluster, given a distribution of initial positions and velocities.In both sections, we endeavour to keep the lengthier calculations to appendices, while quoting the physical results in the main text.In §5 we summarise the results and highlight their implications.In this paper, we consider strictly the equal mass case, although many of the formulaederived here extend immediately to the unequal case (in particular, most of §2).
In the final stages of this work, Atallah et al. ( 2024) published a numerical study of the three-body binary formation.They analysed the results of three-body numerical experiments with positive energy, which provides us with direct data to test our analytic predictions.We therefore structured the calculation of the binary-formation rate in §3 to allow for a simple comparison with their results.As the results of Atallah et al. (2024) are integrated over a distribution of energies and angular momenta, the calculations of §2 are not immediately comparable with them, (but as tested via their integrals in §3).

RESONANT UNBOUND THREE-BODY SCATTERINGS
Consider three bodies coming from infinity, with masses m1, m2, and m3, and total energy E > 0 and angular momentum J.For point particles, the possible scattering outcomes are: (i) three unbound stars, or (ii) a bound binary and a single unbound star.In this section, we focus on the probability for the outcome (ii), for a given set of conserved quantities, i.e. for given values of E, and the total angular momentum (in the centre-of-mass frame), J.To do so, we have to describe the allowed region of phase-space (in § §2.1 and 2.2).The resultant orbital parameter distribution follows, from which we derive scaling relations with the spatial cut-off (see below) and the hard-soft-boundary.The phase-space volume for a positive-energy threebody system is infinite, so we must require that the three bodies interact inside a region of size R0, which will be described in §2.1 below.For now, suffice it to say that if they don't, then for R0 → ∞, the interaction is not chaotic, so indeed these should be excluded.
If the outcome is a bound binary and a single star, then we label the remnant binary parameters with a sub-script b, and those of the orbit of the single star about the binary's centre of mass with a sub-script s.However, we still do that if there are three unbound stars, simply reserving the subscript b for the pair with the lowest two-body energy, and we still refer to this pair as a 'binary', albeit unbound.
Then, in principle, as long as the outcome is highly sensitive to initial conditions, one can still apply the theory of approximate solutions to the three-body problem, similar to situations with negative total energy.While there are no prolonged democratic resonances -phases during which the particles are in approximate energy equipartition (Heggie & Hut 1993) -in a positive-energy three-body interaction, we contend that the question of binary formation still depends very sensitively on the initial conditions, because the energy exchanges required to form binaries are as large as the total energy, and the general three-body problem is chaotic.Therefore, while one, single scattering problem does not cover all the available phase space, an ensemble of them would, which is the more relevant case for astrophysical applications.The results of this paper can be viewed as a calculation of the distribution of the ergodic subset of three-body binary formation.
Under these assumptions, we find (Ginat & Perets 2021a) that the outcome probability distribution f (E b , S|E, J) for the binary's orbital parameters (those of the outer orbit are automatically given by conservation laws) is where Es = E − E b , etc., and θmax is a function of its arguments, explicitly defined in Stone & Leigh (2019, supplementary information; denoted lmax there).The quantity Rs will be defined in §2.2 below.The second factor of θmax is the usual one (Stone & Leigh 2019), which comes from integrating over the mean anomaly of the outer orbit, explicitly accounting for the finite distance the ejected star can have from the centre of mass of the binary as the triple leaves the strong-interaction region.The first factor of θmax has a similar origin: if E b > 0, the mean anomaly of the binary is also likewise restricted, by R0: essentially, we require that the separation of the binary be smaller than R0.As θmax(R, E, L) ∼ R −3/2 E 3/2 , the first factor of θmax cancels the divergence at E b = 0, while the second one cancels the divergence of E −3/2 s at Es = 0.The parameter Rs is defined as in equation ( 10) of Ginat & Perets (2021a) for the E < 0 case, and in §2.2 below for E > 0.

Cut-Offs
The only thing left to specify is the dependence of R0 on the conserved quantities.We envisage two possibilities for R0: first, one may approximate R0 by the 'virial' radius where we have defined This captures the intuition, that the distance of the closest approach of the triple should be such that the potential energy is of the same order of magnitude as the kinetic energy, at least there.Otherwise, a large energy exchangenecessary for the formation of hard binaries -is unlikely.
Alternatively, the other possibility is to relax this requirement, and simply demand that the distance of the closest triple approach be less than the system's Hill radius, i.e. that at least there, the interaction is a three-body interaction (to leading order), and that the other cluster stars can be neglected.This amounts to where d ≡ n −1/3 , with n being the cluster's number density.
As mentioned above, d may be defined as the binary's Hill radius in its motion in the cluster, by equating the tidal force that the latter exerts on it with its own gravity, in which case d ∝ M 1/3 .For our purposes, what matters is that d is independent of the triple's energy or angular momentum.
The possible definitions are While case 2 might be more justifiable as more generic, the first, more restrictive case is required for significant energy exchanges, if a hard binary is to be produced.We will discuss both choices for R0 in this paper.

Additional Phase-Space Constraints
Contrary to the negative-energy three-body problem, it is strictly impossible for a positive-energy triple to exist in a democratic resonance for an extended time.This implies that there could be triple configurations that cannot dynamically arrive at a binary-single system, despite having the right conserved quantities.Hence, one must impose additional constraints, to ensure that a triple encounter would be able to form a binary.The fundamental assumption of this work is that after these extra ones are imposed, the system is able to reach all the allowed binary-single configurations.Though the above-mentioned various assumptions are not trivial, the successful comparison with the few-body results supports their validity.While for the negative energy three-body problem, in the case that a hard binary forms, the single star can be ejected at a distance rs ≤ R− from the centre of mass of the remnant binary, with (Ginat & Perets 2021a) where we use β = 1.5 as in Ginat & Perets (2023).However, in our case, there are two complications: (i) the remnant binary might not be hard, i.e. it might have a b ≫ Gµ b µM m b |E| , and (ii) it could be unbound completely.
For case (ii), there is no binary formation, which means that there is no restriction on the ejection of star s at all, and in that case, rs ≤ R0 is the only condition.Besides, for marginally bound binaries (case (ii)), the requirement rs ≤ R− does not make sense any more, because it stemmed from a requirement of a hierarchy of energies (Ginat & Perets 2021a), which becomes invalid.Instead, we extend the unbound upper limit to these binaries.For 0 , we therefore also only require rs ≤ R0.We choose ε = 1/10, but have verified that the results remain unchanged up to ε = 1/100.
In summary, the cut-off on the separation rs between the single star and the centre-of-mass of the other two, when the interaction is deemed to have terminated, is rs ≤ Rs, where It is this Rs that enters the θmax in equation (1).

Binary Formation For Fixed Conserved Quantities
When evaluated using techniques similar to those of Ginat & Perets (2021a), the marginal energy distribution of the remnant binary we find is shown in figure 1.
As only the second θmax depends on the inclination, one can derive a marginal eccentricity-energy distribution, given that a bound binary forms, in much the same manner as Ginat & Perets (2023).The result is where ∆ψ and Ap are defined in that work (equation ( 4) there and the text around it).The marginal eccentricity distribution is plotted in figure 2. Due to the extra factor of θmax(R0, E b , S), the result is highly non-thermal, even for large values of the total angular momentum J. Intuitively, the reason that this distribution is super-thermal, is because of the additional factor of θmax(R0, E b , S): this introduces a constraint on the pericentre -namely that it be smaller than R0, which is of course easier to fulfil when the eccentricity is high; thus, there is more phase-space volume associated with high eccentricity that the binary can occupy.The joint energy-eccentricity distribution is plotted in figure 3, for different values of the parameter One can see immediately that the softest binaries -those with large E/ |E b | -tend to be the most eccentric.If the system is in a cluster with temperature T , another dimensionless parameter we will require is ζ, defined by The total probability that a binary form is just the integrand of equation ( 1) over E b < 0. This is plotted in figure 4, as a function of the only non-dimensional combination of conserved quantities available, for equal masses, i.e. of κ, for both options for R0 considered here, for a fixed d.One can see that when J becomes too large, the probability of hard binary formation plummets to 0.
An inspection of figure 4 shows that p bin changes very little with κ, from 0.1257 to 0.2353 in the range plotted.This is true, of course, only for low values of κ, because κ is bounded from above by Sundman's inequality (Sundman 1913), which effectively implies that there exists an upper bound κmax above which p bin = 0. We give explicit expressions for κmax in appendix A. When evaluating the rate in §3 below, we therefore approximate it as a constant, viz.
with p bin being the mean value or p bin (κ), and Θ is Heavi- side's function.p bin still depends on R0, as will be explained below.It is 0.1574 for R0 = 150GM µ/E in the top panel of figure 4.

Hard-Binary Formation
The probability of the formation of a hard binary is This is closely related to the quantity For a fixed value of R0 = d0, p0 depends only on κ, but f bd scales like (energy) −3 , so its integral scales like (energy) −2 ; additionally, κ ∝ E, and therefore We test this scaling in appendix B, where we also explain why in the limit where GM 2 R 0 E ≪ 1. (Let us stress that this is not true for GM 2 R 0 E of order unity.)This, and the above, imply that the entire dependence of p hard on the parameters at hand is This scaling behaviour allows us to calculate p hard for one set of values of the parameters R0 and kBT , tabulate it and use that when computing the hard binary formation probability at any value of the parameters, and any κ.This will simplify the rate calculation in the next section considerably.

BINARY FORMATION RATE CALCULATION
So far, we have explored the binary formation probability for fixed total energy and angular momentum.In a realistic cluster, of course, these vary from one interaction to another, and a more relevant quantity is the average value of p bin , over an ensemble of triples, with a certain distribution of initial conditions.The purpose of this section is to find this probability, which we call P bin , and to characterise it.The rate calculation follows Aarseth & Heggie (1976) (see also Goodman & Hut 1993;Ivanova et al. 2005).The rate per star, at which a star of m1 becomes (on average) involved in a 3-body interaction leading to the formation of a binary with component masses m1 and m2, is given by The integration domain for R and r depends on the initial condition distribution.Here, all stars are assumed to have a Gaussian velocity distribution with temperature T , and the triple's centre-of-mass motion has already been integrated out.The Jacobi coordinates are and V and v are their respective associated velocities, where xn is the position of star n.
The rate per unit volume (instead of per star) is simply equation (17) multiplied by n1, i.e.
By the law of total probability, Γv factories into a product of the total interaction rate, Γ3UB, and the probability that a binary forms, given a three-body interaction, P bin .The former evaluates to (Atallah et al. 2024 where it is assumed that the bodies have random positions, uniformly distributed in a sphere of radius R0, and a Gaussian velocity distribution, with a 1-dimensional dispersion σ.
Then, Γv = Γ3UBP bin , where As mentioned above, the binary formation probability, p bin , given a triple interaction, is only a function of κ and R0 in the equal-mass case, but may also depend on the massratios if they are unequal.At any rate, we calculate it directly by integrating equation (1) over the allowed parameter space (for a bound remnant binary), summing over the three possibilities for star s, and normalising.
A perhaps more interesting quantity is obtained by replacing p bin by the probability for creating a hard binary, p hard , which gives the rate of formation of hard 3-body binaries, which we denote by P hard .The non-trivial shape of p hard (κ) (see figure 4) requires a more complicated calculation; both calculations, for P hard and P bin are performed in appendix C. We restrict the calculations to the ζ ≪ 1 régime, which is the physically relevant one.We now move on to describe the integration domain and derive P bin , P hard , the rate, and the orbital parameter distribution.

Occurrence of an Interaction
As alluded to above in §2, for a binary to form, one of the stars must pass sufficiently close to the centre of mass of the other two, to exchange enough energy with them to make them bound.This is not encapsulated in the κ dependence of p bin per se and must be imposed explicitly when evaluating the integrals for P bin and P hard because while the outcome does depend sensitively on the initial conditions, it does not only depend on them via κ.
We can derive a necessary condition for this to happen, based on the minimum requirement of a strong enough perturbation in the impulsive limit (cf.Aarseth & Heggie 1976).In the time-reversed configuration, where a single star ionises a binary, it is required, at a minimum, that the star would be able to give the binary enough energy to unbind it.That is, in the impulse approximation (which is a lower limit) this should be multiplied by √ 2 because both stars are affected.Hence, We therefore restrict R as follows: firstly, R ≤ R0, and additionally, where here E i s = µsV 2 /2 and E i b = µ b v 2 /2.The first case is the same as ( 24); the second case is less stringent because, for the formation of a soft binary, we require merely that one of the stars impart an impulse on the other two, which is as large as their relative momentum at (hyperbolic) pericentre.Whether subsequent interactions in the same encounter then form a binary is determined by p bin .

Computation of The Formation Probability
We are now in a position to evaluate the integrals in equation ( 22), both for P bin and for the analogous P hard .This is done in appendix C, where we also show analytically that, as functions of ζ, as ζ → 0. In particular, this implies that as R0 → ∞, both Γ and Γv tend to a constant, for hard binaries, but grow as R 3 0 if one wishes to include soft binaries, too.We show the full results of calculating P bin and P hard as functions of ζ, for the equal mass case in figure 5. We plot them as functions of a dimension-less parameter χ1 ∝ 1/ζ, which is defined in Atallah et al. (2024), for ease of comparison.In the equal mass case, χ1 = 27/(2ζ).We see good agreement between our results and those of Atallah et al. (2024).1 Surprisingly, the asymptotic χ1 dependence is not reached until very large values of χ1.The two curves, P hard,1 and P hard,2 are plotted when using the first of conditions (25) or the second for hard binaries, respectively (equations (C14) and (C18)).
It is also instructive to compare the total hard binary  that is, the contribution of ergodic hard binary formation to the total rate is roughly 70%, or Γv = 0.51 G 5 m 5 n 3 σ 9 . (29)

Orbital Parameter Probabilities
One can play the same exercise, and compute the total formation probability of binaries with given values of the semimajor axis and eccentricity, (a b , e b ), given that a binary does form.This is a convex combination of functions of the form (8), with E and J sampled from the distribution of initial conditions described above.That is, where we have made the R0 dependence of f bd explicit.We evaluate this 6-dimensional integral using Monte Carlo integration.The result is shown in figure 6.Also shown in that figure is the line e b = 1 − R0/a b , below which there is zero probability of forming binaries.This effectively forces the soft binaries, with a b ≥ R0, to be very eccentric.To compare with figure 7 of Atallah et al. (2024) we also plot the cumulative distribution P (≤ e b |a b ).This is calculated by integrating the probability density in the left panel of figure 6 over eccentricity up to e b , and normalising, at each value of a b .The result is shown in the right panel of figure 6.One can see that the two colour maps agree qualitatively, in both a b ≪ R0 and a b ≫ R0 limits, although our P (≤ e b |a b ) has a ridge which is not present in the simulation results of Atallah et al. (2024); at present, we do not know what gives rise to this ridge.

APPLICATION
To summarise, let us describe briefly how to apply the results of this paper to modelling of binary populations in a cluster with a given temperature and density: first, one should calculate ζ with equation (10).Then, compute P bin with equation (C22) and P hard,1,2 with equations (C14) and (C18).
The rate is then Γ3UB times the probability.For P hard,1 , if ζ ≪ 1 the hard binary formation rate is directly given by ( 29).However, one should bear in mind that the orbital parameters a b and e b of the binaries that form are highly correlated, and their distribution is also highly dependent on ζ.The conditional probability, P (e b |a b ) does not appear to depend strongly on ζ from figure 6, and is thus independent of the environment in the astrophysically realistic limit ζ ≪ 1.In creating a realistic distribution of three-body binaries one must, therefore, draw them from the joint distribution P (a b , e b ) in figure 6.For these one has to use equation ( 8) and the procedure at the end of §3, or equivalently read these off figure 6.

DISCUSSION AND SUMMARY
This paper consisted of two main parts: the first is concerned with the three-body problem, giving a probability distribution for the outcome of a three-body scattering experiment, for fixed values of the conserved quantities, i.e.E, J and the masses.The second part of the paper was an application of the first to a cluster -a system with a certain distribution of conserved quantities (in the equal mass case).This distribution was a Gaussian velocity distribution, while initial positions were assumed to be uniform in a sphere of radius R0.
Everything done here assumed some kind of ergodicity or phase-space mixing.The main argument made to support this was that while indeed there are no democratic resonances for an extended period of time, the problem is still chaotic, and whether a binary forms and its final characteristics, do depend very sensitively on the initial conditions because the binding energies required for a binary are much smaller than E. Alternatively, at least for the negativeenergy problem, we know that detailed balance arguments produce the same probability distribution as the one obtained from ergodic mixing in the available phase-space volume (Ginat & Perets 2021a).Here, though, we do not find the same result as Goodman & Hut (1993), which was based on detailed balance, primarily because that work did not account explicitly for angular momentum conservation, but also because in the negative-energy case, the balance was between bound triples (in a democratic resonance) and binaries and singles, while here there are no democratic resonances.
However, this can also be viewed from the opposite perspective, where the rates that we calculated here are those for three-body binary formation via an ergodic channel.They were slightly lower than those computed from numerical simulations (cf.figure 5), possibly due to non-ergodic three-body binary formation also playing some rôle, or due to inexact matching of the initial condition distribution.But the agreement, albeit imperfect -and the correct scalings with ζ -act as an a posteriori validation of the ergodicity assumption.
We found that soft binaries were formed at a much higher rate and that their eccentricity distribution tended to be highly non-thermal (this arose primarily from the restriction e b ≥ 1 − R0/a b ), while hard binaries had a thermal distribution, except at low angular momenta.When κ ≈ 0, the zero angular momentum case for the three-dimensional three-body problem was already known to be very similar to the two-dimensional problem (Valtonen & Karttunen 2006;Stone & Leigh 2019;Parischewsky et al. 2021;Samsing et al. 2022;Ginat & Perets 2023;Trani et al. 2024;Fabj & Samsing 2024), which also yields super-thermal eccentricity distributions.That soft binaries are highly preferred (which is corroborated by Goodman & Hut 1993;Atallah et al. 2024) suggests that these binaries are important in clusters, despite being mostly ignored in the literature,2 which is why the findings of this paper are not in disagreement with those of, e.g., Geller et al. (2019), who observed a sub-thermal distribution of binaries in a cluster Monte Carlo simulation but with only hard binaries tracked.Indeed, the minuscule minority of soft binaries that do survive encounters with single stars, and become hard, might not be negligible in comparison with the number of the three-body binaries that form hard at the outset.Moreover, if any gas is present, it is suggested in Rozner & Perets (2024) that dynamical friction from it could counteract the softening effect, and shield some soft binaries from destruction.
That three-body binary formation is a natural mechanism for generating eccentric wide binaries might bear on the recent conundrum of eccentric wide binaries detected in Gaia (Tokovinin 2020;Hamilton 2022;Hwang et al. 2022a,b;Hamilton & Modak 2023;Modak & Hamilton 2023), although those binaries are not globular cluster ones.
One of the key questions raised by Atallah et al. (2024) is the scaling of P hard with ζ (equivalently χ1).In a typical scattering experiment (e.g. if gravity was a short-distance force), one would expect that the rate of binary formation converges as the maximum impact parameter tends to infinity.This would require P hard ∼ ζ 5 , which is what we find if we use the first of equations ( 25) to calculate the binary formation rate (P hard,1 in figure 5).However, Atallah et al. (2024) fit a shallower, ζ 9/2 scaling to their results, which incidentally is what we find if we also use the second condition of (25) also for hard binaries (P hard,2 in figure 5); this would imply that the total binary formation rate would tend to infinity as R0 → ∞.One should observe, though, that it is unclear which of P hard,1 and P hard,2 follows the scaling of the simulation results better, and that in both cases, the asymptotic scaling is achieved only at very large values of χ1, for which there aren't any numerical data.Additionally, ζ (equivalently R0) is a physical parameter, which is not infinite in real astrophysical systems -it is very small, but finite, determined by the inter-particle separation (or by the Hill radius at the triple's position in the cluster), and the local velocity dispersion there.
For equal masses this simplifies to

APPENDIX B: SCALING OF THE HARD BINARY FORMATION PROBABILITY WITH PARAMETERS
The two parameters p hard can depend on -in addition to κ -are the ratio E/(kBT ) and the ratio R0E/(GM 2 ) (in the equal mass case).We explained already in §2.4 why so in this appendix we will argue that if R0E/(GM 2 ) ≫ 1, p hard ∝ R −2 0 .This dependence clearly comes from the entire phase-space volume, i.e. the normalisation of f in equation (1).This is readily seen by observing that for hard binaries, the integration boundaries for E b are such that R0/a b ≫ 1, so θmax(R0, E b , S) = 2π always, and nothing seemingly depends on R0.But of course, equation ( 1 16) with a direct integration of equation ( 8) for various values of the cut-off R 0 .We will also henceforth drop the i where it isn't confusing to do so.This transforms P hard into   Observe, that if we had used the second condition in equation ( 25 We also remark that this is the scaling reported by Atallah et al. ( 2024), but that this scaling would imply that the rate does not converge to a constant at arbitrarily large impact parameters.For this reason, we use P hard,1 as P hard in the rest of this paper.

C2 Soft Binaries
We repeat the same procedure as above but only change variables to Delaunay variables in the 'outer' orbit, denoted by s.This results in

Figure 1 .
Figure 1.The marginal energy distribution of the remnant binary, for equal masses.Top: R 0 = d, with d = 150GM µ/E.Bottom: the case R 0 = min {R E , d}.

Figure 2 .
Figure 2. The marginal eccentricity distribution of the remnant binary, given that the remnant binary is bound, for equal masses.This is calculated by integrating equation (8) over the domain E b < 0 (and normalising).Top: R 0 = d, with d = 150GM µ/E; middle: the case R 0 = min {R E , d}; bottom: same as the top panel, but for hard binaries only, i.e. those with |E b | ≥ E.

Figure 3 .
Figure 3.The joint energy-eccentricity distribution of the remnant binary, i.e. equation (8), given that the remnant binary is bound, for equal masses.The plots in the top row show the distribution for different values of the parameter κ = J 2 |E| G 2 M 2 µ 3 .The plots in the bottom row are the probability distributions for the dimension-less angular momentum s ≡ 1 − e 2 b , given z.The black line is the line e b = 1 − R 0 /a b -values of s above that threshold cannot form.

Figure 4 .
Figure 4.The binary formation probability, as a function of the angular momentum, made dimensionless, for equal masses.Top: the case R 0 = d, with d = 150 GM µ/E; bottom: the case R 0 = R E .As hard binaries require significant energy exchanges to form, essentially all hard binaries form for R 0 ≤ R E , so that their relative fraction is lower when R 0 can be larger -because only soft binaries form at R E ≤ R 0 ≤ d.

Figure 5 .
Figure 5.The binary formation probability, as a function of χ 1 = 27/(2ζ), for equal masses.This assumes a density n = 10 5 pc −3 , m = 1 M ⊙ , and a velocity dispersion σ = 10km s −1 .Data from the numerical experiments of Atallah et al. (2024, figure 3) are over-plotted for comparison., and re-scaled according to their equation (20).Details are in the text and appendix C.
formation rate to Γv, to the results of previous works: Goodman & Hut (1993) (cf.Heggie & Hut 2003) obtain, for the equal

Figure 6 .
Figure 6.Left: The joint probability distribution of the semi-major axis and the eccentricity of three-body binaries in a cluster.Right: the cumulative eccentricity distribution, given a semi-major axis.These are evaluated for equal masses at ζ = 0.09 (top) and ζ = 9 × 10 −3 (bottom).The black line is e b = 1 − R 0 /a b .While the joint distributions depend on ζ, the conditional eccentricity distribution does not seem to do so, in the ζ ≪ 1 limit.See text for more details.

Figure B2 .
Figure B2.A comparison of the scaling in equation (16) with a direct integration of equation (8) for various values of the cut-off R 0 .

Figure B3 .
Figure B3.A plot of the integral of p hard over κ, for R 0 E = 150GM µ.See text for details.

∝
is the equation that we integrate numerically for the plot in figure5, specifically, the line denoted by P hard,1 there.Let us now show that this yields a ζ 5 ∝ χ ζ 6 .In the limit of small ζ, i.e. large R0, the only other dependence on R0 in equation (C14) can come from θ b ap θ s ap .These are obtained from the constraints in equation (25), which imply that in fact θ s ap is independent of R0, while the argument of θ b ap , (recall that this is the Delaunay angle of a hyperbolic trajectory, which is unbounded).Together, we find P hard,1 ∼ R ), instead of the first, we would have had θ s ap ∼ √ r, which is dominated by the large r régime, i.e. θ s ap ∼ √ R0.This extra power of 1/2 would have resulted in an overall P hard ∝ R −9/2 0 .This is what is plotted in figure 5 as P hard,2 , viz.