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

Daniel Marín Pina1 and Mark Gieles1,2
1Departament de Física Quàntica i Astrofísica, Institut de Ciències del
Cosmos, Universitat de Barcelona, Martí i Franquès 1, E-08028 Barcelona, Spain
2ICREA, Pg. Lluís Companys 23, E08010 Barcelona, Spain E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### 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 ($\sim 5$), especially for low-$N$ clusters ($\sim
100$), or in the case of two-mass models, low number of black holes. We argue
that the presence of multiple binaries is suppressed by a high rate of binary-
binary interactions, which efficiently ionise one of the binaries involved.
These also lead to triple formation and potentially gravitational wave (GW)
captures, which may provide an explanation for the recently reported high
efficiency of in-cluster mergers in models of low-mass clusters ($\lesssim
10^{5}\,{\rm M}_{\rm\odot})$.

###### keywords:

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

††pubyear: 2023††pagerange: Demographics of three-body binary black holes in
star clusters: implications for gravitational waves–Demographics of three-body
binary black holes in star clusters: implications for gravitational waves

##  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 et al., 2021) – of which the vast majority are binary
black holes (BBH) – which requires us to tackle the question of the origin of
these sources.

Dynamical interactions in the dense cores of stellar clusters have been put
forward (Portegies Zwart & McMillan, 2000) as one of the most likely formation
channels for BBH mergers (Rodriguez et al., 2022; Kremer et al., 2020;
Chattopadhyay et al., 2022; Antonini & Gieles, 2020; Kumamoto et al., 2020;
Banerjee, 2021a; Mapelli et al., 2022). Current GW constraints are compatible
with a majority of massive mergers having formed dynamically (Antonini et al.,
2023). The predictions for the production of BBH mergers in this channel hinge
on the population of BBHs, which points towards the importance of fully
characterising the demographics of these systems. Because the majority of
massive stars form in binaries or higher order multiples (Sana et al., 2012;
Moe & Di Stefano, 2017), a large fraction of BHs end up with a companion after
stellar evolution (hereafter ‘primordial BBHs’). It has been suggested,
however, that in massive clusters the binaries that result in BBH mergers are
dominated by the three-body111We 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 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 et al., 2023). It
is, therefore, important to study the demographics of three-body binaries to
understand dynamical mergers.

Goodman (1984) showed that in steady post-collapse evolution, the number of
three-body binaries, $N_{\rm b}$, depends on the total number of stars in the
cluster, $N$, as $N_{\rm b}\propto 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 organised 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\gtrsim 8000$ stars (Bettwieser & Sugimoto, 1984;
Breeden et al., 1994) and in two-mass models where the number of BHs is
$N_{\rm\mathrm{BH}}\gtrsim 2000$ (Breen & Heggie, 2012). Both of these
assumptions allow us to express relations between several cluster properties,
which are discussed below and summarised in Table 1 and Table 2.

###  2.1 Single-mass clusters

Table 1: Scaling laws for single-mass cluster properties in the steady post core-collapse evolutionary phase. The coefficients $\alpha_{\rm 1}$ and $\alpha_{\rm 2}$ are found from the $N$-body simulations. Value | Scaling law | Coefficients
---|---|---
Velocity dispersion in the core | $\sigma_{\rm c}^{2}=\alpha_{\rm 1}\frac{2}{15}\frac{GmN}{r_{\rm h}}$ | $\alpha_{\rm 1}\simeq 1.4$
Half-mass to core radius ratio | $\frac{r_{\rm h}}{r_{\rm c}}=\alpha_{\rm 2}N^{2/3}$ | $\alpha_{\rm 2}\simeq 0.13$
Net creation rate of binaries | $\Gamma_{\rm b}=0.75\frac{G^{5}m^{5}n_{\rm c}^{3}}{\sigma_{\rm c}^{9}}$ |
Binary hardening rate | $\dot{z}=3.8\frac{G^{2}m^{2}n_{\rm c}}{\sigma^{3}_{\rm c}}$ |
Core radius | $\frac{4\pi G}{9}\rho_{\rm 0}r_{\rm c}^{2}=\sigma_{\rm c}^{2}$ |

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

$\displaystyle\sigma^{2}=\frac{1}{3}\langle
v^{2}\rangle=\frac{1}{3}\frac{GM}{2r_{\rm
v}}\simeq\frac{2}{15}\frac{GNm}{r_{\rm h}}.$ (1)

where $G$ is the gravitational constant, $m$ is the mass of the stars and
$\langle v^{2}\rangle$ is their mean-square velocity. The central dispersion
can be expressed in terms of $\sigma$ by defining a numerical coefficient
$\alpha_{\rm 1}$

$\displaystyle\sigma_{\rm c}^{2}=\alpha_{\rm 1}\sigma^{2}.$ (2)

