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.

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 ﬁnd 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, ∼ 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 R ∼ 12(4 . 4)(1 . 2) yr − 1 Gpc 3 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 ﬁnd that surveys with a limiting magnitude m bol = 25 can detect ∼ 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.


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 compo-Contact e-mail:manuel.arcasedda@gssi.itnent 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 mgap = 40−150 M 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 χBH ∼ 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 χBH 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) where q < 1 is the binary mass ratio, χ1,2 are the two component spin vectors and L is the binary orbital angular momentum.Observations of BBH mergers suggest that χ eff may increase at increasing the binary merger mass ratio, although some merging binaries exhibit a negative value of χ 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 2018Banerjee , 2021Banerjee , 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 semianalytic 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 M pc −3 ) and massive (> 100, 000 M ) 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 ρ = (1.2 × 10 4 − 1.5 × 10 7 ) M 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 code 1 .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.

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) × 10 5 , half-mass radius RHM = (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.
Each cluster is modelled according to a King (1966) profile with adimensional potential well W0 = 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 mZAMS = (0.08−150) M .Stars in primordial binaries are paired depending on their mass, with stars heavier than > 5 M paired according to 1 https://github.com/nbody6ppgpu/Nbody6PPGPU-beijinga flat mass-ratio distribution, and lighter stars paired randomly.Binary eccentricities are distributed according to a thermal distribution, P (e)de = e 2 de, 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 gal = 1.78×10 11 M .All Dragon-II clusters are placed on a circular orbit around this galaxy model at a distance of R 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 mBH 15 M , there are no BHs with a mass above 30 M , or the simulated time exceeds at least one relaxation time.As a result, the simulation time in Dragon-II models spans a range Tsim = 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.
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 follow-ing 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 αCE − λCE scheme, which enables us to regulate the fraction of orbital energy injected into the envelope (αCE) and to scale the binding energy of the envelope by a factor λCE.In this work, we adopt αCE = 3 and λ = 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 mBH,max = 40.5 M (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 fc (Rizzuto et al. 2021(Rizzuto et al. , 2022)), which we set to fc = 0.5.We refer the reader to Rizzuto et al. (2021Rizzuto et al. ( , 2022) ) for a discussion about the impact of fc on BH evolution.

Modelling the final stages of compact object binary mergers
The dynamics of relativistic binaries is followed via the orbitaverage 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, 2020a,b, 2021;Rizzuto et al. 2021Rizzuto et al. , 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 dec = 2kGm bin /c 2 = kR 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 M ), star cluster masses (< 10 6 M ) and half-mass radii (0.1 − 3 pc) explored in this work, it is easy to show that the binary-single interaction timescale t2−1 = (nσΣ) −1 -with n the cluster density, σ the cluster velocity dispersion, and Σ the binary cross section -is generally > 10 8 larger than the binary inspiral timescale, tinsp ∝ a 4 /(m1m2m bin ) (Peters & Mathews 1963).Moreover, the typical merger time for a binary with mass m bin < 200 M and separation a dec is generally tinsp < 100 yr, i.e. much smaller than the cluster crossing time, tstep ∼ 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) 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 χ = 0.5 and dispersion σχ = 0.2, and Maxwellian distribution with dispersion σχ = 0.2.
In this work, whenever spins are taken into account during the simulation we assume a Gaussian distribution with χ = 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: Here, η ≡ qBBH/(1 + qBBH) 2 is the symmetric mass ratio, Ξ ≡ 2( S2 + q 2 BBH S1)/(1 + qBBH) 2 , and the subscripts ⊥ and mark the perpendicular and parallel directions of the BH spin vector ( S) with respect to the direction of the binary angular momentum.We assume A = 1.2 × 10 4 km s −1 , B = −0.93,H = 6.9 × 10 3 km s −1 , and ξ = 145 • González et al.  (2007); Lousto & Zlochower (2008), V11 = 3677.76km s −1 , and VA,B,C = (2.481,1.793, 1.507) × 10 3 km s −1 .The vector ∆ is defined as ∆ The angle between the direction of the infall at merger and the in-plane component of ∆, i.e. φ∆, is drawn from a uniform distribution, while φ1 = 0 − 2π, 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 + 1th 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 M 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 M .

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 RHM = 0.8pc and N = 300k, 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 M .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

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 (mBH < 40.5 M ), 12 involve two mass-gap BHs, and 21 involve one BH below the gap and a massgap BH.Six BH-BH mergers involve a primary with a mass mBH,1 = (5.4−7.1)M and a companion with mass mBH,2 = (2.55 − 3.6) M , 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 fc = 0.5 of the star mass is accreted onto the BH during an accretion event (for a discussion about the impact of fc 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 clus-  ter 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, IMBHs2 , 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 mBH,1 > 30 M and a companion with a mass in the range mBH,2 = (20 − 50) M .Conversely, mergers occurring outside the cluster -or ejected mergers -are characterised by a mass-ratio q > 0.6 and a primary mass typically mBH,1 < 40 M .
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 mBH,1 = (20 − 50) M and mBH,2 = (20 − 40) M , 3) mergers with one component previously in a primordial binary are characterised by a heavy primary, mBH,1 > 40 M , and a heavy companion, mBH,2 > 20 M .A similar trend is observed in recent Nbody simulations tailored to relatively light star clusters, i.e. with mass < 8, 000 M (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.
Dynamics and binary stellar evolution deeply impact also the properties of stellar-size mergers.For example, "dynami-  cal" 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 m1 = (5 − 40) M , with the mass-ratio weakly increasing at increasing the primary mass: note that for m1 15 M 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 m1 > 40.5 M .In this case, the mass ratio decreases with the primary mass as expected from the mass function limit, with companion masses in the range m2 = (30 − 50) M .
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 m1,2 > 30 M 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 mBH,1 20 M 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 mBH,1 > 20 M , 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, 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).

