© The Authors 2025

1 Introduction

The detection of gravitational waves (GW) from the merger of binary black holes (BBH) by the LIGO-Virgo-Kagra (LVK) interferometers is now a routine phenomenon. More than a hundred BBH mergers have been reported in the first three observing runs (Abbott et al. 2023a) and the ongoing fourth observing run. Despite the rapid increase of available data, the origin of the GW sources is still unclear.

The two classical formation pathways for BBH mergers are the isolated channel and the dynamical channel. In the former, the BBH forms from a binary of massive stars in which the system evolves with negligible interaction with its surroundings and the orbit shrinks by either common envelope evolution (e.g. Belczynski et al. 2007; Postnov & Yungelson 2014; Belczynski et al. 2016) or stable mass transfer (e.g. Inayoshi et al. 2017; Gallegos-Garcia et al. 2021; van Son et al. 2022a,b). Other related pathways are chemically homogeneous evolution and Population III binaries (see Mandel & Broekgaarden 2022, Sect. 4.1 for a review on the topic). In the dynamical channel, the formation and evolution of the binary depend on the interactions of the BBH with other black holes (BHs) and stars in dense stellar systems. These include active galactic nuclei (Bartos et al. 2017; Tagawa et al. 2020; Samsing et al. 2022;

Trani et al. 2024b), open clusters (Banerjee 2017, 2018a,b; Rastello et al. 2019), and globular clusters (GC) (Portegies Zwart & McMillan 2000; Rodriguez et al. 2016; Antonini & Gieles 2020; Trani et al. 2021; Antonini et al. 2023; Dall’Amico et al. 2024).

A ‘smoking-gun’ signal that differentiates the isolated and dynamical channels is the presence of residual eccentricity in the frequency band of the LVK interferometers (see also Antonini et al. 2025 for an analysis using the effective and precessing spin parameters χ_eff and χ_p). Once the binary is assembled, the backreaction from GW emission very efficiently circularises the orbit, usually before the binary becomes observable at LVK frequencies (Peters 1964). Therefore, for a BBH merger to have a non-negligible eccentricity, the BBH itself must be assembled very close in time to the moment of merger¹ (see e.g. Fig. 9 of Spera et al. 2022). This is only attainable in the dynamical channel, if the two BHs merge after an almost head-on approach.

In dense clusters, mergers due to the random encounter of two previously unbound BHs (single-single interaction) are rare, and the rate for this process is negligible (Samsing et al. 2020).

The leading mechanisms for the production of eccentric GWs are the merger of two BHs in an interaction involving a BBH and an unbound third BH (binary-single interaction) or two BBHs (binary-binary interaction). These interactions can have chaotic intermediate states (Hut & Bahcall 1983; Sigurdsson & Phinney 1993; Trani et al. 2024a) and thus the probability of having a close encounter of two BHs is higher than in single-single interactions.

In the models of Zevin et al. (2019) and Kremer et al. (2020), binary-single interactions are ∼10 times more common than binary-binary interactions. However, both these channels are found to roughly contribute the same to the merger rate, as their models predict that each binary-binary interaction is more likely to trigger a merger compared to a binary-single interaction.

It is not known whether the BBHs involved in binary-binary actions in the aforementioned simulations formed dynamically (so-called ‘three-body’ BBHs), or formed from massive star binaries (hereafter, primordial BBHs). If they are primordial BBHs, any prediction for the GW merger rate and eccentricity distribution is sensitive to uncertainties in binary properties at birth and binary evolution, such as common envelope efficiency. However, three-body BBHs form from previously unbound BHs via dynamical interactions (Heggie 1975; Aarseth & Heggie 1976; Goodman & Hut 1993; Atallah et al. 2024) and can merge in subsequent interactions. In this channel, the properties of the BBHs are determined only by the cluster dynamics, and therefore are largely independent of the uncertainties in the binary stellar evolution models. Previous studies have shown that three-body BBHs dominate the production of GW mergers in clusters (e.g. Hong et al. 2018). Dynamical formation tends to happen in the central, denser regions of the cluster and results in wide BBHs, most of which are easily disrupted due to interactions with single BHs. The remaining stable BBHs are likely to interact with one another, where the outcome of these BBH − BBH interactions is usually the ionisation of the wider binary. An equilibrium is achieved when a cluster only has one dynamically active BBH, as further dynamically assembled BBHs are rapidly ionised in BBH-BBH interactions (Marín Pina & Gieles 2024). When this happens, the BBH − BBH interaction rate roughly equals the BBH formation rate.

The above process enhances the BBH − BBH interaction rate, which could lead to a higher GW capture rate than predicted by (fast) models that only consider binary-single interactions. Furthermore, according to Goodman (1984) the number of dynamical binaries (which is the result of a balance between binary formation and disruption) is expected to be higher in clusters with a low number of BHs, because the central dynamics are determined by the properties of the BH sub-system (Breen & Heggie 2013). This points towards the higher relative importance of BBH-BBH interactions in low-mass clusters or massive clusters that have ejected most of their BHs. Indeed, different approaches for modelling BBH mergers in clusters present a discrepancy at low initial cluster masses in their ratio of incluster to ejected mergers. For example, Antonini et al. (2019) predicts that, for the typical GC ($$\sim {10}^{5-6}{\mathrm{M}}_{\odot }$$), nearly half of all mergers are in-cluster; whereas for the N-body models of lowermass clusters of Banerjee (2021) this ratio is $$\gtrsim 80\%$$. A recent work (Barber & Antonini 2025) argues that almost all mergers of non-primordial binaries happen in-cluster. GW captures in BBH − BBH interactions, which happen only in-cluster and are included in the N-body models but not the fast (approximate) models, may explain the differences found.

In this paper, we link the origin of BBH − BBH interactions in GCs to dynamical BBH formation and study their outcomes, including GWs with measurable eccentricity and dynamical triple formation. The paper is organised as follows: in Section 2 we enumerate our initial conditions and summarise the methods used to perform the simulations. In Section 3, we present a model for the population of BBHs in star clusters and their BBH-BBH interactions. In Section 4, we put forward a model for the GW outcomes of these interactions. In Section 5, we study dynamical triple formation due to BBH − BBH interactions and GW mergers due to von Zeipel-Lidov-Kozai oscillations in these triples. Our results are discussed in Section 6, and the conclusions are summarised in Section 7.

2 Methods

2.1 Summary of CMC

We use the public suite of GC models run with Cluster Monte Carlo (CMC). CMC is a highly efficient code designed to accurately simulate the time evolution of dense collisional star clusters. It uses a Hénon-like approach (Hénon 1971), where the position of stars is randomly sampled from their orbit around the cluster and the relaxation from multiple weak twobody encounters is treated as a single effective encounter (Joshi et al. 2000). This shortens the computing time and allows for an exploration of the parameter space of GC initial conditions. However, compared to earlier approaches, the integration is done on a star-by-star basis, which allows it to recover short-range physics such as stellar evolution and, more importantly for this paper, strong few-body interactions.

In CMC, strong interactions are sampled from pairs of neighbouring stars and/or binaries according to their interaction rate. This interaction rate is computed from their relative cross section, defined as the area covered within the impact parameter that leads to a pericentre passage of r_p if the binaries were taken to be point masses. For binary-single encounters, they fix $$r_{\mathrm{p}}=Xa$$ (with a the semimajor axis of the binary and X a numerical constant), whereas for binary-binary encounters they use $$r_{\mathrm{p}}=X(a_0+a_1)$$ (with a₀, a₁ the semimajor axes of the binaries). The authors find that using a prefactor of X = 2 captures almost all the relevant energy-generating binary interactions (Fregeau & Rasio 2007). For this paper, we define ‘strong’ interactions as those with X ≤ 2. After sampling, these strong interactions are numerically resolved using FEWBODY (Fregeau 2012). Further details on the CMC code can be found on Rodriguez et al. (2022).

We use the catalogue of 148 CMC models of Kremer et al. (2020), which is constructed by varying the number of stars N, viral radius r_v, metallicity Z, and galactocentric distance R_G. The choice of initial conditions represents a sampling of localUniverse GCs, and is $$N\in [2,4,8,16]\times {10}^5,r_{\mathrm{v}}\in [0.5,1,2,4]\mathrm{p}\mathrm{c}$$, $$Z\in [0.0002,0.002,0.02]$$, and $$R_{\mathrm{G}}\in [2,8,20]\mathrm{k}\mathrm{p}\mathrm{c}$$. This accounts for 144 models, plus four extra models with $$N=3.2\times {10}^6$$, $$R_{\mathrm{G}}=20\mathrm{k}\mathrm{p}\mathrm{c},Z\in [0.0002,0.02]$$, and $$r_{\mathrm{v}}\in [1,2]\mathrm{p}\mathrm{c}$$.

