# The Dragon-II simulations - III. Compact binary mergers in clusters with up
to 1 million stars: mass, spin, eccentricity, merger rate and pair instability
supernovae rate.

Manuel Arca Sedda1,2,3, Albrecht W. H. Kamlah4,3, Rainer Spurzem3,5,6,
Francesco Paolo Rizzuto7, Mirek Giersz9, Thorsten Naab7, Peter Berczik3,10,11
1 Gran Sasso Science Institute (GSSI), 67100 L’Aquila, Italy
2 Physics and Astronomy Department Galileo Galilei, University of Padova,
Vicolo dell’Osservatorio 3, I–35122, Padova, Italy
3 Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität
Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany
4 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
5 National Astronomical Observatories and Key Laboratory of Computational
Astrophysics, Chinese Academy of Sciences, 20A Datun Rd.,Chaoyang District,
100101, Beijing, China
6 Kavli Institute for Astronomy and Astrophysics, Peking University, Yiheyuan
Lu 5, Haidian Qu, 100871, Beijing, China
7 Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2,
FI-00014, Helsinki, Finland
8 Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85740,
Garching, Germany
9 Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, ul.
Bartycka 18, 00-716 Warsaw, Poland
10 Konkoly Observatory, Research Centre for Astronomy and Earth Sciences,
Eötvös Loránd Research Network (ELKH), MTA Centre of Excellence, Konkoly Thege
Miklós út 15-17, 1121 Budapest, Hungary
11 Main Astronomical Observatory, National Academy of Sciences of Ukraine, 27
Akademika Zabolotnoho St., 03680, Kyiv, Ukraine
Contact<EMAIL_ADDRESS>
###### Abstract

Compact binary mergers forming in star clusters may exhibit distinctive
features that can be used to identify them among observed gravitational-wave
(GW) sources. Such features likely depend on the host cluster structure and
the physics of massive star evolution. Here, we dissect the population of
compact binary mergers in the Dragon-II simulation database, a suite of 19
direct $N$-body models representing dense star clusters with up to $10^{6}$
stars and $<33\%$ of stars in primordial binaries. We find a substantial
population of black hole binary (BBH) mergers, some of them involving an
intermediate-mass BH (IMBH), and a handful mergers involving a stellar BH and
either a neutron star (NS) or a white dwarf (WD). Primordial binary mergers,
$\sim 30\%$ of the whole population, dominate ejected mergers. Dynamical
mergers, instead, dominate the population of in-cluster mergers and are
systematically heavier than primordial ones. Around $20\%$ of Dragon-II
mergers are eccentric in the LISA band and $5\%$ in the LIGO band. We infer a
mean cosmic merger rate of $\mathcal{R}\sim 12(4.4)(1.2)$ yr-1 Gpc3 for BBHs,
NS-BH, and WD-BH binary mergers, respectively, and discuss the prospects for
multimessenger detection of WD-BH binaries with LISA. We model the rate of
pair-instability supernovae (PISNe) in star clusters and find that surveys
with a limiting magnitude $m_{\rm bol}=25$ can detect $\sim 1-15$ yr-1 PISNe.
Comparing these estimates with future observations could help to pin down the
impact of massive star evolution on the mass spectrum of compact stellar
objects in star clusters.

###### keywords:

methods: numerical – galaxies: star clusters: general – stars: general, black
holes

††pubyear: 2022††pagerange: The Dragon-II simulations - III. Compact binary
mergers in clusters with up to 1 million stars: mass, spin, eccentricity,
merger rate and pair instability supernovae rate.–The Dragon-II simulations -
III. Compact binary mergers in clusters with up to 1 million stars: mass,
spin, eccentricity, merger rate and pair instability supernovae rate.

##  1 Introduction

In less than a decade, the LIGO-Virgo-Kagra (LVK) collaboration discovered 76
confident gravitational-wave (GW) sources associated to merging stellar black
holes (BHs) and neutron stars (NSs) (The LIGO Scientific Collaboration et al.,
2021b). This number raises up to 90 if one considers the population of events
with a probability to have an astrophysical origin $>0.5$ (The LIGO Scientific
Collaboration et al., 2021a), and it is destined to further increase by the
end of the fourth observation run. Measurable quantities like component
masses, spins, or the orbital eccentricity, and the merger rate of different
types of compact binary mergers can represent the keys to identify the
signatures of different formation channels (Arca Sedda & Benacquista, 2019;
Arca Sedda et al., 2020b; Zevin et al., 2020; Arca Sedda et al., 2023;
Bouffanais et al., 2021; Mapelli et al., 2022). From the theoretical
standpoint, there is a plethora of mechanisms proposed to explain the
formation of compact binary mergers, like isolated binary evolution
(Belczynski et al., 2002; Dominik et al., 2012; Belczynski et al., 2016;
Giacobbo et al., 2018; Spera et al., 2019), dynamical pairing in dense star
clusters (Miller & Hamilton, 2002; Downing et al., 2010; Rodriguez et al.,
2016; Askar et al., 2017; Banerjee, 2018; Di Carlo et al., 2019; Rizzuto et
al., 2022), formation in AGN disks (McKernan et al., 2012; Stone et al., 2017;
McKernan et al., 2018; Tagawa et al., 2020), secular dynamics involving three
compact objects (Antonini et al., 2018; Vigna-Gómez et al., 2021) or a binary
orbiting a supermassive black hole (Antonini & Perets, 2012; Hoang et al.,
2018; Fragione et al., 2019; Arca Sedda, 2020b), and primordial BH evolution
(Carr & Hawking, 1974; Carr et al., 2016; Sasaki et al., 2016). The majority
of the aforementioned mechanisms relies on the assumption that compact objects
are the relic of massive stars, and therefore they suffer the uncertainties
affecting stellar evolution.

For example, the insurgence of pair instability supernova (PISN) and
pulsational pair instability supernova (PPISN) mechanisms can carve in the BH
mass spectrum the so-called upper-mass gap, a region extending in the range
$m_{\rm gap}=40-150{\rm~{}M}_{\odot}$ where no remnants are expected. The
boundaries of the gap are highly uncertain and depend on many poorly
constrained quantities, like stellar rotation, rate of nuclear reactions,
stellar evolution model (Woosley & Heger, 2021; Vink et al., 2021; Stevenson
et al., 2019; Farmer et al., 2019; Costa et al., 2021; Iorio et al., 2022).
The presence of several upper mass-gap BH candidates in the LVK source
catalogue poses the question about the origin of these BHs. Stellar mergers,
star-BH interactions, and repeated BH mergers represent possible pathways to
overcome (P)PISN (e.g. Spera et al., 2019; Banerjee, 2022; Costa et al., 2022;
Ballone et al., 2022) and produce merging compact objects in dense star
clusters (e.g. Rodriguez et al., 2018; Di Carlo et al., 2020a; Kremer et al.,
2020; Di Carlo et al., 2021; Rizzuto et al., 2021; Arca-Sedda et al., 2021c;
Rizzuto et al., 2022).

Spins could carry crucial information on the BH formation scenario and help
placing constraints on the evolution of massive stars, but little is known
about the distribution of stellar BH natal spins. Observations of merging BHs
indicate that the spin distribution follows a Maxwellian distribution, with a
peak around $\chi_{\rm BH}\sim 0.2-0.5$ (The LIGO Scientific Collaboration et
al., 2021b). However, stellar BHs detected in low-mass X-ray binaries (LMXBs)
are characterised by spins broadly distributed in the whole allowed range
(Fragos & McClintock, 2015), whilst those in high-mass X-ray binaries (HMXBs)
involve BHs almost maximally spinning (see e.g. Qin et al., 2019; Reynolds,
2021). Despite these differences may suffer observation biases, they may
represent peculiarities of different evolutionary pathways. Efficient angular
momentum transport driven by magnetic stars could trigger the formation of BHs
with natal spins as small as $\chi_{\rm BH}\lesssim 0.01$, a mechanism
proposed for BHs forming from single stars and in binaries with a negligible
mass transfer among the components (Fuller & Ma, 2019). Significant mass
transfer, instead, has been proposed to produce BHs with spin in a broad range
in LMXB, even for BHs spinless at birth (Fragos & McClintock, 2015), and
nearly extremal BHs in HMXBs (Qin et al., 2019; Gallegos-Garcia et al., 2022).
Common envelope evolution in massive stellar binaries can lead to merging BBHs
consisting in a nearly non-rotating BH (Qin et al., 2018; Bavera et al.,
2020), although this strongly depends on the stellar evolution adopted
(Belczynski et al., 2020), and a BH companion with a spin spanning the whole
allowed range of values (Qin et al., 2018; Bavera et al., 2020; Belczynski et
al., 2020).

Amplitude aside, also the alignment of the spin vectors among themselves and
with the binary angular momentum can affect both the waveform, the final
merger remnant mass and spin, and the recoil kick (e.g. see Equation 10). From
an "observational" perspective, measuring the spin components is intrinsically
hard and their directions generally vary owing to precession, thus the spin of
observed mergers can be characterised through the so-called effective spin
parameter (Racine, 2008; Santamaría et al., 2010; Ajith et al., 2011)

$\chi_{\rm eff}=\frac{\vec{\chi}_{1}+q\vec{\chi}_{2}}{1+q}\cdot\vec{L},$ (1)

where $q<1$ is the binary mass ratio, $\vec{\chi}_{1,2}$ are the two component
spin vectors and $\vec{L}$ is the binary orbital angular momentum.
Observations of BBH mergers suggest that $\chi_{\rm eff}$ may increase at
increasing the binary merger mass ratio, although some merging binaries
exhibit a negative value of $\chi_{\rm eff}$ (The LIGO Scientific
Collaboration et al., 2021b), a feature generally associated with dynamical
sources.

The orbital eccentricity at merger could represent another distinguishing
feature of compact binary mergers, as dynamical interactions could trigger the
formation of fairly eccentric ($>0.1$) sources contrarily to mergers forming
from isolated binaries (see e.g. Nishizawa et al., 2016). It has been recently
claimed that up to four LVK sources may be eccentric (Romero-Shaw et al.,
2019, 2020; Gayathri et al., 2022; Romero-Shaw et al., 2022), although the
effects of eccentricity and precession can lead to degeneracies in GW data
analysis, making the eccentricity a poorly constrained quantity (see e.g.
Romero-Shaw et al., 2023).

Alongside GWs, the detection of (P)PISNe can represent a key piece to
understand the final stages of massive stars’ life. So far, only a few, most
of which controversial, PISN and PPISN candidates have been observed in the
last two decades (Gal-Yam et al., 2009; Gomez et al., 2019; Woosley & Smith,
2022). The rarity of PISNe observations sets an intrinsic limit on the
frequency of PISNe in star clusters, a quantity poorly constrained in
theoretical and numerical models.

Dynamical interactions among stars in dense and massive star clusters can
trigger both the formation of merging binaries and the development of PISNe,
either from single massive stars or from stellar merger products. Young and
intermediate-age star clusters are particularly interesting environments where
these sources can form, because they are still in their dynamical youth, when
cluster mass-loss and expansion did not yet affected substantially the cluster
structure and the interaction rate among stars is maximal. There is vast
literature investigating the formation and evolution of merging BHs in star
clusters via different techniques, e.g. direct $N$-body simulations (Banerjee
et al., 2010; Downing et al., 2010; Di Carlo et al., 2019; Rastello et al.,
2020; Arca Sedda, 2020a; Di Carlo et al., 2021; Banerjee, 2018, 2021, 2022;
Wang et al., 2022; Chattopadhyay et al., 2022; Rizzuto et al., 2022), Monte
Carlo simulations (Rodriguez et al., 2016; Askar et al., 2017; Kremer et al.,
2019; Rodriguez et al., 2019; Ye et al., 2020; Kremer et al., 2020;
Maliszewski et al., 2022), and semi-analytic tools (Fragione & Kocsis, 2018;
Antonini et al., 2019; Arca Sedda & Benacquista, 2019; Arca Sedda et al.,
2020b; Antonini & Gieles, 2020b; González et al., 2021; Mapelli et al., 2021;
Arca Sedda et al., 2023; Mapelli et al., 2022; Antonini et al., 2022; Kritos
et al., 2022). However, there is lack of direct $N$-body simulations of
particularly dense ($>10^{5}~{}{\rm~{}M}_{\odot}$ pc-3) and massive
($>100,000{\rm~{}M}_{\odot}$) star clusters, owing to the computational cost
required to simulate such systems. Exploring this range of mass and densities
with $N$-body models can complement the already existing simulations and can
offer a term of comparison to Monte Carlo simulations (see Figure 1 in Arca
Sedda et al 2023a, hereafter AS-I).

In this work, which represents the third of a series, we present results from
the Dragon-II star cluster database, a suite of 19 direct $N$-body simulations
of young and intermediate-age star clusters comprised of up to 1 million stars
and up to $33\%$ of stars initially in binaries, characterised by typical
densities $\rho=(1.2\times 10^{4}-1.5\times 10^{7}){\rm~{}M}_{\odot}$ pc-3.

In our previous papers, we focused on the general properties of our cluster
models and their compact object populations (paper AS-I) and the processes
that regulate the formation and growth of IMBHs (Arca Sedda et al 2023b,
hereafter AS-II).

Here, we dissect the properties of BH-BH, BH-NS, and BH-WD mergers developing
in the Dragon-II star cluster database, a suite of 19 direct $N$-body
simulations of star clusters comprised of up to 1 million stars and up to
$33\%$ of stars initially in binaries (details about these models are
discussed in our companion paper AS-I), performed with the Nbody6++GPU
code111https://github.com/nbody6ppgpu/Nbody6PPGPU-beijing. The paper is
organised as follows: in Section 2 we briefly summarise the main features of
our models; Section 3 discusses the main properties of compact binary mergers
in our models, focusing on the component masses and mass ratios, the
eccentricity at merger, and the possible signatures that can identify their
formation history; in Section 4 we explore the impact of BH natal spins onto
the global properties of the population, and we adopt a cosmologically
motivated framework to infer the compact binary merger rate, the detection
perspectives for future low-frequency GW detections, and the frequency rate
and detection perspectives in magnitude limited surveys of PISNe; Section 5
summarises the main results of this work.

##  2 Numerical methods

###  2.1 The Dragon-II clusters

The Dragon-II simulation database consists of 19 star cluster models
characterised by an initial number of stars $N=(1.2,~{}3,~{}6,~{}10)\times
10^{5}$, half-mass radius $R_{\rm HM}=(0.48,~{}0.80,~{}1.76)$ pc, and an
initial binary fraction $f_{b}=0.05-0.2$. In the following, we briefly
summarise the main properties of Dragon-II clusters, referring the interested
readers to our companion paper AS-I for more details on the run properties.

To initialise the Dragon-II clusters we exploit the McLuster tool (Küpper et
al., 2011; Kamlah et al., 2022; Leveque et al., 2022, Leveque in prep.).

Each cluster is modelled according to a King (1966) profile with adimensional
potential well $W_{0}=6$. We adopt an initial metallicity $Z=0.0005$, typical
of several clusters possibly hosting a dense sub-system of compact objects or
an IMBH, like NGC3201 or NGC6254 (Askar et al., 2018; Arca Sedda et al., 2018;
Weatherford et al., 2020).

Star masses are drawn according to a Kroupa (2001) initial mass function
limited in the range $m_{\rm ZAMS}=(0.08-150){\rm~{}M}_{\odot}$. Stars in
primordial binaries are paired depending on their mass, with stars heavier
than $>5{\rm~{}M}_{\odot}$ paired according to a flat mass-ratio distribution,
and lighter stars paired randomly. Binary eccentricities are distributed
according to a thermal distribution, $P(e){\rm d}e=e^{2}{\rm d}e$, while
initial semimajor axes are assigned according to a distribution flat in
logarithmic values limited between the sum of stars’ radii and a maximum value
of $50$ AU.

The host galaxy potential is modelled through a Keplerian potential assuming a
total mass of $M_{\rm gal}=1.78\times 10^{11}{\rm~{}M}_{\odot}$. All Dragon-II
clusters are placed on a circular orbit around this galaxy model at a distance
of $R_{\rm clu}=13.3$ kpc. The adopted galaxy mass and orbital radius lead to
a value of the circular velocity compatible with what is observed in the Milky
Way.