Delay times
The delay time of Dragon-II mergers (tGW), 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 tGW (0.5, 1.5, 10) Gyr.However, when the delay time is normalised to the initial half-mass relaxation time (t rlx ) of the cluster, the overall tGW nicely distribute around a main peak located at tGW/t rlx 8 − 30.The exact location of the peak depends on the definition of t rlx .For the sake of clarity, in the plot we use three different expressions of t rlx taken from Gatto et al. ( 2021 The three peaks that appear in the tGW distribution find a clear explanation looking at the tGW/t rlx distribution.In fact, the first peak at tGW = 500 Myr corresponds to mergers happening in simulations with t 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.

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 suffi-ciently 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 subpopulation 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 mBH1,2 = (7.1 + 2.55) M 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 (∼ 3 km/s) that triggers the eccentricity increase.
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)  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 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, ex- plaining 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. e = tmerg 0 edt/ tmerg 0 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) The step-like behaviour of the e − fp 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 e > 0.7 across the whole 0.01-100 Hz frequency spectrum, whilst ejected mergers preserve a moderate eccentricity e < 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 de- tectors more suited to observe GW sources merging outside their parent cluster.

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 = 120k, RHM = 0.47pc, the dynamical BH mergers develop in a simulation with N = 300k, RHM = 0.47pc, and the one forming from a primordial binary in a simulation with N = 120k, RHM = 0.8pc.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.