The advantage of using a Hénon-like approach to simulating GCs, as opposed to a direct N-body integration, is that it allows us to perform a Monte Carlo sampling of the strong interactions. As described above, each strong encounter is considered as a scattering with fixed velocity, impact parameter, masses, semimajor axes, eccentricities and stellar type and radius. These parameters are shown schematically in Fig. 1. The other 8 scattering parameters, which characterise the binaries’ phases and orientations (for each binary, the inclination, argument of the periapsis, true anomaly and longitude of the ascending node), are randomly sampled. This allows us to re-simulate all the interactions with different parameters to increase our statistics without the need to perform costly simulations of entire GCs. We expand on this in the following section.

Fig. 1 Schematic diagram of a BBH − BBH interaction with some of its relevant parameters.

2.2 Summary of TSUNAMI

For each set of BBH − BBH interactions found in CMC we ran ten different scattering experiments, sampling the phases and orientations of the binaries. This sampling was carried out using the FEWBODY code from Fregeau (2012); however, we chose not to run the scattering simulations with FEWBODY, as it only includes the Newtonian and 2.5PN term. Instead, we used TSUNAMI (Trani et al. 2019), a state-of-the-art direct few-body integrator with post-Newtonian (PN) treatment of the equations of motion up to the 3.5PN term (inclusive), which is needed to accurately compute the dynamics before we identify a GW merger. TSUNAMI uses a logarithmic Hamiltonian and time-transformed leapfrog algorithm (Mikkola & Tanikawa 1999a,b) combined with a Bulirsch-Stoer extrapolation algorithm (Stoer & Bulirsch 1980). Each scattering was evolved in time until the interaction was resolved according to the criteria of Appendix A. Specifically, we assumed that a capture happens when the distance between two BHs is smaller than 10 times the sum of their Schwarzschild radii; we discuss this choice in Sect. 6.1. During an interaction, some configurations may lead to arbitrarily long-lived unstable intermediate states. In order to account for these cases, we label any scattering experiment that reached a simulation time of $$t>{10}^6{\mathrm{A}\mathrm{U}}^{3/2}{G}^{-1/2}{\mathrm{M}}_{\odot }^{-1/2}=0.16\mathrm{M}\mathrm{y}\mathrm{r}$$ as unresolved.

3 BBH-BBH interactions

In this section, we characterise the population of BBHs in star clusters and the mechanism that drives them to interact.

3.1 Three-body BBHs

The population of BBHs in a GC is not a simple fiducial evolution of the population of massive stellar binaries. After an initial settling phase, the cluster reaches mass segregation, where the heaviest components - the BHs, after stellar evolution - populate the cluster’s core. Per Hénon (1961), the energetic budget of a cluster is provided by the core, and thus the global properties of the cluster dictate the few-body BH dynamics in the central region. Specifically, the energy that fuels the cluster evolution is supplied by BBHs in the core, where a binary acts as a reservoir of negative binding energy, thus supplying (positive) energy to the rest of the cluster via interactions. According to

Heggie (1975); Goodman (1984), after core collapse there is a state of steady evolution where the energy demands of the cluster are balanced by the creation and hardening of binaries. When including BHs, it is the BBH formation in the BH subcluster that provides energy to the BH subcluster, which is in thermal contact with the stars (Breen & Heggie 2013). These binaries are formed dynamically from the encounter of at least three unbound BHs (Heggie 1975; Tanikawa et al. 2012; Atallah et al. 2024; Ginat & Perets 2024), so they are usually referred to as three-body binaries. When formed, these BBHs have large semimajor axes, a, with a distribution that formally diverges as $$a\to \infty$$ (Heggie 1975; Retterer 1980). This is because in a cluster, any two very separated BHs with sufficiently low relative velocity may be considered as bound. We avoid this issue by only considering ‘hard’ BBHs, which have a (positive) binding energy $$x=G{m}_1{m}_2/\left(2a\right)$$ that is greater than $$m_{\mathrm{B}\mathrm{H}}{\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}^2$$, with m_BH and σ_c,BH the mean mass and velocity dispersion of BHs in the cluster core, respectively. The hardness of a binary can be quantified via its hardness ratio $$z\equiv x/\left(m_{\mathrm{B}\mathrm{H}}{\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}^2\right)$$, so that hard binaries have $$z\ge 1$$. We only consider hard binaries because those that are less bound (i.e. ‘soft’, $$z<1$$) are easily dissociated by interactions with other unbound BHs. Opposite to this, Heggie (1975) showed that hard BBHs tend to become more bound (harden) via these interactions, as a result imparting some kinetic energy to the third body and the centre of mass of the binary, which provides energy to the cluster. These BBHs end their life when their recoil kick after an interaction is enough to eject them from the cluster, when their semimajor axis shrinks enough such that they merge before the next interaction takes place, or when they undergo a GW capture during an interaction. In Section 6.2, we discuss the effect of including soft BBHs interactions to the GW capture rate.

The population of BBHs in a GC can be divided in three separate categories: ‘primordial’², formed from the isolated evolution of a binary of massive stars without dynamical effects; ‘three-body’ binaries, where the component BHs were originally unbound stars and the binary is assembled dynamically; and ‘exchange’ binaries, that were formed dynamically but at least one of the component BHs was originally in a primordial binary.

It is reasonable to assume that a high primordial BBH fraction leads to a high BBH − BBH interaction rate. However, it is not explicitly shown in studies (Banerjee 2021; Kremer et al. 2020; Zevin et al. 2019) that BBH − BBH interactions are indeed the result of primordial BBH. We first analyse the nature of the BBHs involved in BBH − BBH interactions. Specifically, this is done by checking whether the IDs of the BHs involved in the interaction were in binaries in the first snapshot of the simulation.

In our set of CMC simulations, we find that less than 0.1% of the strong³ BBH − BBH interactions are among two primordial BBHs. Most of the interactions happen between two three-body binaries (84%), with 15% of interactions involving at least one exchange binary and ∼1% involving two exchanged binaries. Therefore, we find that three-body BBHs dominate in BBH − BBH interactions. This can be explained by the continuous formation of BBHs near the hard-soft boundary and their immediate disruption in BBH-BBH interactions (Marín Pina & Gieles 2024).

Therefore, in order to study the contribution of binary-binary interactions to the BBH merger rate in general and the eccentric GW production in particular, we focus on the dynamically assembled three-body BBHs.

3.2 BBH − BBH interaction rate

The goal of this section is to understand the cluster processes that lead to close BBH-BBH interactions. Following Marín Pina & Gieles (2024), we assume that there is only one three-body BBH at any given moment in time, that continuously interacts with - and ionises - wide, recently formed three-body BBHs. This model predicts that most interactions among three-body BBHs will have large semimajor axis ratios, as these new three-body BBHs form closer to the hard-soft boundary, while the long-lived stable binary has undergone a process of hardening via binarysingle interactions.

This is indeed what we observe in the CMC simulations, where we find a median semimajor axis ratio, α, (defined as α =$$a_{\mathrm{\text{max}}}/a_{\mathrm{\text{min}}}$$, with a_max and a_min the semimajor axis of the largest and smallest binary of the interaction, respectively) of $$\alpha \sim 50$$. The complete distribution, shown in Fig. 2, can be explained by the interaction of a binary near the hard-soft boundary ($$z\sim 1$$, in the simulations we find a median $$z\simeq 8$$ for the BBH with a = a_max) with a stable, long-lived binary ($$z\gg 1$$, in the simulations we find a median $$z\simeq 600$$ for the BBH with $$a=a_{\mathrm{\text{min}}}$$). We now assume that the distribution of energies of stable binaries follows some power-law $$p\left(z\right)\propto z^{-b}$$. From the models of Marín Pina & Gieles (2024, Fig. 14), the index of the exponent is $$0.9\lesssim b\le$$ 1.4. Then, assuming that a_max is very large, the distribution of α must follow the same power-law dependence, $$p\left(\alpha \right)\propto {\alpha }^{-b}$$. This is indeed what is seen in Fig. 2, where the distribution can be described by a power law with exponent b = 1.1.

In order to quantify the relative importance of BBH − BBH interactions, we will compare the interaction rate of the stable BBH with other $$\mathrm{B}\mathrm{B}\mathrm{H}\mathrm{s},{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{b}}$$, to the rate of interaction of this same central binary with single $$\mathrm{B}\mathrm{H}\mathrm{s},{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{s}}$$. As the BBH − BBH interaction rate is limited by dynamical BBH formation, it is roughly equal to the binary formation rate per unit volume and binding energy, C, which is (Goodman & Hut 1993)

$$C\left(x\right)=\frac{2{\pi }^2}{3\sqrt{3}}\frac{n_{\mathrm{B}\mathrm{H}}^3{G}^5{m}_{\mathrm{B}\mathrm{H}}^4}{{\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}^{11}}\frac{z^{-9/2}}{(1+(5z{)}^{-1})(1+e^{-z})},$$(1)