The resulting tidal radius is much larger than the cluster half-mass radius.
Therefore, Dragon-II models are initially underfilling their Roche lobe, which
implies that the initial impact of the host galaxy potential is negligible.

All simulations are terminated when either the mean BH mass falls below
$\langle m_{\rm BH}\rangle\lesssim 15{\rm~{}M}_{\odot}$, there are no BHs with
a mass above $30{\rm~{}M}_{\odot}$, or the simulated time exceeds at least one
relaxation time. As a result, the simulation time in Dragon-II models spans a
range $T_{\rm sim}=0.1-2.3$ Gyr, corresponding to $0.8-80$ times the initial
half-mass relaxation time (see also Table 1).

Over the simulated time, we find a nice overlap (see also Figure 2 in paper
AS-I) between the evolution of Dragon-II clusters’ mass and half-mass radius
and observed properties of young and intermediate-age massive clusters in the
Milky Way (Portegies Zwart et al., 2010), the Magellanic clouds (Gatto et al.,
2021), and other galaxies in the local Universe like Henize 2-10 (Nguyen et
al., 2014) or M83 (Ryon et al., 2015). In this sense, Dragon-II models can
represent one possible evolutionary pathways of (relatively) young massive
clusters.

###  2.2 The Nbody6++GPU code

Dragon-II simulations have been performed with the Nbody6++GPU code (Wang et
al., 2015), a state-of-the-art direct $N$-body integrator that runs on high
performance computing hardware equipped with graphic-processing-units (GPUs,
Spurzem, 1999; Nitadori & Aarseth, 2012; Wang et al., 2015). The code is part
of the famous NBODY code series that was pioneered almost sixty years ago by
Sverre Aarseth (Aarseth et al., 1974; Spurzem, 1999; Aarseth, 1999, 2003;
Aarseth et al., 2008; Nitadori & Aarseth, 2012; Wang et al., 2015; Kamlah et
al., 2022).

The code implements a 4th-order Hermite integrator scheme with adaptive time-
step based on the Ahmad-Cohen scheme for neighbours (Ahmad & Cohen, 1973), and
implements a treatment for close encounters and few-body dynamics via the
Kustaanheimo-Stiefel regularisation (Stiefel & Kustaanheimo, 1965) and chain
regularisation (Mikkola & Aarseth, 1993).

Stellar evolution in Nbody6++GPU is based on an upgraded version of the
population synthesis code BSE (Hurley et al., 2002). The main features of this
state-of-the-art version, named BSE++, are described in detail in Kamlah et
al. (2022) (but see also Banerjee et al., 2020). We adopt the so-called level
B of stellar evolution (see Kamlah et al., 2022), whose main characteristics
are: delayed supernova (SN) scheme (Fryer et al., 2012), pair- and pulsation
pair-instability supernova (PISN and PPISN) treated following Belczynski et
al. (2016); Banerjee (2021), fallback prescription for NS/BH natal kicks, and
metallicity-dependend winds for massive stars Vink et al. (2001); Belczynski
et al. (2010). We refer the reader to Kamlah et al. (2022) and paper AS-I for
further details.

The common envelope phase in binaries is modelled through the widely known
$\alpha_{\rm CE}-\lambda_{\rm CE}$ scheme, which enables us to regulate the
fraction of orbital energy injected into the envelope ($\alpha_{\rm CE}$) and
to scale the binding energy of the envelope by a factor $\lambda_{\rm CE}$. In
this work, we adopt $\alpha_{\rm CE}=3$ and $\lambda=0.5$ (Giacobbo & Mapelli,
2018; Kamlah et al., 2022).

The adopted stellar evolution recipes imply that the stellar BH mass-spectrum
in Dragon-II clusters is limited to $m_{\rm BH,max}=40.5{\rm~{}M}_{\odot}$
(Belczynski et al., 2016), unless BHs form from stellar mergers or star-BH
interactions. In the latter case, Nbody6++GPU parametrises the amount of mass
accreted in a strong star-BH interaction or collision via an accretion
parameter $f_{c}$ (Rizzuto et al., 2021, 2022), which we set to $f_{c}=0.5$.
We refer the reader to Rizzuto et al. (2021, 2022) for a discussion about the
impact of $f_{c}$ on BH evolution.

####  2.2.1 Modelling the final stages of compact object binary mergers

The dynamics of relativistic binaries is followed via the orbit-average
formalism (Peters, 1964), which enables us to follow the evolution of compact
binaries and their coalescence inside the cluster, similarly to previous works
(see e.g. Di Carlo et al., 2019; Di Carlo et al., 2020a, b, 2021; Rizzuto et
al., 2021, 2022; Rastello et al., 2021; Torniamenti et al., 2022).

In its current implementation, Nbody6++GPU follows the dynamics of
relativistic binaries also if they are part of triples (see e.g. Rizzuto et
al., 2021) and multiple systems, as well as if they form via hyperbolic
interactions.

However, the BBH evolution is not followed down to the merger, rather the
binary is decoupled from dynamics and promtply merged when the BBH pericentre
falls below a critical value, which we set to $10^{2}$ Schwarzschild radii,
i.e. $a_{\rm dec}=2kGm_{\rm bin}/c^{2}=kR_{\rm Sch}$ with $k=100$.

Adopting such limiting separation ensures that the binary is unlikely to
undergo any further interaction with surrounding stars before merging.
Considering the range of binary masses ($1-300{\rm~{}M}_{\odot}$), star
cluster masses ($<10^{6}{\rm~{}M}_{\odot}$) and half-mass radii ($0.1-3$ pc)
explored in this work, it is easy to show that the binary—single interaction
timescale $t_{2-1}=(n\sigma\Sigma)^{-1}$ – with $n$ the cluster density,
$\sigma$ the cluster velocity dispersion, and $\Sigma$ the binary cross
section – is generally $>10^{8}$ larger than the binary inspiral timescale,
$t_{\rm insp}\propto a^{4}/(m_{1}m_{2}m_{\rm bin})$ (Peters & Mathews, 1963).
Moreover, the typical merger time for a binary with mass $m_{\rm
bin}<200{\rm~{}M}_{\odot}$ and separation $a_{\rm dec}$ is generally $t_{\rm
insp}<100$ yr, i.e. much smaller than the cluster crossing time, $t_{\rm
step}\sim 10^{5}$ yr.

Therefore, our procedure ensures reliabile results while reducing the
computational effort required to simulate the evolution of a binary with an
orbital period of minutes or hours.

The pre-merger stages of the merging binary orbits are reconstructed by
retrieving the orbital parameters at decoupling and integrating the orbit via
the Peters (1964) equations:

$\displaystyle\frac{{\rm d}a}{{\rm d}t}$ $\displaystyle=$
$\displaystyle-\frac{64}{5}\beta(m_{1},m_{2})\frac{F(e)}{a^{3}},$ (2)
$\displaystyle\frac{{\rm d}e}{{\rm d}t}$ $\displaystyle=$
$\displaystyle-\frac{304}{15}\beta(m_{1},m_{2})\frac{eG(e)}{a^{4}},$ (3)

with

$\displaystyle F(e)$ $\displaystyle=$
$\displaystyle(1-e^{2})^{-7/2}\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right);$
(4) $\displaystyle\beta(m_{1},m_{2})$ $\displaystyle=$
$\displaystyle(G^{3}/c^{5})m_{1}m_{2}(m_{1}+m_{2});$ (5) $\displaystyle G(e)$
$\displaystyle=$
$\displaystyle(1-e^{2})^{-5/2}\left(1+\frac{121}{304}e^{2}\right).$ (6)

Along with the orbital evolution we calculate the associated GW strain and
frequency (see e.g. Peters & Mathews, 1963; Kocsis et al., 2012; Arca Sedda et
al., 2021a).

Natal spins of stellar BHs can be assigned according to different
distribution, three of which are based on physical stellar model, namely the
“Geneva”, “MESA”, and “Fuller” models (see e.g. Kamlah et al., 2022), and four
are rather generic, namely zero-spins, uniform spin distribution, Gaussian
spin distribution with mean value $\chi=0.5$ and dispersion
$\sigma_{\chi}=0.2$, and Maxwellian distribution with dispersion
$\sigma_{\chi}=0.2$.

In this work, whenever spins are taken into account during the simulation we
assume a Gaussian distribution with $\chi=0.5$ for stellar BHs, whilst for
IMBHs we decide on a case by case basis, depending on the IMBH formation
scenario (see paper AS-II).

Compact binary merger products are assigned a final mass and spin calculated
via numerical relativity fitting formulae (Jiménez-Forteza et al., 2017; Arca
Sedda et al., 2020b) and a relativistic recoil, generated by asymmetric GW
emission (Campanelli et al., 2007; Lousto & Zlochower, 2008; Lousto et al.,
2012), expressed via the following relation:

$\displaystyle\vec{v}_{\rm GW}=$ $\displaystyle
v_{m}\hat{e}_{\bot,1}+v_{\bot}(\cos\xi\hat{e}_{\bot,1}+\sin\xi\hat{e}_{\bot,2})+v_{\parallel}\hat{e}_{\parallel},$
(7) $\displaystyle v_{m}=$ $\displaystyle A\eta^{2}\sqrt{1-4\eta}(1+B\eta),$
(8) $\displaystyle v_{\bot}=$
$\displaystyle\displaystyle{\frac{H\eta^{2}}{1+q_{\rm
BBH}}}\left(S_{2,\parallel}-q_{\rm BBH}S_{1,\parallel}\right),$ (9)
$\displaystyle v_{\parallel}=$
$\displaystyle\displaystyle{\frac{16\eta^{2}}{1+q_{\rm
BBH}}}\left[V_{11}+V_{A}\Xi_{\parallel}+V_{B}\Xi_{\parallel}^{2}+V_{C}\Xi_{\parallel}^{3}\right]\times$
(10) $\displaystyle\times\left|\vec{S}_{2,\bot}-q_{\rm
BBH}\vec{S}_{1,\bot}\right|\cos(\phi_{\Delta}-\phi_{1}).$

Here, $\eta\equiv q_{\rm BBH}/(1+q_{\rm BBH})^{2}$ is the symmetric mass
ratio, $\vec{\Xi}\equiv 2(\vec{S}_{2}+q_{\rm BBH}^{2}\vec{S}_{1})/(1+q_{\rm
BBH})^{2}$, and the subscripts $\bot$ and $\parallel$ mark the perpendicular
and parallel directions of the BH spin vector ($\vec{S}$) with respect to the
direction of the binary angular momentum. We assume $A=1.2\times 10^{4}$ km
s-1, $B=-0.93$, $H=6.9\times 10^{3}$ km s-1, and $\xi=145^{\circ}$ González et
al. (2007); Lousto & Zlochower (2008), $V_{11}=3677.76$ km s-1, and
$V_{A,B,C}=(2.481,1.793,1.507)\times 10^{3}$ km s-1. The vector $\vec{\Delta}$
is defined as $\vec{\Delta}\equiv(M_{a}+M_{b})^{2}(\vec{S}_{b}-q_{\rm
BBH}\vec{S}_{a})/(1+q_{\rm BBH})$. The angle between the direction of the
infall at merger and the in-plane component of $\vec{\Delta}$, i.e.
$\phi_{\Delta}$, is drawn from a uniform distribution, while
$\phi_{1}=0-2\pi$, which represents the phase of the binary, is extracted
between the two limiting values according to a uniform distribution. In
Nbody6++GPU , the user can decide to set the GW recoil to zero or to a fixed
value, or to calculate it self-consistently via Eqs. 7 and 10, in which case
the kick is assigned to the remnant and the resulting energy correction is
included in a similar way as it is done for natal BH kicks.

As described in detail in paper AS-II, in this paper series we adopt a
simplified approach to investigate the impact of GW recoil in the simulations,
owing to the fact that the relatively small sample of mergers does not enable
us to filter out the inevitable stochastic effect of the BH spin directions
and amplitudes on the kick amplitude. The approach consists of three steps.
First, we run all simulations assuming no GW recoil. Second, for each merger
event in each simulation we evaluate the GR recoil assuming different
distribution for BHs natal spins and we determine whether the remnant is
likely to be retained or not in the cluster. Third, if a BH undergoes $n$
mergers in a simulation with zero GW kick, we restart the simulation shortly
before the $n$-th merger event and enable GW kicks assuming a spin for the
merging components that depends on the BH formation history. This enables us
to verify whether the remnant can be retained in the cluster and eventually
merge again in a $n+1$th merger generation. In paper AS-II, we have shown that
this approach permits us to highlight the fact that even when GW kicks are not
taken into account, Newtonian dynamics is sufficient to eject all BH remnants
from the parent cluster via strong binary-single encounters.

We note that none of the mergers with component masses $<100{\rm~{}M}_{\odot}$
undergo multiple mergers. This suggests that even adopting zero GW recoils may
have a negligible impact on the formation of compact binary mergers with mass
$<100{\rm~{}M}_{\odot}$.

##  3 Results

###  3.1 The population of black hole binary mergers in Dragon-II clusters

In this section we describe the main results of our simulations, focusing on
the population of compact binary mergers. Table 1 summarizes the main
properties of Dragon-II clusters and their compact objects.

The population of BHs formed in Dragon-II models and, in general, in star
clusters likely suffers both the effects of single and binary stellar
evolution and stellar dynamics. To highlight this aspect we show in Figure 1,
for the models with $R_{\rm HM}=0.8$pc and $N=300$k, the so-called initial to
final mass relation (IFMR) that links the masses of compact objects and their
stellar progenitors. The plot is dissected into BHs with a progenitor
initially single or in a primordial binary system.

The population of BHs forming from single stars generally follows the
expectations of the adopted stellar evolution recipes (e.g. see Figure B1 in
Kamlah et al., 2022). Deviation from the general trend owes to initially
single stars that got caught in a pair and underwent mass-transfer. The IFMR
of BHs formed from stars in primordial binaries is more complex, being
characterised, for example, by BHs in the upper mass-gap with masses in the
range $40.5-80{\rm~{}M}_{\odot}$. This highlights the crucial role of binary
stellar evolution and dynamics in sculpting the population of BHs in star
clusters (e.g., see also Figure 2 in Di Carlo et al., 2019).

Figure 1: Zero age main sequence mass and final mass of stellar BH progenitors. We dissect the population into progenitors that are initially single (top panel, red circles) and those that are in a primordial binary (bottom panel, light blue triangles). $N_{*}$ | $M_{c}$ | $R_{h}$ | $f_{b}$ | $N_{\rm sim}$ | $T_{\rm rlx}$ | $T_{\rm seg}$ | $T_{\rm sim}$ | $N_{\rm GW,in}$ | $N_{\rm GW,out}$ | $M_{\rm max}$ | $M_{\rm max,fin}$ | $N_{>30}$ | $N_{>40}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$10^{3}$ | $10^{5}{\rm~{}M}_{\odot}$ | pc |  |  | Myr | Myr | Myr |  |  | ${\rm~{}M}_{\odot}$ | ${\rm~{}M}_{\odot}$ |  |
120 | 0.7 | 1.75 | 0.05 | 2 | 99 | 2.1 | 2379 | 2326 | 0 | 2 | 2 | 0 | 64 | 76 | 25 | 34 | 0 | 2 | 0 | 0
300 | 1.8 | 1.75 | 0.05 | 2 | 142 | 2.7 | 1196 | 1422 | 0 | 2 | 2 | 2 | 69 | 77 | 40 | 40 | 13 | 13 | 5 | 1
1000 | 5.9 | 1.75 | 0.05 | 2 | 233 | 3.4 | 207 | 194 | 1 | 1 | 4 | 4 | 81 | 146 | 52 | 70 | 149 | 169 | 72 | 85
120 | 0.7 | 1.75 | 0.2 | 2 | 99 | 2.1 | 1710 | 1540 | 2 | 2 | 0 | 2 | 232 | 81 | 38 | 28 | 2 | 0 | 0 | 0
300 | 1.7 | 1.75 | 0.2 | 2 | 142 | 2.7 | 519 | 793 | 1 | 0 | 7 | 5 | 92 | 77 | 65 | 47 | 26 | 26 | 8 | 14
600 | 3.5 | 1.75 | 0.2 | 2 | 189 | 3.4 | 205 | 126 | 0 | 0 | 2 | 5 | 87 | 144 | 59 | 84 | 95 | 103 | 45 | 65
120 | 0.7 | 0.80 | 0.2 | 2 | 30 | 0.7 | 1154 | 1201 | 4 | 3 | 4 | 2 | 120 | 132 | 21 | 27 | 0 | 0 | 0 | 0
300 | 1.7 | 0.80 | 0.2 | 2 | 44 | 0.8 | 307 | 309 | 1 | 0 | 1 | 0 | 93 | 107 | 40 | 43 | 15 | 11 | 2 | 2
120 | 0.7 | 0.47 | 0.2 | 2 | 14 | 0.3 | 1149 | 530 | 2 | 2 | 3 | 1 | 350 | 92 | 50 | 30 | 1 | 0 | 1 | 0
300 | 1.7 | 0.47 | 0.2 | 1 | 20 | 0.4 | 148 | - | 4 | - | 3 | - | 245 | - | 48 | - | 22 | - | 9 | -
Table 1: Col. 1-4: initial number of stars, cluster mass and half-mass radius,
primordial binary fraction. Col. 5: number of indipendent realisations. Col.
6-7: initial half-mass relaxation and segregation time. Col. 8: simulated
time. Col. 9-10: number of compact object mergers inside the cluster. Col. 11:
maximum BH mass during the simulation. Col. 12: maximum BH mass at the end of
the simulation. Col. 13-14: number of BHs with a mass $m_{\rm
BH}>30{\rm~{}M}_{\odot}$ or $>40{\rm~{}M}_{\odot}$ at the last simulation
snapshot.

