Demographics of three-body binary black holes in star clusters: implications for gravitational waves

To explain both the dynamics of a globular cluster and its production of gravitational waves from coalescing binary black holes, it is necessary to understand its population of dynamically-formed (or, `three-body') binaries. We provide a theoretical understanding of this population, benchmarked by direct $N$-body models. We find that $N$-body models of clusters on average have only one three-body binary at any given time. This is different from theoretical expectations and models of binary populations, which predict a larger number of binaries ($\sim 5$), especially for low-$N$ clusters ($\sim 100$), or in the case of two-mass models, low number of black holes. We argue that the presence of multiple binaries is suppressed by a high rate of binary-binary interactions, which efficiently ionise one of the binaries involved. These also lead to triple formation and potentially gravitational wave (GW) captures, which may provide an explanation for the recently reported high efficiency of in-cluster mergers in models of low-mass clusters ($\lesssim 10^5\,{\rm M}_\odot)$.


INTRODUCTION
The first detection of a gravitational wave (GW) signal (Abbott et al. 2016) was one of the key milestones for the advent of the multimessenger era of astrophysics.This has allowed the direct study of, among others, dark systems such as black hole (BH) binaries.To date, 90 GW signals from mergers of binary compact objects have been detected (Abbott et al. 2021;The LIGO Scientific Collaboration et al. 2021) -of which the vast majority are binary black holes (BBH) -which requires us to tackle the question of the origin of these sources.
Dynamical interactions in the dense cores of stellar clusters have been put forward (Portegies Zwart & McMillan 2000) as one of the most likely formation channels for BBH mergers (Rodriguez et al. 2022;Kremer et al. 2020;Chattopadhyay et al. 2022;Antonini & Gieles 2020;Kumamoto et al. 2020;Banerjee 2021a;Mapelli et al. 2022).Current GW constraints are compatible with a majority of massive mergers having formed dynamically (Antonini et al. 2023).The predictions for the production of BBH mergers in this channel hinge on the population of BBHs, which points towards the importance of fully characterising the demographics of these systems.Because the majority of massive stars form in binaries or higher order multiples (Sana et al. 2012;Moe & Di Stefano 2017), a large fraction of BHs end up with a companion after stellar evolution (hereafter 'primordial BBHs').It has been suggested, however, that in massive clusters the binaries that result in BBH mergers are dominated by the three-body 1 BBHs (Chattopadhyay et al. 2022;Torniamenti et al. ★ E-mail: danielmarin@icc.ub.edu 1 We note that Tanikawa, Hut & Makino (2012) showed that their formation tends to involve more than three particles, but for historic reasons we prefer 2022).A possible explanation of this is due to primordial binaries typically having a smaller semi-major axis, and thus a smaller interaction cross-section, than three-body binaries (Barber et al. 2023).It is, therefore, important to study the demographics of three-body binaries to understand dynamical mergers.Goodman (1984) showed that in steady post-collapse evolution, the number of three-body binaries,  b , depends on the total number of stars in the cluster, , as  b ∝  −1/3 .These models predict that a handful of binaries is expected in single-mass clusters of a few hundred stars.However, -body models of such clusters find instead that there is typically only one binary present (Giersz & Heggie 1994b).The aim of this paper is to understand the demographics of the population of three-body binaries in single-mass and two-mass star clusters, where the latter are intended to understand the behaviour of clusters with stellar-mass BHs.We do this by revisiting the model for steady post-collapse evolution and by comparing this to results from direct -body calculations.
The paper is organised as follows: in Section 2 we revisit the theory for the population of three-body binaries.In Section 3 we present the -body simulations, of single-mass and two-mass models.In Section 4 we compare the predictions and the -body models.The cause of the aforementioned discrepancy is discussed in Section 5.In Section 6 we discuss its implications for GW observations.to refer to dynamically-formed binaries from single bodies as three-body binaries In this section we revisit the models for steady post-collapse evolution of both single-mass (Goodman 1984) and two-mass (Breen & Heggie 2013) clusters.
Throughout this paper we consider star clusters in the post-corecollapse evolutionary phase, when binary formation and hardening is ongoing and there is a balance between the energy production by binaries in the core and energy transport by two-body relaxation (Hénon 1975).Furthermore, we assume that this post-collapse evolution can be approximated by a steady evolution of the cluster parameters, i.e. we ignore gravothermal oscillations that occur in single-mass models with  ≳ 8000 stars (Bettwieser & Sugimoto 1984;Breeden et al. 1994) and in two-mass models where the number of BHs is  BH ≳ 2000 (Breen & Heggie 2012).Both of these assumptions allow us to express relations between several cluster properties, which are discussed below and summarised in Table 1 and Table 2.

Single-mass clusters
In post-collapse, we expect the cluster to be nearly isothermal within the half-mass radius  h , so that the one-dimensional velocity dispersion is similar within the core and within  h .For a cluster with mass  and energy  in virial equilibrium with a virial radius  v = −  2 /(4) and positions and velocities sampled from a Plummer (1911) model, we have  h ≃ 0.8 v and the one-dimensional velocity dispersion, , then obeys where  is the gravitational constant,  is the mass of the stars and ⟨ 2 ⟩ is their mean-square velocity.The central dispersion can be expressed in terms of  by defining a numerical coefficient  1 Next, we define the relation between  h and the core radius,  c .The corresponding equation can be obtained for the steady post collapse evolution via Hénon's principle (Hénon 1975) by noting that the energy flow through  h needs to be provided by the core.The exact form of the  h / c ratio requires one to specify the energy-generating mechanism.This derivation has already been carried out in Heggie & Hut (2003, p. 265) and Breen & Heggie (2013) (although the former did not keep the numerical pre-factors) when the bulk of the energy is set to come from binary hardening via three-body processes.By assuming that the mean mass, the Coulomb logarithm (ln Λ ≃ ln(0.11),Giersz & Heggie 1994a) in the denominator of the central relaxation timescale, and the dimensionless ratio of the central potential (| 0 |) to the central velocity dispersion, | 0 |/ 2 c , all are constant, they obtain where  2 is a constant.If the assumptions above are relaxed, the final ratio would include a dependence on (| 0 |/ 2 c )/ln Λ but, since both these quantities scale weakly with , we are justified in taking  2 as a constant.In Section 4 we determine  2 and a possible additional -dependence from -body simulations.
We complement the above relations with the relation for the number of stars in the cluster core,  c .Starting from the definition of the core radius (Heggie & Hut 2003, p. 71 We then use the fact that, for the isothermal model and King (1966) models with high concentration, the central density can be expressed in the average density within the core as  0 ≃ 1.9 c , allowing us to obtain a measure of the mass contained within the core radius.Then, the number of components (singles + binaries) within the core,  c , is 2 The  i parameters in this theoretical prediction are of order unity and their values and possible -dependences will be determined in Section 4 by fits to -body models.
The binaries that we are concerned with are those that do not break up after an encounter with a typical single star of the cluster, meaning that their binding energy  =  2 /(2) is greater than  h ≃  2 c (where  is the semimajor axis of the binary).We refer to these binaries with the usual nomenclature of hard binaries, as opposed to soft binaries which have  <  h , and  >  h =  2 /(2 h ).Heggie (1975) showed that hard binaries become on average more bound (that is, they harden) after an interaction with a third star.The fractional increase of  is on average Δ ≃ 0.4 for resonant encounters 3 (Spitzer 1987) and the released energy results in a velocity kick of the single star and the centre-of-mass of the binary.Since each encounter gives the binary a momentum kick that scales with , at some moment  is sufficiently large such that the subsequent velocity kick of the binary is above the escape velocity from the centre of the cluster.From conservation of energy and linear momentum and assuming equal masses, it can be shown that the upper limit of the binding energy of three-body binaries is 1984), which corresponds to a semimajor axis of  ej =  2 /(2 ej ) .In order to obtain the total number of binaries, we will compute the distribution of (hard) binaries per unit volume and per unit  ≡ /( 2 c ) and integrate it from This approach results in a scaling law for the number of hard binaries (Goodman 1984) that we will reproduce with a detailed analysis of the numerical pre-factors.
The first model for the energy distribution of three-body binaries in post-collapse evolution was presented by Retterer (1980) and was later refined by Goodman & Hut (1993).We will re-derive the binary distribution using the Goodman & Hut pre-factors, because they are based on a large set of scattering experiments by Heggie & Hut (1993).We begin by considering the net creation rate of binaries per unit volume (Goodman & Hut 1993) and the hardening rate (Heggie & Hut 1993) defined as the rate of energy generation of the binaries Since this hardening rate is independent of  itself, the steady postcollapse distribution of hard binaries is uniform in  between the minimum and maximum binding energy.The distribution per unit volume,  (), is Assuming that  ej ≫  h and that the relevant binary physics happens within the volume  c of the cluster core, the total number of binaries is 5 The above equation can be simplified by using the scaling laws in Table 1 to obtain 6 The ratio | 0 | / 2 c is dimensionless and thus only depends on .In Section 4.1 we quantify this dependence and show that it is very weak.By neglecting it, Goodman (1984) found that the scaling of the number of binaries with the number of stars in the cluster could be approximated as  b ∝  −1/3 .

Two-mass clusters
In the post-collapse evolution of a cluster, the most massive stars have already ended their evolution, so the stellar population can be approximated by two distinct mass bins: that of heavy stellar-mass BHs with masses of  BH ≃ 20 M ⊙ and that of light stars with  ★ ≃ 0.5 M ⊙ .Due to the large difference between these bins, we will assume that all BHs have the same mass  BH and all stars have 5 Compared to Goodman (1984), we use an updated value of the numerical prefactor, so their prediction for the number of binaries is roughly three times larger than ours. 6So far we have assumed that  ≃  s +  b , with  s the number of single stars, and that the mean mass within the core,  c , is equal to the mean mass of the cluster, that is,  ≃  c .If we explicitly keep the / c dependence in the above equations, we obtain that the result for  h / c (equation 3) and  b (equation 9) need to be multiplied by / c and (/ c ) 3 , respectively.The result for  c (equation 5) would remain unchanged.the same mass  ★ ≪  BH .Each of these mass bins contribute a total of  BH and  ★ to the total cluster mass, respectively.We will continue referring to the heavy and light components as BHs and stars, respectively, but we note that we are considering Newtonian gravity in the point-mass limit.
To describe the two-mass cluster, we need two extra parameters with respect to the single-mass case: the mass ratio,  ≡  BH /, which sets what fraction of the cluster mass is in BHs, and the stellar mass ratio,  ≡ ⟨⟩ / BH , which sets how massive the BHs are with respect to the average mass of stars and BHs, ⟨⟩, in the cluster.Both of these quantities are defined such that ,  ∈ (0, 1).Some authors (for example, Breen & Heggie 2013) give alternative definitions of these parameters with respect to the mass of stars,  ★ and  ★ , instead of total and average mass.Both of these definitions are interchangeable, and converge in the limit of a small number of BHs.We prefer the definition given above as it simplifies the equations that follow.
The clusters of interest will be those in which the light stars form the bulk of the cluster mass 7 ,  ★ ≫  BH .This allows us to compare our results to the two-mass model of Breen & Heggie (2013).Their models assume that equipartition between stars and BHs can not be achieved, which applies if (Spitzer 1969) After core collapse, the BHs in a two-mass cluster behave as a denser BH subsystem embedded in a larger star cluster.The central region of the cluster will then be populated mainly by BHs, whose evolution is powered by BBHs.This energy is then transferred outwards in the BH subsystem until it is deposited in the cluster of stars surrounding the BHs.In order to obtain the equations relating the cluster properties, we will make use of the model of Breen & Heggie (2013).These scaling laws, which are the analogues of equations (1-5) for two-mass clusters, are summarised in Table 2.The results are equivalent to what we would have found if we considered the BH sub-cluster in isolation, with the only key difference being the central potential term, which here includes the contribution of both the stars and the BHs.We reproduce the equations here, for completeness 7 Although this is a standard assumption, there is a sparkling interest in low-density clusters near dissolution where most of the mass is in BHs (see for example Banerjee & Kroupa 2011 and the discussion about Palomar 5 in Gieles et al. 2021).Nevertheless, we will not focus on such systems.
As in Section 2.1, the number of binaries can be computed from the integral (equation 9) of their binding energy distribution.To compute the number  BBH of BBHs, we will use the form in equation ( 8) but changing the variables to their counterparts in the BH sub-cluster, e.g. ↦ →  BH .This recovers the same form of the two-mass binary binding energy distribution from Retterer (1980), i.e.
BBH ≃ 3.0 which we can express as In parallel to the result for the single-mass case, we obtain the prediction that the number of BBHs scales as  BBH ∝  −1/3 BH , although in this case the contribution of the | 0 | / 2 c,BH term is larger.In the following section we will test the above equations against -body models.

DESCRIPTION OF THE 𝑁-BODY MODELS
We simulate a set of single-mass and two-mass clusters and compare their global properties and the number of three-body binaries to the theory presented in the previous section.To isolate the dynamics of binary formation and disruption, we neglect primordial binaries, stellar and binary star evolution, post-Newtonian physics, and the Galactic tidal field.We sample initial positions and velocities from a Plummer (1911) model using the McLuster implementation from Küpper et al. (2011).The exact choice of initial positions and velocities is not important for this study, as we will only consider the evolution beyond core collapse, after which the information about the initial conditions is erased (Lynden-Bell & Eggleton 1980).Our simulations are run in Hénon units, i.e.  =  0 = −4 0 = 1 (Hénon 1971), where  0 is the initial cluster mass and  0 is the initial energy of the cluster.For our assumption of virial equilibrium, the initial virial radius is We use the direct -body code PeTar (Wang et al. 2020b), which is a high-performance code used for the modelling of large- collisional systems.It is based on a hybrid method: a Barnes-Hut (Barnes & Hut 1986) tree method for long-range forces and a few-body integrator for short-range interactions.The integrator for close interactions, sdar (Wang et al. 2020a), combines an explicit Hermite integrator for weakly perturbed binaries and a slow-down method for more compact subgroups.
The centre of the cluster as well as  c are defined following Casertano & Hut (1985); McMillan et al. (1990), with  h determined from the bound stars in a reference frame where the density centre of the cluster is in the origin.We consider binaries as pairs of close stars whose binding energy is greater than  2 c -which is proportional to the average kinetic energy of single stars in the core -and that are bound to the cluster.

Single-mass models
We have run a family of single-mass models with logarithmicallyspaced , in the range 128 ≤  ≤ 8192, where each model is run several times for statistical significance.The input  and number of runs are shown in Table 3. Within the aforementioned simplifications, single-mass models are scale-free, and thus  is the only parameter that defines the post-collapse evolution of the cluster.The simulations are run through core collapse, which happens at  cc ≃ 16 rh,0 for Plummer models (Cohn 1980), with  rh,0 being the initial relaxation timescale (Spitzer & Hart 1971), evaluated at the start of the simulation, where  = 1 and ln Λ is the Coulomb logarithm with Λ ≃ 0.11 for single-mass models (Giersz & Heggie 1994a).The data for the post-collapse evolution is then taken starting from  0 = 1.1 cc up until  f =  0 + 20 rh,0 , at constant intervals of 4 -body times.Values found for the cluster properties are averaged over all snapshots between the times  0 and  f , with their corresponding error bars being the standard deviation among different runs.The fittings are done with linear regressions, where we take the uncertainties into account by weighting the data points with the inverse squared value of their errors.

Two-mass models
For the two-mass -body models, we take our data from  0 = 1.1 cc to  f =  0 + 2 rh,0 .We use a lower number of  rh,0 than in the singlemass models because in the two-mass models the binary phenomena happen within the BH subcluster and thus the relevant timescale is shorter.The impact of the adopted simulation time is discussed in Section 6.1.For the core-collapse timescale of two-mass clusters, we use the result of Kim & Lee (1997) as well as their expression for  rh,0 (equation ( 18), but setting Λ = 0.4).The values of the three parameters  BH ,  and  for the set of simulations that we have run is summarised in Table 4.We set  ∈ [5 × 10 4 , 1 × 10 5 , 2 × 10 5 ] and store the snapshots every 8 -body times.For the mass ratio  we adopt  ∈ [1/50, 1/25], to approximate metal-poor (log 10 (/ ⊙ ) ≃ −1.5) and metal-rich (log 10 (/ ⊙ ) ≃ −0.5) clusters, respectively.For the mass fraction , we have chosen  ∈ [0.025, 0.05, 0.1].Note that, for all combinations of these parameters, our clusters avoid gravothermal oscillations (Breen & Heggie 2012) and satisfy the above criteria for  1) and equation ( 2)) as a function of  for the single-mass  -body models (crosses with error bars).Fits for a constant  1 and an  -dependent  ′ 1 are shown with grey and red lines, respectively.
Spitzer instability (Spitzer 1969), such that the evolution should be well-described by the theory of Breen & Heggie (2013).

Single-mass clusters
In this section we calibrate and compare the theoretical predictions of the single-mass model of Section 2.1 to the -body models of Section 3.1.We also allow the theory to include slight additional -dependences in the  1 ,  2 parameters (that we label as  ′ 1 and  ′ 2 , respectively).In Fig. 1, we show the ratio  1 =  2 c / 2 as a function of .We find that a constant  1 (equation 2) is a good approximation only for  ≲ 10 3 , and  ′ 1 increases slightly at larger Eq.20 (theory) Eq. 21 (theory + corr.)N -body data Using  i independent of  (gray) we are able to roughly match the  -body models, whereas using the  -dependant  ′ i coefficients (red) yields a very accurate match.
.A scaling  ′ 1 = 1.1(/10 2 ) 0.06 describes the -body data well.This increase is likely due to the increase of  h / c with , leading to a higher central density and dispersion, relative to the global values, for clusters with larger .
The  h / c ratio is shown Fig. 2. The scaling of equation ( 3) is a good approximation for constant  2 for  ≃ 10 3 − 10 4 , whereas it deviates slightly at lower values of , which, as discussed in Section 2.1, may be due to the weak -dependence of the Coulomb logarithm in the denominator of the central relaxation timescale and in the | 0 | / 2 c term.A scaling of  ′ 2 = 0.19(/10 2 ) −0.09 provides a good fit to the -body data.
From these scaling laws, we are able to accurately predict the number of stars in the core (Fig. 3).The prediction for  c with constant  i (equation 5) has the same form than the one given in Goodman (1984), which yields a rough estimate of  c from the -body models within <50% over two orders of magnitude in .Using the -dependent expressions for  ′ 1 and  ′ 2 yields which accurately describes the -body data.Furthermore, in Fig. 4 we show the ratio | 0 |/ 2 c as a function of , which can be approximated by the power-law relation Now all ingredients of the expression for the number of binaries have been discussed, it is time to combine them and compare to  b we find in the simulations.We obtain Which can be corrected by the deviations from the scaling laws above to yield such that a cluster with ∼ 100 stars is expected to have ∼ 3 binaries.
In Fig. 5 we show both the prediction for  b and the -body results.The prediction (equation 24) approximately reproduces the number of binaries in the core at large  ( ≳ 10 3 ), whereas in the smallest clusters  b is over-predicted by an order of magnitude.It could be reasonable to consider binaries outside of the core, as the interactions with single stars increase the apocentre of their orbits as they evolve in their lifecycle.When we consider all binaries within  h , the number of binaries in the -body models is higher, and increases slightly with , that is, the opposite -dependence.The simulations are compatible with the clusters only having a single central binary at any point in time, which spends a fraction of its lifetime outside the core.

Two-mass clusters
In this section we compare the theoretical prediction of the twomass model of Section 2.2 to the -body models of Section 3.2.The velocity dispersion within the BH sub-cluster's core is shown in Fig. 6, where we find reasonable agreement for a constant  1,BH (equation 12), but a slight decline of  1,BH with  BH is preferred, equivalent to the single-mass case.For the ratio  h,BH / c,BH in the BH sub-cluster (Fig. 7), we find that the scaling law of equation ( 13) is a valid approximation to the data.Similarly to the single-mass case, the deviation at low  BH can be attributed to the weak  BH dependency of the Coulomb logarithm.However, in the two-mass case the deviation is larger due to the smaller number of BHs in the core.
From these relations, we can predict the number of BHs in the sub-cluster core (Fig. 8).As in the single-mass case, constant  i,BH (equation 15) give a rough approximation to  c,BH , with the discrepancy being largest at lower  BH .Using the corrected which accurately reproduces the results from the -body models for  BH ≳ 50, although it slightly overpredicts the results for the sub-clusters with the fewest black holes in the core ( c,BH ≲ 7).
The above scaling laws are complemented by a fit to | 0 | / 2 c,BH (Fig. 9).From the three parameters of the model, we find that this value only correlates strongly with .The best-fit power law is Although this equation may be taken to imply a very deep central potential, this is not the case.The apparently large value is an artefact of expressing the potential in units of the central velocity dispersion Using the calibrations for the scaling laws (equations 12-15), we can evaluate the prediction for the number of BBHs of equation ( 17 This can be corrected with the deviations from the scaling laws as Although the two-mass steady post-collapse model is able to describe the other properties of the -body model, it over-predicts the number of BBHs (Fig. 10).This result is very similar to the singlemass case, although in this case the discrepancy is larger.In general, the discrepancy is most notable in models with smaller  and smaller  BH .

THE CAUSE OF THE DEARTH OF BINARIES: BINARY-BINARY INTERACTIONS
In the previous sections we find that the theoretical predictions for the population of binaries are unable to describe the number of threebody binaries in the -body runs, even if all the other elements of the model do match the simulations.This discrepancy had already been observed in Giersz & Heggie (1994b), but no cause was definitively identified.In this section we search for an explanation for the discrepancy between the theory and the -body models.Because binary formation is well understood (Heggie 1975), we should search for a missing ingredient in the model for binary evolution.The models of Retterer (1980) and Goodman & Hut (1993)  binaries and single stars, and binary-binary interactions are ignored with the argument that they are rare.However, they start to play a role if the binary fraction in the core is above some critical value.The most likely outcome of a binary-binary interaction is that at least one of the binaries is ionised (∼ 95% of the outcomes for the scattering experiments in Antognini & Thompson 2016), which can explain the dearth of binaries.Because the total number of binaries in the -body models is O (1), ignoring interactions among binaries is a justified assumption for clusters with  c ≫ 1.However, as we have shown in Fig. 3,  c ≲ 20 for  ≲ 10 3 , whilst the theory predicts  b ≳ 2. In fact, the predicted binary fraction in the core is  b =  b / c ≃ 0.6(/10 2 ) −0.9 .Therefore, if the model for the binary energy distribution in the steady post-collapse model is correct, then binary-binary interactions would be important, so we need to understand their effect on shaping the energy distribution of binaries.We start by estimating the critical binary fraction in the core, above which binary-binary interactions become more frequent than binarysingle interactions.The interaction cross-section for a binary with a single goes as Σ bs ∝  p  tot / 2 bs , where  bs is the relative velocity,  p ≃  is the relevant minimum distance (with  the semimajor axis of the binary), and  tot = 3 is the total mass of the interacting binary and single.Assuming equipartition, the typical relative velocity among binaries,  bb , is a factor of √︁ 3/4 lower than among binaries and singles ( bb = √︁ 3/4 bs ).Then, taking  tot = 4 and  p = 2 for binary-binary interactions, we estimate that the cross section for these is Σ bb ≃ 3.6Σ bs .The ratio of the rate of binary-binary interactions to the rate of binary-single interactions, Γ bb /Γ bs , is This estimate shows that binary-binary interactions are expected to be more frequent than binary-single interactions for8  b ≳ 0.3 in the core.Although this may seem high, the core contains of order 10 stars, so this criterion is already met if 3 binaries are created.This is approximately the number of hard binaries predicted by Goodman & Hut (1993) for clusters  ≃ 100, and we note that their models predict an even higher formation rate of soft binaries which have a larger cross section for interactions and can therefore also contribute to binary ionisation in a way that is not captured by the Goodman & Hut (1993) model.We will therefore test whether binary-binary interactions can reduce the number of binaries seen in the -body models.

Single-mass models
As stated above, the most likely outcome of a binary-binary interaction is that at least one of the binaries is ionised.If the timescale for this ionisation mechanism is shorter than the time needed to form new binaries, then the predicted uniform binding energy distribution (equation 8) is never realised.Instead, one would find a decline in binaries at higher , because they are destroyed at some point during their life-cycle.Furthermore, this deviation should be larger at smaller , where the theory predicts the highest binary fraction.In Fig. 11 we show the histogram of binary binding energies within the core during the steady post-collapse evolution of three cluster models.As predicted, we observe that the predicted uniform distribution is not realised, supporting the idea that binaries are destroyed before they reach the maximum binding energy.
The qualitative difference of binary-binary interactions with respect to binary-single interactions is the possibility of formation of stable triple systems (Zevin et al. 2019).Although the formation of stable triples via binary-single interactions is not energetically forbidden, the probability of such process is zero (Heggie & Hut 2003, p. 211).Thus, the existence of dynamically-formed stable triple systems is a necessary (but not sufficient, Section 6.2) condition to confirm the importance of binary-binary interactions in a cluster.The formation of triples can easily be measured in an -body model, where we identify triples as bound states of three stars that verify the stability criterion from Mardling & Aarseth (2001).The results for the average number of stable, dynamically-formed triples  t in the single-mass model is shown in Fig. 12, where we see that not only such triples are present, but their number is an inverse function of the number of stars of the cluster.
In order to gauge the impact of binary-binary interactions in the binary population, we will estimate the relative importance of the triple formation rate -assuming the theoretically expected binary distribution function in equation ( 8) -with respect to the net binary formation rate Γ b (equation 6).For triple formation in binary-binary interactions to be negligible, this ratio should be much smaller than one.We can estimate the triple formation rate with the corresponding cross sections.We will use the results of Antognini & Thompson (2016), where the authors considered scattering experiments of equal-mass,  circular binaries in the Newtonian point-particle limit.Their quoted result for the cross section of stable triple formation is The cross section depends on the semimajor axes  1 ,  2 of the incoming binaries and the ratio v of their relative velocity to their critical velocity, v = / √︁ ( 1 +  2 )/( 1  2 ).The total rate of 10 2 10 3 10 4 Avg. number of components (N ) Triple vs. binary formation rates (Γ t,bb /Γ b ) Figure 13.Ratio of the triple formation rate Γ t,bb to the binary formation rate Γ b .This ratio is ≫ 1 for  ≲ 10 3 , which confirms the importance of binarybinary interactions in that  regime.The ratio scales as Γ t,bb /Γ b ∝  −0.8 .triple formation in binary-binary encounters per unit volume can be obtained in the 'n-sigma-v' formalism by the integration of these variables using, respectively, the distribution of  (equation 8) and the distribution for the relative velocities.If we assume that the distribution of velocities in the core is Maxwellian with a dispersion of  c , then the dispersion of relative velocities among singles is √ 2 c .If we then take binaries to be twice as massive and assume equipartition between singles and binaries, we obtain that the relative velocities of binaries also follow a Maxwellian distribution with a dispersion of  c (  Maxwell () ∝  2 / 3 c exp(−0.5 2 / 2 c )).Therefore, the triple formation rate is The ratio Γ t,bb /Γ b is then expressed only as a function of  by using the scaling relations in Table 1 and shown in Fig. 13.The ratio is greater than one for  ≲ 2000 and declining approximately as  −0.8 .This supports the idea that binary-binary interactions are destroying three-body binaries faster than they form, and that this effect is especially relevant in lower-mass clusters.In turn, this gives a theoretical explanation for the observed presence of a single threebody binary in -body simulations, as every time a new binary forms it is rapidly destroyed in a binary-binary interaction.

Two-mass models
In section 5.1 we showed that, for single-mass models, binary-binary interactions are non-negligible even when no primordial binaries are present.We will now argue that this is a general prediction that also applies to the population of BBHs in two-mass models.
As before, we show in Fig. 15 that a uniform BBH binding energy distribution is not realised, with the discrepancy being greatest at lower  BH and higher .Contrary to the single-mass case, here the models with lower  BH have higher  ej because of their smaller  (see equation 29).Furthermore, we show that BBH-BBH interactions happen by looking at the formation of stable triple BH systems.Indeed, in Fig. 14 we can see that stable BH triple formation happens for all clusters, independently of their number size, BH mass fraction or mass ratio.In the two-mass case, the values of  BH correspond to the smaller end of the range of  in the single-mass case, so we do not observe the decrease in the number of triples that is seen at large  in the single-mass models.The presence of stable triples in -body models without primordial binaries had already been noticed before, e.g.Aarseth (2012); Banerjee (2018).Although stable with respect to their internal kinematics, the triples that form in our models are rather soft and can easily be destroyed via interactions with unbound singles.
In parallel to our argument in the previous section, we evaluate the rate of BBH-BBH interactions.We use the same form than in equation ( 33), although we will now use the scaling relations in equations (12-15) to obtain the ratio of rates as a function of  BH and .This ratio of rates is shown in Fig. 16.As can be seen, the ratio is much larger than one for all values of the number of BHs, mass ratio or mass fraction, which explains the discrepancy between the observed and predicted number of binaries in Fig. 10.We can therefore conclude that BBH-BBH interactions are a key ingredient to shape the distribution of three-body BBHs.In the next section we discuss possible caveats in our arguments and alternative interpretations.

Validity of the assumptions
In order for the binaries to populate the binary binding energy distribution in the absence of binary-binary interactions (Fig. 11 and Fig. 15), the properties of the cluster must evolve slowly with time.The theory (dotted lines) does not match the observed distribution (full lines), with the discrepancy being highest at lower  BH ,  and higher .
The criterion for this is that the timescale at which the cluster expands,  exp , is larger than the binary lifecycle,  bin ≃ ( ej −  h )/  ≃  ej / .We can define the timescale for expansion similarly to  bin by assuming that  exp ≃  h /  h .Using Hénon's principle (Hénon 1975), we can relate the energy production of the core to its global properties as with  ≃ 0.1 (Hénon 1965;Alexander & Gieles 2012).If we assume that the cluster energy scales as  ∝   2 / h and that the mass loss is negligible, we obtain This can be compared to the binary lifecycle, which we find from the scaling laws of Section 4.1 We also give the value of this timescale in terms of the dynamical time  dyn by using that  rh / dyn ≃ 0.1/( log Λ) (Heggie & Hut 2003).
The binary timescale is approximately a constant number of -body times, weakly decreasing with .This somewhat counter-intuitive result is a consequence of the larger central concentration of high- models, leading to a faster evolution of the binaries.Alternatively, if we had considered homologous evolution ( c ∝  h and  c ∝ ), this timescale would increase strongly with ,  bin ∝  2  dyn .
For the single-mass case, we find  exp >  bin for  ≳ 800, so we are justified in taking the slow evolution approximation in larger clusters.For smaller clusters, the expansion leads to a decrease in the binding energy distribution, as (equation 8) This implies that, for a fixed , the expansion of the cluster would lead to a lowering of the binding energy distribution with time, so we would find a pile-up of binaries at large , because these were formed earlier when  () was higher.We actually see a lack of binaries at large , making the absence of binaries even more important for  ≲ 800.
Another timescale to consider is the duration of our simulations,  f −  0 .If this timescale is much shorter than the binary lifetime  bin , the binaries would not have enough time to populate the binding energy distribution, even in the absence of binary-binary interactions.The duration of the simulations is  f −  0 = 20 rh,0 .Because the cluster expands, the actual number of elapsed relaxation times is less than 20.We can estimate the instantaneous relaxation timescale by assuming that (Hénon 1965)  rh () =  rh,0  <  cc (/ cc ) rh,0  >  cc (38) So we obtain which we can compare to equation ( 36) to show that  f −  0 >  bin for  ≳ 380, so this condition is fulfilled in all models but the least massive ones.A similar calculation for the two-mass models yields (ignoring the variation of the Coloumb logarithms) 40) which is larger than unity for models with  BH ≳ 150 (450) for  = 1/25 (1/50) and  = 0.025.This shows that our two-mass models may not have been evolved for long enough to be comfortably in the regime where the binding energy distribution is well populated and so may explain why we do not see an obvious decrease in the number of stable triples with  BH in Fig. 14 as in the single-mass case.To check whether there is evolution of the binary energy distribution, we compared  () in the first half to  () in the second half of the simulation interval.The two distributions are statistically indistinguishable for all  BH , suggesting that the simulation time is long enough for  () to have reached its steady post-collapse shape.
Yet another possible caveat is the convergence of the data points to the true underlying distribution.For low- models, we have many runs and therefore the convergence is attained, but for the most massive models, we have only a handful of runs (see Table 3).In this case, the binary lifetime is longer than the -body simulation timescale and thus convergence is attained by combining multiple snapshots from the few runs.

Alternative pathways towards stable triple formation:
'binary-single-single triple formation' Binary-binary interactions are not the only mechanism towards stable triple formation.Alternatively, one could consider triple formation in a binary-single-single interaction, in the same way as three-body binary formation, but with one component being a pre-existing binary.To estimate the rate of hard binary formation in three-body encounters, Heggie & Hut (2003) consider the rate at which two unbound stars approach closer than  h and multiply that rate by the probability that a third star is present within the volume, which is roughly  c  3 h .We can estimate the rate of triple formation due to binary-single-single interactions, Γ t,bss , by a similar argument, although performing the change  c ↦ →  b and  h ↦ → 3 h (so that the outer orbit of the triple is about ∼ 3 h ).Then, the ratio of the rate of triple production over binary formation is roughly , where we use the scaling  c () from equation (20) and, as seen in the body models, we assume that there is a single binary in the core.Therefore, the production of triples due to binary-single-single interactions could be of similar importance as binary-binary interactions.This applies to soft triples with large outer semimajor axis of 3 h .The formation rate of 'hard' triples is a factor of 27 lower.The dependence in Γ t,bb /Γ b is steeper ( −0.8 , see Figs. 13 and 16) so we conclude that binary-single-single triple formation may occur, but that triple formation via binary-binary encounters dominates for the smallest- clusters.In addition, we showed in Fig. 16 that triple formation in binary-binary encounters is more important in two-mass models, whereas the above scaling for binary-single-single triple formation should be the same for two-mass models, but with  replaced by  BH .

The effect of galactic tides
So far we have neglected the effect of the galactic tidal field.Here, we discuss how our results would be affected by the inclusion of galactic tides.We start by considering the flow of energy through the cluster's half-mass radius, which is governed by the coefficient  (equation 35).It is shown in Gieles et al. (2011) that the value of  in tidally-limited clusters and isolated clusters is equal to within 20%, and that the first half of evolution is similar in tidally-limited and isolated clusters.Furthermore, the authors estimate that two thirds of Milky Way globular clusters are in this first expansion phase.Per Hénon's principle, this energy flow is balanced by the production of energy in the core, and therefore the density profile within a few  h and the overall results of this work are independent of the underlying host galaxy potential.The only difference that the tidal stripping introduces is the possibility of evaporation pathways where the mass loss rate of stars is higher than the mass loss rate of BHs, such that the cluster evolves to a 100% BH cluster ( → 1, see Banerjee & Kroupa 2011;Gieles et al. 2021).As stated in Section 2.2, we do not consider high- clusters, although we roughly expect them to behave like the single-mass models in the limit  → 1.Nevertheless, observations point towards most clusters in the Milky Way having a small BH mass fraction of less than 1% (Dickson et al. 2023).

Implications for GW astronomy
In this work, we did not include post-Newtonian terms; however, here we discuss how the high rate of binary-binary encounters (Section 5) and the resulting stable triples could lead to GW inspirals by captures.These captures happen when two BHs have a close approach, radiate away orbital energy and merge in a short timescale due to the emission of GWs (Samsing 2018).According to the scattering experiments by Zevin et al. (2019), binary-binary interactions are roughly five times more likely than binary-single interactions to produce such mergers due to their more complex resonant intermediate states.
The binary-binary interactions described in this paper have at least one binary near the hard-soft boundary.During a resonant intermediate state, two of the BHs can enter a highly eccentric orbit that ends in a GW capture.This orbit, therefore, starts with a large semimajor axis ( ∼  h ), which requires a nearly radial orbit ( ∼ 1) to trigger a GW capture.We can use the capture criterion by Samsing (2018) to quantify what the eccentricity, , must be To put it in context, for a 10 5 M ⊙ cluster with  h = 1 pc, a capture between two 10 M ⊙ BHs would start at ∼ 60 AU and 1 −  ∼ 10 −6 .As the binary inspirals, this eccentricity may be radiated away.To estimate this, we find the peak GW frequency  from Wen (2003), (1 + ) 1.1954  (1 −  2 ) 3/2 , (42) which together with the Peters' equations (Peters 1964) allows us to compute the eccentricity when the binary enters the LIGO-Virgo-KAGRA frequency band (  ≃ 10 Hz).What we find is that, due to the large initial value of , the eccentricity is radiated away before the GW emission becomes observable and thus the waveform appears nearly circular.For the above values of ,  h ,  and  BH , we find a merger that has circularised to  10Hz ∼ 10 −2 , which is currently not detectable (Lower et al. 2018), but it is predicted that third generation GW detectors should be able to find this population of BBH mergers, if it exists.Furthermore, if we assume that after each resonant intermediate state the eccentricity is sampled from a thermal distribution (Heggie 1975), about ∼ 63% (∼ 0%) of captures that start with an initial semimajor axis of  h ( ej ) have a remaining  10Hz < 0.1 at 10 Hz.The sampling of  has a chance of producing mergers that form in the LIGO-Virgo-KAGRA band,  > 10 Hz, with a probability of ∼ 12% (∼ 46%).
An alternative path to mergers via binary-binary interactions is the formation of stable triples.Stable triples with sufficiently high eccentricity can have large jumps in the eccentricity of the inner orbit (Lidov-Kozai cycles, per Lidov 1962;Kozai 1962).As above, extreme values of the eccentricity may lead to GW captures.Per the same argument, the involved semimajor axes may be sufficiently large that the eccentricity may be of the order of  10Hz ∼ 10 −2 when the binary enters the LIGO-Virgo-KAGRA band.Modelling with post-Newtonian terms is required to quantify the importance of mergers via these two channels.
The short timescale of the above processes implies that the merger happens at a location near the multiple-body interaction.This constitutes an explanation for the finding that the majority of mergers occur in-cluster -as opposed to ejected binaries -in -body models of low-mass clusters ( ≲ 10 5 M ⊙ , Rastello et al. 2019;Banerjee 2021b;Barber et al. 2023).These models have  ≲ 1.5 × 10 5 , and for their initial mass function and metallicity  ≃ 0.04 and  ≃ 1/50 such that  BH ≲ 100, which is similar to the values considered in our two-mass -body models.As we explored different  and , our results can account for the different metallicities and other parameters in the more realistic models.The high fraction of in-cluster mergers found in these -body models are in contrast to fast models that assume a single active binary, such as cBHBd (Antonini et al. 2023), that predicts an in-cluster merger fraction of ∼ 40%.Although their assumption of a single, hard binary in the core at any given time is supported by our -body models, we here show that an important ingredient is missing in these fast models that would increase the contribution of dynamically assembled BBH mergers in relatively low-mass star clusters.Banerjee (priv. communication) indeed finds in-cluster mergers that form with a high  (≳ 100 AU) and extremely radial orbit, which nearly circularise at 10 Hz.These mergers may be a signature of binary-binary interactions at the hard-soft boundary.These binaries contribute to the eccentricity distribution by increasing the expected rate of mergers at lower eccentricities ( ≳ 10 −2 ).Such mildly eccentric mergers could be disentangled from quasi-circular mergers in future GW detectors, such as the Einstein Telescope (Maggiore et al. 2020) and Cosmic Explorer (Evans et al. 2021).In conclusion, we present an additional step towards a complete prediction for the rate of dynamically-formed BBH mergers that can be detected in current and upcoming GW experiments.
Figure 2. Ratio  h / c (equation 3) as a function of  for the single-mass  -body models (crosses with error bars).Fits for a constant  2 and an dependent  ′ 2 are shown with grey and red lines, respectively.

Figure 3 .
Figure3.Number of stars in the cluster core as a function of  for the single-mass  -body models (crosses with error bars).Using  i independent of  (gray) we are able to roughly match the  -body models, whereas using the  -dependant  ′ i coefficients (red) yields a very accurate match.

Figure 6 .
Figure 6.Ratio  2 c,BH / 2 BH (equation 12) as a function of  BH for the twomass  -body models (crosses with error bars).Fits for a constant  1,BH and an  -dependent  ′ 1,BH are shown with grey and red lines, respectively.

Figure 11 .
Figure11.Distribution of the binding energies of the binaries as a function of  for the single-mass  -body models.The theory (dotted lines) does not match the observed distribution (full lines), with the discrepancy being highest at lower  and higher .

Figure 12 .
Figure 12.Number of dynamically-formed stable triples  t as a function of  (crosses with error bars) in the single-mass  -body model.The non-zero value of  t implies the presence of binary-binary interactions.

Figure 14 .
Figure 14.Number of dynamically-formed stable BH triples  t,BH as a function of  BH and  (crosses with error bars) in the two-mass  -body model.The non-null value of  t,BH implies the presence of BBH-BBH interactions.

Figure 15 .
Figure 15.Distribution of the binding energies of the BBHs as a function of  BH and  for the two-mass  -body models with  = 50000 and  = 1/50.The theory (dotted lines) does not match the observed distribution (full lines), with the discrepancy being highest at lower  BH ,  and higher .

Figure 16 .
Figure 16.Dimensionless ratio of the BH triple formation rate Γ t,bb to the BBH formation rate Γ b as a function of the number of BHs  BH (top panel)and the mass fraction  (bottom panel).This ratio is much greater than one, which points towards the non-negligible relevance of BBH-BBH interactions when modelling the three-body BBH population.The ratio scales approximately as Γ t,bb /Γ b ∝  −1.2  −1.2 BH .

Table 2 .
Scaling laws for two-mass cluster properties in the steady post corecollapse evolutionary phase.The coefficients  1,BBH and  2,BBH are found from the  -body simulations.

Table 3 .
Values of the input parameters in the  -body runs of the single-mass clusters.Each row in the table represents a different cluster model, which is run multiple times (as shown in the last column) to obtain sufficient statistical significance.

Table 4 .
Values of the input parameters in the  -body runs of the two-mass clusters.Each row in the table represents a different cluster model, which is run multiple times (as shown in the last column) to obtain sufficient statistical significance.
Ratio  h,BH / c,BH (equation 13) as a function of  BH for the twomass  -body models (crosses with error bars).Fits for a constant  2,BH and an  -dependent  ′ 2,BH are shown with grey and red lines, respectively.
assume that the binding energies of binaries evolve because of interactions between Figure8.Number of BHs in the sub-cluster core as a function of  BH for the two-mass  -body models (crosses with error bars).Using  i,BH independent of  BH (gray) we are able to roughly match the  -body models, whereas using the  BH -dependant  ′ i,BH coefficients (red) yields a good match.