with n_BH the central number density of the BHs. Note that this is different to the net formation rate $$\bar{C}$$ discussed in Goodman & Hut (1993, see also Atallah et al. 2024; Ginat & Perets 2024), which is the result obtained after considering binary ionisation due to binary-single interactions, as well as hardening and softening through the hard-soft boundary. We assume that these processes happen within the core volume, V_c, so the binary-binary interaction rate is

$${\mathrm{\Gamma }}_{\mathrm{b}\mathrm{b}}=m_{\mathrm{B}\mathrm{H}}{\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}^2{V}_{\mathrm{c}}{\int }_1^{\infty }C\left(z\right)\mathrm{d}z.$$(2)

For the binary-single interaction rate, we will use the $$n-\mathrm{\Sigma }-v$$ approach (as in e.g. Hills & Day 1976)

$${\mathrm{\Gamma }}_{\mathrm{b}\mathrm{s}}=n_{\mathrm{B}\mathrm{H}}⟨{\mathrm{\Sigma }}_{\mathrm{b}\mathrm{s}}v⟩=n_{\mathrm{B}\mathrm{H}}{\int }_0^{\infty }{\mathrm{\Sigma }}_{\mathrm{b}\mathrm{s}}\left(v\right)f_{\mathrm{\text{Maxwell}}}\left(v\right)v\mathrm{d}v,$$(3)

where v is the relative velocity between single BH and BBH, and we assume equipartition between binaries and singles and a Maxwellian velocity distribution, $$f_{\mathrm{\text{Maxwell}}}\left(v\right)\propto$$$$v^2\mathrm{exp }\left(-v^2/\left(3{\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}^2\right)\right)$$, where $${\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}$$ is the (one-dimensional) central velocity dispersion of the BHs. We use the cross-section given by (Heggie & Hut 2003, p. 217)

$${\mathrm{\Sigma }}_{\mathrm{b}\mathrm{s}}=\frac{6\pi G{m}_{\mathrm{B}\mathrm{H}}a}{v^2}.$$(4)

As discussed above, we assume that the stable binary has a binding energy ∼50 times larger than the hard-soft boundary. We can write the dependency on $${\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}$$ in Equations (2) and (3) in terms of central density by using the definition of the core radius of the BH subsystem, $$\left(4\pi G/9\right){\rho }_0{r}_{\mathrm{c},\mathrm{B}\mathrm{H}}^2={\sigma }_{\mathrm{c},\mathrm{B}\mathrm{H}}^2$$ (Heggie & Hut 2003, p. 71). Finally, taking into account that for the isothermal and King (1966) models with high concentration, the central density can be expressed in the average density within the core as $${\rho }_0\simeq 1.9{\rho }_{\mathrm{c},\mathrm{B}\mathrm{H}}$$, we can express the ratio of rates as

$$\frac{{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{b}}}{{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{s}}}=\frac{5.7}{N_{\mathrm{c},\mathrm{B}\mathrm{H}}},$$(5)

where N_c,BH is the number of BHs in the core of the BH subcluster. For the two-mass model of Breen & Heggie (2013), this quantity depends only on the number of BHs in the cluster (Marín Pina & Gieles 2024, Equation (25)), so that

$$\frac{{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{b}}}{{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{s}}}\simeq 0.3{\left(\frac{N_{\mathrm{B}\mathrm{H}}}{{10}^2}\right)}^{-1/3}.$$(6)

This result implies that clusters with fewer BHs tend to have more binary-binary interactions with respect to binary-single interactions. In Fig. 3 we compare this analytic result to what we extract from the CMC models and we find good agreement. We can compare this to the result that we would obtain if we were to assume that the interactions are among primordial binaries. If we consider a BH binary fraction of f_b, the ratio of rates would be (Marín Pina & Gieles 2024)

$$\frac{{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{b}}}{{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{s}}}=3.1f_{\mathrm{b}}.$$(7)

For a reference time of ∼100Myr after cluster formation, all our cluster simulations have f_b below the initial binary fraction of 5%, likely due to BBH disruption at birth. In Fig. 3 we show the interaction rate ratio for the limiting case of $$f_{\mathrm{b}}=5\%$$. This prediction does not recover the N_BH dependence and thus cannot explain BBH-BBH interactions. We conclude that the high BBH − BBH interaction rate in the CMC models can only be explained by three-body binaries. In the following sections, we study the outcomes of these interactions: GW captures and triple formation.

Fig. 2 Probability density function of the ratio of the semimajor axes in BBH − BBH interactions, as calculated from the CMC models.

Fig. 3 Ratio of the BBH − BBH interaction rate and the binary-single interaction rate as a function of the cluster’s number of BHs at the time the interactions took place. The blue dots correspond to the median values obtained from the CMC models (contour is the 64% region). In green, the analytical prediction of Eq. (6), which considers the interaction of three-body BBHs ; yellow is the naïve prediction of Eq. (7), which only considers the interactions of primordial binaries.

4 GW captures in BBH-BBH interactions

In this section, we study the GW captures in BBH-BBH interactions. We begin by deriving an analytical approximation to the probability of merging (Sect. 4.1), which we compare to the set of our BBH-BBH scatterings sampled from realistic star cluster models (Sect. 4.2). In Sect. 4.3, we define the astrophysical weighting of our cluster models. In Sects. 4.4 and 4.5, we obtain predictions for the merger rate and GW eccentricity distribution, as well as redshift dependence.

4.1 Analytical approximation

In order to do a full characterisation of the population of GW captures in BBH − BBH populations, it is necessary to first understand the mechanism that drives the captures in each interaction. In this section, we will present a simple analytical model for these captures.

The gravitational four-body problem is chaotic, so in general we expect small differences in the initial conditions to grow exponentially. Therefore, our goal in this section is to provide distributions of outcomes for fixed initial parameters. We work in the gravitational focusing limit, i.e. $$v_{\infty }\ll v_{\mathrm{\text{crit}}}$$, with v_∞ the relative velocity at infinity and v_crit the critical velocity for which the total energy in the centre-of-mass frame of the four-body system is zero. In the gravitational focusing limit, we can set the relative velocity at infinity to zero. This approximation is justified by the cluster simulations, in which only ∼3% of the BBH-BBH interactions have $$v_{\infty }/v_{\mathrm{c}}>1$$, and over half of them have $$v_{\infty }/v_{\mathrm{c}}<0.1$$. In general, the masses of the BHs involved in BBH − BBH interactions are comparable, so we will work in the equal mass case in this Section. Then, the mass can be rescaled so the only parameters are α and a_max.

BBH-BBH scatterings can be separated in two different regimes. Scatterings with $$\alpha \gtrsim 1$$ are likely to be resonant with many intermediate states (Zevin et al. 2019), and thus have a high probability of resulting in a GW capture. During interactions with $$\alpha \gg 1$$, however, the tight binary may pass through the wider binary unaffected, or only have an encounter with one of its BHs. Therefore, interactions in this second regime are less likely to lead to GW captures of the wider binary, and the only possibility for GW mergers are those where a_min is sufficiently small that the inner binary was already close to merger before the interaction. We expect the change in these two regimes to happen at some α_crit. This can be measured from the scattering simulations by noticing that all mergers at $$\alpha >{\alpha }_{\mathrm{\text{crit}}}$$ happen between the two BHs in the inner binary, whereas mergers at $$\alpha <{\alpha }_{\mathrm{\text{crit}}}$$ may happen between any two BHs. This is explored in the following section.

In the case of resonant BBH − BBH scatterings, their intermediate states can be classified in two qualitatively different categories: ‘democratic states’, where the distance between the four BHs is comparable, and ‘hierarchical states’, where the distance between two of the BHs is much shorter than the other distances and the system behaves like an (unstable) hierarchical multiple. Democratic states are highly chaotic, so any information about previous orbital properties is erased. This allows us to model a resonant BBH − BBH interaction as a series of hierarchical states, each of which with a semimajor axis and eccentricity sampled from a certain distribution. This approach has been developed by Samsing et al. (2014) for the binary-single case. We begin by explaining its key arguments, and later extend the approach to the four-body case.

In the framework of Samsing et al. (2014), each binary-single interaction has a mean number of hierarchical intermediate states (IMS) of $$N_{\mathrm{\text{IMS}}}=20$$ (Samsing et al. 2014; Rando Forastier et al. 2025). Each interaction is considered to end in a GW capture if at any given hierarchical intermediate state the eccentricity is such that, in a single pericentre passage, the energy radiated via GW emission is equal to the binding energy of the original binary. To leading order (Hansen 1972), the critical eccentricity for the binary-single case, $$e_{\mathrm{\text{crit}},3}$$, is

$$e_{\mathrm{\text{crit}},3}\simeq 1-2.96{\left(\frac{Gm}{a{c}^2}\right)}^{5/7}.$$(8)