####  3.1.1 Component masses and formation channels

The population of compact binary mergers in Dragon-II consists in 75 BH-BH, 2
NS-BH, and 1 WD-BH. Among BH-BH mergers, 45 involve two BHs below the PPISN
maximum mass ($m_{\rm BH}<40.5{\rm~{}M}_{\odot}$), 12 involve two mass-gap
BHs, and 21 involve one BH below the gap and a mass-gap BH. Six BH-BH mergers
involve a primary with a mass $m_{\rm BH,1}=(5.4-7.1){\rm~{}M}_{\odot}$ and a
companion with mass $m_{\rm BH,2}=(2.55-3.6){\rm~{}M}_{\odot}$, i.e. just
above the threshold separating NSs and BHs in our models. All these low-mass
mergers are in primordial binaries. As discussed in paper AS-II, the BHs in
the upper-mass gap mostly form in a star-BH accretion event, either by purely
dynamical interactions or stellar evolution. We stress that, throughout our
models, we assume that a fraction $f_{c}=0.5$ of the star mass is accreted
onto the BH during an accretion event (for a discussion about the impact of
$f_{c}$ on simulations, see Rizzuto et al., 2021).

When GW recoil are “switched off” 4 mergers involve a second or third
generation BH, i.e. which underwent one or two previous mergers. The inclusion
of GW recoil reduces the number of total mergers to 74. For a detailed
discussion about the impact of GW recoil, see paper AS-II.

Figure 2 shows the component masses and mass ratio of Dragon-II mergers and of
mergers observed during the first three LVK observation campaign, collected in
the so-called GWTC-3 catalogue (The LIGO Scientific Collaboration et al.,
2021b, a). The plot includes all mergers occurring in the cluster or outside
the cluter after being ejected via dynamical interactions considering zero GW
kicks.

This plot illustrates the wealth of information hid in the Dragon-II star
clusters: we find mergers in the upper-mass gap, IMBHs222In this work we set a
mass treshold of $M_{\rm IMBH,min}=100{\rm~{}M}_{\odot}$ to discern between
BHs and IMBHs., repeated mergers, and in a handful cases also BHs merging with
either a NS or a WD. Interestingly, we find that mergers occurring inside the
cluster are characterised by a primary with mass $m_{\rm
BH,1}>30{\rm~{}M}_{\odot}$ and a companion with a mass in the range $m_{\rm
BH,2}=(20-50){\rm~{}M}_{\odot}$. Conversely, mergers occurring outside the
cluster — or ejected mergers — are characterised by a mass-ratio $q>0.6$ and a
primary mass typically $m_{\rm BH,1}<40{\rm~{}M}_{\odot}$.

Figure 2: Primary mass (x-axis) and mass-ratio (y-axis) of BH mergers in
Dragon-II simulations, occurring inside the cluster (points) or after
dynamical ejection (diamonds). The colour map identifies the mass of the
secondary. The dotted lines corresponds to a companion mass of $m_{\rm
BH,2}=3,~{}5,~{}10,~{}30,~{}50,~{}100{\rm~{}M}_{\odot}$. Mergers with a
primary mass on the right side of the black dashed line produce an IMBH. The
red shaded area roughly identifies the upper-mass gap region. The grey area
identify mergers containing a BH and another type of compact object: either a
BH in the putative lower mass-gap, a NS or a WD. Shaded grey points represent
observed BH mergers from the GWTC-3 catalogue (The LIGO Scientific
Collaboration et al., 2021b).

The number of mergers occurring inside the cluster (31) is comparable to that
of binaries that merge after being ejected from the cluster (47), thus
suggesting that in-cluster mergers can made-up the $40\%$ of the total merger
population. Among all of them, 27 are from primordial binaries (3 inside, 24
ejected), whilst 51 (28 inside, 23 ejected) are from dynamical binaries.

Figure 3 shows the primary and companion mass of mergers originated from
primordial, dynamical, or mixed binaries, with the latter identifying binary
mergers in which at least one component was originally in a primordial binary.
The plot exhibits some interesting features: 1) mergers from primordial
binaries tend to have nearly equal-mass components, 2) purely dynamical
mergers have masses that occupy a tight region of the plane with $m_{\rm
BH,1}=(20-50){\rm~{}M}_{\odot}$ and $m_{\rm BH,2}=(20-40){\rm~{}M}_{\odot}$,
3) mergers with one component previously in a primordial binary are
characterised by a heavy primary, $m_{\rm BH,1}>40{\rm~{}M}_{\odot}$, and a
heavy companion, $m_{\rm BH,2}>20{\rm~{}M}_{\odot}$. A similar trend is
observed in recent $N$-body simulations tailored to relatively light star
clusters, i.e. with mass $<8,000{\rm~{}M}_{\odot}$ (Torniamenti et al., 2022).

As deeply discussed in paper AS-II, the crucial role of primordial binary
dynamics is highlighted by the fact that all the IMBHs in Dragon-II clusters
but one have an ancestor that was member of a primordial binary, regardless of
the IMBH formation scenario.

Figure 3: Masses of the primary and companion merging BHs from primordial
binaries (yellow diamonds), dynamical binaries (red circles), and binaries
with at least one component being a former primordial binary member (blue
squares).

Dynamics and binary stellar evolution deeply impact also the properties of
stellar-size mergers. For example, “dynamical” and “primordial” mergers occupy
two well separated regions of the primary mass - mass ratio plane.

The vast majority of primordial binary mergers occupy a region delimited by
$q>0.6$ and $m_{1}=(5-40){\rm~{}M}_{\odot}$, with the mass-ratio weakly
increasing at increasing the primary mass: note that for $m_{1}\lesssim
15{\rm~{}M}_{\odot}$ mergers have mass-ratio $q=0.6-1$, whilst mergers with a
heavier primary have mass ratio $q>0.85$.

Dynamical mergers, instead, form in the right hand-side of Figure 3, generally
at $m_{1}>40.5{\rm~{}M}_{\odot}$. In this case, the mass ratio decreases with
the primary mass as expected from the mass function limit, with companion
masses in the range $m_{2}=(30-50){\rm~{}M}_{\odot}$.

We can identify three relatively well separated regions: low BH masses
($m_{\rm BH,1}<15{\rm~{}M}_{\odot}$) and widely distributed mass ratio
($q=0.6-1$) dominated by primordial binary mergers, BH masses in the range
$m_{\rm BH,1}=(15-40.5){\rm~{}M}_{\odot}$ and high mass ratios ($q>0.9$)
dominated by primordial binary mergers, and heavy BH primaries ($m_{\rm
BH,1}>40.5{\rm~{}M}_{\odot}$) with relatively massive companions
($m_{2}=30-50{\rm~{}M}_{\odot}$) dominated by dynamical mergers.

In Dragon-II clusters, most binaries merging outside the cluster originate
from primordial binaries and their ejection is typically triggered by the BH
natal kick. However, all ejected mergers with component masses
$m_{1,2}>30{\rm~{}M}_{\odot}$ have a dynamical origin, owing to the adopted
stellar evolution recipes.

We note that, given the limited simulation time, the population of mergers in
Dragon-II clusters may lack some element that could form later in the cluster
life, beyond several relaxation times. These late mergers would unavoidably
have a dynamical origin, or at most a "mixed" origin, because all BHs formed
in primordial binaries undergo a binary exchange or have been ejected in
Dragon-II clusters over the simulated time. Moreover, late mergers will likely
have smaller masses compared to those shown in Figure 3. This is mostly due to
the BH-burning process, by which the average BH mass decreases over time (see
e.g. paper AS-I; Rodriguez et al., 2015; Chatterjee et al., 2017). As a
consequence, some BH mergers forming at late time may have properties similar
to the primordial binary mergers shown in Figure 3.

Figure 4 shows the mass distribution of the primary BH in Dragon-II mergers,
dissected into in-cluster/ejected mergers and primordial/dynamical ones.
Ejected binaries dominate the $m_{\rm BH,1}\lesssim 20{\rm~{}M}_{\odot}$ mass
range, whilst at larger primary masses their number and distribution is
similar to that of in-cluster mergers. Dynamical mergers completely dominate
the population of mergers with $m_{\rm BH,1}>20{\rm~{}M}_{\odot}$, while
primordial mergers dominate the population of lighter mergers. Noteworthy, we
see that the primary mass distribution for Dragon-II mergers nicely overlap
with the sample of mergers in the GWTC-3 catalogue, i.e. the catalogue of BBH
mergers detected by the LVK collaboration (The LIGO Scientific Collaboration
et al., 2021a). However, a thorough comparison between modelled and observed
mergers would require to take into account observation biases (see e.g.
Fishbach & Holz, 2017; Arca Sedda et al., 2020b; Arca Sedda et al., 2023). For
this reason, we also overlay to our data the cosmic BH mass distribution
inferred from GW detections.

Comparing models and observations can be crucial to assess the impact of
different formation channels on the population of BH-BH mergers (see e.g. Arca
Sedda & Benacquista, 2019; Arca Sedda et al., 2020b; Bavera et al., 2020;
Zevin et al., 2021; Arca Sedda et al., 2023; Mapelli et al., 2022). Our
Dragon-II models suggest, for example, that BH mergers developing in star
clusters could produce a substantial amount of mergers from primordial
binaries. The progenitor binary could, in some cases, suffer the impact of
dynamical interactions which may alter their orbital parameters. Nonetheless,
in most cases BH mergers from primordial binaries could represent "isolated
binary merger impostors", because they have properties typical of merging
binaries developing within the isolated formation scenario but form in a
dynamical environment. Taking into account the impact of these sources with a
sort of mixed formation channel is crucial to correctly quantify the role of
different formation channels in determining the shape of the mass distribution
of detected merging BHs (see also Arca Sedda et al., 2023).

Moreover, Dragon-II models highlights the role of dynamics in determining the
formation of BH mergers with masses inside, and beyond, the mass-gap,
supporting and complementing previous works on the topic based either on
smaller, or lower-density, $N$-body cluster models and Monte Carlo simulations
(e.g. Kremer et al., 2020; Di Carlo et al., 2020a; González et al., 2021;
Rizzuto et al., 2022; Banerjee, 2022).

Figure 4: Mass distribution of primary BH mergers from the GWTC-3 catalogue
(The LIGO Scientific Collaboration et al., 2021a) (orange straight line) and
from Dragon-II simulations (filled light blue steps). The BH mass distribution
inferred from LVK data is overlaid to the histograms (black line), with the
shaded area ancompassing the $90\%$ credible level. Top panel: Simulated
mergers are dissected into those occurring inside the cluster (dotted black
line) and after dynamical ejection (dashed black line). Bottom panel:
Simulated mergers are dissected into those forming from primordial binaries
(dotted black line) or via dynamical interactions (dashed black line).

####  3.1.2 Delay times

The delay time of Dragon-II mergers ($t_{\rm GW}$), defined as the time
elapsed from the beginning of the simulation to the binary merger, is rather
peculiar. As show in Figure 5, it exhibits three peaks at $t_{{\rm
GW}}\simeq(0.5,~{}1.5,~{}10)$ Gyr. However, when the delay time is normalised
to the initial half-mass relaxation time ($t_{\rm rlx}$) of the cluster, the
overall $t_{\rm GW}$ nicely distribute around a main peak located at $t_{\rm
GW}/t_{\rm rlx}\simeq 8-30$. The exact location of the peak depends on the
definition of $t_{\rm rlx}$. For the sake of clarity, in the plot we use three
different expressions of $t_{\rm rlx}$ taken from Gatto et al. (2021) (GR21),
Antonini & Rasio (2016) (AR16), or Rizzuto et al. (2021) (RN21).

The three peaks that appear in the $t_{\rm GW}$ distribution find a clear
explanation looking at the $t_{\rm GW}/t_{\rm rlx}$ distribution. In fact, the
first peak at $t_{\rm GW}=500$ Myr corresponds to mergers happening in
simulations with $t_{\rm rlx}=50-100$ Myr, whilst the second peak corresponds
to mergers occurring in clusters with a longer relaxation time (see Table 1).
This interesting features suggests, on the one hand, that the delay time
depends intrinsically on the cluster initial properties, as they determine the
relaxation time, and, on the other hand, that dynamical processes operate in a
similar way over a relaxation time regardless of the cluster structure. The
third peak, instead, corresponds to ejected binaries that merge outside the
cluster, which are mostly products of primordial binaries ejected via SN
explosion during the formation of one of the BHs in the pair.

Figure 5: Top panel: delay time distribution for Dragon-II mergers. Bottom
panel: as in top panel, but the time is normalised to the initial cluster
relaxation time calculated following Gatto et al. (2021) (red filled steps),
Binney & Tremaine (2008) (red dashed steps), or Spitzer (1987) (red dotted
steps).

####  3.1.3 Eccentricities

One intriguing question that arose since the first detection of GWs is whether
it is possible to untangle different formation channels in the observed
population of BH mergers. Among all the parameters at play, the orbital
eccentricity could represent the key to answer this question. Broadly
speaking, in fact, most BH mergers forming via binary stellar evolution are
expected to feature a negligible eccentricity close to merger, either because
the BBH progenitor undergoes common envelope, which shrinks and circularise
the orbit, or because the BBH separation is initially so large that GW
emission circularise the orbit before the merger. Binaries developing in star
clusters, instead, can form with high eccentricity and sufficiently small
separation that the merger occurs on a timescale shorter than GW
circularisation.

At the lowest-order level, binaries merging in galactic fields, often called
isolated binary mergers, are expected to be circular GW sources, whilst at
least some of those developing in star clusters and galactic nuclei, named
dynamical mergers, can preserve a significant eccentricity (i.e. $e>0.1$) when
entering the typical frequency bands of GW detectors.

This simplistic division between isolated and dynamical binaries does not take
into account several layers of complication. For example, it is well known
that star clusters and stellar nurseries may contain a large fraction of
binaries, especially among the population of massive stars, where the
percentage of paired stars attains values as large as $50-100\%$ (Sana et al.,
2012; Moe & Di Stefano, 2017). If primordial binaries evolve on a timescale
shorter than the typical timescale of dynamical interactions, star cluster
could harbor a sub-population of compact binary mergers with properties pretty
similar to those forming in galactic fields, e.g. low eccentricities or
peculiar component masses and mass-ratios.

