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

ABSTRACT

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 (∼5), especially for low-N clusters (∼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 ionize one of the binaries involved. These also lead to triple formation and potentially gravitational wave 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)$|⁠.

black hole physics, gravitational waves, binaries: general, stars: kinematics and dynamics

1 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 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 (Antonini & Gieles 2020; Kremer et al. 2020; Kumamoto, Fujii & Tanikawa 2020; Banerjee 2021; Chattopadhyay et al. 2022; Mapelli et al. 2022; Rodriguez 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 characterizing 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¹ BBHs (Chattopadhyay et al. 2022; Torniamenti et al. 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, Chattopadhyay & Antonini 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, N_b, depends on the total number of stars in the cluster, N, as N_b ∝ N^−1/3. These models predict that a handful of binaries is expected in single-mass clusters of a few hundred stars. However, N-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 N-body calculations.

The paper is organized as follows: in Section 2 we revisit the theory for the population of three-body binaries. In Section 3, we present the N-body simulations, of single-mass and two-mass models. In Section 4, we compare the predictions and the N-body models. The cause of the aforementioned discrepancy is discussed in Section 5. In Section 6, we discuss its implications for GW observations.

2 THE NUMBER OF THREE-BODY BINARIES IN STEADY POST-COLLAPSE EVOLUTION

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-core-collapse 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 N ≳ 8000 stars (Bettwieser & Sugimoto 1984; Breeden, Cohn & Hut 1994) and in two-mass models where the number of BHs is N_BH ≳ 2000 (Breen & Heggie 2012). Both of these assumptions allow us to express relations between several cluster properties, which are discussed below and summarized in Tables 1 and 2.

2.1 Single-mass clusters

In post-collapse, we expect the cluster to be nearly isothermal within the half-mass radius r_h, so that the one-dimensional velocity dispersion is similar within the core and within r_h. For a cluster with mass M and energy E in virial equilibrium with a virial radius r_{|$\mathrm{ v}$|} = −GM²/(4E) and positions and velocities sampled from a Plummer (1911) model, we have r_h ≃ 0.8r_v and the one-dimensional velocity dispersion, σ, then obeys

$$\begin{eqnarray} \sigma ^2 = \frac{1}{3}\langle v^2\rangle =\frac{1}{3}\frac{GM}{2r_\mathrm{ v}}\simeq \frac{2}{15}\frac{G N m}{r_\mathrm{ h}}. \end{eqnarray}$$

(1)

where G is the gravitational constant, m is the mass of the stars, and 〈v²〉 is their mean-square velocity. The central dispersion can be expressed in terms of σ by defining a numerical coefficient α₁

$$\begin{eqnarray} \sigma _\mathrm{ c}^2=\alpha _1\sigma ^2. \end{eqnarray}$$

(2)

Next, we define the relation between r_h and the core radius, r_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 r_h needs to be provided by the core. The exact form of the r_h/r_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.11N), Giersz & Heggie 1994a] in the denominator of the central relaxation time-scale, and the dimensionless ratio of the central potential (|ϕ₀|) to the central velocity dispersion, |$|\phi _0|/\sigma _\mathrm{ c}^2$|⁠, all are constant, they obtain

$$\begin{eqnarray} r_\mathrm{ h}/r_\mathrm{ c} \simeq \alpha _2 N^{2/3}, \end{eqnarray}$$

(3)

where α₂ is a constant. If the assumptions above are relaxed, the final ratio would include a dependence on |$(|\phi _0|/\sigma _\mathrm{ c}^2)/\ln \Lambda$| but, since both these quantities scale weakly with N, we are justified in taking α₂ as a constant. In Section 4, we determine α₂ and a possible additional N-dependence from N-body simulations.

We complement the above relations with the relation for the number of stars in the cluster core, N_c. Starting from the definition of the core radius (Heggie & Hut 2003, p. 71)

$$\begin{eqnarray} \frac{4\pi G}{9}\rho _0r_\mathrm{ c}^2 = \sigma _\mathrm{ c}^2 \end{eqnarray}$$

(4)

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 ρ₀ ≃ 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, N_c, is²

$$\begin{eqnarray} N_\mathrm{ c}=0.98 \frac{\alpha _1}{\alpha _2} \left(\frac{N}{10^2}\right)^{1/3} \end{eqnarray}$$

(5)