The probability of a capture in a single IMS is the fraction of the eccentricity distribution that is above e_crit. Assuming that the eccentricity of any intermediate state is drawn from a thermal distribution, we obtain that

$$p_{\mathrm{\text{IMS}}}={\int }_{e_{\mathrm{\text{crit,}}3}}^{1}2e\mathrm{d}e=1-e_{\mathrm{\text{crit,}}3}^2\simeq 2.96{\left(\frac{Gm}{a{c}^2}\right)}^{5/7}.$$(9)

The probability of merging in a binary-single interaction, $$p_{\mathrm{\text{merge}},3}$$, is equal to the probability of a capture in any of the IMSs, which is

$$p_{\mathrm{\text{merge}},3}=1-{(1-p_{\mathrm{\text{IMS}}})}^{N_{\mathrm{\text{IMS}}}}=1-e_{\mathrm{\text{crit}},3}^{2N_{\mathrm{\text{IMS}}}}.$$(10)

In the BBH − BBH case, we extend this framework by considering a critical eccentricity for a reference binary whose semimajor axis is equal to that of the smallest initial binary, a_min, with component masses equal to the mean mass of all four BHs,$$⟨m⟩$$, i.e.

$$e_{\mathrm{\text{crit}}}\simeq 1-2.96{\left(\frac{G⟨m⟩}{a_{\min }{c}^2}\right)}^{5/7}.$$(11)

We use a_min as the reference semimajor axis because the tighter binary is more likely to merge. Furthermore, we assume that the merger probability has a correction factor f(α) with respect to the binary-single case, such that the probability of merging in a BBH − BBH interaction is

$$p_{\mathrm{\text{merge}}}=f\left(\alpha \right)(1-e_{\mathrm{\text{crit}}}^{2N_{\mathrm{\text{IMS}}}})\simeq f\left(\alpha \right)\times 2N_{\mathrm{\text{IMS}}}(1-e_{\mathrm{\text{crit}}}).$$(12)

This factor accounts for binary-binary interactions being more chaotic than binary-single interactions and thus having more complex resonances. Note that, because $$e_{\mathrm{\text{crit}}}\simeq 1$$, increasing f(α) by a fixed factor is equivalent to increasing N_IMS. Thus, for comparison with the binary-single case, we fix $$N_{\mathrm{\text{IMS}}}=20$$. Since f(α) is purely Newtonian, it can only depend on α, and not the semimajor axes themselves. In the next section, we obtain a fit for f(α) from our BBH-BBH simulations.

4.2 Scattering simulations

We now compare the above analytical approximation for GW captures in BBH-BBH interactions against our detailed TSUNAMI simulations sampled from the CMC cluster models. From the $$\sim 2\times {10}^6$$ scatterings that we performed, 11% were not resolved before the maximum integration time. From the resolved interactions, the most common outcome (53%) is binary preservation⁴, where two BBHs remain after the scattering, followed by single ionisation (26%) and triple formation (21%). Mergers only account for ∼0.1% of all resolved outcomes. For illustrative purposes, we show online the animation of the trajectories of the BHs during two different interactions whose outcomes are GW mergers for $$\alpha \simeq 1$$ and for $$\alpha >1$$.

As predicted in the previous section, we find two different qualitative behaviours in GW mergers depending on their α (see Figs. 4 and 5). Interactions with $$\alpha <{\alpha }_{\mathrm{\text{crit}}}$$ can have complicated resonances and captures can happen between any two BHs, although we find a bias towards the merging of the smaller binary. Combinatorics would predict that, in fully mixed resonant encounters of four BHs, mergers of the smaller binary represent only 1/6 of all captures, with the larger binary representing another 1/6 and the rest, 2/3. However, this would only be the case if the systems are strongly resonant with equal masses and SMA. All our models have $$\alpha >1$$, which leads to a bias of the smaller binary merging, although we do see the expected ratios at $$\alpha \simeq 1$$. Also, about 30−50% of the interactions are non-resonant, meaning that they merge after only one IMS (Trani et al. 2024a; Rando Forastier et al. 2025). This mechanism biases the outcome ratios to overrepresent mergers with the input configuration.

Furthermore, we find that the ratio of merger types is roughly constant for $$3\lesssim \alpha \lesssim 100$$, although the uncertainties allow for a slight increase or decrease. This may be partly explained by a bias towards the merging of the smaller binary due to its shorter merger timescale. A binary in isolation with a SMA of a and a critical eccentricity $$e_{\mathrm{\text{crit}},3}$$ (Equation (8)) has a GW inspiral timescale (Peters 1964)

$$t_{\mathrm{G}\mathrm{W}}\left(a\right)\sim \frac{a^4}{\beta }{(1-e_{\mathrm{\text{crit}},3}^2(a\left)\right)}^{7/2}\propto a^{1.5},$$(13)

where β is a constant that only depends on the masses of the system. Thus, the ratio of the merger timescales of the initial binaries when they are at the critical eccentricity is

$$\frac{t_{\mathrm{G}\mathrm{W}}\left(a_{\max }\right)}{t_{\mathrm{G}\mathrm{W}}\left(a_{\min }\right)}\propto {\alpha }^{1.5}.$$(14)

In comparison, the duration of the interaction t_res can be estimated as N_IMS times the duration of the hierarchical state, each with an orbital period corresponding to a binary of $$a\sim a_{\max }$$, i.e.

$$t_{\mathrm{\text{res}}}\sim N_{\mathrm{\text{IMS}}}\sqrt{\frac{a_{\mathrm{\text{max}}}^3}{Gm}}.$$(15)

For highly eccentric binaries, $$t_{\mathrm{G}\mathrm{W}}\left(a_{\mathrm{\text{min}}}\right)\ll t_{\mathrm{\text{res}}}$$, but $$t_{\mathrm{G}\mathrm{W}}\left(a_{\mathrm{\text{max}}}\right)\gtrsim$$ t_res. In this case, the tight binary can merge in between intermediate states, while the larger binary may be perturbed by the other BHs before it can inspiral. Thus, these different timescales introduce a bias towards the merging of the smaller binary.

On the other side, interactions with $$\alpha >{\alpha }_{\mathrm{\text{crit}}}$$ generally only undergo mergers of the smaller binary. In this limit, the interaction is not truly resonant, so the only possible mergers are due to the smaller binary being so tight that it is driven to coalesce either due to GW emission alone or via a weak perturbation of the eccentricity in a direct interaction⁵ (Heggie & Rasio 1996). We find a transition at $${\alpha }_{\mathrm{\text{crit}}}\sim 500$$; combined with the steep α distribution, most mergers (89%) happen at $$\alpha <{\alpha }_{\mathrm{\text{crit}}}$$.

For encounters with $$\alpha <{\alpha }_{\mathrm{\text{crit}}}$$, we compute f(α) by substituting the merger probability, obtained in the TSUNAMI simulations, in Equation (12). We show the result in Fig. 6, where we plot f(α) (and corresponding uncertainties) in logarithmic bins, which we fit with the following function

$$f\left(\alpha \right)={\lambda }_1{(1+{\left(\frac{\alpha }{{\lambda }_2}\right)}^2)}^{-{\lambda }_3/2}.$$(16)

This function has the desired properties of a flattening at $$\alpha \to 1$$ and a power-law scaling at large α. Specifically, the value of λ₁ fixes the vertical scaling, λ₂ marks the transition from flattening to power-law, and λ₃ fixes the power-law index at large $$\alpha \left(f\right(\alpha )\propto$$ $${\alpha }^{-{\lambda }_3}$$). We perform a non-linear least-squares fit, for which we obtain $${\lambda }_1=3.4,{\lambda }_2=8.6$$, and $${\lambda }_3=1.7$$.

This shows that binary-binary interactions are more likely to trigger a GW capture when their semimajor axes are comparable, and decrease with α. In the case of α = 1, we recover the result of Zevin et al. (2019) where binary-binary interactions are about ∼5 times more likely to trigger a GW capture than binary-single interactions for any value of the semimajor axis. However, for $$\alpha >16$$ (which is about 65% of our interactions), the capture probability is smaller than for an equivalent binarysingle interaction. As a summary, the probability of merging in a BBH − BBH interaction with $$\alpha <500$$ is

$$p_{\mathrm{\text{merge}}}\simeq 0.034{\left(\frac{⟨m⟩}{20{\mathrm{M}}_{\odot }}\right)}^{5/7}{\left(\frac{a_{\mathrm{\text{min}}}}{0.1\mathrm{A}\mathrm{U}}\right)}^{-5/7}{(1+{\left(\frac{\alpha }{8.6}\right)}^2)}^{-0.83}.$$(17)

In the cluster simulations, BBH − BH interactions are about ∼7 times more frequent than BBH-BBH. However, Equations (10) and (17) allow us to compute the probability of capture for all the BBH − BH and BBH − BBH scatterings in the cluster simulations, respectively. We find that each BBH − BBH interaction is about ∼3 times more likely to trigger a merger than a BBH − BH interaction, so the total rate of captures in BBH − BBH interactions is about half that in BBH − BH interactions. This is due to both the contribution of the f(α) term and the different distribution of a_min in BBH-BBH interactions with respect to a in BBH-BH interactions.