With up to $33\%$ of stars initially paired, Dragon-II simulations offer us
the possibility to search for differences between mergers forming entirely via
dynamics and those forming from the evolution of primordial binaries. Figure 6
shows the semimajor axis and eccentricity of all BH-BH mergers in Dragon-II
clusters calculated at the moment of decoupling, i.e. when the GW emission
starts dominating over dynamical perturbations. The plot dissects the Dragon-
II population of BH mergers into those coming from the evolution of primordial
binaries, those assembled purely via dynamical interactions, and those
involving at least one component that was former member of a primordial
binary. The population of dynamical and mixed binaries seem to follow two
different sequences, although the low statistics make hard to understand
whether they actually exist. The population of nearly circular primordial
binaries is evident. These mergers can be considered mimickers of the field
merger population, and constitute the $33\%$ of the whole population of
Dragon-II mergers. Only two of the primordial binaries exhibit a significant
eccentricity and a relatively small separation.

The first is a NS-BH binary, we postpone a discussion about this specific
source to the next subsection.

The second one involves two low-mass BHs, with masses $m_{\rm
BH1,2}=(7.1+2.55){\rm~{}M}_{\odot}$ and eccentricity $e=0.997$. The progenitor
of this merger was a binary that underwent a common envelope phase first,
after which the first BH formed, and later undergo Roche lobe overflow, at the
end of which also the second BH forms and receives a small kick ($\sim 3$
km$/$s) that triggers the eccentricity increase.

Figure 6: Semimajor axis and eccentricity of compact binary merger at
decoupling for primordial (yellow diamonds), dynamical (red points), and mixed
(blue squares) mergers. Labels indicate particularly interesting primordial
mergers with a significant eccentricity at decoupling.

As the binary shrinks and circularises because of GW emission, its frequency
will increase. Therefore, a first step to determine whether a binary merger
can appear eccentric in the sensitivity band of a specific GW detector
requires to compare the binary eccentricity and the corresponding GW
frequency. We show in Figure 7 the characteristic strain - frequency evolution
for all mergers in our sample, assuming that they are located at a redshift
$z=0.05$, i.e. at a luminosity distance of 230 Mpc. To calculate the GW strain
of Dragon-II sources we follows the formalism laid out in Peters & Mathews
(1963); Peters (1964) and the implementation described in Arca Sedda et al.
(2021a) (see Eqs. 30-39). The GW signal from simulated mergers is overlaid to
the sensitivity curves of current and future ground-based and space-based
detectors like LIGO (Acernese et al., 2015), Einstein Telescope (ET Punturo et
al., 2010; Branchesi et al., 2023), DECIGO (Kawamura et al., 2011; Isoyama et
al., 2018; Kawamura et al., 2021), and LISA (Amaro-Seoane et al., 2013,
2022a). The plot highlights how the eccentricity drops as the binary sweeps
across different frequency bands.

Figure 7: Evolution of the binary GW strain as a function of the frequency for
all mergers in Dragon-II simulations, assuming that the sources are located at
a redshift $z=0.05$. The colormap identifies the binary eccentricity along its
orbit. Dotted lines represent sensitivity curves for different detectors, from
left to right: LISA, DECIGO, ET, and LIGO.

The top panel in Figure 8 shows the fraction of mergers with eccentricity
above a given threshold calculated when Dragon-II mergers sweep through five
frequency bands centered in $f_{\rm band}=10^{-3}-10^{-2}-10^{-1}-1-10$ Hz,
i.e. the typical sensitivity bands of space-borne detectors like LISA
($<10^{-2}$ Hz), mid-frequency detectors like DECIGO ($10^{-2}-1$ Hz), and
ground-based detectors like LIGO-Virgo-Kagra or the Einstein Telescope ($>1$
Hz). The plot highlights the fact that around $20-40-5\%$ of all mergers
appear eccentric, i.e. $e>0.1$, while moving through the
$f=10^{-3}-10^{-1}-10^{1}$ Hz frequency bands, respectively. Clearly, the
detectability of these mergers depend on many parameters, among which the
location of the merger and the detector properties. Nonetheless, the plot
makes apparent the importance of future deci-Hz detectors in placing
constraints on the population of eccentric BBH mergers. Moreover, comparing
models with future observations will help to quantify the impact of star
cluster dynamics on the cosmic population of merging BHs.

Noteworthy, the eccentricity carries information about the formation history
of the merger. For example, we find that all mergers with an eccentricity
$e>0.1$ in both the $0.05-1$ Hz and $1-10$ Hz frequency bands occur inside the
cluster. The number of eccentric binaries doubles in the $10^{-2}-1$ Hz
frequency band, but these eccentric binaries appear almost circular while
reaching the ground-based detector band, explaining why it is more likely to
find a Dragon-II merging binary with significant eccentricity while sweeping
through the deci-Hz band.

Any binary merger will spend some time in the detector band before merging. In
order to characterise the evolution of the eccentricity as the binary
inspirals, we calculate the average binary eccentricity weighted with the time
to the inspiral, i.e. $\langle e\rangle=\int_{0}^{t_{\rm
merg}}edt/\int_{0}^{t_{\rm merg}}dt$. Practically, we measure the binary
eccentricity in subsequent time bins from the time of decoupling to the time
of merger and weight it with the remaining time to the merger. This quantity
is shown for all mergers in the bottom panel of Figure 8, along with the
evolution of the eccentricity as a function of the peak frequency (e.g.
O’Leary et al., 2009)

$\displaystyle f_{p}=$ $\displaystyle 0.29{\rm
Hz}\left(\frac{m_{1}+m_{2}}{30{\rm~{}M}_{\odot}}\right)^{1/2}\left(\frac{a}{50{\rm
R}_{\odot}}\right)^{-3/2}\times$ $\displaystyle\times{\rm
ceil}\left[1.15\frac{(1+e)^{1/2}}{(1-e)^{3/2}}\right],$ (11)

The step-like behaviour of the $e-f_{p}$ is due to the ceil function in
Equation 11, which returns the nearest integer larger than the function
argument. The majority of in-cluster mergers clearly show an average
eccentricity $\langle e\rangle>0.7$ across the whole 0.01-100 Hz frequency
spectrum, whilst ejected mergers preserve a moderate eccentricity $\langle
e\rangle<0.4$ in the $f<1$ mHz band. This suggests that GW detectors operating
in different bands can probe different sub-populations of mergers forming in
dense star clusters, with high-frequency detectors being more suited to
observe short-lived, highly eccentric mergers occurring inside star clusters,
and low-frequency detectors more suited to observe GW sources merging outside
their parent cluster.

Figure 8: Top panel: Fraction of mergers with eccentricity above a given
value. Different colors correspond to different frequency bands. Bottom panel:
eccentricity evolution as a function of the GW frequency for all mergers in
Dragon-II models (grey lines) and average value of the eccentricity for
mergers occurring inside the cluster (red points) or after ejection (white
points).

###  3.2 Exotic mergers

Despite the relatively small simulation grid, we also find some exotic
mergers: a dynamical WD-BH and 2 NS-BH mergers, one dynamical and one from a
primordial binary. The three mergers occur in the most dense simulations in
our sample: the WD-BH merger occurs in a simulation with $N=120$k, $R_{\rm
HM}=0.47$pc, the dynamical BH mergers develop in a simulation with $N=300$k,
$R_{\rm HM}=0.47$pc, and the one forming from a primordial binary in a
simulation with $N=120$k, $R_{\rm HM}=0.8$pc. This type of mergers are
particularly rare in star clusters, because dynamical exchanges favor the
replacement of the light component with another BH. Given their rarity, we
discuss in the following the details of the formation and evolution of these
interesting sources.

####  3.2.1 White dwarf - black hole mergers: implications for low-mass X-ray
binaries

The WD-BH binary consists in a BH with mass $m_{\rm BH}=23.1{\rm~{}M}_{\odot}$
and a carbon-oxygen white dwarf (COWD) with mass $m_{\rm
WD}=1.18{\rm~{}M}_{\odot}$. Initially, dynamical interactions pair the BH with
the WD progenitor, a MS star with mass $m_{\rm MS,pro}=4.89{\rm~{}M}_{\odot}$.
The two objects undergo common envelope during the late AGB phase of the
companion, at the end of which the star turns into a WD, after $\sim 105$ Myr.
The resulting WD-BH binary has an "initial" eccentricity of $e=0.2$ and period
of $900$ days. The binary undergoes a series of strong scatterings that cause
a rather chaotic variation of the binary semimajor axis and a systematic
increase of the eccentricity from $e=0.6$ up to $e=0.99994930$ after 135 Myr,
corresponding to $\sim 4$ relaxation times. At this stage, GW emission becomes
sufficiently effective to drive binary coalescence. Figure 9 shows the time
variation of the WD-BH binary semimajor axis and eccentricity before
coalescence.

Figure 9: Time evolution of the semimajor axis (top panel) and eccentricity
(bottom panel) of a dynamical WD-BH merger in Dragon-II clusters.

The WD Roche lobe is larger than the BH innermost stable circular orbit, hence
the WD will likely undergo disruption and start feeding the BH, possibly
evolving into a low-mass X-ray binary. In these regards, it is interesting
noting the observation of a ultracompact X-ray binary in the galactic cluster
47 Tuc (Miller-Jones et al., 2015; Bahramian et al., 2017), likely comprised
of a COWD and a BH (Bahramian et al., 2017), with the BH being probably
heavier than $m_{\rm BH}>9{\rm~{}M}_{\odot}$ (Church et al., 2017). Our model
confirms the possibility to form such type of low-mass X-ray binary via
interactions of stars and BHs in a dense cluster, even in a relatively short
time ($t<200$ Myr). Ultimately, the binary shrinkage driven by GW emission
will disrupt the WD and the mass falling onto the BH could possibly power jets
that can give rise to transients with peak energy $10^{47}-10^{50}$ erg-1 and
duration of a few minutes (Fernández et al., 2019), giving rise to a tidal
disruption event (TDE). Despite this source is the only one undergoing
coalescence, we find a total of 50 WD-BH binaries by the end of the simulation
in all Dragon-II clusters. None of them have orbits such to trigger a TDE
within a Hubble time, unless a strong interaction with some cluster member
pushes the WD on an extremely eccentric orbit. Pushing the orbit to at least
$e>0.9999$($0.99999$) would lead to 1(26) further WD-BH mergers. Note that the
eccentricity value required to trigger a WD TDE may seem extreme, but it is
comparable to the eccentricity achieved by the WD-BH merger, hence testifying
that it is possible to reach such extreme eccentricity values in Dragon-II
clusters.

####  3.2.2 Neutron star - black hole mergers: implications for multimessenger
astronomy

Concerning NS-BH binaries, we find two mergers, one of dynamical origin and
the other forming from the evolution of a primordial binary.

The dynamical NS-BH has a NS with mass $m_{\rm NS}=1.28{\rm~{}M}_{\odot}$ and
a BH with mass $m_{\rm BH}=14.96{\rm~{}M}_{\odot}$. The BH, whose progenitor
had a mass $m_{\rm MS}=26.7{\rm~{}M}_{\odot}$, undergoes a series of chaotic
interactions with a primordial binary containing the NS and its companion,
which eventually leads to the merger. When the binary decouples from the
cluster dynamics, it has a semimajor axis of $a=0.33$ AU and an eccentricity
$e=0.99778817$, corresponding to a GW peak frequency $f_{\rm GW}=0.01$ Hz.

After decoupling, the binary evolution is completely dominated by GW emission
and the variation of its orbital parameters can be described, at first order,
via the Peters (1964) formalism. We find that as the binary sweeps through the
$0.01-0.5-1-10$ Hz GW frequency band the NS-BH merger has a residual
eccentricity of $e_{\rm NSBH}=0.99779-0.9974-0.21-0.02$, thus in principle
measurable with future GW detectors, especially with those operating in the
deci-Hz frequency band. The chirp mass of this merger, $\mathcal{M}_{\rm
chirp}=3.4{\rm~{}M}_{\odot}$, is typical of dynamically assembled NS-BH
mergers (Arca Sedda, 2020a; Rastello et al., 2020; Ye et al., 2020; Arca
Sedda, 2021), but hard to produce with isolated binary evolution (Giacobbo &
Mapelli, 2018; Zevin et al., 2020), although this strongly depends on the
adopted stellar evolution scheme (Broekgaarden et al., 2021).

The primordial NS-BH binary merger, instead, forms from a primordial binary
with initial mass components $m_{1,2}=(26.3+18.7){\rm~{}M}_{\odot}$ and
evolves through a common envelope phase initiated by the primary, which
eventually forms the BH. Shortly after, the binary undergoes a second common
envelope phase and eventually the companion evolves into a NS. Eventually, the
merging binary consists of a BH with mass $m_{\rm BH}=5.6{\rm~{}M}_{\odot}$
and a NS mass $m_{\rm NS}=1.88{\rm~{}M}_{\odot}$. Note that these properties,
are intriguingly similar to GW200115, a GW source detected by the LVK during
the O3 observation campaign, which was characterised by a BH with $m_{\rm
BH}=5.7^{+1.8}_{-2.1}{\rm~{}M}_{\odot}$ and a NS with $m_{\rm
NS}=1.5^{+0.7}_{-0.3}{\rm~{}M}_{\odot}$. When the NS forms, the common
envelope has shrunk the binary from $2.5$ R⊙ to $a=0.6$ R⊙, whilst the natal
kick imparted at formation onto the NS causes an enhancement of the
eccentricity from nearly zero to $e=0.57$. The new orbital parameter as such
that GW emission dominates over dynamics and the binary coalesces in $\sim
7\times 10^{4}$ yr. At decoupling, the binary peak frequency is $f_{{\rm
GW}}\sim 2$ mHz, right in the middle of LISA sensitivity band.

The development of a NS-BH binary merger from a primordial binary in a dense
star cluster highlights the impact of primordial binaries in contributing to
the population of mergers with properties similar to those forming in
isolation, making quite hard untangling their actual origin.

Merging NS-BH binaries are thought to be possible progenitors of several
electromagnetic (EM) transients, like short Gamma Ray Bursts (sGRBs) (e.g. Lee
& Ramirez-Ruiz, 2007) and kilonovae (e.g. Metzger & Berger, 2012). A basic
condition for the possible development of an EM transient is that part of the
NS material remains bound to the BH, forming a disk. The fraction of NS mass
in the disk depends on several quantities, among which the BH-to-NS mass ratio
$m_{\rm BH}/m_{\rm NS}$, the BH spin $\chi$, and the NS compactness $C\equiv
Gm_{\rm NS}/c^{2}R_{\rm NS}$ (Foucart, 2012). As numerical simulations have
shown, in general the larger the $m_{\rm BH}/m_{\rm NS}$ the larger the
minimum spin required for the NS material to form a disk around the BH, and
the larger the spin the larger the amount of matter bound to the BH (Foucart,
2012). Depending on the orbital parameters, the BH tidal field can tear apart
the NS before it enters the BH event horizon, provided that the NS tidal
radius

$R_{\rm tid}=R_{\rm NS}\left(\frac{3m_{\rm BH}}{m_{\rm NS}}\right)^{1/3},$
(12)

exceeds the BH innermost stable circular orbit (ISCO), which for a spinning BH
can be expressed as (Bardeen et al., 1972)

$R_{\rm ISCO}=\frac{Gm_{\rm BH}}{c^{2}}\left[3+Z_{2}-{\rm
sign}(\chi)\left[(3-Z_{1})(3+Z_{1}+2Z_{2})\right]^{1/2}\right],$ (13)

where $Z_{1,2}$ are functions of the BH adimensionless spin parameter $\chi$.
Whilst the condition $R_{\rm tid}/R_{\rm ISCO}<1$ implies that the merger has
no EM emission, the opposite does not ensure the EM counterpart detectability,
as it depends on the geometry of the merger with respect to the observer and
other possible observation biases.