White dwarf -black hole mergers: implications for low-mass X-ray binaries
The WD-BH binary consists in a BH with mass mBH = 23.1 M and a carbon-oxygen white dwarf (COWD) with mass mWD = 1.18 M .Initially, dynamical interactions pair the BH with the WD progenitor, a MS star with mass mMS,pro = 4.89 M .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 ∼ 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 ∼ 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.
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 mBH > 9 M (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.

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 mNS = 1.28 M and a BH with mass mBH = 14.96M .The BH, whose progenitor had a mass mMS = 26.7 M , 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 fGW = 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 eNSBH = 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, M chirp = 3.4 M , 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 m1,2 = (26.3+ 18.7) M 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 mBH = 5.6 M and a NS mass mNS = 1.88 M .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 mBH = 5.7 +1.8  −2.1 M and a NS with mNS = 1.5 +0.7 −0.3 M .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 ∼ 7 × 10 4 yr.At decoupling, the binary peak frequency is fGW ∼ 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 mBH/mNS, the BH spin χ, and the NS compactness C ≡ GmNS/c 2 RNS (Foucart 2012).As numerical simulations have shown, in general the larger the mBH/mNS 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 exceeds the BH innermost stable circular orbit (ISCO), which for a spinning BH can be expressed as (Bardeen et al. 1972) where Z1,2 are functions of the BH adimensionless spin parameter χ.Whilst the condition R tid /RISCO < 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 mBH/mNS = 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 χ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 mBH/mNS = 2.98 and C = 0.27.Even in this case, the spin required for an accretion disk to form is χ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 mBH ∼ 12 M and mNS = 1.6 M (Broekgaarden et al. 2021), which implies a minimum BH spin of χBH 0.8 to permit the formation of a fairly massive (mass > 0.1mNS) disk around the BH (Fryer et al. 2012).

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 (χ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 χ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 χBH = 0.01 and to the other χ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 χBH = 0.5 and dispersion σχ = 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, θi, such that (Arca Sedda et al. 2023) We set n θ = 8, which implies that binaries have a 20(55)% probability to have θ1,2 that differ by less than 5(20)%.Note that n θ = 0 implies the isotropic distribution whilst n θ 1 implies fully aligned spins, i.e. θ1 = θ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 χ eff > 0, and dynamical binaries, for which χ eff ∼ 0. The plots suggest that the STEV model, based on stellar evolution models, leads primordial binaries to have a χ 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 χ eff ∼ 0 and m bin = (40 − 100) M , 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 χ eff values and mBH < 40 M , and one widely distributed around zero χ eff involving massive BHs, mBH > 40 M .
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 M if the Dragon-II BH merger mass is below the upper-mass gap, and in the range 40.5 − 100 M 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 χ eff distribution.Figure 11 shows the distribution of χ 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).

Merger efficiency
As described in the previous section, we have simulated a total Dragon-II mass of Msim = 3.65 × 10 6 M 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 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 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 ∼ 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, GW(RHM), as a function of the average cluster density ρsim , using the following definitions where Msim is the total simulated mass and Nsim is the number of simulations performed for a given value of the halfmass radius, RHM.At fixed value of the binary fraction, this relation is well described by a power-law in the form GW = a( ρHM /1M pc −3 ) b , with a = (0.15 ± 0.07) × 10 −5 and b = 0.25 ± 0.03.
The plot makes clear that increasing the cluster density by two orders of magnitude leads to ∼ 2.5× more mergers.Moreover, it further highlights the role of primordial binaries, showing that clusters with a lower binary fraction have a probability ∼ 50% smaller to develop a merger, at least in the case of RHM = 1.75 pc.

Merger rate for black hole binaries
We define the cosmic merger rate following Santoliquido et al. (2020); Bouffanais et al. ( 2021) where t lb (z) is the lookback time at merger, ψ clus (z ) is the star cluster formation rate when the merging binary formed, ηGW(Z) is the merger efficiency at the metallicity Z, F(z , z, Z) is the number of mergers forming at redshift z 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. 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 i.e. we rescale the stellar star formation rate derived by (Madau & Fragos 2017) by a factor fCFE, 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 fCFE,YC = 0.3 for young clusters (Mapelli et al. 2021) and fCFE,GC = 0.08±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 η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     thus assuming that the number of mergers at redshift z that formed at z is independent on the metallicity distribution.The p(Z, z ) term represents the cosmic fraction of clusters with metallicity in the (Z,Z + dZ) bin at redshift z .We assume that the metallicity follows a log-normal distribution peaked at (Madau & Fragos 2017) with dispersion either σ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 metalrich 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 fCFE = 1, i.e. that all stars form in star clusters, and calculate a merger rate of R = 27 yr −1 Gpc −3 , in broad agreement with the rate inferred for low-mass star clusters (N = 10 2 −5×10 4 M ) (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, fCFE = 0.08 − 0.3.As shown in the plot, we find a simulated merger rate of RGW = (12 ± 7) at redshift z = 0.2.At the same redshift, the BBH merger rate inferred by the LVK is RLVK = 17.9 − 44 yr −1 Gpc −3 (The LIGO Scientific Collaboration et al. 2021a,b).

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: where Mg * is the galaxy stellar mass, δ = 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 where ρ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 rel = 10 9 yr (Harris 2010), or massive and relatively young clusters in the Small Magellanic Cloud (SMC), t rel = 3.2×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 case of NS-BH mergers, instead, the event occurs over a timescale of tGW = (0.04 − 0.5)t rel .The fraction of cluster with an age longer than tGW is fx ∼ 0.94 for clusters in both the Milky Way and the SMC, the resulting frequency rate for NS-BH mergers is RNSBH = (0.13 − 40.7) yr −1 , which implies a volumetric merger rate of RNSBH = (0.027 − 8.7) yr −1 Gpc −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) where we have assumed RWD = 10 4 km.Note that an eccentricity between 0 and 1 would affect fGW 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) 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) where Sc 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-tonoise ratio of (S/N)> 8(15) is D < 1.5 Mpc(0.7 Mpc).
Note that the Andromeda galaxy is ∼ 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 22N (< D) = 2 and find an upper limit to the local merger rate of closeby WD-BH mergers of R WDBH,close < (8.4 × 10 −10 − 5.1 × 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) M .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 mZAMS mPISN = 150 M 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 mPISN.Typical masses for Dragon-II PISN progenitors are in the range (150 − 282) M .
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: where r = −5.02× 104 , s = −2.159,and t = 2.887 (Heger & Woosley 2002;du Buisson et al. 2020), and M 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) which can be converted into an apparent bolometric magnitude via the Pogson's relation being DL 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 He,f = 90.4M as extracted from our models.The value of µ GW bol 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 µ GW bol and a magnitude detection threshold µ lim , we define the fraction of detectable sources as where we adopted σµ = 0.23 .The PISN rate as a function of the redshift can thus be evaluated as: where dV /dz is the comoving volume element and ψ(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 fCFE described in Section 4.2, and that only stars with a metallicity Z 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 µ 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 (SDSS 4 ) or the Palomar Transient Factory (PTF5 ), ongoing, e.g. the Dark Energy Survey (DES 6 ), and future surveys, like the Large Synoptic Survey Telescope (LSST7 ), the Zwicky Transient Facility (ZTF8 ), or the EUCLID mission9 (Moriya et al. 2022).From Figure 14 we see that only future surveys (µ bol 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, µ bol , and the cluster formation efficiency, fCFE, 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 ∼ 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.

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 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 mBH,1 > 30 M and companion in the mBH,2 = 30 − 50 M mass range, whilst ejected mergers involve lighter primaries, mBH,1 < 40 M , 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 M = 3.4 M , 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 χ > 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 ∼ 2 × 10 −5 M −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 ∼ 2.Moreover, the merger efficiency increases with the cluster density following a power-law with slope ∼ 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 R = 5 − 19 (0.027 − 8.7) (3.8 × 10 −4 − 2.3) yr −1 Gpc −3 for BBHs(WD-BH)(NS-BH) mergers, respectively.We predict that, in a 4 yr-long mission, LISA could detect NBBH = 12 ± 7(5 ± 3) BBH mergers (IMRIs) and can identify the WD-BH merger with a signal-to-noise ratio SNR> 8(15) if it occurs within DL < 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 bol = 25 could detect NPISN = 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 ρ = (1.2 × 10 4 − 1.6 × 10 6 ) M 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.

Figure 1 .
Figure1.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).
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 BH > 30 M or > 40 M at the last simulation snapshot.

Figure 2 .
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 BH,2 = 3, 5, 10, 30, 50, 100 M .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).

Figure 3 .
Figure3.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).

Figure 4 .
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).

Figure 5 .
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).

Figure 6 .
Figure6.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.

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

Figure 8 .
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).

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

Figure 10 .
Figure10.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.

Figure 11 .
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).

Figure 12 .
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.

Figure 14 .
Figure 14.PISNe cosmological rate (blue straight line) and detection rate assuming a magnitude limited survey at different threshold values, namely m bol = 17 (red dotted line), 20 (green dasheddotted line), and 25 (orange dashed line), assuming a cluster formation efficiency of f CFE = 1 and assuming that only clusters with a metallicity Z 0.008 can host stars undergoing PISN.

Table 2 .
Volumetric merger rate calculated in the local Universe for different GW sources in Dragon-II clusters.
Figure13.Merger rate density for all simulations in our runs, assuming either metal-poor only (red straight line) or both metalpoor 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 CFE = 0.3 and the lower one corresponding to f CFE = 0.08.

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.