Fig. 4 Fractions of GW mergers in BBH − BBH interactions as a function of the initial α. In green, mergers between the BHs of the smaller BBH ; in yellow, mergers between the BHs of the larger BBH ; in blue, mergers of any other two BHs. Each point represents the mean among multiple runs binned on their value of α, the shaded areas correspond to the standard deviation of the mean.

Fig. 5 GW mergers in BBH − BBH interactions in the initial $$a_{\mathrm{\text{max}}},a_{\mathrm{\text{min}}}$$ parameter space. In green, mergers between the BHs of the smaller BBH ; in yellow, mergers between the BHs of the larger BBH ; in blue, mergers of any other two BHs.

Fig. 6 Correction factor f(α) (defined in Sect. 4.1) as a function of the semimajor axes ratio α. A value of $$f\left(\alpha \right)>1$$ implies that the binarybinary interaction is more likely to trigger a GW merger than a binarysingle interaction where the binary has semimajor axis amin.

4.3 GC population model weighting

The grid of cluster models presented above sufficiently samples the parameter space of current-day Milky Way GCs (Kremer et al. 2020); however, in order to extrapolate our predictions to the GC population in the Universe, we must weight each model according to how common the distribution of its initial parameters is (e.g. in general, there are fewer clusters with high masses than clusters with lower masses). We follow a similar approach as in previous studies on the topic (e.g. Kremer et al. 2020; Antonini et al. 2023; Ye & Fishbach 2024). Firstly, we assume that the distribution of clusters in the Universe follows the quasi-separable distribution

$$\varphi (M_0,r_{\mathrm{v}},R_{\mathrm{G}},Z,z)\equiv {\varphi }_{M_0}\left(M_0\right){\varphi }_{r_{\mathrm{v}}}\left(r_{\mathrm{v}}\right){\varphi }_{R_{\mathrm{G}}}\left(R_{\mathrm{G}}\right){\varphi }_Z(Z,z){\varphi }_z\left(z\right),$$(18)

where each of the ϕ_i is a normalised probability density function. From these, we assume a Schechter-like initial cluster mass function (Schechter 1976)

$${\varphi }_{M_0}\left(M_0\right)\propto M_0^{-2}{e}^{-M_0/M_{\mathrm{s}}},$$(19)

with minimum mass $$M_{\mathrm{\text{min}}}={10}^4{\mathrm{M}}_{\odot }$$, maximum mass $$M_{\mathrm{\text{max}}}=2\times {10}^7{\mathrm{M}}_{\odot }$$, and cutoff mass $$M_{\mathrm{S}}={10}^{6.2}{\mathrm{M}}_{\odot }$$ (Antonini & Gieles 2020). We can use this to compute the number density of GCs formed across cosmic time, $$n_{\mathrm{G}\mathrm{C},0}$$. Assuming that the cluster mass formed per unit volume is $${\rho }_{\mathrm{G}\mathrm{C},0}\simeq 2.4\times {10}^{16}{\mathrm{M}}_{\odot }{\mathrm{G}\mathrm{p}\mathrm{c}}^{-3}$$ (Antonini et al. 2023), we obtain

$$n_{\mathrm{G}\mathrm{C},0}=\frac{{\rho }_{\mathrm{G}\mathrm{C},0}}{{\int }_{M_{\min }}^{M_{\max }}{M}_0{\varphi }_{M_0}\left(M_0\right)\mathrm{d}{M}_0}\simeq 5\times {10}^{11}{\mathrm{G}\mathrm{p}\mathrm{c}}^{-3}.$$(20)

We further weight our models with a log-normal distribution for the metallicity

$${\varphi }_Z(Z,z)\propto \frac{1}{Z}\mathrm{exp }\left(\frac{-{\mathrm{log }}_{10}^2(Z/⟨Z(z)⟩)}{2{\sigma }_Z^2}\right);Z\in [0,1],$$(21)

with parameters $${\sigma }_{\mathrm{Z}}=0.25$$ and $$⟨Z\left(z\right)⟩$$ as given by (Madau & Fragos 2017).

$${\mathrm{log }}_{10}⟨Z\left(z\right)/Z_{\odot }⟩\simeq 0.153-0.074z^{1.34}.$$(22)

We use a normal distribution for the cluster virial radii as in Fishbach & Fragione (2023)

$${\varphi }_{r_{\mathrm{v}}}\left(r_{\mathrm{v}}\right)\propto \mathrm{exp }\left(-\frac{{(r_{\mathrm{v}}-{\mu }_{r_{\mathrm{v}}})}^2}{2{\sigma }_{r_{\mathrm{v}}}^2}\right);r_{\mathrm{v}}\in [0.5,4],$$(23)

with mean $${\mu }_{r_{\mathrm{v}}}=2\mathrm{p}\mathrm{c}$$ and standard deviation $${\sigma }_{r_{\mathrm{v}}}=2\mathrm{p}\mathrm{c}$$. For the cluster formation rate, we follow the procedure in Antonini et al. (2023) and sample the formation redshift of our cluster models, $${\varphi }_z\left(z\right)$$, from the fiducial semi-analytical galaxy formation model of El-Badry et al. (2019), shown in Fig. 8. Finally, we weigh all three Galactocentric radii equally.

4.4 Merger rate

In order to obtain the merger rate, we have to account for initial cluster masses below and above the mass range in our catalogue. We do so by following a procedure similar to Kremer et al. (2020). The number of GW captures in BBH-BBH interactions N_m, shown in Fig. 7, can be described by a power-law dependence on the initial cluster mass. Within the uncertainties, this is independent on the cluster’s initial concentration (also shown in Fig. 7) and metallicity. The data can be well described by

$$N_{\mathrm{m}}={\left(\frac{M_0}{{10}^5{\mathrm{M}}_{\odot }}\right)}^{1/2}.$$(24)

We note that this scales slower than linearly and, therefore, the merger efficiency per unit cluster mass is higher in lowmass clusters. This is a direct consequence of the argument of Sect. 3.2, where we obtained that BBH − BBH interactions were more common in clusters with fewer BHs. If we average N_m with Equation (19), we obtain a mean number of mergers per cluster across the whole initial mass function of $$⟨N_{\mathrm{m}}⟩\simeq 0.6$$. It is important to note, however, that this result depends on the adopted M_min (and to a smaller degree, on M_max).

We use this result to compute the comoving merger rate $$ℛ\left(z\right)$$ by imposing that its integral over time is equal to the (massweighted) number of mergers per cluster, $$⟨N_{\mathrm{m}}⟩$$, times the number density of GCs formed across cosmic time, $$n_{\mathrm{G}\mathrm{C},0}$$, such that

$$⟨N_{\mathrm{m}}⟩n_{\mathrm{G}\mathrm{C},0}={\int }_0^{t_0}ℛ\left(z\right(t\left)\right)\mathrm{d}t,$$(25)

where z(t) is the function that converts from lookback time to redshift and t₀ is the age of the Universe, both assuming the Planck Collaboration VI (2020) cosmology.

The merger rate as a function of redshift for GW captures in BBH − BBH interactions is shown in Fig. 8. These interactions usually happen early in the cluster evolution when the cluster is at its densest. Furthermore, the time from BBH formation to merger is short for GW captures (compared to ejected mergers), so their distribution with redshift is a close tracer of the cluster formation history. This prediction should be independent of the specific details of cluster evolution, so the detection of GWs with non-zero eccentricity can be used to constrain cluster formation models (Romero-Shaw et al. 2021). In our case, we obtain that the peak rate of GW captures in BBH-BBH interactions happens at $$z\simeq 3.7$$. This result can be compared to the redshift distribution of the isolated channel, where the merger rate is predicted to follow the star formation rate (SFR, Madau & Dickinson 2014) with some delay. At $$z\simeq 0$$, the slope of our assumed cluster formation rate (CFR, El-Badry et al. 2019) is similar to that of the SFR, so that the slope of the merger rate as a function of redshift cannot be used to separate the dynamical and the isolated formation channels. However, we predict a peak rate at a higher redshift than the SFR, whereas the isolated channel predicts a peak at $$z\sim 2$$ (Santoliquido et al. 2020). This can be used to separate both channels, although it is a general prediction of cluster dynamics and not specifically of BBH − BBH interactions.

For GW captures in BBH − BBH interactions, we predict a local merger rate of

$$ℛ(z\simeq 0)=0.9{\mathrm{G}\mathrm{p}\mathrm{c}}^{-3}{\mathrm{y}\mathrm{r}}^{-1}.$$(26)