In Dragon-II clusters, the dynamical NS-BH merger is characterised by $m_{\rm
BH}/m_{\rm NS}=11.7$ and compactness $C=0.19$ (assuming a NS radius of $10$
km). As shown in Figures 6-8 of Foucart (2012), the minimum BH spin required
for an accretion disk to form with a mass $10\%$ of the NS mass around such
type of binary is $\chi_{\rm BH}>0.98$. The BH formed in this binary did not
undergo any major interaction with stellar companions that could spin-up it
(Qin et al., 2019; Bavera et al., 2020; Belczynski et al., 2020). Hence, it is
possible that the BH formed with low-spin, according to the Fuller & Ma (2019)
model, hampering the formation of a massive accretion disk around the BH and
minimizing the probability for a EM counterpart to develop. The isolated NS-BH
merger, instead, is characterised by $m_{\rm BH}/m_{\rm NS}=2.98$ and
$C=0.27$. Even in this case, the spin required for an accretion disk to form
is $\chi_{\rm BH}>0.9$. The BH in this binary undergoes a RLO phase, which
could, in principle, spin-up the BH up to extremal values (e.g. Fragos &
McClintock, 2015), although this strongly depends on the stellar evolution
recipes and the binary properties (Qin et al., 2019; Bavera et al., 2020;
Belczynski et al., 2020).

The development of just 2 NS-BH mergers highlights how rare are these type of
objects on the one hand, and make any statistical analysis poor, on the other
hand. Nonetheless, the fact that the NS-BH mergers developed in Dragon-II
clusters seem to be unlikely to feature an EM counterpart supports the idea
that most NS-BH mergers proceed unseen in star clusters (Arca Sedda, 2020a).
For comparison, note that for isolated binaries typically $m_{\rm BH}\sim
12{\rm~{}M}_{\odot}$ and $m_{\rm NS}=1.6{\rm~{}M}_{\odot}$ (Broekgaarden et
al., 2021), which implies a minimum BH spin of $\chi_{\rm BH}\gtrsim 0.8$ to
permit the formation of a fairly massive (mass > $0.1m_{\rm NS}$) disk around
the BH (Fryer et al., 2012).

##  4 Discussion

###  4.1 The impact of natal spins on the properties of stellar black hole
mergers

Spin amplitude and mutual orientation at merger represent two possible
quantities that can help discerning whether a BBH merger results from isolated
stellar binary evolution or stellar dynamics (e.g. Arca Sedda & Benacquista,
2019; Arca Sedda et al., 2020b; Zevin et al., 2021; Arca Sedda et al., 2023;
Mapelli et al., 2022; Banerjee, 2022; Banerjee et al., 2023).

In order to explore the impact of different spin prescriptions on Dragon-II
mergers, we devise two models. The first model (hereafter STEV) assumes that
the spin is intrinsically related to the BH evolutionary pathways. For BHs
formed from single stellar evolution, we assume a negligible spin ($\chi_{\rm
BH}=0.01$) owing to efficient angular momentum transport triggered by the
Tayler-Spruit dynamo (Spruit, 2002; Fuller & Ma, 2019). For upper-mass gap BHs
formed from massive binary evolution we assume that final spins spans the
$\chi_{\rm BH}=0.8-1$ range (Qin et al., 2018; Bavera et al., 2020; Belczynski
et al., 2020; Schrøder et al., 2020). For BHs in primordial binaries, instead,
we assign to one BH a spin value of $\chi_{\rm BH}=0.01$ and to the other
$\chi_{\rm BH}=0.1-1$ (Qin et al., 2018; Bavera et al., 2020).

The second model (GAUS model) assumes, instead, that the spin distribution
follows a Gaussian distribution with mean $\bar{\chi}_{\rm BH}=0.5$ and
dispersion $\sigma_{\chi}=0.2$, regardless the BH past evolution, a case
possibly supported by the population of observed BH-BH mergers (The LIGO
Scientific Collaboration et al., 2021b).

In our analysis, we assume that the spin vectors in dynamical mergers are
isotropically distributed, whilst for primordial mergers we proceeds as
follows. We define an ad-hoc distribution function for the cosine of the angle
between the spin of the $i$-th binary component and the binary angular
momentum, $\theta_{i}$, such that (Arca Sedda et al., 2023)

$P(\cos\theta)=[(\cos\theta+1)/2]^{n_{\theta}+1}.$ (14)

We set $n_{\theta}=8$, which implies that binaries have a $20(55)\%$
probability to have $\theta_{1,2}$ that differ by less than $5(20)\%$. Note
that $n_{\theta}=0$ implies the isotropic distribution whilst $n_{\theta}\gg
1$ implies fully aligned spins, i.e. $\theta_{1}=\theta_{2}$.

For each BBH merger in our sample we select 1,000 values of the spin and spin
directions depending on the aforementioned assumptions, in order to assess
statistically the properties of Dragon-II mergers. The top panels in Figure 10
show the median value and 95th percentile of the effective spin parameter and
remnant BH mass for all BBH mergers in Dragon-II models. As, expected, we can
clearly see a difference between primordial binaries, which have mildly
aligned spins and thus $\chi_{\rm eff}>0$, and dynamical binaries, for which
$\chi_{\rm eff}\sim 0$. The plots suggest that the STEV model, based on
stellar evolution models, leads primordial binaries to have a $\chi_{\rm eff}$
smaller, on average, compared to the GAUSS model. The bottom panels of Figure
10 overlay to a single realisation of the simulated data the observed mergers
from GWTC-3, for comparison’s sake.

Noteworthy, the assumption that BHs form with a negligible spin unless matter
accretion processes are at play (STEV model) leads to a sub-population of
mergers with $\chi_{\rm eff}\sim 0$ and $m_{\rm
bin}=(40-100){\rm~{}M}_{\odot}$, a feature that disappears when a global
Gaussian spin distribution is adopted (GAUS model) as shown in Figure 10.

If BH spins do not strongly depend on stellar evolution processes, but rather
are well described by a general distribution, like a Gaussian, we can identify
two populations in the plot, one with clearly positive $\chi_{\rm eff}$ values
and $m_{\rm BH}<40{\rm~{}M}_{\odot}$, and one widely distributed around zero
$\chi_{\rm eff}$ involving massive BHs, $m_{\rm BH}>40{\rm~{}M}_{\odot}$.

Figure 10: Top panels: median and 95th percentile value of remnant mass and
effective spin parameters of Dragon-II mergers assuming a prescription for BH
natal spin based on stellar evolution models (STEV, left panel) or a Gaussian
distribution (GAUSS, right panel). Bottom panel: effective spin parameter and
remnant mass for one realisation of Dragon-II mergers (red points) and the
observed population of LVK mergers in the GWTC-3 catalog.

In order to improve the poor statistics, we proceed as follows: from the list
of Dragon-II mergers we create an oversampled catalogue by repeating the spin
assignment 100 times and, at each time, selecting a new "mock" BH mass in the
range $2.5-40.5{\rm~{}M}_{\odot}$ if the Dragon-II BH merger mass is below the
upper-mass gap, and in the range $40.5-100{\rm~{}M}_{\odot}$ otherwise. This
way, each real Dragon-II merger will have 100 counterparts with BHs of the
same class (upper mass-gap or not, merger in primordial or dynamical binary),
but enabling to build-up a catalogue sufficiently rich to analyse the overall
$\chi_{\rm eff}$ distribution. Figure 11 shows the distribution of $\chi_{\rm
eff}$ for the augmented sample in STEV and GAUS models. We see that the STEV
model follows a narrower distribution compared to the GAUS model, and exhibits
a clear peak around zero owing to the population of BHs formed from single
stars (see also Figure 9 in Arca Sedda et al., 2023).

Figure 11: Effective spin distribution for a population of mergers based on
Dragon-II models assuming model STEV (filled red steps) or GAUS (dashed line).

###  4.2 Compact binary merger rates

####  4.2.1 Merger efficiency

As described in the previous section, we have simulated a total Dragon-II mass
of $M_{\rm sim}=3.65\times 10^{6}{\rm~{}M}_{\odot}$ and find in total 78
mergers when GW recoil is not accounted for, and 74 otherwise. Therefore, the
resulting BH merger efficiency, defined as the ratio between the number of
mergers and the total simulated mass (Ziosi et al., 2014), is given by

$\eta_{\rm GW}=\frac{N_{\rm GW}}{M_{\rm sim}}\simeq(2.0-2.1)\times
10^{-5}{\rm~{}M}_{\odot}^{-1},$ (15)

similar to what inferred for young and open clusters with a similar
metallicity (e.g. Di Carlo et al., 2020b; Santoliquido et al., 2020; Rastello
et al., 2021). Note that given the limited simulation time our estimate could
represent a lower limit to the total merger efficiency in massive young and
intermediate-age clusters. Nonetheless, we note that as the cluster loses mass
and expands, the binary formation rate and binary–single interaction rate will
sharply decrease until the point in which it will be unlikely for tight
binaries to form and merge within a Hubble time.

Interestingly, at fixed value of the half-mass radius, the merger efficiency
changes sensibly with the initial binary fraction, being

$\eta_{\rm GW,fb}=\begin{cases}2.3\times
10^{-5}{\rm~{}M}_{\odot}^{-1}&f_{b}=0.20,\\\ 1.2\times
10^{-5}{\rm~{}M}_{\odot}^{-1}&f_{b}=0.05.\\\ \end{cases}$ (16)

This highlights the role of primordial binaries in determining the formation
of merging compact objects. For comparison, note that the merger efficiency
derived in Rastello et al. (2021) is based on star cluster models containing
$\sim 40\%$ of stars in primordial binaries.

To further explore the impact of cluster properties on the merger efficiency,
we show in Figure 12 the average merger efficiency per cluster, $\epsilon_{\rm
GW}(R_{\rm HM})$, as a function of the average cluster density
$\langle\rho_{\rm sim}\rangle$, using the following definitions

$\displaystyle\epsilon_{\rm GW}(R_{\rm HM})$ $\displaystyle=$
$\displaystyle\frac{N_{\rm GW}}{(M_{\rm sim}/N_{\rm sim})},$ (17)
$\displaystyle\langle\rho_{\rm sim}\rangle$ $\displaystyle=$
$\displaystyle\frac{M_{\rm sim}}{N_{\rm sim}R_{\rm HM}^{3}}$ (18)

where $M_{\rm sim}$ is the total simulated mass and $N_{\rm sim}$ is the
number of simulations performed for a given value of the half-mass radius,
$R_{\rm HM}$. At fixed value of the binary fraction, this relation is well
described by a power-law in the form $\epsilon_{\rm GW}=a(\langle\rho_{\rm
HM}\rangle/1{\rm M}_{\odot}{\rm pc}^{-3})^{b}$, with $a=(0.15\pm 0.07)\times
10^{-5}$ and $b=0.25\pm 0.03$.

The plot makes clear that increasing the cluster density by two orders of
magnitude leads to $\sim 2.5\times$ more mergers. Moreover, it further
highlights the role of primordial binaries, showing that clusters with a lower
binary fraction have a probability $\sim 50\%$ smaller to develop a merger, at
least in the case of $R_{\rm HM}=1.75$ pc.

Figure 12: Merger efficiency as a function of the average cluster density.
Different symbols correspond to different values of the cluster half-mass
radius. Open circles mark simulations with the lowest binary fraction.

####  4.2.2 Merger rate for black hole binaries

We define the cosmic merger rate following Santoliquido et al. (2020);
Bouffanais et al. (2021)

$\displaystyle\mathcal{R}(z)=$ $\displaystyle\frac{{\rm d}}{{\rm d}t_{\rm
lb}(z)}\int_{0}^{z_{\rm max}}\psi_{\rm clus}(z^{\prime})\frac{{\rm d}t_{\rm
lb}(z)}{{\rm d}z^{\prime}}{\rm d}z^{\prime}$ $\displaystyle\times\int_{Z_{\rm
min}}^{Z_{\rm max}}\eta_{\rm GW}(Z)\mathcal{F}(z^{\prime},z,Z){\rm d}Z,$ (19)

where $t_{\rm lb}(z)$ is the lookback time at merger, $\psi_{\rm
clus}(z^{\prime})$ is the star cluster formation rate when the merging binary
formed, $\eta_{\rm GW}(Z)$ is the merger efficiency at the metallicity $Z$,
$\mathcal{F}(z^{\prime},z,Z)$ is the number of mergers forming at redshift
$z^{\prime}$ and merging at redshift $z$ in environments with metallicity $Z$.
The adoption of Equation 19 enables us to compare Dragon-II simulation results
with those obtained for low-mass star clusters (Santoliquido et al., 2020).
Note that this procedure does not take into account possible effects related
to the initial cluster mass function, which could indeed have an impact on the
overall merger rate (see e.g. Antonini & Gieles, 2020a). Nonetheless, the
similarity between the merger efficiency derived from Dragon-II simulations
and that obtained by Santoliquido et al. (2020) for low-mass clusters suggests
that it is possible to safely utilise the merger efficiency as a proxy of the
overall number of mergers per unit mass in the whole range of possible cluster
masses. This choice, although representing an approximation, permits us to
avoid the inclusion of a cluster mass function in Equation 19 and all the
related uncertainties, like the cluster mass function boundaries and
functional form.

We adopt a cosmic star cluster formation rate in the form

$\psi_{\rm clus}(z)=\frac{0.01(1+z)^{2.6}f_{\rm
CFE}}{1+\left[(1+z)/{3.2}\right]^{6.2}}~{}{\rm~{}M}_{\odot}{\rm yr}^{-1}{\rm
Mpc}^{-3},$ (20)

i.e. we rescale the stellar star formation rate derived by (Madau & Fragos,
2017) by a factor $f_{\rm CFE}$, which represents the cluster formation
efficiency, i.e. the fraction of star formation that goes into bound clusters.
Although uncertain, observations and models suggest that the cluster formation
efficiency (CFE) can be as large as $f_{\rm CFE,YC}=0.3$ for young clusters
(Mapelli et al., 2021) and $f_{\rm CFE,GC}=0.08\pm 0.03$ (Bastian, 2008) for
globular clusters, regardless of the star formation history. In the following,
we adopt both young and globular cluster CFE values to constrain the BBH
merger rate in our simulations.

For dynamical mergers, it has been shown that the merger efficiency $\eta_{\rm
GW}(Z)$ remains almost constant in the range $Z<10^{-3}$, and decreases
roughly by an order of magnitude at solar values (Di Carlo et al., 2020a, b;
Rastello et al., 2021). Since our models have all the same metallicity,
$Z=0.005$, to infer the merger rate we assume that the merger efficiency is
constant at $Z<0.005$ and reduces by 10 times at larger metallicities (see
e.g. Figure 1 in Santoliquido et al., 2020). Moreover, we factorise the
function $F(z,z^{\prime},Z)=p(Z,z^{\prime})N(z,z^{\prime})$, thus assuming
that the number of mergers at redshift $z$ that formed at $z^{\prime}$ is
independent on the metallicity distribution. The $p(Z,z^{\prime})$ term
represents the cosmic fraction of clusters with metallicity in the
($Z$,$Z+dZ$) bin at redshift $z^{\prime}$. We assume that the metallicity
follows a log-normal distribution peaked at (Madau & Fragos, 2017)

${\rm
Log}\left\langle\frac{Z(z)}{{\rm~{}Z}_{\odot}}\right\rangle=0.153-0.074z^{1.34},$
(21)

with dispersion either $\sigma_{Z}=0.2-0.5-0.8$ (Schrøder et al., 2020;
Bouffanais et al., 2021; Bavera et al., 2021; Zevin et al., 2021).

Since all Dragon-II models have the same metallicity, to infer the simulated
merger rate we integrate Equation 19 under two assumptions, one conservative
and one optimistic. In the conservative case, we consider only clusters with a
metallicity $Z<0.005$ and assume that they have a similar merger rate
efficiency (Santoliquido et al., 2020). In the optimistic case, instead, we
include in the integration also clusters with metallicity larger than the
simulated one, reducing for metal-rich clusters the simulated merger
efficiency by a factor 10, as expected from low-$N$ simulations of young
clusters (Di Carlo et al., 2019).

To compare with similar estimates in the literature, we first set $f_{\rm
CFE}=1$, i.e. that all stars form in star clusters, and calculate a merger
rate of $\mathcal{R}=27{\rm~{}yr}^{-1}{\rm~{}Gpc}^{-3}$, in broad agreement
with the rate inferred for low-mass star clusters ($N=10^{2}-5\times
10^{4}{\rm~{}M}_{\odot}$) (Di Carlo et al., 2020b; Rastello et al., 2021;
Santoliquido et al., 2020) and semi-analytic models of young and globular
clusters (see e.g. Mapelli et al., 2021).