Next, we define the relation between $r_{\rm h}$ and the core radius, $r_{\rm
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_{\rm h}$ needs to be provided by the core. The exact form of the
$r_{\rm h}/r_{\rm 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\Lambda\simeq\ln(0.11N)$, Giersz & Heggie 1994a) in the
denominator of the central relaxation timescale, and the dimensionless ratio
of the central potential ($|\phi_{\rm 0}|$) to the central velocity
dispersion, $|\phi_{\rm 0}|/\sigma_{\rm c}^{2}$, all are constant, they obtain

$\displaystyle r_{\rm h}/r_{\rm c}\simeq\alpha_{\rm 2}N^{2/3},$ (3)

where $\alpha_{\rm 2}$ is a constant. If the assumptions above are relaxed,
the final ratio would include a dependence on $(|\phi_{\rm 0}|/\sigma_{\rm
c}^{2})/\ln\Lambda$ but, since both these quantities scale weakly with $N$, we
are justified in taking $\alpha_{\rm 2}$ as a constant. In Section 4 we
determine $\alpha_{\rm 2}$ 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_{\rm c}$. Starting from the definition of the core radius
(Heggie & Hut, 2003, p. 71)

$\displaystyle\frac{4\pi G}{9}\rho_{\rm 0}r_{\rm c}^{2}=\sigma_{\rm c}^{2}$
(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 $\rho_{\rm 0}\simeq 1.9\rho_{\rm 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_{\rm c}$,
is222Unless 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.

$\displaystyle N_{\rm c}=0.98\frac{\alpha_{\rm 1}}{\alpha_{\rm
2}}\left(\frac{N}{10^{2}}\right)^{1/3}$ (5)

The $\alpha_{\rm\it 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^{2}/(2a)$ is greater than $x_{\rm h}\simeq m\sigma_{\rm
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_{\rm h}$, and $a>a_{\rm h}=Gm^{2}/(2x_{\rm 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 $\Delta\simeq 0.4$ for resonant encounters333The
fractional energy change per interaction is lower when considering all
interactions, but we use the result for resonant encounters because $\Delta$
is not well-defined for non-resonant encounters (Hut, 1993). (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_{\rm ej}=\frac{6}{\Delta}m\left|\phi_{\rm
0}\right|=15m\left|\phi_{\rm 0}\right|$ (Goodman, 1984), which corresponds to
a semimajor axis of $a_{\rm ej}=Gm^{2}/(2x_{\rm 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_{\rm c}^{2})$ and integrate
it from $z_{\rm h}=1$ to $z_{\rm ej}=15\left|\phi_{\rm 0}\right|/\sigma_{\rm
c}^{2}$, where $\left|\phi_{\rm 0}\right|/\sigma_{\rm c}^{2}$ is a
dimensionless central potential, which is a measure of the (instantaneous)
central concentration444This quantity is reminiscent of the dimensionless
central potential $W_{\rm 0}$ of King’s model (King, 1966), and for $W_{\rm
0}\gtrsim 5$ indeed $W_{\rm 0}\simeq|\phi_{\rm 0}|/\sigma_{\rm c}^{2}$.. 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)

$\displaystyle\Gamma_{\rm b}=0.75\frac{G^{5}m^{5}n_{\rm c}^{3}}{\sigma_{\rm
c}^{9}},$ (6)

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

$\displaystyle\dot{z}=3.8\frac{G^{2}m^{2}n_{\rm c}}{\sigma^{3}_{\rm c}}.$ (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

$\displaystyle f(z)=\frac{\Gamma_{\rm b}}{\dot{z}}=0.20\frac{G^{3}m^{3}n_{\rm
c}^{2}}{\sigma_{\rm c}^{6}}.$ (8)

Assuming that $z_{\rm ej}\gg z_{\rm h}$ and that the relevant binary physics
happens within the volume $V_{\rm c}$ of the cluster core, the total number of
binaries is555Compared 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.

$\displaystyle N_{\rm b}$ $\displaystyle=\int^{z_{\rm ej}}_{\rm z_{\rm
h}}f(z)V_{\rm c}\mathrm{d}z\simeq 3.0\frac{G^{3}m^{3}n_{\rm
c}^{2}}{\sigma_{\rm c}^{8}}V_{\rm c}\left|\phi_{\rm 0}\right|.$ (9)

The above equation can be simplified by using the scaling laws in Table 1 to
obtain666So far we have assumed that $N\simeq N_{\rm s}+N_{\rm b}$, with
$N_{\rm s}$ the number of single stars, and that the mean mass within the
core, $m_{\rm c}$, is equal to the mean mass of the cluster, that is, $m\simeq
m_{\rm c}$. If we explicitly keep the $m/m_{\rm c}$ dependence in the above
equations, we obtain that the result for $r_{\rm h}/r_{\rm c}$ (equation 3)
and $N_{\rm b}$ (equation 9) need to be multiplied by $m/m_{\rm c}$ and
$(m/m_{\rm c})^{3}$, respectively. The result for $N_{\rm c}$ (equation 5)
would remain unchanged.

$\displaystyle N_{\rm b}=2.8\frac{\alpha_{\rm 2}}{\alpha_{\rm
1}}\frac{\left|\phi_{\rm 0}\right|}{\sigma_{\rm
c}^{2}}\left(\frac{N}{10^{2}}\right)^{-1/3}.$ (10)

The ratio $\left|\phi_{\rm 0}\right|/\sigma_{\rm 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_{\rm b}\propto N^{-1/3}$.

###  2.2 Two-mass clusters

Table 2: Scaling laws for two-mass cluster properties in the steady post core-collapse evolutionary phase. The coefficients $\alpha_{\rm 1,BBH}$ and $\alpha_{\rm 2,BBH}$ are found from the $N$-body simulations. Value | Scaling law | Coefficients
---|---|---
Velocity dispersion in the core | $\sigma_{\rm c,\mathrm{BH}}^{2}=\alpha_{\rm 1,\mathrm{BH}}\frac{2}{15}\frac{GM_{\rm\mathrm{BH}}}{r_{\rm h,\mathrm{BH}}}$ | $\alpha_{\rm 1,\mathrm{BH}}\simeq 1.5$
Half-mass to core radius ratio | $\frac{r_{\rm h,\mathrm{BH}}}{r_{\rm c,\mathrm{BH}}}=\alpha_{\rm 2,\mathrm{BH}}N_{\rm\mathrm{BH}}^{2/3}$ | $\alpha_{\rm 2,\mathrm{BH}}\simeq 0.088$
Net creation rate of binaries | $\Gamma_{\rm b}=0.75\frac{G^{5}m_{\rm\mathrm{BH}}^{5}n_{\rm c,\mathrm{BH}}^{3}}{\sigma_{\rm c,\mathrm{BH}}^{9}}$ |
Binary hardening rate | $\dot{z}=3.8\frac{G^{2}m_{\rm\mathrm{BH}}^{2}n_{\rm c,\mathrm{BH}}}{\sigma^{3}_{\rm c,\mathrm{BH}}}$ |
Core radius | $\frac{4\pi G}{9}\rho_{\rm 0,\mathrm{BH}}r_{\rm c,\mathrm{BH}}^{2}=\sigma_{\rm c,\mathrm{BH}}^{2}$ |

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_{\rm\mathrm{BH}}\simeq 20\,{\rm M}_{\rm\odot}$ and that of light stars with
$m_{\rm\star}\simeq 0.5\,{\rm M}_{\rm\odot}$. Due to the large difference
between these bins, we will assume that all BHs have the same mass
$m_{\rm\mathrm{BH}}$ and all stars have the same mass $m_{\rm\star}\ll
m_{\rm\mathrm{BH}}$. Each of these mass bins contribute a total of
$M_{\rm\mathrm{BH}}$ and $M_{\rm\star}$ 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, $\mu\equiv M_{\rm\mathrm{BH}}/M$, which
sets what fraction of the cluster mass is in BHs, and the stellar mass ratio,
$q\equiv\left\langle m\right\rangle/m_{\rm\mathrm{BH}}$, which sets how
massive the BHs are with respect to the average mass of stars and BHs,
$\left\langle m\right\rangle$, in the cluster. Both of these quantities are
defined such that $\mu,q\in(0,1)$. Some authors (for example, Breen & Heggie,
2013) give alternative definitions of these parameters with respect to the
mass of stars, $M_{\rm\star}$ and $m_{\rm\star}$, 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 mass777Although 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., $M_{\rm\star}\gg M_{\rm\mathrm{BH}}$. This allows us to compare
our results to the two-mass model of Breen & Heggie (2013). Their models
assume that equipartition between stars and BHs can not be achieved, which
applies if (Spitzer, 1969)

$\displaystyle\frac{M_{\rm\mathrm{BH}}}{M_{\rm\star}}\left(\frac{m_{\rm\mathrm{BH}}}{m_{\rm\star}}\right)^{3/2}>0.16.$
(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 summarised in Table 2. The results are equivalent
to what we would have found if we considered the BH sub-cluster in isolation,
with the only key difference being the central potential term, which here
includes the contribution of both the stars and the BHs. We reproduce the
equations here, for completeness

$\displaystyle\sigma_{\rm c,\mathrm{BH}}^{2}=\alpha_{\rm
1,\mathrm{BH}}\sigma^{2}_{\rm\mathrm{BH}}=\alpha_{\rm
1,\mathrm{BH}}\frac{2}{15}\frac{GM_{\rm\mathrm{BH}}}{r_{\rm h,\mathrm{BH}}}$
(12) $\displaystyle\frac{r_{\rm h,\mathrm{BH}}}{r_{\rm
c,\mathrm{BH}}}=\alpha_{\rm 2,\mathrm{BH}}N_{\rm\mathrm{BH}}^{2/3}$ (13)
$\displaystyle\frac{4\pi}{9}G\rho_{\rm 0,\mathrm{BH}}r_{\rm
c,\mathrm{BH}}^{2}=\sigma_{\rm c,\mathrm{BH}}^{2}$ (14) $\displaystyle N_{\rm
c,\mathrm{BH}}=0.98\frac{\alpha_{\rm 1,\mathrm{BH}}}{\alpha_{\rm
2,\mathrm{BH}}}\left(\frac{N_{\rm\mathrm{BH}}}{10^{2}}\right)^{1/3}$ (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_{\rm\mathrm{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\mapsto m_{\rm\mathrm{BH}}$. This recovers the same form of the two-mass
binary binding energy distribution from Retterer (1980), i.e.

$\displaystyle N_{\rm\mathrm{BBH}}\simeq
3.0\frac{G^{3}m_{\rm\mathrm{BH}}^{3}n_{\rm c,\mathrm{BH}}^{2}}{\sigma_{\rm
c,\mathrm{BH}}^{8}}V_{\rm c,\mathrm{BH}}\left|\phi_{\rm 0}\right|$ (16)

which we can express as

$\displaystyle N_{\rm\mathrm{BBH}}=2.8\frac{\alpha_{\rm
2,\mathrm{BH}}}{\alpha_{\rm 1,\mathrm{BH}}}\frac{\left|\phi_{\rm
0}\right|}{\sigma_{\rm
c,\mathrm{BH}}^{2}}\left(\frac{N_{\rm\mathrm{BH}}}{10^{2}}\right)^{-1/3}$ (17)

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

##  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_{\rm 0}=-4E_{\rm 0}=1$ (Hénon,
1971), where $M_{\rm 0}$ is the initial cluster mass and $E_{\rm 0}$ is the
initial energy of the cluster. For our assumption of virial equilibrium, the
initial virial radius is $r_{\rm v,0}=-GM_{\rm 0}^{2}/(2W_{\rm 0})=1$, where
$W_{\rm 0}=-1/2$ is the gravitational energy. The corresponding unit of time
is $\tau_{\rm dyn}=(r_{\rm v,0}^{3}/GM_{\rm 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 et al., 2020a), combines an
explicit Hermite integrator for weakly perturbed binaries and a slow-down
method for more compact subgroups.

The centre of the cluster as well as $r_{\rm c}$ are defined following
Casertano & Hut (1985); McMillan et al. (1990), with $r_{\rm 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_{\rm 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\leq N\leq 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 $\tau_{\rm cc}\simeq 16\tau_{\rm rh,0}$ for Plummer models (Cohn,
1980), with $\tau_{\rm rh,0}$ being the initial relaxation timescale (Spitzer
& Hart, 1971),

$\displaystyle\tau_{\rm rh}=0.138\frac{N}{\psi\ln\Lambda}\sqrt{\frac{r_{\rm
h}^{3}}{GM}},$ (18)

evaluated at the start of the simulation, where $\psi=1$ and $\ln\Lambda$ is
the Coulomb logarithm with $\Lambda\simeq 0.11N$ for single-mass models
(Giersz & Heggie, 1994a). The data for the post-collapse evolution is then
taken starting from $t_{\rm 0}=1.1\tau_{\rm cc}$ up until $t_{\rm f}=t_{\rm
0}+20\tau_{\rm 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_{\rm 0}$ and $t_{\rm 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$ | $t_{\rm f}/\tau_{\rm 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_{\rm
0}=1.1\tau_{\rm cc}$ to $t_{\rm f}=t_{\rm 0}+2\tau_{\rm rh,0}$. We use a lower
number of $\tau_{\rm 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 timescale is shorter. The impact of the adopted simulation time is
discussed in Section 6.1. For the core-collapse timescale of two-mass
clusters, we use the result of Kim & Lee (1997)

$\displaystyle\tau_{\rm
cc}=5.3\left(\frac{m_{\rm\star}}{m_{\rm\mathrm{BH}}}\right)^{2}\left(\frac{N_{\rm\star}}{N_{\rm\mathrm{BH}}}\right)^{1/2}\tau_{\rm
rh,0},$ (19)

as well as their expression for $\tau_{\rm rh,0}$ (equation (18), but setting
$\Lambda=0.4N$). The values of the three parameters $N_{\rm\mathrm{BH}}$, $q$
and $\mu$ for the set of simulations that we have run is summarised in Table
4. We set $N\in[5\times 10^{4},1\times 10^{5},2\times 10^{5}]$ and store the
snapshots every 8 $N$-body times. For the mass ratio $q$ we adopt
$q\in[1/50,1/25]$, to approximate metal-poor ($\log_{\rm
10}(Z/Z_{\rm\odot})\simeq-1.5$) and metal-rich ($\log_{\rm
10}(Z/Z_{\rm\odot})\simeq-0.5$) clusters, respectively. For the mass fraction
$\mu$, we have chosen $\mu\in[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$ | $\mu$ | $N_{\rm\mathrm{BH}}$ | Num. of runs
---|---|---|---|---
$5\times 10^{4}$ | 1/50 | 0.025 | 25 | 4
$5\times 10^{4}$ | 1/25 | 0.025 | 50 | 4
$5\times 10^{4}$ | 1/50 | 0.05 | 50 | 4
$5\times 10^{4}$ | 1/25 | 0.05 | 100 | 4
$5\times 10^{4}$ | 1/50 | 0.1 | 100 | 4
$5\times 10^{4}$ | 1/25 | 0.1 | 200 | 4
$1\times 10^{5}$ | 1/50 | 0.025 | 50 | 2
$1\times 10^{5}$ | 1/25 | 0.025 | 100 | 2
$1\times 10^{5}$ | 1/50 | 0.05 | 100 | 2
$1\times 10^{5}$ | 1/25 | 0.05 | 200 | 2
$1\times 10^{5}$ | 1/50 | 0.1 | 200 | 2
$1\times 10^{5}$ | 1/25 | 0.1 | 400 | 2
$2\times 10^{5}$ | 1/50 | 0.025 | 100 | 1
$2\times 10^{5}$ | 1/25 | 0.025 | 200 | 1
$2\times 10^{5}$ | 1/50 | 0.05 | 200 | 1
$2\times 10^{5}$ | 1/25 | 0.05 | 400 | 1
$2\times 10^{5}$ | 1/50 | 0.1 | 400 | 1
$2\times 10^{5}$ | 1/25 | 0.1 | 800 | 1

##  4 Comparing the steady post-collapse model to $N$-BODY simulations

###  4.1 Single-mass clusters

Figure 1: Ratio $\sigma_{\rm c}^{2}/\sigma^{2}$ (equation (1) and equation
(2)) as a function of $N$ for the single-mass $N$-body models (crosses with
error bars). Fits for a constant $\alpha_{\rm 1}$ and an $N$-dependent
$\alpha^{\prime}_{\rm 1}$ are shown with grey and red lines, respectively.
Figure 2: Ratio $r_{\rm h}/r_{\rm c}$ (equation 3) as a function of $N$ for
the single-mass $N$-body models (crosses with error bars). Fits for a constant
$\alpha_{\rm 2}$ and an $N$-dependent $\alpha^{\prime}_{\rm 2}$ are shown with
grey and red lines, respectively. 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 $\alpha_{\rm\it i}$ independent of $N$ (gray) we are able
to roughly match the $N$-body models, whereas using the $N$-dependant
$\alpha^{\prime}_{\rm\it i}$ coefficients (red) yields a very accurate match.
Figure 4: Ratio $\left|\phi_{\rm 0}\right|/\sigma_{\rm c}^{2}$ (equation 22)
as a function of $N$ for the single-mass $N$-body models (crosses with error
bars). 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.

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
$\alpha_{\rm 1}$, $\alpha_{\rm 2}$ parameters (that we label as
$\alpha^{\prime}_{\rm 1}$ and $\alpha^{\prime}_{\rm 2}$, respectively). In
Fig. 1, we show the ratio $\alpha_{\rm 1}=\sigma_{\rm c}^{2}/\sigma^{2}$ as a
function of $N$. We find that a constant $\alpha_{\rm 1}$ (equation 2) is a
good approximation only for $N\lesssim 10^{3}$, and $\alpha^{\prime}_{\rm 1}$
increases slightly at larger $N$. A scaling $\alpha^{\prime}_{\rm
1}=1.1(N/10^{2})^{0.06}$ describes the $N$-body data well. This increase is
likely due to the increase of $r_{\rm h}/r_{\rm c}$ with $N$, leading to a
higher central density and dispersion, relative to the global values, for
clusters with larger $N$.

The $r_{\rm h}/r_{\rm c}$ ratio is shown Fig. 2. The scaling of equation (3)
is a good approximation for constant $\alpha_{\rm 2}$ for $N\simeq
10^{3}-10^{4}$, 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 timescale and in the
$\left|\phi_{\rm 0}\right|/\sigma_{\rm c}^{2}$ term. A scaling of
$\alpha^{\prime}_{\rm 2}=0.19(N/10^{2})^{-0.09}$ provides a good fit to the
$N$-body data.

From these scaling laws, we are able to accurately predict the number of stars
in the core (Fig. 3). The prediction for $N_{\rm c}$ with constant
$\alpha_{\rm\it i}$ (equation 5) has the same form than the one given in
Goodman (1984),

$\displaystyle N_{\rm c}=10\left(\frac{\alpha_{\rm
1}}{1.3}\right)\left(\frac{\alpha_{\rm
2}}{0.13}\right)^{-1}\left(\frac{N}{10^{2}}\right)^{1/3},$ (20)

which yields a rough estimate of $N_{\rm c}$ from the $N$-body models within
<50% over two orders of magnitude in $N$. Using the $N$-dependent expressions
for $\alpha^{\prime}_{\rm 1}$ and $\alpha^{\prime}_{\rm 2}$ yields

$\displaystyle N_{\rm
c}^{\prime}=0.58\left(\frac{N}{10^{2}}\right)^{0.17}N_{\rm c}$ (21)

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

$\displaystyle\frac{\left|\phi_{\rm 0}\right|}{\sigma_{\rm c}^{2}}\simeq
7.2\left(\frac{N}{10^{2}}\right)^{0.1}.$ (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_{\rm b}$ we find in
the simulations. We obtain

$\displaystyle N_{\rm b}=1.9\left(\frac{\alpha_{\rm
2}}{0.13}\right)\left(\frac{\alpha_{\rm
1}}{1.3}\right)^{-1}\left(\frac{\left|\phi_{\rm 0}\right|/\sigma_{\rm
c}^{2}}{7.2}\right)\left(\frac{N}{10^{2}}\right)^{-1/3}$ (23)

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

$\displaystyle N^{\prime}_{\rm
b}=1.7\left(\frac{N}{10^{2}}\right)^{-0.17}N_{\rm b},$ (24)

such that a cluster with $\sim 100$ stars is expected to have $\sim 3$
binaries.

In Fig. 5 we show both the prediction for $N_{\rm b}$ and the $N$-body
results. The prediction (equation 24) approximately reproduces the number of
binaries in the core at large $N$ ($N\gtrsim 10^{3}$), whereas in the smallest
clusters $N_{\rm 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_{\rm 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.

###  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 $\alpha_{\rm 1,\mathrm{BH}}$ (equation 12), but a
slight decline of $\alpha_{\rm 1,BH}$ with $N_{\rm\mathrm{BH}}$ is preferred,
equivalent to the single-mass case. For the ratio $r_{\rm
h,\mathrm{BH}}/r_{\rm c,\mathrm{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_{\rm\mathrm{BH}}$
can be attributed to the weak $N_{\rm\mathrm{BH}}$ dependency of the Coulomb
logarithm. However, in the two-mass case the deviation is larger due to the
smaller number of BHs in the core.

From these relations, we can predict the number of BHs in the sub-cluster core
(Fig. 8). As in the single-mass case, constant $\alpha_{\rm{\it
i},\mathrm{BH}}$ (equation 15) give a rough approximation to $N_{\rm
c,\mathrm{BH}}$,

$\displaystyle N_{\rm c,\mathrm{BH}}=16\left(\frac{\alpha_{\rm
1,\mathrm{BH}}}{1.5}\right)\left(\frac{\alpha_{\rm
2,\mathrm{BH}}}{0.088}\right)^{-1}\left(\frac{N_{\rm\mathrm{BH}}}{10^{2}}\right)^{1/3}$
(25)

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

$\displaystyle N^{\prime}_{\rm
c,\mathrm{BH}}=0.76\left(\frac{N_{\rm\mathrm{BH}}}{10^{2}}\right)^{0.23}N_{\rm
c,\mathrm{BH}},$ (26)

which accurately reproduces the results from the $N$-body models for
$N_{\rm\mathrm{BH}}\gtrsim 50$, although it slightly overpredicts the results
for the sub-clusters with the fewest black holes in the core ($N_{\rm
c,\mathrm{BH}}\lesssim 7$).

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

$\displaystyle\frac{\left|\phi_{\rm 0}\right|}{\sigma_{\rm
c,\mathrm{BH}}^{2}}\simeq 35\left(\frac{\mu}{0.025}\right)^{-0.62}.$ (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

$\displaystyle\frac{\left|\phi_{\rm 0}\right|}{\sigma_{\rm c}^{2}}\simeq
4\left(\frac{\mu}{0.025}\right)^{-0.22}\left(\frac{q}{1/25}\right)^{0.4}.$
(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

$\displaystyle N_{\rm\mathrm{BBH}}=6\left(\frac{\alpha_{\rm
2,\mathrm{BH}}}{0.088}\right)\left(\frac{\alpha_{\rm
1,\mathrm{BH}}}{1.5}\right)^{-1}\left(\frac{\left|\phi_{\rm
0}\right|/\sigma_{\rm
c,\mathrm{BH}}^{2}}{35}\right)\left(\frac{\mu}{0.025}\right)^{0.62}\left(\frac{N_{\rm\mathrm{BH}}}{10^{2}}\right)^{-1/3}.$
(29)

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

$\displaystyle
N^{\prime}_{\rm\mathrm{BBH}}=1.3\left(\frac{N_{\rm\mathrm{BH}}}{10^{2}}\right)^{-0.23}N_{\rm\mathrm{BBH}}.$
(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 $\mu$ and smaller $N_{\rm\mathrm{BH}}$.

Figure 6: Ratio $\sigma_{\rm c,\mathrm{BH}}^{2}/\sigma_{\rm\mathrm{BH}}^{2}$
(equation 12) as a function of $N_{\rm\mathrm{BH}}$ for the two-mass $N$-body
models (crosses with error bars). Fits for a constant $\alpha_{\rm
1,\mathrm{BH}}$ and an $N$-dependent $\alpha^{\prime}_{\rm 1,\mathrm{BH}}$ are
shown with grey and red lines, respectively.  Figure 7: Ratio $r_{\rm
h,\mathrm{BH}}/r_{\rm c,\mathrm{BH}}$ (equation 13) as a function of
$N_{\rm\mathrm{BH}}$ for the two-mass $N$-body models (crosses with error
bars). Fits for a constant $\alpha_{\rm 2,\mathrm{BH}}$ and an $N$-dependent
$\alpha^{\prime}_{\rm 2,\mathrm{BH}}$ are shown with grey and red lines,
respectively.  Figure 8: Number of BHs in the sub-cluster core as a function
of $N_{\rm\mathrm{BH}}$ for the two-mass $N$-body models (crosses with error
bars). Using $\alpha_{\rm{\it i},\mathrm{BH}}$ independent of
$N_{\rm\mathrm{BH}}$ (gray) we are able to roughly match the $N$-body models,
whereas using the $N_{\rm\mathrm{BH}}$-dependant $\alpha^{\prime}_{\rm{\it
i},\mathrm{BH}}$ coefficients (red) yields a good match. Figure 9: Ratio
$\left|\phi_{\rm 0}\right|/\sigma_{\rm c,\mathrm{BH}}^{2}$ (equation 27) as a
function of $\mu$ for the two-mass $N$-body models (crosses with error bars).
Figure 10: Number of BBHs as a function of $N_{\rm\mathrm{BH}}$ and $\mu$ 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 interactions

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 ionised ($\sim
95\%$ of the outcomes for the scattering experiments in Antognini & Thompson
2016), which can explain the dearth of binaries. Because the total number of
binaries in the $N$-body models is $\mathcal{O}(1)$, ignoring interactions
among binaries is a justified assumption for clusters with $N_{\rm c}\gg 1$.
However, as we have shown in Fig. 3, $N_{\rm c}\lesssim 20$ for $N\lesssim
10^{3}$, whilst the theory predicts $N_{\rm b}\gtrsim 2$. In fact, the
predicted binary fraction in the core is $f_{\rm b}=N_{\rm b}/N_{\rm c}\simeq
0.6(N/10^{2})^{-0.9}$. Therefore, if the model for the binary energy
distribution in the steady post-collapse model is correct, then binary-binary
interactions would be important, so we need to understand their effect on
shaping the energy distribution of binaries.

We start by estimating the critical binary fraction in the core, above which
binary-binary interactions become more frequent than binary-single
interactions. The interaction cross-section for a binary with a single goes as
$\Sigma_{\rm bs}\propto Gr_{\rm p}m_{\rm tot}/v_{\rm bs}^{2}$, where $v_{\rm
bs}$ is the relative velocity, $r_{\rm p}\simeq a$ is the relevant minimum
distance (with $a$ the semimajor axis of the binary), and $m_{\rm tot}=3m$ is
the total mass of the interacting binary and single. Assuming equipartition,
the typical relative velocity among binaries, $v_{\rm bb}$, is a factor of
$\sqrt{3/4}$ lower than among binaries and singles ($v_{\rm
bb}=\sqrt{3/4}v_{\rm bs}$). Then, taking $m_{\rm tot}=4m$ and $r_{\rm p}=2a$
for binary-binary interactions, we estimate that the cross section for these
is $\Sigma_{\rm bb}\simeq 3.6\Sigma_{\rm bs}$. The ratio of the rate of
binary-binary interactions to the rate of binary-single interactions,
$\Gamma_{\rm bb}/\Gamma_{\rm bs}$, is

$\displaystyle\frac{\Gamma_{\rm bb}}{\Gamma_{\rm bs}}\simeq\frac{N_{\rm
b}^{2}\Sigma_{\rm bb}v_{\rm bb}}{N_{\rm b}N_{\rm s}\Sigma_{\rm bs}v_{\rm
bs}}\simeq 3.1f_{\rm b}$ (31)

This estimate shows that binary-binary interactions are expected to be more
frequent than binary-single interactions for888This derivation assumes a
cluster with no primordial binaries or triples. See Leigh & Sills (2011) for a
discussion of the general case. $f_{\rm b}\gtrsim 0.3$ in the core. Although
this may seem high, the core contains of order 10 stars, so this criterion is
already met if 3 binaries are created. This is approximately the number of
hard binaries predicted by Goodman & Hut (1993) for clusters $N\simeq 100$,
and we note that their models predict an even higher formation rate of soft
binaries which have a larger cross section for interactions and can therefore
also contribute to binary ionisation in a way that is not captured by the
Goodman & Hut (1993) model. We will therefore test whether binary-binary
interactions can reduce the number of binaries seen in the $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 ionised. If the timescale for this
ionisation mechanism is shorter than the time needed to form new binaries,
then the predicted uniform binding energy distribution (equation 8) is never
realised. Instead, one would find a decline in binaries at higher $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 realised, 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_{\rm 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_{\rm t}$ as a
function of $N$ (crosses with error bars) in the single-mass $N$-body model.
The non-zero value of $N_{\rm 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 $\Gamma_{\rm 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

$\displaystyle\Sigma(a_{\rm 1},a_{\rm 2},v)=\pi\left(a_{\rm 1}^{2}+a_{\rm
2}^{2}\right)\left(\frac{\hat{v}^{2}}{0.567}+\frac{\hat{v}^{6}}{0.0279}\right)^{-1}$
(32)

The cross section depends on the semimajor axes $a_{\rm 1}$, $a_{\rm 2}$ of
the incoming binaries and the ratio $\hat{v}$ of their relative velocity to
their critical velocity, $\hat{v}=v/\sqrt{Gm(a_{\rm 1}+a_{\rm 2})/(a_{\rm
1}a_{\rm 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 $\sigma_{\rm c}$, then the dispersion of relative velocities
among singles is $\sqrt{2}\sigma_{\rm 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 $\sigma_{\rm c}$ ($f_{\rm Maxwell}(v)\propto
v^{2}/\sigma_{\rm c}^{3}\exp(-0.5v^{2}/\sigma_{\rm c}^{2}$)). Therefore, the
triple formation rate is

$\displaystyle\Gamma_{\rm t,bb}=\iiint\frac{f(z_{\rm 1})f(z_{\rm 2})}{2}f_{\rm
Maxwell}(v)\Sigma(z_{\rm 1},z_{\rm 2},v)v\mathrm{d}z_{\rm 1}\mathrm{d}z_{\rm
2}\mathrm{d}v$ (33)

The ratio $\Gamma_{\rm t,bb}/\Gamma_{\rm 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\lesssim 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 $\Gamma_{\rm t,bb}$ to the
binary formation rate $\Gamma_{\rm b}$. This ratio is $\gg 1$ for $N\lesssim
10^{3}$, which confirms the importance of binary-binary interactions in that
$N$ regime. The ratio scales as $\Gamma_{\rm t,bb}/\Gamma_{\rm b}\propto
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. 15 that a uniform BBH binding energy distribution
is not realised, with the discrepancy being greatest at lower
$N_{\rm\mathrm{BH}}$ and higher $z$. Contrary to the single-mass case, here
the models with lower $N_{\rm BH}$ have higher $z_{\rm ej}$ because of their
smaller $\mu$ (see equation 29). Furthermore, we show that BBH-BBH
interactions happen by looking at the formation of stable triple BH systems.
Indeed, in Fig. 14 we can see that stable BH triple formation happens for all
clusters, independently of their number size, BH mass fraction or mass ratio.
In the two-mass case, the values of $N_{\rm\mathrm{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);
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: Number of dynamically-formed stable BH triples $N_{\rm
t,\mathrm{BH}}$ as a function of $N_{\rm\mathrm{BH}}$ and $\mu$ (crosses with
error bars) in the two-mass $N$-body model. The non-null value of $N_{\rm
t,\mathrm{BH}}$ implies the presence of BBH-BBH interactions. Figure 15:
Distribution of the binding energies of the BBHs as a function of
$N_{\rm\mathrm{BH}}$ and $\mu$ for the two-mass $N$-body models with $N=50000$
and $q=1/50$. The theory (dotted lines) does not match the observed
distribution (full lines), with the discrepancy being highest at lower
$N_{\rm\mathrm{BH}}$, $\mu$ and higher $z$.

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_{\rm\mathrm{BH}}$ and $\mu$. 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 $\Gamma_{\rm
t,bb}$ to the BBH formation rate $\Gamma_{\rm b}$ as a function of the number
of BHs $N_{\rm\mathrm{BH}}$ (top panel) and the mass fraction $\mu$ (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_{\rm t,bb}/\Gamma_{\rm
b}\propto\mu^{-1.2}N_{\rm\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 (Fig. 11 and Fig. 15), the
properties of the cluster must evolve slowly with time. The criterion for this
is that the timescale at which the cluster expands, $\tau_{\rm exp}$, is
larger than the binary lifecycle, $\tau_{\rm bin}\simeq(z_{\rm ej}-z_{\rm
h})/\dot{z}\simeq z_{\rm ej}/\dot{z}$. We can define the timescale for
expansion similarly to $\tau_{\rm bin}$ by assuming that $\tau_{\rm exp}\simeq
r_{\rm h}/\dot{r}_{\rm h}$. Using Hénon’s principle (Hénon, 1975), we can
relate the energy production of the core to its global properties as

$\displaystyle\dot{E}\simeq\zeta\frac{|E|}{\tau_{\rm rh}}$ (34)

with $\zeta\simeq 0.1$ (Hénon, 1965; Alexander & Gieles, 2012). If we assume
that the cluster energy scales as $E\propto GM^{2}/r_{\rm h}$ and that the
mass loss is negligible, we obtain

$\displaystyle\tau_{\rm exp}=\frac{\left|E\right|}{\dot{E}_{\rm
h}}\simeq\frac{\tau_{\rm rh}}{\zeta}.$ (35)

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

$\displaystyle\tau_{\rm bin}$ $\displaystyle\simeq\frac{z_{\rm
ej}}{\dot{z}}\simeq
300\left(\frac{\psi\log\Lambda}{10}\right)\left(\frac{N}{10^{2}}\right)^{-1.2}\tau_{\rm
rh},$ $\displaystyle\simeq 300\left(\frac{N}{10^{2}}\right)^{-0.2}\tau_{\rm
dyn}.$ (36)

We also give the value of this timescale in terms of the dynamical time
$\tau_{\rm dyn}$ by using that $\tau_{\rm rh}/\tau_{\rm dyn}\simeq
0.1N/(\psi\log\Lambda)$ (Heggie & Hut, 2003). The binary timescale 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_{\rm
c}\propto r_{\rm h}$ and $N_{\rm c}\propto N$), this timescale would increase
strongly with $N$, $\tau_{\rm bin}\propto N^{2}\tau_{\rm dyn}$.

For the single-mass case, we find $\tau_{\rm exp}>\tau_{\rm bin}$ for
$N\gtrsim 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)

$\displaystyle f(z)\propto\frac{m^{3}n_{\rm c}^{2}}{\sigma_{\rm
c}^{6}}\propto\frac{r_{\rm h}^{3}N_{\rm c}^{2}}{N^{3}r_{\rm c}^{6}}\propto
N^{5/3}r_{\rm h}^{-3}$ (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\lesssim 800$.

Another timescale to consider is the duration of our simulations, $t_{\rm
f}-t_{\rm 0}$. If this timescale is much shorter than the binary lifetime
$\tau_{\rm 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_{\rm f}-t_{\rm
0}=20\tau_{\rm rh,0}$. Because the cluster expands, the actual number of
elapsed relaxation times is less than 20. We can estimate the instantaneous
relaxation timescale by assuming that (Hénon, 1965)

$\displaystyle\tau_{\rm rh}(t)=\begin{cases}\tau_{\rm rh,0}&t<\tau_{\rm cc}\\\
(t/\tau_{\rm cc})\tau_{\rm rh,0}&t>\tau_{\rm cc}\end{cases}$ (38)

So we obtain

$\displaystyle t_{\rm f}-t_{\rm 0}=\langle\tau_{\rm rh}\rangle\int^{t_{\rm
f}}_{\rm t_{\rm 0}}\frac{\mathrm{d}t}{\tau_{\rm rh}(t)}=12\langle\tau_{\rm
rh}\rangle$ (39)

which we can compare to equation (6.1) to show that $t_{\rm f}-t_{\rm
0}>\tau_{\rm bin}$ for $N\gtrsim 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)

$\displaystyle\frac{t_{\rm f}-t_{\rm 0}}{\tau_{\rm bin}}\simeq
0.3\left(\frac{N_{\rm
BH}}{10^{2}}\right)^{0.8}\left(\frac{q}{1/50}\right)^{1.1}\left(\frac{\mu}{0.025}\right)^{-1.3}$
(40)

which is larger than unity for models with $N_{\rm BH}\gtrsim 150\ (450)$ for
$q=1/25\ (1/50)$ and $\mu=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_{\rm BH}$ in
Fig. 14 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_{\rm 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 timescale 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_{\rm h}$ and multiply
that rate by the probability that a third star is present within the volume,
which is roughly $n_{\rm c}a_{\rm h}^{3}$. We can estimate the rate of triple
formation due to binary-single-single interactions, $\Gamma_{\rm t,bss}$, by a
similar argument, although performing the change $n_{\rm c}\mapsto n_{\rm b}$
and $a_{\rm h}\mapsto 3a_{\rm h}$ (so that the outer orbit of the triple is
about $\sim 3a_{\rm h}$). Then, the ratio of the rate of triple production
over binary formation is roughly $\Gamma_{\rm t,bss}/\Gamma_{\rm b}\sim n_{\rm
b}(3a_{\rm h})^{3}/(n_{\rm c}a_{\rm h}^{3})=27/N_{\rm c}\simeq
3(N/10^{2})^{-1/3}$, where we use the scaling $N_{\rm 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_{\rm h}$. The formation rate of ‘hard’ triples is a factor of 27 lower.
The $N$-dependence in $\Gamma_{\rm t,bb}/\Gamma_{\rm 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_{\rm 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 $\zeta$ (equation 35). It is
shown in Gieles et al. (2011) that the value of $\zeta$ in tidally-limited
clusters and isolated clusters is equal to within 20%, and that the first half
of evolution is similar in tidally-limited and isolated clusters. Furthermore,
the authors estimate that two thirds of Milky Way globular clusters are in
this first expansion phase. Per Hénon’s principle, this energy flow is
balanced by the production of energy in the core, and therefore the density
profile within a few $r_{\rm h}$ and the overall results of this work are
independent of the underlying host galaxy potential. The only difference that
the tidal stripping introduces is the possibility of evaporation pathways
where the mass loss rate of stars is higher than the mass loss rate of BHs,
such that the cluster evolves to a 100% BH cluster ($\mu\rightarrow 1$, see
Banerjee & Kroupa, 2011; Gieles et al., 2021). As stated in Section 2.2, we do
not consider high-$\mu$ clusters, although we roughly expect them to behave
like the single-mass models in the limit $\mu\rightarrow 1$. Nevertheless,
observations point towards most clusters in the Milky Way having a small BH
mass fraction of less than 1% (Dickson et al., 2023).

###  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 timescale due to the emission of GWs (Samsing,
2018). According to the scattering experiments by Zevin et al. (2019), binary-
binary interactions are roughly five times more likely than binary-single
interactions to produce such mergers due to their more complex resonant
intermediate states.

The binary-binary interactions described in this paper have at least one
binary near the hard-soft boundary. During a resonant intermediate state, two
of the BHs can enter a highly eccentric orbit that ends in a GW capture. This
orbit, therefore, starts with a large semimajor axis ($a\sim a_{\rm h}$),
which requires a nearly radial orbit ($e\sim 1$) to trigger a GW capture. We
can use the capture criterion by Samsing (2018) to quantify what the
eccentricity, $e$, must be

$\displaystyle 1-e\lesssim\left(\frac{4.6Gm}{ac^{2}}\right)^{5/7}.$ (41)

To put it in context, for a $10^{5}\text{ M}_{\rm\odot}$ cluster with $r_{\rm
h}=1\text{ pc}$, a capture between two $10\text{ M}_{\rm\odot}$ BHs would
start at $\sim 60\text{ AU}$ and $1-e\sim 10^{-6}$. As the binary inspirals,
this eccentricity may be radiated away. To estimate this, we find the peak GW
frequency $f$ from Wen (2003),

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

which together with the Peters’ equations (Peters, 1964) allows us to compute
the eccentricity when the binary enters the LIGO-Virgo-KAGRA frequency band
($f\simeq 10\text{ 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_{\rm h}$, $a$ and $m_{\rm\mathrm{BH}}$, we find a merger that has
circularised to $e_{\rm 10Hz}\sim 10^{-2}$, which is currently not detectable
(Lower et al., 2018), but it is predicted that third generation GW detectors
should be able to find this population of BBH mergers, if it exists.
Furthermore, if we assume that after each resonant intermediate state the
eccentricity is sampled from a thermal distribution (Heggie, 1975), about
$\sim 63\%$ ($\sim 0\%$) of captures that start with an initial semimajor axis
of $a_{\rm h}$ ($a_{\rm ej}$) have a remaining $e_{\rm 10Hz}<0.1$ at 10 Hz.
The sampling of $e$ has a chance of producing mergers that form in the LIGO-
Virgo-KAGRA band, $f>10\text{ Hz}$, with a probability of $\sim 12\%$ ($\sim
46\%$).

An alternative path to mergers via binary-binary interactions is the formation
of stable triples. Stable triples with sufficiently high eccentricity can have
large jumps in the eccentricity of the inner orbit (Lidov-Kozai cycles, per
Lidov 1962; Kozai 1962). As above, extreme values of the eccentricity may lead
to GW captures. Per the same argument, the involved semimajor axes may be
sufficiently large that the eccentricity may be of the order of $e_{\rm
10Hz}\sim 10^{-2}$ when the binary enters the LIGO-Virgo-KAGRA band. Modelling
with post-Newtonian terms is required to quantify the importance of mergers
via these two channels.

The short timescale of the above processes implies that the merger happens at
a location near the multiple-body interaction. This constitutes an explanation
for the finding that the majority of mergers occur in-cluster – as opposed to
ejected binaries – in $N$-body models of low-mass clusters ($M\lesssim
10^{5}\,{\rm M}_{\rm\odot}$, Rastello et al. 2019; Banerjee 2021b; Barber et
al. 2023). These models have $N\lesssim 1.5\times 10^{5}$, and for their
initial mass function and metallicity $\mu\simeq 0.04$ and $q\simeq 1/50$ such
that $N_{\rm BH}\lesssim 100$, which is similar to the values considered in
our two-mass $N$-body models. As we explored different $q$ and $\mu$, 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\%$. 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
circularise at 10 Hz. These mergers may be a signature of binary-binary
interactions at the hard-soft boundary. These binaries contribute to the
eccentricity distribution by increasing the expected rate of mergers at lower
eccentricities ($e\gtrsim 10^{-2}$). Such mildly eccentric mergers could be
disentangled from quasi-circular mergers in future GW detectors, such as the
Einstein Telescope (Maggiore et al., 2020) and Cosmic Explorer (Evans et al.,
2021). In conclusion, we present an additional step towards a complete
prediction for the rate of dynamically-formed BBH mergers that can be detected
in current and upcoming GW experiments.

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

## References

  * Aarseth (2012) Aarseth S. J., 2012, MNRAS, 422, 841
  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2021) Abbott R., et al., 2021, Physical Review X, 11, 021053
  * Alexander & Gieles (2012) Alexander P. E. R., Gieles M., 2012, MNRAS, 422, 3415
  * Antognini & Thompson (2016) Antognini J. M. O., Thompson T. A., 2016, MNRAS, 456, 4219
  * Antonini & Gieles (2020) Antonini F., Gieles M., 2020, Phys. Rev. D, 102, 123016
  * Antonini et al. (2023) Antonini F., Gieles M., Dosopoulou F., Chattopadhyay D., 2023, MNRAS, 522, 466
  * Banerjee (2018) Banerjee S., 2018, MNRAS, 481, 5123
  * Banerjee (2021a) Banerjee S., 2021a, MNRAS, 503, 3371
  * Banerjee (2021b) Banerjee S., 2021b, MNRAS, 503, 3371
  * Banerjee & Kroupa (2011) Banerjee S., Kroupa P., 2011, ApJ, 741, L12
  * Barber et al. (2023) Barber J., Chattopadhyay D., Antonini F., 2023, arXiv e-prints, p. arXiv:2310.02376
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
  * Bettwieser & Sugimoto (1984) Bettwieser E., Sugimoto D., 1984, MNRAS, 208, 493
  * Breeden et al. (1994) Breeden J. L., Cohn H. N., Hut P., 1994, ApJ, 421, 195
  * Breen & Heggie (2012) Breen P. G., Heggie D. C., 2012, MNRAS, 420, 309
  * Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779
  * Casertano & Hut (1985) Casertano S., Hut P., 1985, ApJ, 298, 80
  * Chattopadhyay et al. (2022) Chattopadhyay D., Hurley J., Stevenson S., Raidani A., 2022, MNRAS, 513, 4527
  * Cohn (1980) Cohn H., 1980, ApJ, 242, 765
  * Dickson et al. (2023) Dickson N., Smith P. J., Hénault-Brunet V., Gieles M., Baumgardt H., 2023, arXiv e-prints, p. arXiv:2308.13037
  * Evans et al. (2021) Evans M., et al., 2021, arXiv e-prints, p. arXiv:2109.09882
  * Gieles et al. (2011) Gieles M., Heggie D. C., Zhao H., 2011, MNRAS, 413, 2509
  * Gieles et al. (2021) Gieles M., Erkal D., Antonini F., Balbinot E., Peñarrubia J., 2021, Nature Astronomy, 5, 957
  * Giersz & Heggie (1994a) Giersz M., Heggie D. C., 1994a, MNRAS, 268, 257
  * Giersz & Heggie (1994b) Giersz M., Heggie D. C., 1994b, MNRAS, 270, 298
  * Goodman (1984) Goodman J., 1984, ApJ, 280, 298
  * Goodman & Hut (1993) Goodman J., Hut P., 1993, ApJ, 403, 271
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Heggie & Hut (1993) Heggie D. C., Hut P., 1993, ApJS, 85, 347
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press
  * Hénon (1965) Hénon M., 1965, Annales d’Astrophysique, 28, 62; English translation: arXiv:1103.3498
  * Hénon (1971) Hénon M. H., 1971, Ap&SS, 14, 151
  * Hénon (1975) Hénon M., 1975, in A. Hayli ed., Proc. IAU Symp. 69, Dynamics of the Solar Systems. Reidel, Dordrecht, p. 133.
  * Hut (1993) Hut P., 1993, ApJ, 403, 256
  * Kim & Lee (1997) Kim S. S., Lee H. M., 1997, Journal of Korean Astronomical Society, 30, 115
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Kremer et al. (2020) Kremer K., et al., 2020, ApJ, 903, 45
  * Kumamoto et al. (2020) Kumamoto J., Fujii M. S., Tanikawa A., 2020, MNRAS, 495, 4268
  * Küpper et al. (2011) Küpper A. H. W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300
  * Leigh & Sills (2011) Leigh N., Sills A., 2011, MNRAS, 410, 2370
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Lower et al. (2018) Lower M. E., Thrane E., Lasky P. D., Smith R., 2018, Phys. Rev. D, 98, 083028
  * Lynden-Bell & Eggleton (1980) Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483
  * Maggiore et al. (2020) Maggiore M., et al., 2020, J. Cosmology Astropart. Phys., 2020, 050
  * Mapelli et al. (2022) Mapelli M., Bouffanais Y., Santoliquido F., Arca Sedda M., Artale M. C., 2022, MNRAS, 511, 5797
  * Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  * McMillan et al. (1990) McMillan S., Hut P., Makino J., 1990, ApJ, 362, 522
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Peters (1964) Peters P. C., 1964, Physical Review, 136, 1224
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Rastello et al. (2019) Rastello S., Amaro-Seoane P., Arca-Sedda M., Capuzzo-Dolcetta R., Fragione G., Tosta e Melo I., 2019, MNRAS, 483, 1233
  * Retterer (1980) Retterer J. M., 1980, AJ, 85, 249
  * Rodriguez et al. (2022) Rodriguez C. L., et al., 2022, ApJS, 258, 22
  * Samsing (2018) Samsing J., 2018, Phys. Rev. D, 97, 103014
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Spitzer (1969) Spitzer L. J., 1969, ApJ, 158, L139
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters. Princeton University Press
  * Spitzer & Hart (1971) Spitzer Lyman J., Hart M. H., 1971, ApJ, 164, 399
  * Tanikawa et al. (2012) Tanikawa A., Hut P., Makino J., 2012, New Astron., 17, 272
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration et al., 2021, arXiv e-prints, p. arXiv:2111.03606
  * Torniamenti et al. (2022) Torniamenti S., Rastello S., Mapelli M., Di Carlo U. N., Ballone A., Pasquato M., 2022, MNRAS, 517, 2953
  * Wang et al. (2020a) Wang L., Nitadori K., Makino J., 2020a, MNRAS, 493, 3398
  * Wang et al. (2020b) Wang L., Iwasawa M., Nitadori K., Makino J., 2020b, MNRAS, 497, 536
  * Wen (2003) Wen L., 2003, ApJ, 598, 419
  * Zevin et al. (2019) Zevin M., Samsing J., Rodriguez C., Haster C.-J., Ramirez-Ruiz E., 2019, ApJ, 871, 91