Assuming a total rate of BBH mergers of $$23.9{\mathrm{G}\mathrm{p}\mathrm{c}}^{-3}{\mathrm{y}\mathrm{r}}^{-1}$$ (Abbott et al. 2023b), we estimate that the BBH − BBH channel contributes to about ∼4%.

Fig. 7 Number of mergers per cluster due to GW captures in BBH − BBH interactions as a function of initial cluster mass, separated by initial rv (crosses) and averaged over initial rv (points).

Fig. 8 In blue, merger rate due to captures in BBH − BBH interactions as a function of redshift. In green, the same rate but only considering eccentric captures ($$e>0.1$$, measured at $$f_{22}=10\mathrm{H}\mathrm{z}$$ as explained in Sect. 4.5). In yellow, the cluster formation rate of El-Badry et al. (2019), in the same units as the merger rate.

4.5 Eccentricity distribution

To identify the GWs produced in a BBH-BBH interaction, it will be sufficient to study the dynamics of the two merging BHs as if they were isolated, because their relative separation is orders of magnitude smaller than the separation to the other two BHs in the interaction (but see Samsing et al. 2024, for the possibility of detecting phase shifts in the GW signal due to a nearby companion). Even with this simplification, the inclusion of PN terms causes the orbits to deviate from a Keplerian parametrisation, and therefore the eccentricity is ill-defined. There exist multiple ways of defining eccentricity in General Relativity in such a way that it converges to the Newtonian value in the weak-field regime (Shaikh et al. 2023). In PN theory, we can use a quasi-Keplerian parametrisation with three different values of the eccentricity: temporal e_t, radial e_r, and angular e_ϕ; which are only equal in the Newtonian regime. We discuss the nomenclature for these eccentricities in Appendix B. In this section we report on e_t; by convention, we will simply call it e, although keeping the above discussion in mind.

A further complication of this quasi-Keplerian parametrisation is the inclusion of non-conservative PN terms, which start at 2.5 PN. These terms describe the decay and circularisation of the orbit via GW emission, which implies that the eccentricity is no longer a constant of motion. The production of GWs is very efficient at circularising the binary, so it becomes necessary to define a point in the orbit at which to report the eccentricity. To overcome this issue, we report the value of the eccentricity when the frequency of the 22 -mode, f₂₂, is equal to 10 Hz. This is usually taken to be the moment at which a quasi-circular binary enters the LIGO-Virgo-Kagra frequency band. The choice of reference frequency (f₂₂) is motivated by recommendations in the literature that aim to standardise the definition of eccentricity (Ramos-Buades et al. 2022; Shaikh et al. 2023; Vijaykumar et al. 2024). Our definition matches that of Shaikh et al. (2023), used in GW data analysis (Vijaykumar et al. 2024). However, previous work in the field uses the peak GW frequency (Wen 2003) as their reference value. Both these definitions agree only at small e (Vijaykumar et al. 2024). We compute f₂₂ as

$$f_{22}\equiv 2f_{\mathrm{\text{orb}}}=\frac{\mathrm{\Phi }n}{\pi }$$(27)

where $$f_{\mathrm{\text{orb}}},\mathrm{\Phi }$$ and n are the orbital frequency, angle of return to the periapsis and mean motion, respectively. The last two quantities are computed in 3PN approximation as Memmesheimer et al. (2004, their equations $$25\mathrm{c},\mathrm{k}$$).

In order to obtain the eccentricity distribution of captures in BBH-BBH interactions as observable by an Earth-based GW detector, we use the list of mergers found in the above scattering experiments. We associate a weight with each merger, as described in Sect. 4.3, to take into account the distribution of cluster parameters across the Universe, with the total rate normalised as in Equation (25). The results are shown in Fig. 9.

For comparison, we also plot the simple analytical model of Sect. 4.1. During a four-body hierarchical intermediate state, the inner binary has a semimajor axis of a_IMS. If we model this IMS as a point-like inner binary with two orbiting BHs, due to conservation of energy we obtain a range for a_IMS of

$$a_{\mathrm{\text{IMS}}}\in [1,5]{(\frac{1}{a_{\mathrm{\text{min}}}}+\frac{1}{a_{\mathrm{\text{max}}}})}^{-1}.$$(28)

We obtain the predicted eccentricity distribution by sampling a thermal distribution of eccentricities above e_crit for each BBH − BBH scattering. Then, we evolve these eccentricities (starting at $$a=a_{\mathrm{\text{IMS}}}$$) using 2.5PN orbit-averaged equations (Peters 1964) until the GW emission reaches $$f_{22}=10\mathrm{H}\mathrm{z}$$. In Fig. 9 we show these results, weighted both according to their p_merge (Equation (17)) and their parameter distribution (Sect. 4.3).

As can be seen, the distribution of eccentricities is roughly flat in $${\mathrm{log }}_{10}e$$ between -2.5 and 0, with ∼10% of mergers having $$e>0.1$$; and the simple model can roughly reproduce the simulations. Our results for the eccentricity distribution can be used directly as a prior for Numerical Relativity studies in eccentric waveforms.

The above results assume that $$z\simeq 0$$, so that f₂₂ is the same in the source frame and the detector frame. If we fix the eccentricity to $$f_{22}=10\mathrm{H}\mathrm{z}$$ in the detector frame, we are actually measuring the eccentricity at higher frequencies in the source frame. The effect of this is to shift the eccentricity distribution towards lower values of e and reduce the fraction of eccentric mergers for higher z. In Fig. 8, we show the rate of mergers with $$e>0.1$$ at $$f_{22}=10\mathrm{H}\mathrm{z}$$ in the detector frame, as a function of redshift.

Note that, since our definition of eccentricity is different from the commonly used Wen (2003), there is not a clear comparison between our results and e.g. Antonini & Gieles (2020). However, both definitions agree at low e (Vijaykumar et al. 2024), so that the fraction of mergers with $$e>0.1$$ is roughly independent of the eccentricity definition. We find that ∼10% of our mergers have $$e>0.1$$, whereas Antonini & Gieles (2020) find that about half of their GW captures in binary-single scatterings have $$e>0.1$$. A similar fraction of eccentric mergers in binary-single scatterings was found in Trani et al. (2024a). This discrepancy can be attributed to both the different methods to integrate the orbit as well as the intrinsic differences between binary-single and binary-binary scatterings.

Fig. 9 Merger rate due to captures in BBH − BBH interactions as a function of the eccentricity (measured at $$f_{22}=10\mathrm{H}\mathrm{z}$$, as explained in Sect. 4.5). In blue, the results obtained from the TSUNAMI simulations. In green and yellow, the analytical model of Sections 4.1 and 4.5, obtained from the CMC data, and computed from the minimum and maximum aIMS (Equation (28)), respectively. For clarity, the area between these lines is coloured in light grey.

5 Triple formation and von Zeipel-Lidov-Kozai oscillations

In the Newtonian point-mass approximation, the possibility that a binary-single interaction leads to the formation of a stable triple is zero (Heggie & Hut 2003, p. 211). This can be seen heuristically by considering time reversibility: an isolated threebody system that is stable in the future needs to be stable in the past due to the symmetry of the Newtonian equations of motion under time reversal. Therefore, binary-binary interactions are the dominant dynamical process for stable triple formation. In this section, we explore dynamically assembled triples using our scattering simulations. As explained in Table A.1, we classify as triples any endstate that verifies the Mardling & Aarseth (2001) stability criterion (but see Lalande & Trani 2022; Hayashi et al. 2022, 2023, for more precise dynamical stability criteria).

In Fig. 10, we show the semimajor axes of the inner $$\left(a_{\mathrm{i}}\right)$$ and outer (a_o) binaries of the dynamically assembled BH triples in our models. Generally, these triples have large a_o, with a median semimajor axis ratio between outer and inner binary of ∼120, and a mode of ∼14. The distribution of eccentricities, also shown in Fig. 10, is close to thermal for both the inner (e_i) and outer (e_o) binaries. However, there is a lack of high-eccentricity outer binaries; this is because such triples are unstable (Mardling & Aarseth 2001). A similar distribution for dynamically formed triples in low-mass star clusters was found in Trani et al. (2022).

Since dynamically assembled triples generally have large a₀, they are easily destroyed or destabilised in triple-single interactions. We can estimate the timescale as $$t_{\mathrm{t}\mathrm{s}}=1/\left(n{\mathrm{\Sigma }}_{\mathrm{t}\mathrm{s}}v\right)$$. Here, n is the number density of stars around the triple, that we assume equal to the mean number density in the cluster core; v their relative velocity, which we assume equal to the root mean square density-weighted velocity at the cluster centre; and $${\mathrm{\Sigma }}_{\mathrm{\text{ts}}}$$ the cross-section of resonant triple-single interactions, as given by (Antognini & Thompson 2016, below Equation (26))