A more reliable estimate of the merger rate is shown in Figure 13 for both the
conservative and optimistic cases, and assuming different values of the
cluster formation efficiency, $f_{\rm CFE}=0.08-0.3$. As shown in the plot, we
find a simulated merger rate of $\mathcal{R}_{\rm GW}=(12\pm 7)$ at redshift
$z=0.2$. At the same redshift, the BBH merger rate inferred by the LVK is
$\mathcal{R}_{\rm LVK}=17.9-44{\rm~{}yr}^{-1}{\rm~{}Gpc}^{-3}$ (The LIGO
Scientific Collaboration et al., 2021a, b).

Figure 13: Merger rate density for all simulations in our runs, assuming either metal-poor only (red straight line) or both metal-poor and rich (black dashed line) clusters. The shaded area embrace two limiting values of the cluster formation efficiency, with the upper limit corresponding to $f_{\rm CFE}=0.3$ and the lower one corresponding to $f_{\rm CFE}=0.08$. Source type | $\mathcal{R}_{\rm loc}$
---|---
| ${\rm~{}yr}^{-1}{\rm~{}Gpc}^{-3}$
BBH | $5-19$
NS-BH | $0.027-8.7$
WD-BH | $3.8\times 10^{-4}-2.3$
Table 2: Volumetric merger rate calculated in the local Universe for different
GW sources in Dragon-II clusters.

####  4.2.3 Merger rate for exotic mergers

In Dragon-II simulations we find 3 elusive mergers: one WD-BH and two NS-BH
mergers. Despite they are evidently too scarce to allow a statistical
treatment, we can exploit them to attempt a rough, order of magnitude,
estimate of the merger rates for these two classes of GW sources assembled in
star clusters as:

$R_{\rm xBH}(<D)=\frac{N_{\rm x}}{M_{\rm sim}}f_{x}\delta M_{g*}N(<D)t_{\rm
rel}^{-1},$ (22)

where $M_{g*}$ is the galaxy stellar mass, $\delta=0.001-0.01$ is the fraction
of galaxy mass made up by star clusters (Spitler & Forbes, 2009; Georgiev et
al., 2010; Harris et al., 2013; Arca-Sedda & Capuzzo-Dolcetta, 2014; Gnedin et
al., 2014; Webb & Leigh, 2015), $f_{x}$ is the fraction of clusters with a
given property (e.g. age within a certain range), $t_{\rm rel}$ is the cluster
relaxation time, and $N(<D)$ is the number of MW equivalent galaxies within a
given cosmological distance $D$ (Abadie et al., 2010)

$N(<D)=4\pi/3(2.26)^{-3}\left(\frac{D}{{\rm
Mpc}}\right)^{3}\frac{\rho_{g}}{{\rm Mpc^{-3}}},$ (23)

where $\rho_{g}=0.0116$ Mpc-3 is the galaxy number density in the local
Universe (Kopparapu et al., 2008). Moreover, we consider typical relaxation
times of either globular clusters, $t_{\rm rel}=10^{9}$ yr (Harris, 2010), or
massive and relatively young clusters in the Small Magellanic Cloud (SMC),
$t_{\rm rel}=3.2\times 10^{7}$ yr (Gatto et al., 2021).

Note that the relaxation time of Galactic clusters is inferred from their
present time properties. Depending on the amount of mass lost and the level of
cluster expansion, it could be possible that the initial relaxation time was
relatively shorter and therefore the number of dynamically old globular
clusters is larger than what we see at present. In these regards, note that
the relaxation time of SMC clusters, which are generally younger than a few
Gyr, is sensibly smaller compared to Milky Way globulars, possibly because
relaxation processes did not have time to sufficiently influence the cluster
dynamics.

In the following calculations, we consider Milky Way-like galaxies only,
$M_{g*}=6\times 10^{10}{\rm~{}M}_{\odot}$ (Licquia & Newman, 2015) located
within $D=1$ Gpc. In the Milky Way, there are only $\sim 4$ out of 155
globular clusters with an age larger than 1 relaxation time, whilst around
half of clusters in the SMC satisfy this requirement, thus $f_{x}\sim
0.025-0.5$.

This implies a frequency rate for WD-BH mergers in the local Universe of
$R_{\rm WDBH}=(1.8\times 10^{-3}-10.8)$ yr-1, corresponding to a volumetric
merger rate $\mathcal{R}_{\rm WDBH}=RV_{\rm com}^{-1}(1{\rm~{}Gpc})=(3.8\times
10^{-4}-2.3){\rm~{}yr}^{-1}{\rm~{}Gpc}^{-3}$.

In the case of NS-BH mergers, instead, the event occurs over a timescale of
$t_{\rm GW}=(0.04-0.5)t_{\rm rel}$. The fraction of cluster with an age longer
than $t_{\rm GW}$ is $f_{x}\sim 0.94$ for clusters in both the Milky Way and
the SMC, the resulting frequency rate for NS-BH mergers is $R_{\rm
NSBH}=(0.13-40.7)$ yr-1, which implies a volumetric merger rate of
$\mathcal{R}_{\rm NSBH}=(0.027-8.7){\rm~{}yr}^{-1}{\rm~{}Gpc}^{-3}$.

###  4.3 Multimessenger sources: prospects for LISA detection

Over the next decade, the network of ground-based detectors will be
complemented by LISA, possibly the first space-borne low-frequency detector.
LISA will be able to possibly catch the GW echoes of merging stellar BHs,
IMBHs, and nearby WD and NS binaries. While we postpone a detailed discussion
about BBHs forming in young massive clusters detectable with LISA to a
forthcoming paper, we focus in the following on the handful exotic mergers
that develop in our Dragon-II models.

Let us consider the case of a WD-BH merger. We have shown in Section 3.2 that
such a source could appear as an X-ray binary and give rise to a TDE once the
WD approaches too closely the BH. Assuming that the binary evolves solely via
GW emission, and adopting the Peters & Mathews (1963) formalism to evolve the
binary until the merger, we find that around 6 months prior to the merger, the
WD will overfill the Roche lobe and start the X-ray binary phase. At
disruption, the frequency of the associated GW emission is given by (Hils &
Bender, 1995; Kobayashi et al., 2004; Rosswog et al., 2009; Dai & Blandford,
2013; Fragione et al., 2020)

$\displaystyle f_{\rm GW}\simeq$ $\displaystyle 0.09{\rm
Hz}\left(1+\frac{M_{\rm WD}}{M_{\rm BH}}\right)\times$
$\displaystyle\times\left(\frac{M_{\rm
WD}}{0.6{\rm~{}M}_{\odot}}\right)^{1/2}\left(\frac{R_{\rm WD}}{10^{4}{\rm
km}}\right)^{-3/2}=0.13{\rm Hz},$ (24)

where we have assumed $R_{\rm WD}=10^{4}$ km. Note that an eccentricity
between 0 and 1 would affect $f_{\rm GW}$ by less than $20\%$ (Fragione et
al., 2020). The amplitude of the emitted signal at disruption will be (Robson
et al., 2019; Fragione et al., 2020)

$\displaystyle h_{c}\simeq$ $\displaystyle 2\times 10^{-20}\left(\frac{T_{\rm
obs}}{4{\rm yr}}\right)^{1/2}\left(\frac{D_{L}}{10{\rm
Mpc}}\right)^{-1}\left(\frac{M_{\rm
BH}}{10{\rm~{}M}_{\odot}}\right)^{0.66}\times$
$\displaystyle\times\left(\frac{M_{\rm
WD}}{0.6{\rm~{}M}_{\odot}}\right)^{1.58}\left(\frac{R_{\rm WD}}{10^{4}{\rm
km}}\right)^{-1.75}\simeq 10^{-19}.$ (25)

Since the WD will disrupt completely as crossing its Roche limit, the
associated GW emission will appear as a burst (Rosswog et al., 2009). For such
source, the corresponding signal-to-noise ratio (S/N) for LISA can be written
as (see e.g. Robson et al., 2019)

$\left({\rm S}/{\rm
N}\right)=f^{2/3}\frac{h_{c}}{S_{c}}=1.2\left(\frac{D_{L}}{10{\rm
Mpc}}\right)^{-1},$ (26)

where $S_{c}$ is the detector sensitivity curve in terms of characteristic
strain (Robson et al., 2019) and we have exploited the intrinsic dependence on
the measurable GW strain and the source luminosity distance $D$.

If the merger occurs inside the Milky Way, i.e. at $D<0.1$ Mpc, it would
appear as a loud source in LISA, with (S/N)$>120$. More in general, the
maximum distance at which LISA could detect such merger with a minimum signal-
to-noise ratio of (S/N)$>8(15)$ is $D<1.5{\rm~{}Mpc}$($0.7{\rm~{}Mpc}$).

Note that the Andromeda galaxy is $\sim 0.7-0.8$ Mpc away from us, therefore
to roughly estimate the probability for a closeby WD-BH merger we can replace
in Equation 22 $N(<D)=2$ and find an upper limit to the local merger rate of
closeby WD-BH mergers of $R_{\rm WDBH,close}<(8.4\times 10^{-10}-5.1\times
10^{-6})$ yr-1.

###  4.4 The pair-instability supernova rate for massive star clusters:
perspectives for detection via magnitude limited surveys

The onset of IMBH formation and the development of BBH mergers depend
intrinsically on the cluster radius and initial density, the amount of stars
initially in a binary, and the stellar evolution recipes adopted – e.g. BH
matter accretion efficiency, the physics of PISNe and PPISNe.

In these regards, the fact that PISNe are rare events for which a smoking gun
has not been observed yet (e.g. Woosley et al., 2007; Gal-Yam et al., 2009;
Terreran et al., 2017; Gomez et al., 2019; Woosley & Smith, 2022), offers us
the possibility to use this physical process as a diagnostic quantity in
Dragon-II models. In practice, we can infer the PISN rate in Dragon-II
simulations and compare such rate with current observation limits to explore
whether our simulations produce unrealistically large PISN frequency rates. As
described in paper AS-II, in Dragon-II models PISNe develop either in single
stars or in stellar merger products, provided that their core Helium reaches a
mass in the range $(64-130){\rm~{}M}_{\odot}$. This offer us a unique
possibility to explore the impact of PISNe in star clusters, taking
simultaneously into account the impact of stellar mergers in the overall
population of PISN progenitors. According to the adopted stellar evolution, in
a simple stellar population only stars heavier than $m_{\rm ZAMS}\geq m_{\rm
PISN}=150{\rm~{}M}_{\odot}$ could undergo a PISN event, i.e. larger than the
maximum stellar mass adopted for the initial mass function. Instead, in
Dragon-II models we find 23 stars that undergo supernova. All these stars are
either in a primordial binary or are captured in a binary before the explosion
and undergo one or more stellar merger and accretion events that bring the
star mass above $m_{\rm PISN}$. Typical masses for Dragon-II PISN progenitors
are in the range $(150-282){\rm~{}M}_{\odot}$.

The simulated PISN efficiency can be defined similarly to the compact object
merger rate, i.e. $\eta_{\rm PISN}=N_{\rm PISN}/M_{\rm sim}=6.2\times
10^{-6}{\rm~{}M}_{\odot}^{-1}$.

To calculate the PISN rate, we follow the approach adopted by du Buisson et
al. (2020). Firstly, we assume that the Ni mass of the massive star that goes
off as a PISN can be calculated via the following equation:

${\rm Log}\left(M_{\rm Ni}/{\rm~{}M}_{\odot}\right)=r\left(M_{\rm
He,f}/{\rm~{}M}_{\odot}\right)^{s}+t,$ (27)

where $r=-5.02\times 10^{4}$, $s=-2.159$, and $t=2.887$ (Heger & Woosley,
2002; du Buisson et al., 2020), and $M_{\rm He,f}$ is the final mass of the
star He core. The Ni mass is used to infer the peak bolometric magnitude
exploiting an Arnett-like relation (Arnett, 1982)

$\Upsilon_{\rm bol,Ni}^{{\rm GW}}=-19.2-2.5{\rm Log}\left(M_{\rm
Ni}/0.6{\rm~{}M}_{\odot}\right),$ (28)

which can be converted into an apparent bolometric magnitude via the Pogson’s
relation

$\mu_{\rm bol}^{{\rm GW}}=\Upsilon_{\rm bol,Ni}^{{\rm GW}}+5{\rm
Log}(D_{L}/10{\rm pc}),$ (29)

being $D_{L}$ the luminosity distance. To simplify the calculations, we adopt
for the He mass, which is the main ingredient to calculate the Ni mass, an
average value of $M_{\rm He,f}=90.4{\rm~{}M}_{\odot}$ as extracted from our
models.

The value of $\mu_{\rm bol}^{{\rm GW}}$ is used to calculate whether a PISN
can be detected in a magnitude limited survey. Assuming to have a population
of PISNe with apparent magnitude distributed according to a Gaussian around
$\mu_{\rm bol}^{{\rm GW}}$ and a magnitude detection threshold $\mu_{\rm
lim}$, we define the fraction of detectable sources as

$f_{\rm GSS}=0.5\left[1+{\rm erf}\left(\frac{\mu_{\rm lim}-\mu_{\rm bol}^{{\rm
GW}}}{\sqrt{2}\sigma_{\mu}}\right)\right],$ (30)

where we adopted $\sigma_{\mu}=0.2$333We verified that varying $\sigma_{\mu}$
in the range $0.1-0.3$ has little effect on our results..

The PISN rate as a function of the redshift can thus be evaluated as:

$\mathcal{R}_{\rm PISN}(z)=\int_{z_{1}}^{z_{2}}\frac{dV}{dz}\psi(z)\eta_{\rm
PISN}f_{\rm GSS}(z)f_{\rm Z}(z)dz,$ (31)