The α_i parameters in this theoretical prediction are of order unity and their values and possible N-dependences will be determined in Section 4 by fits to N-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 x = Gm²/(2a) is greater than |$x_\mathrm{ h}\simeq m\sigma _\mathrm{ c}^2$| (where a 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 x < x_h, and a > a_h = Gm²/(2x_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 x is on average Δ ≃ 0.4 for resonant encounters³ (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 x, at some moment x 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 |$x_{ej}=\frac{6}{\Delta }m \left| \phi _0 \right|= 15 m \left| \phi _0 \right|$| (Goodman 1984), which corresponds to a semimajor axis of a_ej = Gm²/(2x_ej). In order to obtain the total number of binaries, we will compute the distribution of (hard) binaries per unit volume and per unit |$z\equiv x/(m \sigma _\mathrm{ c}^2)$| and integrate it from z_h = 1 to |$z_{ej}=15 \left| \phi _0 \right|/\sigma _\mathrm{ c}^2$|⁠, where |$\left| \phi _0 \right|/\sigma _\mathrm{ c}^2$| is a dimensionless central potential, which is a measure of the (instantaneous) central concentration⁴ 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)

$$\begin{eqnarray} \Gamma _\mathrm{ b}=0.75\frac{G^5m^5 n_\mathrm{ c}^3}{\sigma _\mathrm{ c}^9}, \end{eqnarray}$$

(6)

and the hardening rate (Heggie & Hut 1993) defined as the rate of energy generation of the binaries

$$\begin{eqnarray} \dot{z}=3.8\frac{G^2 m^2 n_\mathrm{ c}}{\sigma ^3_\mathrm{ c}}. \end{eqnarray}$$

(7)

Since this hardening rate is independent of z itself, the steady post-collapse distribution of hard binaries is uniform in z between the minimum and maximum binding energy. The distribution per unit volume, f(z), is

$$\begin{eqnarray} f(z)=\frac{\Gamma _\mathrm{ b}}{\dot{z}}=0.20 \frac{G^3m^3\mathrm{ n}_c^2}{\sigma _\mathrm{ c}^6}. \end{eqnarray}$$

(8)

Assuming that z_ej ≫ z_h and that the relevant binary physics happens within the volume V_c of the cluster core, the total number of binaries is⁵

$$\begin{eqnarray} N_\mathrm{ b}=\int ^{z_{ej}}_{z_\mathrm{ h}} f(z) V_\mathrm{ c} \mathrm{d}z\simeq 3.0\frac{G^3m^3n_\mathrm{ c}^2}{\sigma _\mathrm{ c}^8} V_\mathrm{ c} \left| \phi _0 \right|. \end{eqnarray}$$

(9)

The above equation can be simplified by using the scaling laws in Table 1 to obtain⁶

$$\begin{eqnarray} N_\mathrm{ b} = 2.8\frac{\alpha _2}{\alpha _1} \frac{\left| \phi _0 \right|}{\sigma _\mathrm{ c}^2}\left(\frac{N}{10^2}\right)^{-1/3}. \end{eqnarray}$$

(10)

The ratio |$\left| \phi _0 \right|/\sigma _\mathrm{ c}^2$| is dimensionless and thus only depends on N. 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 N_b ∝ N^−1/3.

Table 1.

Scaling laws for single-mass cluster properties in the steady post core-collapse evolutionary phase. The coefficients α₁ and α₂ are found from the N-body simulations.

Value	Scaling law	Coefficients
Velocity dispersion in the core	|$\sigma _\mathrm{ c}^2=\alpha _1\frac{2}{15}\frac{GmN}{r_\mathrm{ h}}$|	α1 ≃ 1.4
Half-mass to core radius ratio	|$\frac{r_\mathrm{ h}}{r_\mathrm{ c}}=\alpha _2 N^{2/3}$|	α2 ≃ 0.13
Net creation rate of binaries	|$\Gamma _\mathrm{ b}=0.75\frac{G^5m^5 n_\mathrm{ c}^3}{\sigma _\mathrm{ c}^9}$|
Binary hardening rate	|$\dot{z}=3.8\frac{G^2 m^2 n_\mathrm{ c}}{\sigma ^3_\mathrm{ c}}$|
Core radius	|$\frac{4\pi G}{9}\rho _0r_\mathrm{ c}^2 = \sigma _\mathrm{ c}^2$|

2.2 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 |$m_\mathrm{BH}\simeq 20\, {\rm M}_\odot$| and that of light stars with |$m_\star \simeq 0.5\, {\rm M}_\odot$|⁠. Due to the large difference between these bins, we will assume that all BHs have the same mass m_BH and all stars have the same mass m_⋆ ≪ m_BH. Each of these mass bins contribute a total of M_BH and M_⋆ 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, μ ≡ M_BH/M, which sets what fraction of the cluster mass is in BHs, and the stellar-mass ratio, q ≡ 〈m〉/m_BH, which sets how massive the BHs are with respect to the average mass of stars and BHs, 〈m〉, in the cluster. Both of these quantities are defined such that μ, q ∈ (0, 1). Some authors (for example, Breen & Heggie 2013) give alternative definitions of these parameters with respect to the mass of stars, M_⋆ and m_⋆, 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⁷, M_⋆ ≫ M_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 cannot be achieved, which applies if (Spitzer 1969)

$$\begin{eqnarray} \frac{M_\mathrm{BH}}{M_\star }\left(\frac{m_\mathrm{BH}}{m_\star }\right)^{3/2}\gt 0.16. \end{eqnarray}$$

(11)

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 summarized 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

$$\begin{eqnarray} \sigma _{\mathrm{ c}, \mathrm{BH}}^2=\alpha _{1, \mathrm{BH}} \sigma ^2_\mathrm{BH}= \alpha _{1, \mathrm{BH}}\frac{2}{15} \frac{G M_\mathrm{BH}}{r_{\mathrm{ h}, \mathrm{BH}}} \end{eqnarray}$$

(12)

$$\begin{eqnarray} \frac{r_{\mathrm{ h}, \mathrm{BH}}}{r_{\mathrm{ c}, \mathrm{BH}}}=\alpha _{2, \mathrm{BH}} N_\mathrm{BH}^{2/3} \end{eqnarray}$$

(13)

$$\begin{eqnarray} \frac{4 \pi }{9} G \rho _{0, \mathrm{BH}} r_{\mathrm{ c}, \mathrm{BH}}^2 = \sigma _{\mathrm{ c}, \mathrm{BH}}^2 \end{eqnarray}$$

(14)

$$\begin{eqnarray} N_{\mathrm{ c}, \mathrm{BH}}=0.98 \frac{\alpha _{1,\mathrm{BH}}}{\alpha _{2,\mathrm{BH}}} \left(\frac{N_\mathrm{BH}}{10^2}\right)^{1/3} \end{eqnarray}$$

(15)

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 N_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. m↦m_BH. This recovers the same form of the two-mass binary binding energy distribution from Retterer (1980), i.e.

$$\begin{eqnarray} N_{\mathrm{BBH}}\simeq 3.0\frac{G^3m_\mathrm{BH}^3n_{\mathrm{ c}, \mathrm{BH}}^2}{\sigma _{\mathrm{ c}, \mathrm{BH}}^8} V_{\mathrm{ c}, \mathrm{BH}} \left| \phi _0 \right| \end{eqnarray}$$

(16)

which we can express as

$$\begin{eqnarray} N_{\mathrm{BBH}} = 2.8\frac{\alpha _{2,\mathrm{BH}}}{\alpha _{1,\mathrm{BH}}} \frac{\left| \phi _0 \right|}{\sigma _{\mathrm{ c}, \mathrm{BH}}^2} \left(\frac{N_\mathrm{BH}}{10^2}\right)^{-1/3} \end{eqnarray}$$

(17)

In parallel to the result for the single-mass case, we obtain the prediction that the number of BBHs scales as |$N_{\mathrm{BBH}}\propto N_\mathrm{BH}^{-1/3}$|⁠, although in this case the contribution of the |$\left| \phi _0 \right|/\sigma _{\mathrm{ c}, \mathrm{BH}}^2$| term is larger. In the following section we will test the above equations against N-body models.

Table 2.

Scaling laws for two-mass cluster properties in the steady post core-collapse evolutionary phase. The coefficients α_1,BBH and α_2,BBH are found from the N-body simulations.

Value	Scaling law	Coefficients
Velocity dispersion in the core	|$\sigma _{\mathrm{ c}, \mathrm{BH}}^2=\alpha _{1, \mathrm{BH}}\frac{2}{15} \frac{G M_\mathrm{BH}}{r_{\mathrm{ h}, \mathrm{BH}}}$|	α1, BH ≃ 1.5
Half-mass to core radius ratio	|$\frac{r_{\mathrm{ h}, \mathrm{BH}}}{r_{\mathrm{ c}, \mathrm{BH}}}=\alpha _{2, \mathrm{BH}} N_\mathrm{BH}^{2/3}$|	α2, BH ≃ 0.088
Net creation rate of binaries	|$\Gamma _\mathrm{ b}=0.75\frac{G^5m_\mathrm{BH}^5 n_{\mathrm{ c}, \mathrm{BH}}^3}{\sigma _{\mathrm{ c}, \mathrm{BH}}^9}$|
Binary hardening rate	|$\dot{z}=3.8\frac{G^2 m_\mathrm{BH}^2 n_{\mathrm{ c}, \mathrm{BH}}}{\sigma ^3_{\mathrm{ c}, \mathrm{BH}}}$|
Core radius	|$\frac{4 \pi G}{9} \rho _{0, \mathrm{BH}} r_{\mathrm{ c}, \mathrm{BH}}^2 = \sigma _{\mathrm{ c}, \mathrm{BH}}^2$|

3 DESCRIPTION OF THE N-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. G = M₀ = −4E₀ = 1 (Hénon 1971), where M₀ is the initial cluster mass and E₀ is the initial energy of the cluster. For our assumption of virial equilibrium, the initial virial radius is |$r_{\mathrm{ v},0}=-GM_0^2/(2W_0)=1$|⁠, where W₀ = −1/2 is the gravitational energy. The corresponding unit of time is |$\tau _{\mathrm{ dyn}}=(r_{\mathrm{ v},0}^3/GM_0)^{1/2}$|⁠.

We use the direct N-body code petar (Wang et al. 2020b), which is a high-performance code used for the modelling of large-N 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, Nitadori & Makino 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 r_c are defined following Casertano & Hut (1985) and McMillan, Hut & Makino (1990), with r_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 |$m\sigma _\mathrm{ c}^2$| – which is proportional to the average kinetic energy of single stars in the core – and that are bound to the cluster.

3.1 Single-mass models

We have run a family of single-mass models with logarithmically spaced N, in the range 128 ≤ N ≤ 8192, where each model is run several times for statistical significance. The input N and number of runs are shown in Table 3. Within the aforementioned simplifications, single-mass models are scale-free, and thus N 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 time-scale (Spitzer & Hart 1971),

$$\begin{eqnarray} \tau _{\mathrm{ rh}}=0.138 \frac{N}{\psi \ln \Lambda } \sqrt{\frac{r_\mathrm{ h}^3}{GM}}, \end{eqnarray}$$

(18)

evaluated at the start of the simulation, where ψ = 1 and ln Λ is the Coulomb logarithm with Λ ≃ 0.11N for single-mass models (Giersz & Heggie 1994a). The data for the post-collapse evolution is then taken starting from t₀ = 1.1τ_cc up until t_f = t₀ + 20τ_rh,0, at constant intervals of 4 N-body times. Values found for the cluster properties are averaged over all snapshots between the times t₀ and t_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.

Table 3.

Values of the input parameters in the N-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.

N	tf/τdyn	Num. of runs
128	143	64
256	227	32
512	376	16
1024	642	8
2048	1120	4
4096	1986	4
8192	3567	4

3.2 Two-mass models

For the two-mass N-body models, we take our data from t₀ = 1.1τ_cc to t_f = t₀ + 2τ_rh,0. We use a lower number of τ_rh,0 than in the single-mass models because in the two-mass models the binary phenomena happen within the BH subcluster and thus the relevant time-scale is shorter. The impact of the adopted simulation time is discussed in Section 6.1. For the core-collapse time-scale of two-mass clusters, we use the result of Kim & Lee (1997)

$$\begin{eqnarray} \tau _{\mathrm{ cc}}=5.3 \left(\frac{m_\star }{m_\mathrm{BH}}\right)^2\left(\frac{N_\star }{N_\mathrm{BH}}\right)^{1/2}\tau _{\mathrm{ rh},0}, \end{eqnarray}$$

(19)

as well as their expression for τ_rh,0 (equation 18, but setting Λ = 0.4N). The values of the three parameters N_BH, q, and μ for the set of simulations that we have run is summarized in Table 4. We set N ∈ [5 × 10⁴, 1 × 10⁵, 2 × 10⁵] and store the snapshots every eight N-body times. For the mass ratio q we adopt q ∈ [1/50, 1/25], to approximate metal-poor (log₁₀(Z/Z_⊙) ≃ −1.5) and metal-rich (log₁₀(Z/Z_⊙) ≃ −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 Spitzer instability (Spitzer 1969), such that the evolution should be well-described by the theory of Breen & Heggie (2013).

Table 4.

Values of the input parameters in the N-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.

N	q	μ	NBH	Num. of runs
5 × 104	1/50	0.025	25	4
5 × 104	1/25	0.025	50	4
5 × 104	1/50	0.05	50	4
5 × 104	1/25	0.05	100	4
5 × 104	1/50	0.1	100	4
5 × 104	1/25	0.1	200	4
1 × 105	1/50	0.025	50	2
1 × 105	1/25	0.025	100	2
1 × 105	1/50	0.05	100	2
1 × 105	1/25	0.05	200	2
1 × 105	1/50	0.1	200	2
1 × 105	1/25	0.1	400	2
2 × 105	1/50	0.025	100	1
2 × 105	1/25	0.025	200	1
2 × 105	1/50	0.05	200	1
2 × 105	1/25	0.05	400	1
2 × 105	1/50	0.1	400	1
2 × 105	1/25	0.1	800	1

4 COMPARING THE STEADY POST-COLLAPSE MODEL TO N-BODY SIMULATIONS

4.1 Single-mass clusters

In this section, we calibrate and compare the theoretical predictions of the single-mass model of Section 2.1 to the N-body models of Section 3.1. We also allow the theory to include slight additional N-dependences in the α₁, α₂ parameters (that we label as |$\alpha ^{\prime }_1$| and |$\alpha ^{\prime }_2$|⁠, respectively). In Fig. 1, we show the ratio |$\alpha _1=\sigma _\mathrm{ c}^2/\sigma ^2$| as a function of N. We find that a constant α₁ equation (2) is a good approximation only for N ≲ 10³, and |$\alpha ^{\prime }_1$| increases slightly at larger N. A scaling |$\alpha ^{\prime }_1 = 1.1 (N/10^2)^{0.06}$| describes the N-body data well. This increase is likely due to the increase of r_h/r_c with N, leading to a higher central density and dispersion, relative to the global values, for clusters with larger N.

Figure 1.

Ratio |$\sigma _\mathrm{ c}^2/\sigma ^2$| (equations 1 and 2) as a function of N for the single-mass N-body models (crosses with error bars). Fits for a constant α₁ and an N-dependent |$\alpha ^{\prime }_1$| are shown with grey and red lines, respectively.

The r_h/r_c ratio is shown Fig. 2. The scaling of equation (3) is a good approximation for constant α₂ for N ≃ 10³–10⁴, whereas it deviates slightly at lower values of N, which, as discussed in Section 2.1, may be due to the weak N-dependence of the Coulomb logarithm in the denominator of the central relaxation time-scale and in the |$\left| \phi _0 \right|/\sigma _\mathrm{ c}^2$| term. A scaling of |$\alpha ^{\prime }_2 = 0.19 (N/10^2)^{-0.09}$| provides a good fit to the N-body data.

Figure 2.

Ratio r_h/r_c equation (3) as a function of N for the single-mass N-body models (crosses with error bars). Fits for a constant α₂ and an N-dependent |$\alpha ^{\prime }_2$| are shown with grey and red lines, respectively.

From these scaling laws, we are able to accurately predict the number of stars in the core (Fig. 3). The prediction for N_c with constant α_i equation (5) has the same form than the one given in Goodman (1984),

$$\begin{eqnarray} N_\mathrm{ c}=10 \left(\frac{\alpha _1}{1.3}\right) \left(\frac{\alpha _2}{0.13}\right)^{-1} \left(\frac{N}{10^2}\right)^{1/3}, \end{eqnarray}$$

(20)

which yields a rough estimate of N_c from the N-body models within <50 per cent over two orders of magnitude in N. Using the N-dependent expressions for |$\alpha ^{\prime }_1$| and |$\alpha ^{\prime }_2$| yields

$$\begin{eqnarray} N_\mathrm{ c}^{\prime }=0.58\left(\frac{N}{10^2}\right)^{0.17} N_\mathrm{ c} \end{eqnarray}$$

(21)

which accurately describes the N-body data. Furthermore, in Fig. 4 we show the ratio |$|\phi _0|/\sigma _\mathrm{ c}^2$| as a function of N, which can be approximated by the power-law relation

$$\begin{eqnarray} \frac{\left| \phi _0 \right|}{\sigma _\mathrm{ c}^2}\simeq 7.2 \left(\frac{N}{10^2}\right)^{0.1}. \end{eqnarray}$$

(22)

Now all ingredients of the expression for the number of binaries have been discussed, it is time to combine them and compare to N_b we find in the simulations. We obtain

$$\begin{eqnarray} N_\mathrm{ b} = 1.9\left(\frac{\alpha _2}{0.13}\right)\left(\frac{\alpha _1}{1.3}\right)^{-1}\left(\frac{\left| \phi _0 \right|/\sigma _\mathrm{ c}^2}{7.2}\right) \left(\frac{N}{10^2}\right)^{-1/3} \end{eqnarray}$$

(23)

Which can be corrected by the deviations from the scaling laws above to yield

$$\begin{eqnarray} N^{\prime }_\mathrm{ b} = 1.7 \left(\frac{N}{10^2}\right)^{-0.17} N_\mathrm{ b}, \end{eqnarray}$$

(24)

such that a cluster with ∼100 stars is expected to have ∼3 binaries.

Figure 3.

Number of stars in the cluster core as a function of N for the single-mass N-body models (crosses with error bars). Using α_i independent of N (grey) we are able to roughly match the N-body models, whereas using the N-dependant |$\alpha ^{\prime }_{\it i}$| coefficients (red) yields a very accurate match.

Figure 4.

Ratio |$\left| \phi _0 \right|/\sigma _\mathrm{ c}^2$| equation (22) as a function of N for the single-mass N-body models (crosses with error bars).

In Fig. 5, we show both the prediction for N_b and the N-body results. The prediction equation (24) approximately reproduces the number of binaries in the core at large N (N ≳ 10³), whereas in the smallest clusters N_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 r_h, the number of binaries in the N-body models is higher, and increases slightly with N, that is, the opposite N-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.

Figure 5.

Number of binaries as a function of N for the single-mass N-body models (crosses with error bars). The theory over-predicts the number of binaries.

4.2 Two-mass clusters

In this section, we compare the theoretical prediction of the two-mass model of Section 2.2 to the N-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 N_BH is preferred, equivalent to the single-mass case. For the ratio r_h,BH/r_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 N_BH can be attributed to the weak N_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.

Figure 6.

Ratio |$\sigma _{\mathrm{ c}, \mathrm{BH}}^2/\sigma _\mathrm{BH}^2$| equation (12) as a function of N_BH for the two-mass N-body models (crosses with error bars). Fits for a constant α_1,BH and an N-dependent |$\alpha ^{\prime }_{1, \mathrm{BH}}$| are shown with grey and red lines, respectively.

Figure 7.

Ratio r_h,BH/r_c,BH equation (13) as a function of N_BH for the two-mass N-body models (crosses with error bars). Fits for a constant α_2,BH and an N-dependent |$\alpha ^{\prime }_{2, \mathrm{BH}}$| are shown with grey and red lines, respectively.

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 N_c,BH,

$$\begin{eqnarray} N_{\mathrm{ c}, \mathrm{BH}} = 16 \left(\frac{\alpha _{1, \mathrm{BH}}}{1.5}\right) \left(\frac{\alpha _{2, \mathrm{BH}}}{0.088}\right)^{-1} \left(\frac{N_\mathrm{BH}}{10^2}\right)^{1/3} \end{eqnarray}$$

(25)

with the discrepancy being largest at lower N_BH. Using the corrected |$\alpha ^{\prime }_{{\it i}, \mathrm{BH}}$| yields

$$\begin{eqnarray} N^{\prime }_{\mathrm{ c}, \mathrm{BH}} = 0.76\left(\frac{N_\mathrm{BH}}{10^2}\right)^{0.23} N_{\mathrm{ c}, \mathrm{BH}}, \end{eqnarray}$$

(26)

which accurately reproduces the results from the N-body models for N_BH ≳ 50, although it slightly overpredicts the results for the sub-clusters with the fewest black holes in the core (N_c,BH ≲ 7).

Figure 8.

Number of BHs in the sub-cluster core as a function of N_BH for the two-mass N-body models (crosses with error bars). Using α_i,BH independent of N_BH (grey) we are able to roughly match the N-body models, whereas using the N_BH-dependant |$\alpha ^{\prime }_{{\it i}, \mathrm{BH}}$| coefficients (red) yields a good match.

The above scaling laws are complemented by a fit to |$\left| \phi _0 \right|/\sigma _{\mathrm{ c}, \mathrm{BH}}^2$| (Fig. 9). From the three parameters of the model, we find that this value only correlates strongly with μ. The best-fitting power law is

$$\begin{eqnarray} \frac{\left| \phi _0 \right|}{\sigma _{\mathrm{ c}, \mathrm{BH}}^2}\simeq 35 \left(\frac{\mu }{0.025}\right)^{-0.62}. \end{eqnarray}$$

(27)

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 of the BHs. Using Breen & Heggie (2013, equation 4), we can write it in terms of the central (total) velocity dispersion, yielding a much smaller value

$$\begin{eqnarray} \frac{\left| \phi _0 \right|}{\sigma _\mathrm{ c}^2}\simeq 4 \left(\frac{\mu }{0.025}\right)^{-0.22} \left(\frac{q}{1/25}\right)^{0.4}. \end{eqnarray}$$

(28)

Using the calibrations for the scaling laws (equations 12–15), we can evaluate the prediction for the number of BBHs of equation (17) to obtain

$$\begin{eqnarray} N_{\mathrm{BBH}} &=& 6 \left(\frac{\alpha _{2, \mathrm{BH}}}{0.088}\right)\left(\frac{\alpha _{1, \mathrm{BH}}}{1.5}\right)^{-1}\left(\frac{\left| \phi _0 \right|/\sigma _{\mathrm{ c}, \mathrm{BH}}^2}{35}\right)\\ && \times \left(\frac{\mu }{0.025}\right)^{0.62} \left(\frac{N_\mathrm{BH}}{10^2}\right)^{-1/3}. \end{eqnarray}$$

(29)

This can be corrected with the deviations from the scaling laws as

$$\begin{eqnarray} N^{\prime }_{\mathrm{BBH}}=1.3 \left(\frac{N_\mathrm{BH}}{10^2}\right)^{-0.23} N_{\mathrm{BBH}}. \end{eqnarray}$$

(30)

Although the two-mass steady post-collapse model is able to describe the other properties of the N-body model, it over-predicts the number of BBHs (Fig. 10). This result is very similar to the single-mass case, although in this case the discrepancy is larger. In general, the discrepancy is most notable in models with smaller μ and smaller N_BH.

Figure 9.

Ratio |$\left| \phi _0 \right|/\sigma _{\mathrm{ c}, \mathrm{BH}}^2$| equation (27) as a function of μ for the two-mass N-body models (crosses with error bars).

Figure 10.

Number of BBHs as a function of N_BH and μ for the two-mass N-body models (crosses with error bars). The number of BBHs is roughly constant and the theory vastly overpredicts the number of binaries.

5 THE CAUSE OF THE DEARTH OF BINARIES: BINARY–BINARY INTERACTION

In the previous sections, we find that the theoretical predictions for the population of binaries are unable to describe the number of three-body binaries in the N-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 N-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) assume that the binding energies of binaries evolve because of interactions between 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 ionized (⁠|$\sim 95~{{\ \rm per\ cent}}$| 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 N-body models is |$\mathcal {O}(1)$|⁠, ignoring interactions among binaries is a justified assumption for clusters with N_c ≫ 1. However, as we have shown in Fig. 3, N_c ≲ 20 for N ≲ 10³, whilst the theory predicts N_b ≳ 2. In fact, the predicted binary fraction in the core is f_b = N_b/N_c ≃ 0.6(N/10²)^−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 binary–single interactions. The interaction cross-section for a binary with a single goes as |$\Sigma _{\mathrm{ bs}}\propto G r_\mathrm{ p} m_{\mathrm{ tot}}/v_{\mathrm{ bs}}^2$|⁠, where v_bs is the relative velocity, r_p ≃ a is the relevant minimum distance (with a the semimajor axis of the binary), and m_tot = 3m is the total mass of the interacting binary and single. Assuming equipartition, the typical relative velocity among binaries, v_bb, is a factor of |$\sqrt{3/4}$| lower than among binaries and singles (⁠|$v_{\mathrm{ bb}}=\sqrt{3/4}v_{\mathrm{ bs}}$|⁠). Then, taking m_tot = 4m and r_p = 2a 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

$$\begin{eqnarray} \frac{\Gamma _{\mathrm{ bb}}}{\Gamma _{\mathrm{ bs}}}\simeq \frac{N_\mathrm{ b}^2\Sigma _{\mathrm{ bb}}v_{\mathrm{ bb}}}{N_\mathrm{ b}N_\mathrm{ s}\Sigma _{\mathrm{ bs}}v_{\mathrm{ bs}}}\simeq 3.1 f_\mathrm{ b} \end{eqnarray}$$

(31)

This estimate shows that binary–binary interactions are expected to be more frequent than binary–single interactions for⁸ f_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 three binaries are created. This is approximately the number of hard binaries predicted by Goodman & Hut (1993) for clusters N ≃ 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 ionization 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 N-body models.

5.1 Single-mass models

As stated above, the most likely outcome of a binary–binary interaction is that at least one of the binaries is ionized. If the time-scale for this ionization mechanism is shorter than the time needed to form new binaries, then the predicted uniform binding energy distribution equation (8) is never realized. Instead, one would find a decline in binaries at higher z, because they are destroyed at some point during their life-cycle. Furthermore, this deviation should be larger at smaller N, 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 realized, supporting the idea that binaries are destroyed before they reach the maximum binding energy.

Figure 11.

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

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 N-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 N_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.

Figure 12.

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

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

$$\begin{eqnarray} \Sigma (a_1, a_2, v) = \pi \left(a_1^2+a_2^2\right)\left(\frac{\hat{v}^2}{0.567}+\frac{\hat{v}^6}{0.0279}\right)^{-1} \end{eqnarray}$$

(32)

The cross-section depends on the semimajor axes a₁, a₂ of the incoming binaries and the ratio |$\hat{v}$| of their relative velocity to their critical velocity, |$\hat{v}=v/\sqrt{Gm(a_1+a_2)/(a_1a_2)}$|⁠. The total rate of 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 z 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 |$\sqrt{2}\sigma _\mathrm{ 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 (⁠|$f_{\mathrm{ Maxwell}}(v)\propto v^2/\sigma _\mathrm{ c}^3 \exp (- 0.5v^2/\sigma _\mathrm{ c}^2$|⁠)). Therefore, the triple formation rate is

$$\begin{eqnarray} \Gamma _{\mathrm{ t}, \mathrm{ bb}} = \iiint \frac{f(z_1)f(z_2)}{2}f_{\mathrm{ Maxwell}}(v) \Sigma (z_1, z_2, v) v \mathrm{d}z_1\mathrm{d}z_2 \mathrm{d}v \end{eqnarray}$$

(33)

The ratio Γ_t,bb/Γ_b is then expressed only as a function of N by using the scaling relations in Table 1 and shown in Fig. 13. The ratio is greater than one for N ≲ 2000 and declining approximately as N^−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 three-body binary in N-body simulations, as every time a new binary forms it is rapidly destroyed in a binary–binary interaction.

Figure 13.

Ratio of the triple formation rate Γ_t,bb to the binary formation rate Γ_b. This ratio is ≫1 for N ≲ 10³, which confirms the importance of binary–binary interactions in that N regime. The ratio scales as Γ_t,bb/Γ_b ∝ N^−0.8.

5.2 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. 14 that a uniform BBH binding energy distribution is not realized, with the discrepancy being greatest at lower N_BH and higher z. Contrary to the single-mass case, here the models with lower N_BH have higher z_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. 15 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 N_BH correspond to the smaller end of the range of N in the single-mass case, so we do not observe the decrease in the number of triples that is seen at large N in the single-mass models. The presence of stable triples in N-body models without primordial binaries had already been noticed before, e.g. Aarseth (2012) and 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.

Figure 14.

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

Figure 15.

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

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 N_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.

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 N_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 |$\Gamma _{\mathrm{ t, bb}}/\Gamma _\mathrm{ b}\propto \mu ^{-1.2} N_{\mathrm{BH}}^{-1.2}$|⁠.

6 DISCUSSION AND CONCLUSIONS

6.1 Validity of the assumptions

In order for the binaries to populate the binary binding energy distribution in the absence of binary–binary interactions (Figs 11 and 14), the properties of the cluster must evolve slowly with time. The criterion for this is that the time-scale at which the cluster expands, τ_exp, is larger than the binary lifecycle, |$\tau _{\mathrm{ bin}}\simeq (z_{ej}-z_\mathrm{ h}) /\dot{z}\simeq z_{ej}/\dot{z}$|⁠. We can define the time-scale for expansion similarly to τ_bin by assuming that |$\tau _{\mathrm{ exp}}\simeq r_\mathrm{ h}/\dot{r}_\mathrm{ h}$|⁠. Using Hénon’s principle (Hénon 1975), we can relate the energy production of the core to its global properties as

$$\begin{eqnarray} \dot{E} \simeq \zeta \frac{|E|}{\tau _{\mathrm{ rh}}} \end{eqnarray}$$

(34)

with ζ ≃ 0.1 (Hénon 1965; Alexander & Gieles 2012). If we assume that the cluster energy scales as E ∝ GM²/r_h and that the mass-loss is negligible, we obtain

$$\begin{eqnarray} \tau _{\mathrm{ exp}}=\frac{\left| E \right|}{\dot{E}_\mathrm{ h}}\simeq \frac{\tau _{\mathrm{ rh}}}{\zeta }. \end{eqnarray}$$

(35)

This can be compared to the binary lifecycle, which we find from the scaling laws of Section 4.1

$$\begin{eqnarray} \tau _{\mathrm{ bin}}&\simeq \frac{z_{ej}}{\dot{z}}\simeq 300 \left(\frac{\psi \log \Lambda }{10}\right)\left(\frac{N}{10^2}\right)^{-1.2} \tau _{\mathrm{ rh}}, \\ &\simeq 300 \left(\frac{N}{10^2}\right)^{-0.2} \tau _{\mathrm{ dyn}}. \end{eqnarray}$$

(36)

We also give the value of this time-scale in terms of the dynamical time τ_dyn by using that τ_rh/τ_dyn ≃ 0.1N/(ψlog Λ) (Heggie & Hut 2003). The binary time-scale is approximately a constant number of N-body times, weakly decreasing with N. This somewhat counter-intuitive result is a consequence of the larger central concentration of high-N models, leading to a faster evolution of the binaries. Alternatively, if we had considered homologous evolution (r_c ∝ r_h and N_c ∝ N), this time-scale would increase strongly with N, τ_bin ∝ N²τ_dyn.

For the single-mass case, we find τ_exp > τ_bin for N ≳ 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)

$$\begin{eqnarray} f(z)\propto \frac{m^3n_\mathrm{ c}^2}{\sigma _\mathrm{ c}^6}\propto \frac{r_\mathrm{ h}^{3}N_\mathrm{ c}^2}{N^3r_\mathrm{ c}^{6}}\propto N^{5/3}r_\mathrm{ h}^{-3} \end{eqnarray}$$

(37)

This implies that, for a fixed N, 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 z, because these were formed earlier when f(z) was higher. We actually see a lack of binaries at large z, making the absence of binaries even more important for N ≲ 800.

Another time-scale to consider is the duration of our simulations, t_f − t₀. If this time-scale 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 t_f − t₀ = 20τ_rh,0. Because the cluster expands, the actual number of elapsed relaxation times is less than 20. We can estimate the instantaneous relaxation time-scale by assuming that (Hénon 1965)

$$\begin{eqnarray} \tau _{\mathrm{ rh}}(t)=\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\tau _{\mathrm{ rh}, 0}&t\lt \tau _{\mathrm{ cc}}\\ (t/\tau _{\mathrm{ cc}})\tau _{\mathrm{ rh}, 0}&t\gt \tau _{\mathrm{ cc}} \end{array}\right. \end{eqnarray}$$

(38)

So we obtain

$$\begin{eqnarray} t_\mathrm{ f}-t_0=\langle \tau _{\mathrm{ rh}}\rangle \int ^{t_\mathrm{ f}}_{t_0}\frac{\mathrm{d}t}{\tau _{\mathrm{ rh}}(t)}=12 \langle \tau _{\mathrm{ rh}}\rangle \end{eqnarray}$$

(39)

which we can compare to equation (36) to show that t_f − t₀ > τ_bin for N ≳ 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)

$$\begin{eqnarray} \frac{t_\mathrm{ f}-t_0}{\tau _{\mathrm{ bin}}} \simeq 0.3 \left(\frac{N_{\mathrm{ BH}}}{10^2}\right)^{0.8}\left(\frac{q}{1/50}\right)^{1.1} \left(\frac{\mu }{0.025}\right)^{-1.3} \end{eqnarray}$$

(40)

which is larger than unity for models with N_BH ≳ 150 (450) for q = 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 N_BH in Fig. 15 as in the single-mass case. To check whether there is evolution of the binary energy distribution, we compared f(z) in the first half to f(z) in the second half of the simulation interval. The two distributions are statistically indistinguishable for all N_BH, suggesting that the simulation time is long enough for f(z) 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-N 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 N-body simulation time-scale and thus convergence is attained by combining multiple snapshots from the few runs.

6.2 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 a_h and multiply that rate by the probability that a third star is present within the volume, which is roughly |$n_\mathrm{ c}a_\mathrm{ h}^3$|⁠. We can estimate the rate of triple formation due to binary–single–single interactions, Γ_t,bss, by a similar argument, although performing the change n_c↦n_b and a_h↦3a_h (so that the outer orbit of the triple is about ∼3a_h). Then, the ratio of the rate of triple production over binary formation is roughly |$\Gamma _{\mathrm{ t}, \mathrm{ bss}}/\Gamma _{\mathrm{ b}}\sim n_\mathrm{ b} (3a_\mathrm{ h})^3/(n_\mathrm{ c}a_\mathrm{ h}^3)=27/N_\mathrm{ c}\simeq 3(N/10^2)^{-1/3}$|⁠, where we use the scaling N_c(N) from equation (20) and, as seen in the N-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 3a_h. The formation rate of ‘hard’ triples is a factor of 27 lower. The N-dependence in Γ_t,bb/Γ_b is steeper (N^−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-N 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 N replaced by N_BH.

6.3 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, Heggie & Zhao (2011) that the value of ζ in tidally limited clusters and isolated clusters is equal to within 20 per cent, 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 r_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 per cent 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 per cent (Dickson et al. 2023).

6.4 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 time-scale 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 (a ∼ a_h), which requires a nearly radial orbit (e ∼ 1) to trigger a GW capture. We can use the capture criterion by Samsing (2018) to quantify what the eccentricity, e, must be

$$\begin{eqnarray} 1-e\lesssim \left(\frac{4.6 G m}{a c^2}\right)^{5/7}. \end{eqnarray}$$

(41)

To put it in context, for a 10⁵ M_⊙ cluster with r_h = 1pc, a capture between two 10 M_⊙ BHs would start at ∼60 au and 1 − e ∼ 10⁻⁶. As the binary inspirals, this eccentricity may be radiated away. To estimate this, we find the peak GW frequency f from Wen (2003),

$$\begin{eqnarray} f = \frac{1}{\pi }\sqrt{\frac{2Gm_\mathrm{BH}}{a^3}}\frac{(1+e)^{1.1954}}{(1-e^2)^{3/2}}, \end{eqnarray}$$

(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 (f ≃ 10 Hz). What we find is that, due to the large initial value of a, the eccentricity is radiated away before the GW emission becomes observable and thus the waveform appears nearly circular. For the above values of M, r_h, a, and m_BH, we find a merger that has circularized to e_10 Hz ∼ 10⁻², 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 |$\sim 63~{{\ \rm per\ cent}}$| (⁠|$\sim 0~{{\ \rm per\ cent}}$|⁠) of captures that start with an initial semimajor axis of a_h (a_ej) have a remaining e_10 Hz < 0.1 at 10 Hz. The sampling of e has a chance of producing mergers that form in the LIGO-Virgo-KAGRA band, f > 10Hz, with a probability of |$\sim 12~{{\ \rm per\ cent}}$| (⁠|$\sim 46~{{\ \rm per\ cent}}$|⁠).

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 Kozai 1962; Lidov 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 e_10 Hz ∼ 10⁻² 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 time-scale 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 N-body models of low-mass clusters (⁠|$M\lesssim 10^5\, {\rm M}_\odot$|⁠, Rastello et al. 2019; Banerjee 2021; Barber, Chattopadhyay & Antonini 2023). These models have N ≲ 1.5 × 10⁵, and for their initial mass function and metallicity μ ≃ 0.04 and q ≃ 1/50 such that N_BH ≲ 100, which is similar to the values considered in our two-mass N-body models. As we explored different q 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 N-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 |$\sim 40~{{\ \rm per\ cent}}$|⁠. Although their assumption of a single, hard binary in the core at any given time is supported by our N-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 a (⁠|$\gtrsim 100\, \rm {au}$|⁠) and extremely radial orbit, which nearly circularize 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 (e ≳ 10⁻²). 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.

ACKNOWLEDGEMENTS

The authors thank the referee Nathan Leigh for helpful comments and suggestions. The authors acknowledge financial support from the grants PRE2020-091801, PID2021-125485NB-C22, EUR2020-112157, and 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 grant (AGAUR).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

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 to refer to dynamically formed binaries from single bodies as three-body binaries

2

Unless explicitly stated otherwise, throughout this paper we assume that the binary fraction is sufficiently small such that the mean mass within the cluster is well approximated by the mass of the single stars and that N is equal to the sum of the number of single stars and the number of binaries. The corrections to our results when relaxing this assumption are quantified at the end of this section.

3

The fractional energy change per interaction is lower when considering all interactions, but we use the result for resonant encounters because Δ is not well-defined for non-resonant encounters (Hut 1993).

4

This quantity is reminiscent of the dimensionless central potential W₀ of King’s model (King 1966), and for W₀ ≳ 5 indeed |$W_0\simeq |\phi _0|/\sigma _\mathrm{ c}^2$|⁠.

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.

6

So far we have assumed that N ≃ N_s + N_b, with N_s the number of single stars, and that the mean mass within the core, m_c, is equal to the mean mass of the cluster, that is, m ≃ m_c. If we explicitly keep the m/m_c dependence in the above equations, we obtain that the result for r_h/r_c (equation 3) and N_b (equation 9) need to be multiplied by m/m_c and (m/m_c)³, respectively. The result for N_c (equation 5) would remain unchanged.

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.

8

This derivation assumes a cluster with no primordial binaries or triples. See Leigh & Sills (2011) for a discussion of the general case.