$${\mathrm{\Sigma }}_{\mathrm{t}\mathrm{s}}=\pi a_{\mathrm{o}}^2(\frac{\hat{a}}{{\hat{v}}^{3/2}}+{(\frac{{\hat{v}}^2}{\hat{b}}+\frac{{\hat{v}}^6}{\hat{c}})}^{-1})$$(29)

with parameters⁶ $$\hat{a}=8.89\times {10}^{-3},\hat{b}=0.275$$, and $$\hat{c}=0.044$$, and the velocity written as $$\hat{v}=v/v_{\mathrm{\text{crit}}}$$.

After a triple is dynamically assembled, the inner binary may change its orbital inclination and eccentricity due to the interaction with the outer component, an effect that is known as von Zeipel-Lidov-Kozai (ZLK) oscillations (von Zeipel 1910; Lidov 1962; Kozai 1962). If the eccentricity reaches a sufficiently high value, the inner binary may be driven to merge due to GW emission at pericentre. We will now estimate the relative contribution of this channel to the observable GW rate.

We begin by assuming that the ZLK effect is only relevant for triples with a sufficiently large inclination of the outer to the inner orbit, i₀. For triples with $$\left|\cos i_0\right|<3/5$$, the eccentricity of the inner binary can reach a value of $$e_{\mathrm{\text{max}}}^{\mathrm{\text{DA}}}$$, given by Naoz (2016); Mangipudi et al. (2022)

$$e_{\mathrm{\text{max}}}^{\mathrm{D}\mathrm{A}}=\sqrt{1-\frac{5}{3}{\mathrm{cos }}^2{i}_0}.$$(30)

The timescale for this process, t_ZLK, is (Liu et al. 2024)

$$t_{\mathrm{Z}\mathrm{L}\mathrm{K}}=2\pi \sqrt{\frac{a_{\mathrm{i}}^3}{G(m_{\mathrm{i},1}+m_{\mathrm{i},2})}}\frac{m_{\mathrm{i},1}+m_{\mathrm{i},2}}{m_{\mathrm{o}}}{\left(\frac{a_{\mathrm{o}}}{a_{\mathrm{i}}}\right)}^3{(1-e_{\mathrm{o}}^2)}^{3/2},$$(31)

where m_i,1 and m_i,2 are the masses of the components of the inner BBH and m_o is the mass of the outer BH. On average, a triple undergoes a number of ZLK cycles, N_ZLK, before being disrupted, with $$N_{\mathrm{\text{ZLK}}}=t_{\mathrm{\text{ts}}}/t_{\mathrm{\text{ZLK}}}$$. We will assume that a merger occurs if the maximum eccentricity reaches a critical value e_crit, which we define similarly to the previous section

$$e_{\mathrm{\text{crit}}}\simeq 1-2.96{\left(\frac{G^5{m}_{\mathrm{i},1}^2{m}_{\mathrm{i},2}^2(m_{\mathrm{i},1}+m_{\mathrm{i},2})}{2a_{\mathrm{i}}^5{c}^{10}}\right)}^{1/7}.$$(32)

Therefore, the inner binary of a triple merges due to the ZLK mechanism if $$\left|\cos i_0\right|<3/5,N_{\mathrm{L}\mathrm{K}}\ge 1$$ and $$e_{\mathrm{\text{max}}}^{\mathrm{D}\mathrm{A}}\ge e_{\mathrm{\text{crit}}}$$. We find that, on average a fraction 0.14% of the dynamically assembled triples in our simulations will have time to merge due to ZLK interactions before they are destabilised by an interaction. Even though there are three orders of magnitude more dynamically assembled triples than GW captures (21% and 0.1% of the resolved BBH − BBH outcomes in our simulations before model weighting, respectively), the GW capture rate is ∼3 times higher than the ZLK induced merger rate in dynamical triples. Then, taking the merger rate of Sect. 4.4, we obtain that for dynamically assembled triples, $$ℛ=0.3{\mathrm{G}\mathrm{p}\mathrm{c}}^{-3}{\mathrm{y}\mathrm{r}}^{-1}$$. This is compatible with previous studies, e.g. Martinez et al. (2020).

We note that this result only applies to triples that have formed in BBH − BBH scatterings, are stable, and then merged due to their secular ZLK evolution. Previous work using N-body models (e.g. Banerjee 2018b) agrees with our result that this channel is subdominant, instead arguing for mergers in dynamically formed unstable triples. This is equivalent to what we define as GW captures during BBH-BBH and BBH-BH encounters, as the interaction is not formally resolved if the triple is unstable.

Fig. 10 Probability density function of the semimajor axes (top) and eccentricity (bottom) for the inner and outer binaries of BH triples dynamically assembled via BBH − BBH interactions.

6 Discussion

6.1 Capture criteria

In this paper, we have assumed that a capture is inevitable when the distance between two BHs is smaller than 10 times the sum of their Schwarzschild radii. A trade-off exists; if we considered a larger cutoff value, we could mislabel hyperbolic encounters as mergers. For a smaller cutoff, we would be integrating the PN equations of motion beyond their validity range. For quasicircular orbits this is mostly fine, as the collision time between two stellar-mass BHs in such an orbit is of the order of seconds. However, when introducing highly eccentric or hyperbolic orbits, the merging time can be arbitrarily long. Some authors have used the orbital semimajor axis as the cutoff criteria, e.g. in the case of Zevin et al. (2019), the authors classify a merger if the semimajor axis is smaller than 10 times the sum of the BHs’ Schwarzschild radii, but as explained above, that allows the periapsis to be arbitrarily small, breaking the PN expansion. Furthermore, there can be hyperbolic encounters of two BHs (Capozziello et al. 2008; De Vittori et al. 2012; Bini & Geralico 2021; Gröbner et al. 2020; García-Bellido & Nesseris 2018) such that their orbit occurs in the strong field regime but their endstate is not a merger (Hopper et al. 2023; Damour et al. 2014; Jaraba & Garcia-Bellido 2021)⁷ The definition of a capture criteria is, then, beyond the scope of the PN expansion, and one must make use of full Numerical Relativity. This is computationally expensive and well beyond the scope of this work, although limited simulations in Numerical Relativity suggest that the energy change at our cutoff distance is such that, for the GW captures in this work, the merger cannot be avoided.

6.2 Soft binaries and wide scatterings

In a cluster, the definition of what counts as strong scattering is somewhat arbitrary. As explained in Sect. 2.1, we consider all interactions that have $$r_{\mathrm{p}}=X(a_0+a_1)$$, with X an arbitrary value that we fix at X = 2. Increasing X would include more scatterings whose most likely outcome would be binary preservation. As such, the outcome ratios in Sect. 4.2 are dependent on the value of X. Nevertheless, we argue that our choice of X is reasonable as it includes almost all resonances and is not too large as to be dominated by weak interactions that only slightly modify the BBH’s orbital parameters (Fregeau & Rasio 2007). However, we miss some of the direct encounters that can lead to a sufficiently large Δe to reach e_crit, while $$\mathrm{\Delta }E\simeq 0$$ (Heggie & Rasio 1996; Rando Forastier et al. 2025; Trani et al. 2024a).

Another limitation of our approach is with respect to soft binaries. These are usually disrupted in binary-single interactions, but when they undergo BBH − BBH scatterings there is a nonzero chance that they induce a GW merger. Per the limitations of the CMC code, our analysis does not include soft binary scattering. However, since interactions of soft and hard binaries have a large α, their contribution to the merger rate can safely be neglected. However, the interaction among two soft binaries may play a role, as their energy distribution increases very sharply with a. For a Saha-like distribution (Heggie 1975; Retterer 1980), the creation rate of soft binaries is $$C\left(x\right)x\propto x^{-7/2}$$. If we assume this to roughly equal the interaction rate Γ, working in the limit of α = 1 we obtain

$$\mathrm{\Gamma }\left(a\right)\simeq C\left(x\right)\left|\frac{\partial x}{\partial a}\right|\propto a^{3/2}.$$(33)

Then, the merger rate is $$ℛ\simeq C\left(x\right)p_{\mathrm{\text{merge}}}\propto a^{11/14}$$ (Equation (9)). This points towards the relevance of soft binary interactions, which may explain the enhanced in-cluster merger rate in the N body models of Barber & Antonini (2025). In a follow-up study, we consider this possibility in more detail.

7 Summary and conclusions

7.1 Summary of key findings

Our work confirms and extends on earlier work on BBH-BBH interactions:

• BBH-BBH interactions are less common than BBH-BH interactions (by a factor of ∼7) but, as they are more likely to trigger a GW capture, their contribution to the merger rate is about half that of BBH − BH interactions. This verifies the results of Zevin et al. (2019) for equal-mass BHs;

• Stable triple formation is a likely outcome of BBH-BBH interactions. These triples generally have very large semimajor axes, so GW captures due to von Zeipel-Lidov-Kozai oscillations in dynamically formed stable triples do not significantly contribute to the merger rate. This is in agreement with previous studies, e.g. Martinez et al. (2020).