where $dV/dz$ is the comoving volume element and $\psi(z)$ is the cosmic star
formation rate, for which we assume the cosmic star formation history in
Equation 20 (Madau & Fragos, 2017) and the same limits for $f_{\rm CFE}$
described in Section 4.2, and that only stars with a metallicity $Z\leq 0.008$
undergo PISNe (Spera & Mapelli, 2017). Figure 14 shows the PISNe rate for the
intrinsic cosmic population and assuming different detection threshold in
magnitude limited surveys, namely $\mu_{\rm bol}=17,~{}20,~{}25$. Note that
these threshold roughly corresponds to the typical maximum detectable
magnitude of already completed, like the Sloan Digital Sky Survey (SDSS444SDSS
home: http://www.sdss.org) or the Palomar Transient Factory (PTF555PTF home:
http://www.ptf.caltech.edu), ongoing, e.g. the Dark Energy Survey (DES666DES
home: http://www.darkenergysurvey.org), and future surveys, like the Large
Synoptic Survey Telescope (LSST777LSST home: http://www.lsst.org), the Zwicky
Transient Facility (ZTF888ZTF home: https://www.ztf.caltech.edu), or the
EUCLID mission999EUCLID home: https://sci.esa.int/web/euclid (Moriya et al.,
2022). From Figure 14 we see that only future surveys ($\mu_{\rm bol}\geq 25$)
will be able to probe the cosmological properties of PISNe, whilst current
surveys could in principle place constraints on PISNe within a redshift
$z<0.3$.

Integrating Equation 31 over the redshift returns the number of detected
sources per year. The possible number of PISN detections per year for
different values of the limiting bolometric magnitude, $\mu_{\rm bol}$, and
the cluster formation efficiency, $f_{\rm CFE}$, is summarized in Table 3.
From the table is clear that the detection of PISNe from star clusters is
still highly unlikely in completed and ongoing surveys, but it could lead to
$\sim 8$ detections per year with the next generation of detectors. Comparing
future PISNe detections with numerical models could have a twofold aim. On the
one hand, it will permit us to shed light on the actual contribution of
massive stars in dense clusters to the overall population of PISNe. On the
other hand, it will provide us with an useful term of comparison to determine
the reliability of cluster simulations.

Figure 14: PISNe cosmological rate (blue straight line) and detection rate assuming a magnitude limited survey at different threshold values, namely $m_{\rm bol}=17$ (red dotted line), $20$ (green dashed-dotted line), and $25$ (orange dashed line), assuming a cluster formation efficiency of $f_{\rm CFE}=1$ and assuming that only clusters with a metallicity $Z\lesssim 0.008$ can host stars undergoing PISN. Table 3: Cosmological and detectable PISNe occurrence rate. Col. 1: cluster formation efficiency. Col. 2-5: number of PISNe per year for the cosmological population, and for magnitude limited surveys with different treshold magnitudes. $f_{\rm CFE}$ | No. of detections per 1 year
---|---
| Cosmological | $m_{\rm bol}=25$ | $m_{\rm bol}=20$ | $m_{\rm bol}=17$
$0.08$ | $24$ | $0.7$ | $0.001$ | $1.8\times 10^{-5}$
$0.3$ | $89$ | $2.6$ | $0.004$ | $6.8\times 10^{-5}$
$1$ | $297$ | $8.8$ | $0.012$ | $2.3\times 10^{-4}$

##  5 Conclusions

In this paper we have presented and discussed the properties of compact binary
mergers and PISNe in the Dragon-II simulations, a suite of direct $N$-body
models representing star clusters with up to 1 million stars and a relatively
large ($10\%-33\%$) binary fraction. Our main results can be summarised as
follows:

  * •

We find a population of 75 BBH, 2 NS-BH, and 1 WD-BH mergers. Among them, 4
BBHs avoid merger when GW recoils are enabled. Mergers occurring inside the
cluster make-up the $\gtrsim 40\%$ of the whole population and are mostly due
to mergers formed via dynamical interactions (dynamical mergers). The
population of ejected mergers, which merge outside the parent cluster, are
equally contributed by mergers formed dynamically and from primordial binaries
(primordial mergers). Typically, in-cluster mergers have primaries with masses
$m_{\rm BH,1}>30{\rm~{}M}_{\odot}$ and companion in the $m_{\rm
BH,2}=30-50{\rm~{}M}_{\odot}$ mass range, whilst ejected mergers involve
lighter primaries, $m_{\rm BH,1}<40{\rm~{}M}_{\odot}$, and are characterised
by fairly large mass ratios, $q>0.6$;

  * •

Mergers forming from primordial binaries are characterised by large mass
ratios and component masses clearly smaller than those formed dynamically.
Among dynamical mergers, the most massive ones are those in which at least one
component had an ancestor in a primordial binary;

  * •

BBH mergers are characterised by a delay time that nicely distribute around a
value of $10-30$ cluster relaxation time. This highlights the fact that the
processes that trigger BBH formation and merger are intrinsically related to
the cluster dynamical evolution;

  * •

The population of mergers forming from dynamical interactions or primordial
binaries is clearly distinguishable from the residual eccentricity of the
binary as it enters in the typical frequency bands of GW detectors, i.e.
$f=0.001-100$ Hz. We find that practically all primordial binaries are
circular at merger, this implying that primordial binaries merge before
dynamics can have an impact on their evolution, whilst around $20-40-5\%$ of
mergers preserve an eccentricity $e>0.1$ when entering the LISA-DECIGO-LIGO
bands. All mergers with $e>0.1$ in the 0.05-1 Hz and 1-10 Hz bands occur
inside the cluster, whilst half of eccentric mergers in the mHz band are
ejected. This hints at the possibility to distinguish the formation history of
a BBH merger from the frequency band in which it is observed;

  * •

We identify three exotic mergers in our sample: a WD-BH binary formed
dynamically and two NS-BH mergers, one formed dynamically and the other from a
primoardial binary. A WD-BH merger that forms after 4 cluster relaxation time
and it is triggered by chaotic interactions that increase the eccentricity up
to an extremal value of $e=0.99994930$. Once the WD approaches sufficiently
close the BH, this type of sources could appear as an ultraluminous X-ray
sources and, ultimately, be a source detectable by LISA if it occurs within
700 kpc from us, i.e. within the distance between the Milky Way and Andromeda.
The dynamical NS-BH binary is characterised by a chirp mass
$\mathcal{M}=3.4{\rm~{}M}_{\odot}$, larger than what predicted by the isolated
stellar evolution scenario, and preserve an eccentricity of $e=0.9974(0.21)$
when crossing a frequency of $f=0.5(1)$ Hz, thus future observations with ET
could help probing the population of closeby, dynamically formed, NS-BH
mergers. The primordial NS-BH binary is not affected by dynamics at all, thus
they can be mistaken for a merger occurring in isolation. This highlights the
importance of star clusters with a large binary fraction as contributors of
the isolated scenario of compact binary mergers. None of the NS-BH mergers are
expected to release EM emission, unless the BHs have a spin $\chi>0.9$;

  * •

We find that comparing the remnant mass and spin of BBH mergers could help
untangling their origin. Using a model based on stellar evolution theories, we
show that primordial binary mergers are characterised by remnant masses
systematically smaller and effective spin parameters systematically larger
than dynamical mergers;

  * •

We derive a BBH merger efficiency of $\sim 2\times
10^{-5}{\rm~{}M}_{\odot}^{-1}$, comparable with the value estimated for low-
mass star clusters. Interestingly, we find that the merger efficiency depends
on the star cluster properties. Decreasing the binary fraction by a factor
$4$, for example, leads to a decrease of the merger efficiency by a factor
$\sim 2$. Moreover, the merger efficiency increases with the cluster density
following a power-law with slope $\sim 0.25$. We adopt a series of
cosmologically motivated assumptions for the cosmic star formation history,
and use them to infer a merger rate density at redshift $z<0.2$ of
$\mathcal{R}=5-19~{}(0.027-8.7)~{}(3.8\times
10^{-4}-2.3){\rm~{}yr}^{-1}{\rm~{}Gpc}^{-3}$ for BBHs(WD-BH)(NS-BH) mergers,
respectively. We predict that, in a 4 yr-long mission, LISA could detect
$N_{\rm BBH}=12\pm 7$($5\pm 3$) BBH mergers (IMRIs) and can identify the WD-BH
merger with a signal-to-noise ratio SNR$>8$($15$) if it occurs within
$D_{L}<1.5$($0.7$) Mpc from us.

  * •

We retrieve the cosmic frequency rate of PISNe, in order to explore the
reliability of our simulations on the one hand, and to make predictions for
PISNe detection from star clusters on the other hand. We find that future
surveys with a limiting magnitude of $m_{\rm bol}=25$ could detect $N_{\rm
PISN}=0.7-8.8$ PISNe per year. Comparing these estimates with future surveys
could help placing constraints on the population of massive stars in dense
star clusters.

The Dragon-II clusters represent a further step forward in the modelling of
young and intermediate-age star clusters, providing the first suite of
simulations that models clusters with both $N>120,000$ stars (up to $10^{6}$),
a high binary fraction (up to $33\%$), and an initial density of
$\rho=(1.2\times 10^{4}-1.6\times 10^{6}){\rm~{}M}_{\odot}$ pc-3. These
simulations complement the vast literature of $N$-body simulations of lower-
mass and lower density star clusters (e.g. Di Carlo et al., 2019; Rizzuto et
al., 2021; Rastello et al., 2021; Banerjee, 2021; Rizzuto et al., 2022;
Banerjee, 2022), and provide the largest catalogue of BH mergers obtained in
direct $N$-body simulations of metal-poor, dense and massive young clusters.

## Acknowledgements

The authors thank the referee for their constructive and helful report. The
authors warmly thank Agostino Leveque for their help and assistance in using
their implementation of the McLuster code, and Giuliano Iorio, Sara Rastello,
and Michela Mapelli for useful comments and discussion.

This work benefited of the support from the Volkswagen Foundation Trilateral
Partnership through project No. 97778 “Dynamical Mechanisms of Accretion in
Galactic Nuclei” and the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) – Project-ID 138713538 – SFB 881 “The Milky Way System”), and by
the COST Action CA16104 “GWverse”. The authors gratefully acknowledge the
Gauss Centre for Supercomputing e.V. for funding this project by providing
computing time through the John von Neumann Institute for Computing (NIC) on
the GCS Supercomputer JUWELS Booster at Jülich Supercomputing Centre (JSC).
Data analysis and part of the runs were conducted on the GRACE-BH HPC
workstation, funded by the European Union’s under the research project GRACE-
BH.

MAS acknowledges funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No.
101025436 (project GRACE-BH, PI: Manuel Arca Sedda).

AWHK is a fellow of the International Max Planck Research School for Astronomy
and Cosmic Physics at the University of Heidelberg (IMPRS-HD). The work of PB
was supported by the Volkswagen Foundation under the special stipend No.
9B870.