Our key novel findings are:

• Most BBH-BBH interactions in GCs occur between threebody binaries and not among the primordial binary population;

• The rate and parameters of BBH − BBH interactions in GCs are largely independent of the primordial binary fraction and assumed binary stellar evolution models, as a consequence of the first point;BBH − BBH interactions tend to involve the dominant dynamically active BBH that interacts and destroys any further formed softer BBH ;

• An important consequence of the previous points is that the efficiency of GW captures depends only on the number of BHs in clusters, and is highest in low-mass clusters;

• Using a GC population model, we find that BBH-BBH interactions produce mergers at a rate of $$ℛ(z\simeq 0)=$$ $$0.9{\mathrm{G}\mathrm{p}\mathrm{c}}^{-3}{\mathrm{y}\mathrm{r}}^{-1}$$, and about 10% of these have $$e>0.1$$;

• This merger rate evolves with redshift according to Fig. 8. Its peak occurs slightly later than the peak of the cluster formation rate;

• The distribution of eccentricities in GW captures (Fig. 9) is roughly flat in $${\mathrm{log }}_{10}e$$ between -2.5 and 0.

7.2 Enhanced merger rate in low-mass clusters

The efficiency of GW captures in BBH-BBH interactions scales inversely with the cluster mass, which implies that it is most relevant for less massive clusters and clusters with fewer BHs. Here, we explore the difference between detailed cluster simulations and approximate models that assume a single dynamically active BBH. We argue that the discrepancy at lower masses can be explained by the missing mechanism of GW captures in BBH − BBH interactions.

In Fig. 11, we show the mean number of mergers per cluster in different mass regimes and for different methods of cluster simulation. For the most massive clusters ($$\gtrsim {10}^5{\mathrm{M}}_{\odot }$$) we plot the Kremer et al. (2020) catalogue generated by the CMC code, whereas for the less massive clusters ($$\lesssim {10}^5{\mathrm{M}}_{\odot }$$) we plot the catalogue of Banerjee (2021) generated by NBODY7 (Aarseth 2003, 2012; Nitadori & Aarseth 2012; Hurley et al. 2000, 2002). We assume that these models are detailed enough in their simulation to accurately reproduce the number of GW captures. Overlaid to these, we show the results of resimulating their initial conditions using the fast code CBHBD (Antonini et al. 2019, 2023). We use the same initial conditions and supernova models, and run the simulations for the same amount of physical time.

As can be seen, the number of mergers per initial cluster mass seems to follow a power-law⁸ $$N_{\mathrm{m}}\propto M_0^{1.4}$$. The CBHBD code can reproduce the number of mergers at higher cluster masses, but underpredicts the number of mergers at $$M_0\lesssim$$ $${10}^5{\mathrm{M}}_{\odot }$$. It appears that adding the effect of GW captures in BBH − BBH interactions bridges the gap with the merger rate of NBODY7 simulations. We plan to include this effect in CBHBD in the future.

Fig. 11 Number of BBH mergers with respect to initial cluster mass for different cluster simulations (crosses with error bars). In red, the number of mergers due to GW captures in GBH-BBH interactions, as in Equation (24). This channel appears to bridge the results of the different codes at lower initial cluster masses.

Data availability

The movies are available at https://www.aanda.org

Acknowledgments

DMP thanks Miguel Martínez and Giacomo Fragione for help in understanding the CMC output format. DMP and MG thank Kyle Kremer and Fabio Antonini for discussions. DMP and MG acknowledge financial support from the grants PRE2020-091801, PID2021-125485NB-C22, EUR2020112157, CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033 (State Agency for Research of the Spanish Ministry of Science and Innovation) and SGR-2021-01069 (AGAUR). AAT acknowledges support from the Horizon Europe research and innovation programs under the Marie Skłodowska-Curie grant agreement no. 101103134.

Appendix A Outcome of the interactions

Appendix B PN definition of eccentricity

As we discussed in Sect. 4.5, in PN theory there are three different values for the eccentricity: temporal e_t, radial e_r, and angular e_ϕ, which are only equal in the Newtonian limit. In this Section, we will give an interpretation to these values, following Memmesheimer et al. (2004). Let us begin by the Kepler equations in Newtonian mechanics

$$R=a(1-e\mathrm{cos }u)$$(B.1)

$$n(t-t_0)=u-e\mathrm{sin }u$$(B.2)

$$\varphi -{\varphi }_0=2\mathrm{arctan }\left[{\left(\frac{1+e}{1-e}\right)}^{1/2}\tan \frac{u}{2}\right]$$(B.3)

where a is the SMA, e is the eccentricity, u is the eccentric anomaly, n is the mean motion, and t is the time. The relative separation between the two bodies is given by the vector ($$R\mathrm{cos }\varphi ,R\mathrm{sin }\varphi ,0$$), and the suffix zero indicates the initial conditions. The above equations can be compared to the quasiKeplerian PN parametrisation of the orbit. For simplicity, in this Section we will only consider 1PN terms, for which the equations read (Damour & Deruelle 1985)

$$R=a_{\mathrm{r}}(1-e_{\mathrm{r}}\cos u)$$(B.4)

$$n(t-t_0)=u-e_{\mathrm{t}}\mathrm{sin }u$$(B.5)

$$\frac{2\pi }{\mathrm{\Phi }}(\varphi -{\varphi }_0)=2\mathrm{arctan }\left[{\left(\frac{1+e_{\varphi }}{1-e_{\varphi }}\right)}^{1/2}\tan \frac{u}{2}\right]$$(B.6)

where 2π/Φ is the angular advance of the periastron per orbit. The comparison between the two sets of equations explains the naming of the different eccentricities, as well as giving an intuition on their convergence in the Newtonian limit.

References

All Tables

All Figures

	Fig. 1 Schematic diagram of a BBH − BBH interaction with some of its relevant parameters.
In the text

	Fig. 2 Probability density function of the ratio of the semimajor axes in BBH − BBH interactions, as calculated from the CMC models.
In the text

Fig. 3 Ratio of the BBH − BBH interaction rate and the binary-single interaction rate as a function of the cluster’s number of BHs at the time the interactions took place. The blue dots correspond to the median values obtained from the CMC models (contour is the 64% region). In green, the analytical prediction of Eq. (6), which considers the interaction of three-body BBHs ; yellow is the naïve prediction of Eq. (7), which only considers the interactions of primordial binaries.

In the text

	Fig. 4 Fractions of GW mergers in BBH − BBH interactions as a function of the initial α. In green, mergers between the BHs of the smaller BBH ; in yellow, mergers between the BHs of the larger BBH ; in blue, mergers of any other two BHs. Each point represents the mean among multiple runs binned on their value of α, the shaded areas correspond to the standard deviation of the mean.
In the text

	Fig. 5 GW mergers in BBH − BBH interactions in the initial $$a_{\mathrm{\text{max}}},a_{\mathrm{\text{min}}}$$ parameter space. In green, mergers between the BHs of the smaller BBH ; in yellow, mergers between the BHs of the larger BBH ; in blue, mergers of any other two BHs.
In the text

	Fig. 6 Correction factor f(α) (defined in Sect. 4.1) as a function of the semimajor axes ratio α. A value of $$f\left(\alpha \right)>1$$ implies that the binarybinary interaction is more likely to trigger a GW merger than a binarysingle interaction where the binary has semimajor axis amin.
In the text

	Fig. 7 Number of mergers per cluster due to GW captures in BBH − BBH interactions as a function of initial cluster mass, separated by initial rv (crosses) and averaged over initial rv (points).
In the text

	Fig. 8 In blue, merger rate due to captures in BBH − BBH interactions as a function of redshift. In green, the same rate but only considering eccentric captures ($$e>0.1$$, measured at $$f_{22}=10\mathrm{H}\mathrm{z}$$ as explained in Sect. 4.5). In yellow, the cluster formation rate of El-Badry et al. (2019), in the same units as the merger rate.
In the text

Fig. 9 Merger rate due to captures in BBH − BBH interactions as a function of the eccentricity (measured at $$f_{22}=10\mathrm{H}\mathrm{z}$$, as explained in Sect. 4.5). In blue, the results obtained from the TSUNAMI simulations. In green and yellow, the analytical model of Sections 4.1 and 4.5, obtained from the CMC data, and computed from the minimum and maximum aIMS (Equation (28)), respectively. For clarity, the area between these lines is coloured in light grey.

In the text

	Fig. 10 Probability density function of the semimajor axes (top) and eccentricity (bottom) for the inner and outer binaries of BH triples dynamically assembled via BBH − BBH interactions.
In the text

	Fig. 11 Number of BBH mergers with respect to initial cluster mass for different cluster simulations (crosses with error bars). In red, the number of mergers due to GW captures in GBH-BBH interactions, as in Equation (24). This channel appears to bridge the results of the different codes at lower initial cluster masses.
In the text