PB acknowledge the support within the grant No. AP14869395 of the Science
Committee of the Ministry of Science and Higher Education of Kazakhstan
("Triune model of Galactic center dynamical evolution on cosmological time
scale"). The work of PB was supported under the special program of the NRF of
Ukraine Leading and Young Scientists Research Support - "Astrophysical
Relativistic Galactic Objects (ARGO): life cycle of active nucleus", No.
2020.02/0346.

RS thanks Max Planck Institute for Astrophysics (Thorsten Naab) for
hospitality during many visits

MG was partially supported by the Polish National Science Center (NCN) through
the grant No. 2021/41/B/ST9/01191.

FPR acknowledges the support by the European Research Council via ERC
Consolidator Grant KETJU (no. 818930).

TN acknowledges the support of the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 -
390783311 of the DFG Cluster of Excellence "ORIGINS”.

## Data Availability

The data from the runs of these simulations and their initial models will be
made available upon reasonable request by the corresponding author. The
Nbody6++GPU code is publicly
available101010https://github.com/nbody6ppgpu/Nbody6PPGPU-beijing. The
McLuster version used in this work will soon be available. A similar version
is described in Leveque et al. (2022).

## References

  * Aarseth (1999) Aarseth S. J., 1999, PASP, 111, 1333
  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations
  * Aarseth et al. (1974) Aarseth S. J., Henon M., Wielen R., 1974, A&A, 37, 183
  * Aarseth et al. (2008) Aarseth S. J., Tout C. A., Mardling R. A., 2008, The Cambridge N-Body Lectures. Vol. 760, doi:10.1007/978-1-4020-8431-7,
  * Abadie et al. (2010) Abadie J., et al., 2010, Classical and Quantum Gravity, 27, 173001
  * Acernese et al. (2015) Acernese F., et al., 2015, Classical and Quantum Gravity, 32, 024001
  * Ahmad & Cohen (1973) Ahmad A., Cohen L., 1973, Journal of Computational Physics, 12, 389
  * Ajith et al. (2011) Ajith P., et al., 2011, Phys. Rev. Lett., 106, 241101
  * Amaro-Seoane et al. (2013) Amaro-Seoane P., et al., 2013, GW Notes, 6, 4
  * Amaro-Seoane et al. (2022a) Amaro-Seoane P., et al., 2022a, arXiv e-prints, p. arXiv:2203.06016
  * Amaro Seoane et al. (2022b) Amaro Seoane P., et al., 2022b, General Relativity and Gravitation, 54, 3
  * Antonini & Gieles (2020a) Antonini F., Gieles M., 2020a, Phys. Rev. D, 102, 123016
  * Antonini & Gieles (2020b) Antonini F., Gieles M., 2020b, MNRAS, 492, 2936
  * Antonini & Perets (2012) Antonini F., Perets H. B., 2012, ApJ, 757, 27
  * Antonini & Rasio (2016) Antonini F., Rasio F. A., 2016, ApJ, 831, 187
  * Antonini et al. (2018) Antonini F., Rodriguez C. L., Petrovich C., Fischer C. L., 2018, MNRAS, 480, L58
  * Antonini et al. (2019) Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008
  * Antonini et al. (2022) Antonini F., Gieles M., Dosopoulou F., Chattopadhyay D., 2022, arXiv e-prints, p. arXiv:2208.01081
  * Arca Sedda (2020a) Arca Sedda M., 2020a, Communications Physics, 3, 43
  * Arca Sedda (2020b) Arca Sedda M., 2020b, ApJ, 891, 47
  * Arca Sedda (2021) Arca Sedda M., 2021, ApJ, 908, L38
  * Arca Sedda & Benacquista (2019) Arca Sedda M., Benacquista M., 2019, MNRAS, 482, 2991
  * Arca-Sedda & Capuzzo-Dolcetta (2014) Arca-Sedda M., Capuzzo-Dolcetta R., 2014, MNRAS, 444, 3738
  * Arca Sedda et al. (2018) Arca Sedda M., Askar A., Giersz M., 2018, MNRAS, 479, 4652
  * Arca Sedda et al. (2020a) Arca Sedda M., et al., 2020a, Classical and Quantum Gravity, 37, 215011
  * Arca Sedda et al. (2020b) Arca Sedda M., Mapelli M., Spera M., Benacquista M., Giacobbo N., 2020b, ApJ, 894, 133
  * Arca Sedda et al. (2021a) Arca Sedda M., Li G., Kocsis B., 2021a, A&A, 650, A189
  * Arca Sedda et al. (2021b) Arca Sedda M., Amaro Seoane P., Chen X., 2021b, A&A, 652, A54
  * Arca-Sedda et al. (2021c) Arca-Sedda M., Rizzuto F. P., Naab T., Ostriker J., Giersz M., Spurzem R., 2021c, ApJ, 920, 128
  * Arca Sedda et al. (2023) Arca Sedda M., Mapelli M., Benacquista M., Spera M., 2023, MNRAS, 520, 5259
  * Arnett (1982) Arnett W. D., 1982, ApJ, 253, 785
  * Askar et al. (2017) Askar A., Szkudlarek M., Gondek-Rosińska D., Giersz M., Bulik T., 2017, MNRAS, 464, L36
  * Askar et al. (2018) Askar A., Arca Sedda M., Giersz M., 2018, MNRAS, 478, 1844
  * Bahramian et al. (2017) Bahramian A., et al., 2017, MNRAS, 467, 2199
  * Ballone et al. (2022) Ballone A., Costa G., Mapelli M., MacLeod M., 2022, arXiv e-prints, p. arXiv:2204.03493
  * Banerjee (2018) Banerjee S., 2018, MNRAS, 473, 909
  * Banerjee (2021) Banerjee S., 2021, MNRAS, 500, 3002
  * Banerjee (2022) Banerjee S., 2022, A&A, 665, A20
  * Banerjee et al. (2010) Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 371
  * Banerjee et al. (2020) Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R., Spurzem R., Wang L., 2020, A&A, 639, A41
  * Banerjee et al. (2023) Banerjee S., Olejak A., Belczynski K., 2023, arXiv e-prints, p. arXiv:2302.10851
  * Bardeen et al. (1972) Bardeen J. M., Press W. H., Teukolsky S. A., 1972, ApJ, 178, 347
  * Bastian (2008) Bastian N., 2008, MNRAS, 390, 759
  * Bavera et al. (2020) Bavera S. S., et al., 2020, A&A, 635, A97
  * Bavera et al. (2021) Bavera S. S., et al., 2021, A&A, 647, A153
  * Belczynski et al. (2002) Belczynski K., Kalogera V., Bulik T., 2002, ApJ, 572, 407
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217
  * Belczynski et al. (2016) Belczynski K., et al., 2016, A&A, 594, A97
  * Belczynski et al. (2020) Belczynski K., et al., 2020, A&A, 636, A104
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Bouffanais et al. (2021) Bouffanais Y., Mapelli M., Santoliquido F., Giacobbo N., Di Carlo U. N., Rastello S., Artale M. C., Iorio G., 2021, MNRAS, 507, 5224
  * Branchesi et al. (2023) Branchesi M., et al., 2023, arXiv e-prints, p. arXiv:2303.15923
  * Broekgaarden et al. (2021) Broekgaarden F. S., et al., 2021, MNRAS, 508, 5028
  * Campanelli et al. (2007) Campanelli M., Lousto C. O., Zlochower Y., Merritt D., 2007, Phys. Rev. Lett., 98, 231102
  * Carr & Hawking (1974) Carr B. J., Hawking S. W., 1974, MNRAS, 168, 399
  * Carr et al. (2016) Carr B., Kühnel F., Sandstad M., 2016, Phys. Rev. D, 94, 083504
  * Chatterjee et al. (2017) Chatterjee S., Rodriguez C. L., Kalogera V., Rasio F. A., 2017, ApJ, 836, L26
  * Chattopadhyay et al. (2022) Chattopadhyay D., Hurley J., Stevenson S., Raidani A., 2022, MNRAS, 513, 4527
  * Church et al. (2017) Church R. P., Strader J., Davies M. B., Bobrick A., 2017, ApJ, 851, L4
  * Costa et al. (2021) Costa G., Bressan A., Mapelli M., Marigo P., Iorio G., Spera M., 2021, MNRAS, 501, 4514
  * Costa et al. (2022) Costa G., Ballone A., Mapelli M., Bressan A., 2022, MNRAS, 516, 1072
  * Dai & Blandford (2013) Dai L., Blandford R., 2013, MNRAS, 434, 2948
  * Di Carlo et al. (2019) Di Carlo U. N., Giacobbo N., Mapelli M., Pasquato M., Spera M., Wang L., Haardt F., 2019, MNRAS, 487, 2947
  * Di Carlo et al. (2020a) Di Carlo U. N., Mapelli M., Bouffanais Y., Giacobbo N., Santoliquido F., Bressan A., Spera M., Haardt F., 2020a, MNRAS, 497, 1043
  * Di Carlo et al. (2020b) Di Carlo U. N., et al., 2020b, MNRAS, 498, 495
  * Di Carlo et al. (2021) Di Carlo U. N., et al., 2021, MNRAS, 507, 5132
  * Dominik et al. (2012) Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2012, ApJ, 759, 52
  * Downing et al. (2010) Downing J. M. B., Benacquista M. J., Giersz M., Spurzem R., 2010, MNRAS, 407, 1946
  * Farmer et al. (2019) Farmer R., Renzo M., de Mink S. E., Marchant P., Justham S., 2019, ApJ, 887, 53
  * Fernández et al. (2019) Fernández R., Margalit B., Metzger B. D., 2019, MNRAS, 488, 259
  * Fishbach & Holz (2017) Fishbach M., Holz D. E., 2017, ApJ, 851, L25
  * Foucart (2012) Foucart F., 2012, Phys. Rev. D, 86, 124007
  * Fragione & Kocsis (2018) Fragione G., Kocsis B., 2018, Phys. Rev. Lett., 121, 161103
  * Fragione et al. (2019) Fragione G., Grishin E., Leigh N. W. C., Perets H. B., Perna R., 2019, MNRAS, 488, 47
  * Fragione et al. (2020) Fragione G., Metzger B. D., Perna R., Leigh N. W. C., Kocsis B., 2020, MNRAS, 495, 1061
  * Fragos & McClintock (2015) Fragos T., McClintock J. E., 2015, ApJ, 800, 17
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Fuller & Ma (2019) Fuller J., Ma L., 2019, ApJ, 881, L1
  * Gal-Yam et al. (2009) Gal-Yam A., et al., 2009, Nature, 462, 624
  * Gallegos-Garcia et al. (2022) Gallegos-Garcia M., Fishbach M., Kalogera V., L Berry C. P., Doctor Z., 2022, ApJ, 938, L19
  * Gatto et al. (2021) Gatto M., et al., 2021, MNRAS, 507, 3312
  * Gayathri et al. (2022) Gayathri V., et al., 2022, Nature Astronomy, 6, 344
  * Georgiev et al. (2010) Georgiev I. Y., Puzia T. H., Goudfrooij P., Hilker M., 2010, MNRAS, 406, 1967
  * Giacobbo & Mapelli (2018) Giacobbo N., Mapelli M., 2018, MNRAS, 480, 2011
  * Giacobbo et al. (2018) Giacobbo N., Mapelli M., Spera M., 2018, MNRAS, 474, 2959
  * Gnedin et al. (2014) Gnedin O. Y., Ostriker J. P., Tremaine S., 2014, ApJ, 785, 71
  * Gomez et al. (2019) Gomez S., et al., 2019, ApJ, 881, 87
  * González et al. (2007) González J. A., Sperhake U., Brügmann B., Hannam M., Husa S., 2007, Phys. Rev. Lett., 98, 091101
  * González et al. (2021) González E., Kremer K., Chatterjee S., Fragione G., Rodriguez C. L., Weatherford N. C., Ye C. S., Rasio F. A., 2021, ApJ, 908, L29
  * Harris (2010) Harris W. E., 2010, arXiv e-prints, p. arXiv:1012.3224
  * Harris et al. (2013) Harris W. E., Harris G. L. H., Alessi M., 2013, ApJ, 772, 82
  * Heger & Woosley (2002) Heger A., Woosley S. E., 2002, ApJ, 567, 532
  * Hils & Bender (1995) Hils D., Bender P. L., 1995, ApJ, 445, L7
  * Hoang et al. (2018) Hoang B.-M., Naoz S., Kocsis B., Rasio F. A., Dosopoulou F., 2018, ApJ, 856, 140
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Iorio et al. (2022) Iorio G., et al., 2022, arXiv e-prints, p. arXiv:2211.11774
  * Isoyama et al. (2018) Isoyama S., Nakano H., Nakamura T., 2018, Progress of Theoretical and Experimental Physics, 2018, 073E01
  * Jiménez-Forteza et al. (2017) Jiménez-Forteza X., Keitel D., Husa S., Hannam M., Khan S., Pürrer M., 2017, Phys. Rev. D, 95, 064024
  * Kamlah et al. (2022) Kamlah A. W. H., et al., 2022, MNRAS, 511, 4060
  * Kawamura et al. (2011) Kawamura S., et al., 2011, Classical and Quantum Gravity, 28, 094011
  * Kawamura et al. (2021) Kawamura S., et al., 2021, Progress of Theoretical and Experimental Physics, 2021, 05A105
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Kobayashi et al. (2004) Kobayashi S., Laguna P., Phinney E. S., Mészáros P., 2004, ApJ, 615, 855
  * Kocsis et al. (2012) Kocsis B., Ray A., Portegies Zwart S., 2012, ApJ, 752, 67
  * Kopparapu et al. (2008) Kopparapu R. K., Hanna C., Kalogera V., O’Shaughnessy R., González G., Brady P. R., Fairhurst S., 2008, ApJ, 675, 1459
  * Kremer et al. (2019) Kremer K., et al., 2019, Phys. Rev. D, 99, 063003
  * Kremer et al. (2020) Kremer K., et al., 2020, ApJ, 903, 45
  * Kritos et al. (2022) Kritos K., Strokov V., Baibhav V., Berti E., 2022, arXiv e-prints, p. arXiv:2210.10055
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Küpper et al. (2011) Küpper A. H. W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300
  * Lee & Ramirez-Ruiz (2007) Lee W. H., Ramirez-Ruiz E., 2007, New Journal of Physics, 9, 17
  * Leveque et al. (2022) Leveque A., Giersz M., Banerjee S., Vesperini E., Hong J., Portegies Zwart S., 2022, MNRAS, 514, 5739
  * Licquia & Newman (2015) Licquia T. C., Newman J. A., 2015, ApJ, 806, 96
  * Lousto & Zlochower (2008) Lousto C. O., Zlochower Y., 2008, Phys. Rev. D, 77, 044028
  * Lousto et al. (2012) Lousto C. O., Zlochower Y., Dotti M., Volonteri M., 2012, Phys. Rev. D, 85, 084015
  * Madau & Fragos (2017) Madau P., Fragos T., 2017, ApJ, 840, 39
  * Maliszewski et al. (2022) Maliszewski K., Giersz M., Gondek-Rosinska D., Askar A., Hypki A., 2022, MNRAS, 514, 5879
  * Mapelli et al. (2021) Mapelli M., Santoliquido F., Bouffanais Y., Arca Sedda M. A., Artale M. C., Ballone A., 2021, Symmetry, 13, 1678
  * Mapelli et al. (2022) Mapelli M., Bouffanais Y., Santoliquido F., Arca Sedda M., Artale M. C., 2022, MNRAS, 511, 5797
  * McKernan et al. (2012) McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MNRAS, 425, 460
  * McKernan et al. (2018) McKernan B., et al., 2018, ApJ, 866, 66
  * Metzger & Berger (2012) Metzger B. D., Berger E., 2012, ApJ, 746, 48
  * Mikkola & Aarseth (1993) Mikkola S., Aarseth S. J., 1993, Celestial Mechanics and Dynamical Astronomy, 57, 439
  * Miller & Hamilton (2002) Miller M. C., Hamilton D. P., 2002, ApJ, 576, 894
  * Miller-Jones et al. (2015) Miller-Jones J. C. A., et al., 2015, MNRAS, 453, 3918
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Moore et al. (2019) Moore C. J., Gerosa D., Klein A., 2019, MNRAS, 488, L94
  * Moriya et al. (2022) Moriya T. J., et al., 2022, A&A, 666, A157
  * Nguyen et al. (2014) Nguyen D. D., Seth A. C., Reines A. E., den Brok M., Sand D., McLeod B., 2014, ApJ, 794, 34
  * Nishizawa et al. (2016) Nishizawa A., Berti E., Klein A., Sesana A., 2016, Phys. Rev. D, 94, 064020
  * Nitadori & Aarseth (2012) Nitadori K., Aarseth S. J., 2012, MNRAS, 424, 545
  * O’Leary et al. (2009) O’Leary R. M., Kocsis B., Loeb A., 2009, MNRAS, 395, 2127
  * Peters (1964) Peters P. C., 1964, Physical Review, 136, 1224
  * Peters & Mathews (1963) Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Punturo et al. (2010) Punturo M., et al., 2010, Classical and Quantum Gravity, 27, 194002
  * Qin et al. (2018) Qin Y., Fragos T., Meynet G., Andrews J., Sørensen M., Song H. F., 2018, A&A, 616, A28
  * Qin et al. (2019) Qin Y., Marchant P., Fragos T., Meynet G., Kalogera V., 2019, ApJ, 870, L18
  * Racine (2008) Racine É., 2008, Phys. Rev. D, 78, 044021
  * Rastello et al. (2020) Rastello S., Mapelli M., Di Carlo U. N., Giacobbo N., Santoliquido F., Spera M., Ballone A., Iorio G., 2020, MNRAS, 497, 1563
  * Rastello et al. (2021) Rastello S., Mapelli M., Di Carlo U. N., Iorio G., Ballone A., Giacobbo N., Santoliquido F., Torniamenti S., 2021, MNRAS, 507, 3612
  * Reynolds (2021) Reynolds C., 2021, in 43rd COSPAR Scientific Assembly. Held 28 January - 4 February. p. 1412 (arXiv:2011.08948)
  * Rizzuto et al. (2021) Rizzuto F. P., et al., 2021, MNRAS, 501, 5257
  * Rizzuto et al. (2022) Rizzuto F. P., Naab T., Spurzem R., Arca-Sedda M., Giersz M., Ostriker J. P., Banerjee S., 2022, MNRAS, 512, 884
  * Robson et al. (2019) Robson T., Cornish N. J., Liu C., 2019, Classical and Quantum Gravity, 36, 105011
  * Rodriguez et al. (2015) Rodriguez C. L., Morscher M., Pattabiraman B., Chatterjee S., Haster C.-J., Rasio F. A., 2015, Phys. Rev. Lett., 115, 051101
  * Rodriguez et al. (2016) Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D, 93, 084029
  * Rodriguez et al. (2018) Rodriguez C. L., Amaro-Seoane P., Chatterjee S., Rasio F. A., 2018, Phys. Rev. Lett., 120, 151101
  * Rodriguez et al. (2019) Rodriguez C. L., Zevin M., Amaro-Seoane P., Chatterjee S., Kremer K., Rasio F. A., Ye C. S., 2019, Phys. Rev. D, 100, 043027
  * Romero-Shaw et al. (2019) Romero-Shaw I. M., Lasky P. D., Thrane E., 2019, MNRAS, 490, 5210
  * Romero-Shaw et al. (2020) Romero-Shaw I., Lasky P. D., Thrane E., Calderón Bustillo J., 2020, ApJ, 903, L5
  * Romero-Shaw et al. (2022) Romero-Shaw I., Lasky P. D., Thrane E., 2022, ApJ, 940, 171
  * Romero-Shaw et al. (2023) Romero-Shaw I. M., Gerosa D., Loutrel N., 2023, MNRAS, 519, 5352
  * Rosswog et al. (2009) Rosswog S., Ramirez-Ruiz E., Hix W. R., 2009, ApJ, 695, 404
  * Ryon et al. (2015) Ryon J. E., et al., 2015, MNRAS, 452, 525
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Santamaría et al. (2010) Santamaría L., et al., 2010, Phys. Rev. D, 82, 064016
  * Santoliquido et al. (2020) Santoliquido F., Mapelli M., Bouffanais Y., Giacobbo N., Di Carlo U. N., Rastello S., Artale M. C., Ballone A., 2020, ApJ, 898, 152
  * Sasaki et al. (2016) Sasaki M., Suyama T., Tanaka T., Yokoyama S., 2016, Phys. Rev. Lett., 117, 061101
  * Schrøder et al. (2020) Schrøder S. L., MacLeod M., Loeb A., Vigna-Gómez A., Mandel I., 2020, ApJ, 892, 13
  * Spera & Mapelli (2017) Spera M., Mapelli M., 2017, MNRAS, 470, 4739
  * Spera et al. (2019) Spera M., Mapelli M., Giacobbo N., Trani A. A., Bressan A., Costa G., 2019, MNRAS, 485, 889
  * Spitler & Forbes (2009) Spitler L. R., Forbes D. A., 2009, MNRAS, 392, L1
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters
  * Spruit (2002) Spruit H. C., 2002, A&A, 381, 923
  * Spurzem (1999) Spurzem R., 1999, Journal of Computational and Applied Mathematics, 109, 407
  * Stevenson et al. (2019) Stevenson S., Sampson M., Powell J., Vigna-Gómez A., Neijssel C. J., Szécsi D., Mandel I., 2019, ApJ, 882, 121
  * Stiefel & Kustaanheimo (1965) Stiefel E., Kustaanheimo P., 1965, Journal für die reine und angewandte Mathematik, 218, 204
  * Stone et al. (2017) Stone N. C., Metzger B. D., Haiman Z., 2017, MNRAS, 464, 946
  * Tagawa et al. (2020) Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 898, 25
  * Terreran et al. (2017) Terreran G., et al., 2017, Nature Astronomy, 1, 713
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration et al., 2021a, arXiv e-prints, p. arXiv:2111.03606
  * The LIGO Scientific Collaboration et al. (2021b) The LIGO Scientific Collaboration et al., 2021b, arXiv e-prints, p. arXiv:2111.03634
  * Torniamenti et al. (2022) Torniamenti S., Rastello S., Mapelli M., Di Carlo U. N., Ballone A., Pasquato M., 2022, MNRAS, 517, 2953
  * Vigna-Gómez et al. (2021) Vigna-Gómez A., Toonen S., Ramirez-Ruiz E., Leigh N. W. C., Riley J., Haster C.-J., 2021, ApJ, 907, L19
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Vink et al. (2021) Vink J. S., Higgins E. R., Sander A. A. C., Sabhahit G. N., 2021, MNRAS, 504, 146
  * Wang et al. (2015) Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P., Kouwenhoven M. B. N., Naab T., 2015, MNRAS, 450, 4070
  * Wang et al. (2022) Wang L., Tanikawa A., Fujii M. S., 2022, MNRAS, 509, 4713
  * Weatherford et al. (2020) Weatherford N. C., Chatterjee S., Kremer K., Rasio F. A., 2020, ApJ, 898, 162
  * Webb & Leigh (2015) Webb J. J., Leigh N. W. C., 2015, MNRAS, 453, 3278
  * Woosley & Heger (2021) Woosley S. E., Heger A., 2021, ApJ, 912, L31
  * Woosley & Smith (2022) Woosley S. E., Smith N., 2022, ApJ, 938, 57
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Ye et al. (2020) Ye C. S., Fong W.-f., Kremer K., Rodriguez C. L., Chatterjee S., Fragione G., Rasio F. A., 2020, ApJ, 888, L10
  * Zevin et al. (2020) Zevin M., Spera M., Berry C. P. L., Kalogera V., 2020, ApJ, 899, L1
  * Zevin et al. (2021) Zevin M., et al., 2021, ApJ, 910, 152
  * Ziosi et al. (2014) Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, MNRAS, 441, 3703
  * du Buisson et al. (2020) du Buisson L., et al., 2020, MNRAS, 499, 5941