The Dragon-II simulations - I. Evolution of single and binary compact objects in star clusters with up to 1 million stars

We present the ﬁrst results of the Dragon-II simulations, a suite of 19 N -body simulations of star clusters with up to 10 6 stars, with up to 33% of them initially paired in binaries. In this work, we describe the main evolution of the clusters and their compact objects (COs). All Dragon-II clusters form in their centre a black hole (BH) subsystem with a density 10 − 100 times larger than the stellar density, with the cluster core containing 50 − 80% of the whole BH population. In all models, the BH average mass steeply decreases as a consequence of BH burning, reaching values (cid:104) m BH (cid:105) < 15 M (cid:12) within 10 − 30 relaxation times. Generally, our clusters retain only BHs lighter than 30 M (cid:12) over 30 relaxation times. Looser clusters retain a higher binary fraction, because in such environments binaries are less likely disrupted by dynamical encounters. We ﬁnd that BH-main sequence star binaries have properties similar to recently observed systems. Double CO binaries (DCOBs) ejected from the cluster exhibit larger mass ratios and heavier primary masses than ejected binaries hosting a single CO (SCOBs). Ejected SCOBs have BH masses m BH = 3 − 20 M (cid:12) , deﬁnitely lower than those in DCOBs ( m BH = 10 − 100 M (cid:12) ).


INTRODUCTION
Massive star clusters in the range (10 4 − 10 6 ) M , like globular clusters or young massive clusters, represent galactic repositories of stellar compact objects, and are ideal laboratories to study the interplay of stellar evolution and dynamics.Several hundreds of stellar black holes (BHs), neutron stars (NSs), and white dwarfs (WDs) are expected to form in a typical massive cluster.
The progress in stellar evolution of massive stars (Woosley et al. 2007;Wang et al. 2015;Giacobbo & Mapelli 2018;Belczynski et al. 2010Belczynski et al. , 2016;;Spera et al. 2019;Vink et al. 2021), partly triggered by the discovery of gravitational-wave (GW) emission by merging BH and NS binaries (The LIGO Scientific Collaboration et al. 2021a,b), has completely changed our understanding of BHs.Stellar models demonstrated that the evolution of single massive stars is significantly influ-enced by the possible development of so-called pair instability supernovae (PISN), which causes the complete disruption of stars that develop an He core with a mass of MHe = 64 − 135 M , and pulsational pair instability supernovae (PPISN), a mechanism that leads to an enhanced mass-loss in stars with a He core mass of MHe = 32 − 64 M .This leads to a maximum stellar BH mass in the range mBH,max = (40−60) M , depending on the theoretical model adopted and the stellar metallicity.Direct consequence of these two processes is the well known upper-mass gap of BHs, a region of the mass-spectrum where no remnants are expected (Woosley et al. 2007).The boundaries of the uppermass gap are highly uncertain and depend on the adopted stellar evolution model and metallicity (Woosley et al. 2007;Wang et al. 2015;Belczynski et al. 2016;Spera & Mapelli 2017;Vink et al. 2021;Costa et al. 2021;Iorio et al. 2022).Only stars with a zero age main sequence mass beyond MZAMS > (200 − 250) M can avoid PISN and, depending on their metallicity, directly collapse to an intermediate-mass BH with little mass loss in the process (see e.g.Spera & Mapelli 2017).Stellar collisions might lead to the formation of BHs in the upper-mass gap (e.g.Spera et al. 2019), thus suggesting that star clusters could be perfect laboratories to form mass-gap BHs (e.g.Di Carlo et al. 2019;Kremer et al. 2020a;Rizzuto et al. 2021;Rastello et al. 2021;Rizzuto et al. 2022;Banerjee 2022), but it is unclear how the stellar merger frequency depends on the cluster initial properties (Rizzuto et al. 2022) or the stellar conditions at merger (Ballone et al. 2023;Costa et al. 2022).
More in general, the formation of a population of compact objects can significantly affect star cluster dynamics.Massive stars and BHs rapidly sink into the cluster centre via mass-segregation, possibly forming a massive subsystem on a core-collapse timescale (e.g.Spitzer 1987;Breen & Heggie 2013;Pavlík et al. 2018;Arca Sedda et al. 2018;Leveque et al. 2022a) which can contract and determine the onset of runaway stellar collisions if the time does not exceed the stellar evolution timescale (e.g.Spitzer 1987;Portegies Zwart & McMillan 2002;Portegies Zwart et al. 2004;Fregeau et al. 2004;Breen & Heggie 2013;Mapelli 2016;Giersz et al. 2015;Arca Sedda et al. 2018;Vergara et al. 2021Vergara et al. , 2022;;Maliszewski et al. 2022).The runaway growth of a massive star can be hampered by the formation of tight binaries that supply energy to the cluster core, cause BH ejection, deplete the cluster's BH reservoir, and eventually kick each other out via super-elastic encounters (Breen & Heggie 2013).
The competing effect of binary energy supply and stellar collisions likely depends on the cluster mass, density, metallicity, the fraction of primordial binaries, the initial mass function and its boundaries, the natal kicks of BHs and NSs, and the compact object mass spectrum.Typically, the exploration of a tiny part of such parameter space is performed with numerical models capable of simultaneously accounting for stellar dynamics and evolution, either via direct N -body (e.g.Wang et al. 2015;Banerjee 2018Banerjee , 2021;;Di Carlo et al. 2019;Rastello et al. 2021;Di Carlo et al. 2021;Chattopadhyay et al. 2022) or Monte Carlo techniques (e.g.Rodriguez et al. 2016;Askar et al. 2017;Rodriguez et al. 2019;Kremer et al. 2020a;Maliszewski et al. 2022).
Direct N -body simulations offer most likely the highest level of accuracy in terms of stellar dynamics modelling, but their computational cost forced the vast majority of works in the literature to focus on star clusters with less than a few × 10 5 stars and/or with a relatively small fraction of primordial binaries (Banerjee 2018(Banerjee , 2021;;Di Carlo et al. 2019;Rastello et al. 2021;Di Carlo et al. 2021), with a few notable exceptions.For example, several works have explored the impact of a large primordial binary fraction, up to 100%, on the dynamics of isotropic (e.g.Heggie et al. 2006;Trenti et al. 2007;Pavlík 2020) and anisotropic (Pavlík & Vesperini 2022a,b) low-mass star cluster models, i.e. with N < 20, 000, with equal-mass stars, and recently in intermediate-mass GCs, i.e.N ∼ 10 5 (Wang et al. 2022).With regards to simulations tailored to represent massive globular clusters, the DRAGON simulations remain the only one that exploited 10 6 particles (Wang et al. 2016).
Since the development of such pioneering simulations, and especially after the discovery of GWs, numerical tools underwent major upgrades in terms of stellar evolution and treatment of relativistic binaries.
In this work, we present the Dragon-II simulation database, a suite of 19 direct N -body simulations performed with the Nbody6++GPU code1 representing star clusters with N = (0.12 − 1) × 10 6 stars, half-mass radius densities in the ρ h = 1.3 × 10 4 − 6.9 × 10 6 M pc −3 range, and a fraction f 2b = 0.10 − 0.33 of stars initially paired in primordial binaries.This work, which is the first one of a series, focuses on the evolution of single and binary BHs and compact objects in massive and dense star clusters, paying particular attention to the relation between the BH population (mass, average BH mass, density) and the cluster properties (mass, radius).Our Dragon-II models explore a portion of the parameter space still uncharted by direct N -body simulations, thus complementing previous works that either rely on Monte Carlo simulations or exploit star cluster models with old stellar evolution recipes or a significantly smaller number of stars.
The paper is organised as follows: Section 2 describes the main properties of the Dragon-II clusters and the improvements integrated in the Nbody6++GPU code; Section 3 presents our main results in terms of overall star cluster evolution (Section 3.1), main properties of single and binary compact objects (Sections 3.2 -3.5), and the possible implementation of N -body outputs into semi-analytic tools (Section 3.6); whilst Section 4 is devoted to summarise the main outcomes of our work.
In the last few years, the code underwent a series of major upgrades related to the treatment of relativistic compact objects (Rizzuto et al. 2021), the implementation of flexible stellar evolution recipes (Kamlah et al. 2022a), and the inclusion of a dedicated treatment for spins (Banerjee et al. 2020;Kamlah et al. 2022a, this work).Here, we expand the possible choices for BH natal spin distribution and implement relativistic recoil for post-merger remnants.In the following, we briefly summarize the features of the code that are most relevant for this work, and discuss the newest upgrades that we implemented into the code and use here for the first time.

Stellar evolution
Nbody6++GPU implements stellar evolution for single and binary stars via the SSE and BSE routines (Hurley et al. 2000(Hurley et al. , 2002)), which we heavily updated to include up-to-date prescriptions for the evolution of massive stars.We refer the reader to Kamlah et al. (2022a) for a comprehensive discussion about the updated stellar evolution encoded in Nbody6++GPU .
In this work, we adopt the level-B of stellar evolution as defined in (Kamlah et al. 2022a, , see their Table A1).This implies that our models take into account the formation of electron-capture supernovae (ECSNe, following Belczynski et al. 2008), the delayed SN scheme (Fryer et al. 2012), and the development of pair-instability (PISN) and pulsational pair instability supernovae (PPISN) (Belczynski et al. 2016).For the formation of compact objects, we adopt mass loss from Belczynski et al. (2010) with additional metallicitydependent correction factors taken from Vink et al. (2001) and a dedicated treatment for mass loss of hot and massive Hrich O/B stars (Vink et al. 2001).The adopted stellar evolution models imply that the maximum BH mass attainable by massive stars with zero-age main-sequence mass < 150 M is mBH,max = 40.5 M (Belczynski et al. 2016;Banerjee et al. 2020;Kamlah et al. 2022a).The BHs falling in the so-called upper mass-gap can still form via stellar collisions, accretion of stellar material onto stellar BHs, and BH-BH mergers, as we discuss in our companion papers.
Natal kicks for NSs forming via ECSNe, accretion induced collapse (AIC), and merger-induced collapse (MIC) are drawn from a Maxwellian distribution with dispersion 3 km/s (see Gessner & Janka 2018;Kamlah et al. 2022a), whilst for all other NSs we adopt a Maxwellian distribution with dispersion 265 km/s (Hobbs et al. 2005).This latter value is adopted also for BHs, but the kick amplitude is reduced by a factor that accounts for the amount of fallback material (Fryer et al. 2012).
For binary stars, we model common envelope evolution via the parametrised αCE − λCE scheme, according to which it is possible 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 a way similar, but not equal, to the one followed by Claeys et al. (2014) (further details about these parameters are discussed in Kamlah et al. 2022a).In this work, we adopt αCE = 3 (Giacobbo & Mapelli 2018;Kamlah et al. 2022a).

Dynamics of compact objects
In particularly dense clusters, stellar interactions can trigger collisions among stars and/or compact objects.The aftermath of such collisions is still a poorly understood process that can crucially affect the formation and evolution of stellar BHs.Whilst the outcome of stellar mergers is better understood, also thanks to recent detailed hydrodynamical simulations coupled with stellar evolution models (Ballone et al. 2023;Costa et al. 2022), it is still unclear how much mass a massive star can accrete onto a stellar BH.Several works have shown that in the case of a star with a mass ∼ (1 − 10) M merging with a stellar BH, there is little accretion as most of the energy is radiated away via jets, although the mechanism is highly uncertain and likely depends on the star structure and evolutionary stage (Guillochon & Ramirez-Ruiz 2013;MacLeod & Ramirez-Ruiz 2015;Cruz-Osorio & Rezzolla 2020;Kremer et al. 2022).Hydrodynamical simulations of star-BH close encounters have shown that up to 70% of the star mass remains bound to the BH, but energy arguments suggest that even a tiny amount of accreted matter, O(10 −3 − 10 −2 M ) would suffice to evaporate the accretion disk and halt the BH growth (Kremer et al. 2022).Nonetheless, recent simulations modelling the common envelope phase of a tight star-BH binary have shown that the BH accretes the stellar core and expels the envelope, a process accompanied by a SN-like transient and spin-up of the BH to nearly extreme values regardless of the initial spin (Schrøder et al. 2020).In multiple main-sequence star collisions, the merger product is expected to have a compact core and a tenuous envelope with densities as low as 10 −10 g cm −3 (Glebbeek et al. 2009).Therefore, if: a) most of the merger product mass is in the core (Glebbeek et al. 2009), and b) the core can efficiently feed the BH (Schrøder et al. 2020), it is reasonable to assume that a BH would accrete a significant fraction of it.
Given the aforementioned uncertainties, in Nbody6++GPU we parametrise the outcome of star-BH collisions via the fraction of star mass accreted onto the BH, fc (Banerjee 2021;Rizzuto et al. 2021).Throughout this paper we adopt fc = 0.5.
Nbody6++GPU also features a treatment for compact binary mergers based on an orbit-averaged formalism (Peters & Mathews 1963;Peters 1964), which enables us to follow the formation and evolution of in-cluster compact binary mergers, a feature implemented in a number of recent works modelling young star clusters (Di Carlo et al. 2019, 2020a,b, 2021;Rizzuto et al. 2021Rizzuto et al. , 2022;;Rastello et al. 2021).
In this work, we present the implementation of three new features of the code: mass and spin of the merger remnant, calculated via numerical relativity fitting formulas (Jiménez-Forteza et al. 2017;Arca Sedda et al. 2020), and the recoil kick imparted by asymmetric GW emission promptly after merging events (Campanelli et al. 2007;Lousto & Zlochower 2008;Lousto et al. 2012).We follow the implementation depicted in our previous works (e.g.Arca Sedda et al. 2020;Arca-Sedda et al. 2021b).

Massive star cluster models with up to one million stars
We generate the 19 Dragon-II star clusters with the updated McLuster software (Küpper et al. 2011), as described in Kamlah et al. (2022a) and Leveque et al. (2022b).
All Dragon-II star clusters are modelled via King (1966) dynamical models with a central dimensionless potential well W0 = 6, and are characterised by three values of the halfmass radius, RHM = 0.47, 0.80, 1.75 pc, four values of the initial number of stars, N = (1.2, 3, 6, 10) × 10 5 , and two values of the primordial binary fraction, as described below.All clusters have the same metallicity Z = 0.0005, a value typical of several clusters proposed to host a dense subsystem of stellar BHs, like NGC3201 or a central intermediate-mass black hole (IMBH), like NGC6254 (see e.g.Askar et al. 2018;Arca Sedda et al. 2018;Weatherford et al. 2020).
All simulations were conducted on the Juwels BOOSTER supercomputer and the GRACE HPC workstation over a ∼ 2 yr timespan.Eventually, the whole Dragon-II database consists of almost 35 Tb of data.
Stellar masses are drawn from the Kroupa (2001) initial mass function limited between m * = 0.08 − 150 M , which implies an initial average stellar mass is m * 0.59 M .
The corresponding initial mass and density scale in Dragon-II clusters are Mc = (0.7 − 5.9) × 10 5 M and densities ρc 1.3 × 10 4 − 6.9 × 10 6 M pc −3 , respectively.All Dragon-II clusters move on a circular orbit at a distance of 13.3 kpc from the centre of a galaxy whose gravitational potential is modelled via a simple Keplerian potential assuming a total galaxy mass of Mg = 1.78 × 10 11 M .As a consequence, our clusters have initially a tidal radius in the range R tid = 67 − 138 pc and they can all be considered as underfilling systems, thus the gravitational field has a smaller impact on the cluster evolution with respect to internal dynamics, at least at the beginning.Dragon-II clusters would underfill their Roche lobe even in the case of a rather extremely eccentric orbit, e.g. e = 0.9.
We assume that a fraction of the total number of stars is initially paired in a primordial binary system.Following in McLuster , we define the binary fraction as the ratio between the number of binaries and the sum of single stars and binaries, f b = n b /(ns + n b ).We set a f b = 0.05 − 0.2 depending on the cluster model as summarized in Table 1.Our simulation grid contains two sets that differ only in f b , thus their comparison could unveil some effects triggered by primordial binary dynamics.Note also our definition of f b implies that the number of stars in binaries over the total is Binaries are initialised assuming the same mass function of single stars and a uniform mass ratio distribution in the range q = 0.1 − 1 for stars heavier than m * > 5 M or random pairing for the lighter ones (Kiminki & Kobulnicky 2012;Sana et al. 2012;Kobulnicky et al. 2014).Following previous works on the same topics, we adopt a thermal distribution of the eccentricity and a semi-major axis distribution flat in logarithmic values, with an upper limit set to 50 AU and a lower limit set by the sum of the stars' radii (Wang et al. 2015;Kamlah et al. 2022a).
In the majority of the cases, for each value of RHM and N we run two simulations with different random seeds to explore possible dependencies on the randomness of the star distribution.The only exception is the case RHM = 0.47 pc and N = 300k stars, which was limited to only one model because of the available computational time.
The simulations are performed until either the average mass of stellar BHs falls below mBH 15 M , no BHs with a mass above 30 M are retained in the cluster, or the simulated time exceeds at least one relaxation time (Spitzer 1987;Binney & Tremaine 2008;Gatto et al. 2021), which can be expressed in the form (Gatto et al. 2021) where γn = 0.11 − 0.4 for a monochromatic mass spectrum (Giersz & Heggie 1996;Binney & Tremaine 2008) but it can be as low as γn = 0.02 for a multi-mass mass spectrum (Giersz & Heggie 1996).These choices result in a physical simulated time ranging between Tsim ∼ 0.1 − 2.3 Gyr and lead to an optimal balance between the computational cost of the simulations and the portion of parameter space that can be explored.Table 1 summarizes the main properties of Dragon-II models.
As sketched in Figure 1, in comparison to the most recent studies based on N -body (e.g.Wang et al. 2015;Banerjee 2018Banerjee , 2021;;Di Carlo et al. 2019;Rastello et al. 2021;  Di Carlo et al. 2021) and Monte Carlo simulations (e.g.Rodriguez et al. 2016;Askar et al. 2017;Rodriguez et al. 2019;Kremer et al. 2020a;Maliszewski et al. 2022), the Dragon-II clusters occupy a region of the N -ρ h plane mostly populated by Monte Carlo simulation grids.This, coupled with the fact that simulations with N > 10 5 stars usually adopt a binary fraction < 20%, makes our Dragon-II simulations an unprecedented grid of models that complements, and further expands, the phase space accessible with direct N -body models.

Star cluster evolution
The Dragon-II clusters were originally devised to explore compact object dynamics, compact binary mergers, and intermediate-mass black hole build-up in dense star clusters, thus they are not meant to be representative of any observed cluster.Nonetheless, it is interesting to compare in Figure 2 the time evolution of the modelled mass and half-mass radius with relatively young, i.e. typical ages 0.1 − 1 Gyr, massive star clusters in the Milky Way (MW), the Small (SMC) and Large Magellanic Cloud (LMC), M31 (Portegies Zwart et al. 2010;Gatto et al. 2021), the Henize 2-10 starburst dwarf galaxy (Nguyen et al. 2014), and the M83 galaxy (Ryon et al. 2015).Over the simulated time, our models overlap with observed clusters, thus indicating that the adopted initial conditions lead to numerical models that can represent one possible evolutionary pathway of some observed clusters.We find that the mass and half-mass radius evolution is well described by the following relations: The values of the fitting parameters, which are summarised in Table 2, are independent of the initial cluster mass, and weakly depend on the initial value of the half-mass radius.This owes to the fact that the mass-segregation time scales with M HM , thus it is mostly affected by the choice of the half-mass radius.
Figure 3 shows the ratio between the final and initial values of RHM as a function of the simulated time, normalised to the initial relaxation time.The plot clearly highlights how the cluster expansion depends only on the dynamical age of the cluster, regardless of the initial cluster mass.By the end of the simulations, our clusters have typically lost ∼ 25 − 50% of their initial mass and their radius has expanded by a factor of 1.5−10, thus implying a reduction of the density at the half-mass radius by up to four orders of magnitude and a reduction of the velocity dispersion of around 1−1.5 times.The drop in density and velocity dispersion crucially affects the rates at which dynamical interactions take place.
A thorough comparison among Dragon-II simulations and the models discussed in the past literature is made hard by the many different assumptions of previous works, like the use of equal-mass stars to represent the cluster, the different binary fraction, the properties of the primordial binary population, the lack of a dedicated treatment to deal with compact binaries, and the use of outdated prescriptions for the evolution of massive stars (mZAMS > 50 M ).
In order to test the new features of the Nbody6++GPU code, we have carried out an extensive comparison of the evolution of star clusters with 110,000 stars in N -body and Monte Carlo simulations in our companion paper (Kamlah et al. 2022a), where we have shown, among other things, that N -body models of the same clusters seem to evolve toward sparser configurations compared to Monte Carlo models with large tidal radii simulated with the MOCCA code.This difference is likely due to the different criteria used to identify escapers in the two methods, which can lead to an early removal of escaping stars in MOCCA simulations compared to Nbody6++GPU .

Stellar and compact object binaries
Mass-segregation of the most massive stars enhances strong dynamical interactions, which can trigger the ejection of the tightest binaries, the ionisation of the loosest ones, and the formation and hardening of new binaries.In the Dragon-II clusters, the processes responsible for the formation and disruption of binaries counterbalance efficiently, determining a slow variation of the overall binary fraction.As shown in Figure 4, the binary fraction decreases by a small fraction, down to f b,f in ∼ 0.16 − 0.18 in models starting with f b = 0.2 and to f b,f in = 0.04−0.05 in models with f b = 0.05.Interestingly, this variation in the binary fraction is similar, within the simulation time, to results obtained for lower-N cluster simulations (see e.g.Heggie et al. 2006).The decrease of the binary fraction is mostly due to the disruption of the softest binaries in the cluster and, for a small fraction (< 5%), to hard binaries that are ejected in strong dynamical interactions.These binaries have typical semi-major axes broadly distributed in the 10 −2 − 5 × 10 2 AU.For the sake of comparison, Figure 5 shows the initial period-mass distribution and mass-ratio of the population of primordial binaries in our models.
Figure 6 shows the distribution of the ratio between the semi-major axis of ejected binaries and the hard-binary separation, both measured at the moment of the ejection, and the ejection velocity distribution for two different simulations.The plot makes clear that the vast majority of ejected binaries are hard and that this population is dominated mostly by binaries with a mass m bin < 2 M .The velocities of the ejected binaries generally remain in the range of 1 − 100 km s −1 , too small compared to the circular velocity of the Galaxy to permit the identification of these escapers as former cluster members.
The upper panel of Figure 7 shows the variation of the fraction of binaries normalised to the total number of stars in a given mass bin and at a given time.Initially, around 35−50% of all stars with a mass above 20 M are initially binary members, with the maximum percentage achieved for stars heavier than 100 M .However, the population of heavy objects is rapidly depleted (note that t/T rlx = 0.22 corresponds in this case to t = 18.8 Myr) owing mostly to stellar/binary evolu- tion, which causes a sharp drop in their number.The maximum stellar mass keeps decreasing over time, whilst a small population of binaries with components in the 5 − 100 M develops -clearly owing to the formation of binaries with one or two BHs.The mass distribution of objects in binary systems, shown in the lower panel of Figure 7, highlights that the number of binaries with at least one component heavier than 10 M is relatively small compared to the total number of objects in binaries.Assuming initially N = 120, 000 stars and f b = 0.2, we see that less than 1, 000 binaries contain a component with a mass m * > 10 M , most of them being former components of a primordial binary.
The progenitors of compact objects, which are the most massive stars and stellar binaries in the cluster, have already sunk into the cluster centre when compact objects form.Therefore, to dissect the properties of compact binaries in Dragon-II clusters, we focus on binaries forming within the cluster half-mass radius, calculated along the cluster evolution.
Figure 8 shows the number of binaries with a WD, NS, or BH as a function of time for all models.
The population of binaries containing at least one WD (dWDs), N dWD , depends on the half-mass radius and binary fraction.At fixed half-mass radius, the number of binaries with a WD significantly decreases at decreasing f b , because most of these binaries are of a primordial origin.In fact, at fixed N stars and RHM, the ratio between the number of dWDs is 4 − 5 times higher in models with f b = 0.2 compared to those with f b = 0.05, thus comparable to the ratio between the initial amount of primordial binaries in one case or the other.At fixed value of f b , instead, the smaller the half-mass radius, the smaller is the number of dWDs.In general, by the end of the simulations we find N dWD 200 − 700 dWDs per cluster.The amount of binaries with a WD monotonically increases over the simulated time, highlighting the competition between WD depletion via dynamical encoun- ters and the formation of new WDs, mostly via binary stellar evolution (see also Figure 4 in Kamlah et al. 2022a).
The evolution of the number of binaries with a NS (dNS) shows two clear peaks at 20 and ∼ 100 Myr.These peaks correspond to the formation of NSs from stars in the highend (the first) and low-end (the second) of the NS progenitor mass range.The drop after each of the peaks is due to NS natal kicks, which cause the ejection of a substantial fraction of NSs from the parent cluster.The width of the peaks is related to the time needed for the NS to leave the cluster, i.e. when their distance from the cluster centre exceeds twice the tidal radius.After the second peak, the number of binaries with a NS decreases in all simulations, regardless of the initial conditions.We find that the largest value of N dNS is reached in the case of RHM = 1.75 pc, f b = 0.2, and N = 600k.At fixed value of RHM and N we find that a larger initial binary fraction leads to a more numerous population of binaries with a NS, around 50% more for models with f b = 0.2.At fixed value of N and f b the number of binaries with a NS increases at increasing values of RHM because in denser clusters it is more likely that massive stellar binaries either are ejected or merge before stellar evolution becomes dominant.
The population of binaries with a BH (dBH), similarly to those with a NS, are characterised by two peaks of forma-  tion, one at around 10 Myr, driven by stellar evolution, and another at later times driven by dynamics.The number of binaries with a BH, N dBH , in the primary peak depends on the initial number of stars -the larger N0 the larger N bBH , whilst the number in the secondary peak depends on both the half-mass radius and binary fraction, although it is hard to discern the effects of different initial conditions in this case.the Dragon-II clusters,

Ejection of single and double compact objects
Over the simulated time, all clusters lose around 20-70 single BHs, depending on the cluster initial conditions, and 10-70 binaries containing either one or two compact objects.Figure 9 shows the mass distribution of ejected single BHs, which is characterised by two peaks, one at mBH ∼ 3 M and another at mBH ∼ 25 M , and a tail that extends up to mBH ∼ 102 M .The first peak is due to the natal kick of NSs and low-mass BHs, with masses in the range mBH = 2.5 − 6 M , and develops in the first 10-50 Myr, whilst the secondary peak is due to dynamical interactions 2 .The population of ejected binaries hardly depends on the cluster initial conditions.Therefore, for the sake of simplicity, we gather the ejected binaries from all simulations to have a statistically larger sample.In the following, we distinguish between binaries containing two compact objects, labelled as DCOB, and those containing one compact object and a star, labelled as SCOB. Figure 10 shows the component mass, semi-major axis, and eccentricity distribution of the ejected binaries in all the Dragon-II clusters.
Around 94% of the ejected binaries are primordial.
A clear difference between double and single compact object binaries arises from these Figures.In total, we find 229 ejected DCOBs of both dynamical (144) and primordial (85) origin.The DCOBs exhibit a similar mass distribution for the primary and the companion, characterised by a plateau in the m1,2 = 2 − 20 M and a clear peak at m1 ∼ 45 M for the primary and m2 ∼ 27 M for the companion.The resulting mass ratio distribution is quite peculiar, with a clear dominance of DCOB with a mass ratio q > 0.6, owing to the tendency of dynamical interactions to pair objects of comparable mass.The eccentricity distribution is dominated by a peak around 0, caused by a sub-population of primordial binaries that underwent the common envelope phase (64.7%), and a nearly flat distribution in the range e = 0.5 − 1.
Additionally, we find 375 ejected SCOBs, the vast majority of which coming from primordial binaries (353) with a small contribution from dynamically assembled systems ( 22).The mass distribution of the compact objects in SCOBs peaks at a value, mCO ∼ 2 − 4 M , in the range of NSs and small BHs, definitely smaller compared to the mass distribution of the stellar companion, which peaks at 10 M , but with a secondary peak at ∼ 0.3 − 0.5 M .The binary mass-ratio distribution of SCOBs clearly differs from DCOBs, showing a peak at q ∼ 0.2 and a decrease toward larger values.The compact object in the SCOBs is mostly a low-mass BH (200) -typically with a mass mBH < 10 M (173) -or a NS (173), and in only two cases a ONeWD (2).The stellar companion is a main-sequence star in the vast majority of the cases (353), followed by core He burning stars (20) (all with a primary weighing < 5 M ), and 2 naked He main-sequence (MS) star.Stellar companions in the MS phase are relatively massive: 18 of them have a mass mMS < 1 M , 245 have a mass in the range 1 < mMS/ M < 10, 74 in the range 10 < mMS M < 20, and just one with a mass mMS = 29 M .All stars in the CHeB phase have a mass in the mCHeB = 5 − 16 M range and are paired with an object lighter than mCO < 5 M , all of them come from primordial binaries.
Focusing on DCOBs, we find a few peculiar and interesting systems.Among all ejected BBHs only 5 merge within a Hubble time, because most BBHs were ejected when the density and velocity dispersion of the cluster had already dropped due to its expansion and mass loss.
In two cases, the ejected BBH contains an IMBH with mass either MIMBH = 120 M or 350 M .In five cases, instead, we find an ejected BBH with a merging time smaller than a Hubble time.Table 3 summarises the number of ejected single and binary BHs, and of BBHs and BH-IMBH binaries that merge within a Hubble time.

Black hole -main sequence star binaries
The sample of known BH-MS star systems has significantly grown over the last few years (e.g.Giesers et al. 2018Giesers et al. , 2019;;Saracino et al. 2022;El-Badry et al. 2022;Shenar et al. 2022;Mahy et al. 2022).Some of the BHs observed in a BH-MS binary appear to reside in star clusters both in the Milky Way (Giesers et al. 2018(Giesers et al. , 2019) ) and the Large Magellanic Cloud (Saracino et al. 2022;Shenar et al. 2022), whilst others appear to be in the Galactic disc (El-Badry et al. 2022).It is an open question whether these BH-MS systems come from primordial or dynamically assembled binaries.In the case of a dynamical origin it is also unknown whether the stellar companion captured the BH or its progenitor.
In these regards, the Dragon-II models offer us a novel way to look for BH-MS binaries in simulated clusters and identify possible properties of BH-MS binaries formed through different channels.Since the Dragon-II cluster database is relatively small and limited to a single metallicity, we cannot perform a comprehensive comparison between observed and simulated BH-MS binaries.Nonetheless, it is il- lustrative to qualitatively compare the properties of BH-MS binaries formed in Dragon-II models and the observed one.For example, Dragon-II models permit us to dissect the population of BH-MS binaries into those forming inside the cluster, some of which have a lifetime much shorter than the cluster life and are disrupted via interactions with other cluster members, or that have been ejected from the cluster.Figure 11 shows the component masses, period, and eccentricity of in-cluster and ejected BH-MS binaries.We assume that in-cluster binaries are those forming at any time over the simulated time, therefore the same binary or one or both components can appear multiple times in the plot.We see that in-cluster binaries are markedly different from ejected binaries.The latter can be divided in two sub-classes.The first sub-class exhibits a short period (P < 0.1 day) and an almost null eccentricity, e ∼ 0. Binaries in this sub-class are characterised by a BH with mass mBH < 10 M and a MS star with a mass in the 2 − 10 M range.They originate from binary evolution, and, in particular, underwent a common envelope phase that shrank the semi-major axis and damped the eccentricity of the binary.The ejection engine of these binaries is a SN explosion.The second sub-class, instead, comprises heavier BHs (mBH = 10 − 100 M ) and lighter MS stars (mMS < 1 M ), and is characterised by eccentricities in the range e = 0.2 − 1, indicating that these binaries come from dynamical interactions sufficiently strong to eject the binary from the cluster.
In-cluster BH-MS binaries can contain BHs and MS stars significantly heavier than the ejected binaries and are characterised by longer periods (P > 10 d) compared to ejected binaries.Most in-cluster binaries with a period P 10 3 d have zero eccentricity, whilst practically all those with a longer period have eccentricity > 0.1 and up to extreme values.
From Figures 11, it is evident that in-cluster binaries exhibit a peculiar distribution in the mBH − mMS, which suggests the existence of two sub-classes.We find that the first class is characterised by a companion with a mass mMS/mBH = k (mBH/1 M ) −1/2 , with k = 2 − 10.Most binaries falling in this class have a period shorter than 100 d, whilst the second class involves binaries with mBH > 10 M and mMS < 5 M .
An even more clear distinction is shown in Figure 12, where the MS-to-BH mass ratio is shown against the orbital period and eccentricity.This plot highlights four interesting peculiarities of in-cluster BH-MS binaries: • the vast majority of binaries with e < 0.1 are primordial.Most of them are characterised by mMS/mBH > 0.3, heavy MS stars mMS > 1 M , and periods below P < 100 d; • primordial binaries with e > 0.1 have larger periods (P = 10 2 − 10 6 d), and similar mass ratio and MS mass as circular primordial binaries; • the vast majority of dynamically formed binaries have e > 0.1 and periods in the range (P = 10 2 − 10 9 d).They are generally characterised by a mass ratio mMS/mBH < 0.3, MS stars with a mass mMS < 10 M and a BH with mass mBH = (10 − 100) M ; • only a handful dynamically formed binaries have e < 0.1, and are all characterised by a period P = 1 − 10 d.
As shown in Figure 12, we find that the longer is the orbital period the larger the binary eccentricity, and almost all binaries with eccentricity e > 0.1 have a period P > 100 d, with a handful exceptions.Most binaries with a period shorter than P < 100 d, instead, are primordial and involve a MS star heavier than mMS > 1 M .
The difference between primordial and dynamical BH-MS binaries is further highlighted in Figure 13, which shows the component masses of these two classes of binaries.From the plot, it is apparent that dynamically assembled binaries dominate the region of the plane with mBH > 10 M and mMS < 10 M .The observed BH-MS binaries have orbital properties quite different from our ejected binaries, especially if we consider the observed period and eccentricity.However, only the quiescent BH candidates in NGC3201 are still associated with a star cluster, whilst the origin of the other binaries is unknown.Two of the six observed binaries (Shenar et al. 2022;Mahy et al. 2022) have component masses compatible with our primordial binaries, one of them (El-Badry et al. 2022) falls in a range where only dynamically assembled binaries are present, and the three sources observed in the Galactic globular cluster NGC3201 have component masses compatible with both in-cluster and ejected binaries.
In our models, the vast majority of ejected binaries have a primordial origin and their small period (P < 0.01 d) owes to mass transfer episodes.The few ejected binaries formed dynamically are characterised by a period P < 1 d, still much shorter than observed values.Wider, and more numerous, ejected binaries could form in substantially looser or lighter star clusters.On the one hand, decreasing the cluster mass or density would enlarge the hard-binary separation and possibly increase the semi-major axis of ejected binaries (Morscher et al. 2015).On the other hand, a smaller cluster mass would correspond to a lower escape velocity and thus it is more likely for binaries to escape the parent cluster.
In principle, MS-MS binaries ejected in the earliest phase of the cluster life could further contribute to the population of BH-MS binaries, but these binaries are removed from our simulations before they can further evolve.Nonetheless, we find that only two ejected MS-MS binaries have at least one component with mass above the threshold for BH formation, i.e. ∼ 18 M , thus ensuring that ejected MS-MS binaries do not contribute to the population of ejected BH-MS binaries.
Among all observed data, the binaries observed in NGC3201 are probably the ones more suited for a comparison with our models, given the metallicity and mass of NGC3201.
From the central and bottom panel of Figure 11, it is apparent that our in-cluster binaries have periods, eccentricities, and BH masses compatible with those observed in NGC3201.The fact that our models do not match well the companion mass may be due to NGC3201's age.In fact, this cluster is relatively old (∼ 11.5 ± 0.4 Gyr VandenBerg et al. 2013), thus its population of binaries has likely been heavily processed over time, and most of its stellar population with super-solar mass already left the MS. Figure 13 favours this interpretation.Note that both the mass of BHs and MS stars in dynamically formed BH-MS binaries tend to be smaller compared to primordial binaries.As the BH-burning process proceeds, the average BH mass will keep decreasing, while stellar evolution processes will deplete the high-end tail of the MS mass distribution, possibly favouring the formation of BH-MS binaries in the region populated by NGC3201 sources.

Black hole subsystem
In all Dragon-II clusters, the segregation time is generally shorter than the stellar evolution timescale of massive stars, therefore massive stars sink to the cluster centre before evolving to BHs.This implies a possible enhancement of the probability for stellar mergers and star-BH collisions.
Given the short segregation times, BHs dominate the dynamics in the cluster core already after a time t = 20 − 40 Myr, making up the 50 − 80% of the mass in the cluster core and around 10% of the mass within the half-mass radius, as shown in Figure 14.Given the amount of mass in BHs enclosed within the core radius, this length scale can be regarded as the BH sub-system scale radius (see e.g.Arca Sedda et al. 2018).
A similar trend of the BH mass fraction inside RHM has been found also in recent simulations performed with the Monte Carlo code MOCCA (Giersz et al. 2019; Wang et al. 2022) and the N -body code PeTar (Wang et al. 2022), which both exploit similar stellar evolution recipes.
Both the primordial binary evolution and the onset of three-body and multiple gravitational scattering favour the formation of binaries containing at least one BH.Figure 15 shows the BH formation efficiency, defined as the ratio between the number of BHs inside the cluster core radius and the initial cluster mass, i.e.BH,BBH = NBH,BBH(< Rc)/M cl,0 .We find that, regardless of the initial cluster mass, half-mass radius, or binary fraction, all models are characterised by BH (0.8 − 2) × 10 −3 M −1 for single BHs and BBH (0.8−2)×10 −4 M −1 for binary BHs.As shown in the right panel of Figure 15, the BH formation efficiency slightly increases with the simulation time, although it is unclear whether this quantity saturates already at tsim/T rlx 10.Note that our definition of BBH implies that a cluster with initial mass 7 × 10 4 (6 × 10 5 ) M contains around 7(60) BHs in a binary system after 10 relaxation times.
It might seem trivial that is independent of the cluster initial conditions, as it suggests that it is just a consequence of the adopted mass function.However, the BH-burning mechanism (e.g.Banerjee et al. 2010;Downing et al. 2011Downing et al. , 2010;;Breen & Heggie 2013;Arca Sedda et al. 2018;Kremer et al. 2020b), by which the most massive BHs pair in binaries that first eject the lighter BHs from the cluster and then get themselves ejected via super-elastic binary-single and binarybinary scatterings, could significantly affect the population of BHs.This does not seem the case in the Dragon-II models.The small spread observed in the BH binary formation efficiency is related to the initial cluster half-mass radius and binary fraction, whilst the weak increase of BBH over time is the result of dynamically formed binaries.
Figures 16-18 show the cluster and BH subsystem density profiles at different times for three cluster models with Looking at the different panels it is possible to identify the signatures of the whole BH burning process as described in Breen & Heggie (2013).Firstly, BHs start forming and interacting, driving the formation of the BH subsystem and its subsequent expansion over a timescale t ∼ T rlx .Secondly, dynamical BH interactions cause the steepening of the BHS density and the contraction of its structure, driven by BH ejections over a time 1 < t/T rlx < 5. Thirdly, the BH subsystem rebounces and expands again, reaching a seemingly stable structure, at least within the simulated time.
Figure 19 shows the BH mass distribution at different times for a model with N = 1.2 × 10 5 stars, RHM = 1.75 pc, and f b = 0.2.This plot shows all BHs inside the cluster at a given time, regardless whether they are components of a binary system or single BHs.For the sake of comparison, we added in the plots the overall BH mass distribution inferred by the LVC (The LIGO Scientific Collaboration et al. 2021b).The plot highlights an initial phase in which the first BHs start to form, some of them falling in the upper-mass gap, but as the evolution proceeds new, lighter, BHs form while the most massive BHs are ejected via binary-single and binary-binary scatterings, as expected in the BH-burning scenario.
Interestingly, our simulations suggest that the evolution of the cluster can naturally lead to the peak around 10 M inferred from GW detections, mostly owing to stellar dynamics that crucially sculpts the BH population.Nonetheless, any comparison among our data, which show all BHs in the cluster, and LVC observations, which are representative of BH mergers, must be taken with a grain of salt.There are other potential explanations for the 10 M peak, like isolated binary stellar evolution (e.g.van Son et al. 2022), impact of primordial binary evolution in star clusters (e.g.Arca Sedda et al. 2021a), metal rich star clusters (e.g.Arca Sedda et al. 2020;Rastello et al. 2021;Chattopadhyay et al. 2022).Hopefully, the new data acquired during the forthcoming four LVC observation run could help pinning down the impact of different processes on the BH mass distribution.
We find that almost all BHs heavier than > 30 M are ejected from the simulated clusters reaching more than ∼ 15 relaxation times.
To further highlight the BH burning process, we reconstruct the time evolution of the average BH mass, mBH , for all BHs enclosed within the half-mass radius.As shown in Figure 20, mBH follows the same trend regardless of the cluster initial condition, namely: i) the most massive BHs form first and the average mass sets close to the peak allowed by the adopted stellar evolution model (35 − 40 M ); ii) more numerous, lighter BHs start to form causing a rapid decrease of the average mass down to 15 − 20 M ; iii) dynamical processes kick in and trigger BH ejection, leading to a secular decrease of the BH average mass down to ∼ 8 − 10 M (see also Giersz et al. 2019).
The similar mBH time evolution observed in different models supports the idea that the BH burning process is substantially due to dynamics.This is further highlighted in Figure 21, which shows the BH average mass as a function of the time normalised to the cluster relaxation time.We find that at a time t > T rlx the average BH mass is well described t/T rlx = 6.9 10 −1 10 0 10 1 t/T rlx = 9.7 Figure 18.Same as in Figure 16, but for one simulation with N = 3 × 10 5 and R HM = 0.47 pc.
by a simple relation: where m BH,rlx = 17.4 ± 0.1 M .Although our models are not meant to be representative of any observed cluster, and although there are certainly many pathways leading to the same final cluster evolutionary stage, our results suggest that old Galactic globular clusters and massive clusters in the Small Magellanic Cloud could be harboring a population of relatively light BHs (see Figure 2).This would explain why observations of BHs in binary systems are generally characterised by masses mBH < 20 M , relatively lighter than the typical value inferred for the population of merging BHs, i.e. mBH,GW 30 M .

Using scaling relations as input for semi-analytic codes
It is well known that N -body simulations of star clusters require generous computational resources to enable an exploration of the phase space and to reach an appreciably long simulated time.The Dragon-II simulations make no exceptions, as they required in total approximately 2.2 million core hours.To overcome this problem, many works have proposed semi-analytic tools specifically devoted to study the evolution of compact objects in the last few years, and especially BH binary mergers (e.g.Arca Sedda et al. 2020;Mapelli et al. 2021;Fragione & Kocsis 2018;Antonini & Gieles 2020b;Antonini et al. 2019;Arca Sedda et al. 2021a;Kritos et al. 2022;Mapelli et al. 2022).
One ingredient missing in some of these fast and accurate codes is a treatment of the co-evolution of the star cluster and the BH population, which may significantly affect the The Dragon-II models could provide important fitting formulas to implement the evolution of under-filling cluster models in such semi-analytic tools.
The overall evolution of Dragon-II star clusters can be described by simple expressions (Equations 6 and 7).If the cluster initial mass and half-mass radius are known, the aforementioned relations enable an accurate description of its evolution, at least in the case of under-filling star cluster models.Moreover, our models offer also insights on the internal evolution of the cluster, providing, for example, details about the mass distribution of ejected single and double compact objects, and the properties of the central black hole subsystem.These ingredients can be easily implemented in semi-analytic tools to obtain a fast and accurate description of compact object dynamics in clusters too massive to be modelled with N -body models.
A simple implementation of the cluster evolution has been already developed by Arca Sedda et al. (2021a) in their B-POP code, showing that the inclusion of cluster mass loss and expansion causes a critical decrease of the probability of high-generation mergers in dense and massive star clusters (see also Antonini & Gieles 2020a).

SUMMARY AND CONCLUSIONS
In this work, we have presented the first results from the Dragon-II star cluster simulations: a suite of 19 direct N -body simulations, performed with the Nbody6++GPU code, modelling the evolution of star clusters with up to 1 million stars and up to 33% of stars initially in a binary, over a timescale of ∼ 0.5 − 2 Gyr.These simulations contain up-to-date stellar evolution models, and for the first time a series of recipes to treat relativistic binaries in terms of merger remnant mass, spin, and post-merger recoil.Our models represent clusters initially under-filling their Roche lobe, and therefore their evolution can be considered quasi-isolated.The Dragon-II models considerably expand the portion of parameter space covered with full N -body simulations, opening the possibility to compare with large-N Monte Carlo models.Clearly, there is a vast number of parameters whose impact on the simulation results remains unclear.For example, adopting a sufficiently large value of the metallicity would imply the impossibility to form IMBHs from stellar collapse.However, we expect that our main conclusions about the properties of the BH population should not be severely affected by cluster metallicity, as they appear to be driven mostly by dynamics.
We find that the amount of primordial binaries seems to poorly affect the overall evolution of the cluster and the evolution of the BH population; however the adopted initial orbital properties could become important when comparing our data with observations, like in the case of BH-MS binaries.For example, a different assumption on the initial mass-ratio distribution could lead to primordial binaries with final BH-MS component masses more similar to the observed one.However, discrepancies among observations and models could arise from a combination of different assumptions, making it hard to pinpoint the main source of uncertainty.
Finally, our simulations model initially underfilling clusters, meaning that the impact of the Galactic field is almost negligible compared to clusters' internal dynamics.This choice enabled us to have a clean view at the impact of stellar interactions on the evolution of the whole cluster and its BH population, and incidentally lead to star cluster models that resemble observed clusters in term of mass and radius.Future simulations adopting filling or overfilling clusters may help understanding whether the evolution of BH subsystems is intrinsically linked to the overall evolution of the host cluster, for example in terms of mass-loss and expansion.The main outcomes of the Dragon-II models can be summarised as follows.
• mass-loss and expansion of Dragon-II clusters is mostly determined by internal dynamics and can be described by simple analytical expressions, with parameters that weakly depend on the initial conditions.The binary fraction varies mildly over the simulated time, within 10 − 15% of its initial values.Nonetheless, stellar evolution and dynamics cause a progressive drop of the fraction of stars in binary systems for primary masses m1 > 2 M [Figures 2-7]; • over a Gyr timescale, Dragon-II clusters contains around 200-700 binaries with at least one WD, whilst the number of binaries with a NS or a BH generally remains below 1-10 and 5-40, respectively.In general, binaries with at least one compact object are more numerous in clusters with a larger initial binary fraction, suggesting that most of these binaries have a primordial origin.Moreover, the denser the cluster is the smaller the number of binaries, owing to energetic dynamical interactions that disrupt binaries more efficiently [Figure 8]; • ejected binaries with one (SCOB) or two (DCOB) compact objects have different properties.DCOBs exhibit masses following a nearly flat distribution around 2 − 20 M and a peak at mBH = 45 M , a peculiar mass-ratio distribution that peaks around q 0.6, and a flat eccentricity distribution in the range e = 0.5 − 1. SCOBs, most of which formed from primordial binaries, typically involve low-mass BHs (mBH = 3 − 10 M ) and fairly massive MS stars (mST = 1 − 10 M ) [Figure 10]; • we find a substantial population of BH-MS binaries in Dragon-II models.Most BH-MS binaries forming inside the cluster have typical BH masses mBH > 10 M , a companion star with mass mMS = 0.7 − 100 M , orbital periods > 10 days, and span the entire eccentricity range.Ejected BH-MS binaries, instead, feature significantly smaller BH masses mBH < 10 M , shorter periods (< 10 days), and are mostly contributed by primordial binaries.We find that the properties of the modelled binaries are compatible with some features of observed BH-MS binaries, especially those observed in the globular cluster NGC3201 ; • in all Dragon-II models, BHs form a long-lived subsystem in the cluster centre already after 0.5 relaxation times, with a typical density 10 − 100 times higher than that of stars.The cluster core radius represents a good proxy of the BH subsystem size, as BHs make up 50 − 80% of the mass enclosed within this radius.We find that the ratio between the number of BHs inside the core radius and the bound cluster mass, which we refer to as formation efficiency, attains values of BH,BBH/ M = 10 −3 (10 −4 ), for single and binary BHs, respectively.This quantity is only mildly dependent on the initial conditions, suggesting that dynamical processes have a relatively minor effect on the overall BH population over the simulation time ; • dynamics in the BH subsystem critically affects the BH mass spectrum, owing to the BH-burning process.The peak of the mass distribution generally shifts from initial values m BH,pk = 25 M down to m BH,pk = 5 − 15, and the average mass steadily decreases after one relaxation time, following an identical evolution regardless of cluster properties .
Our simulations suggest that dynamically old star clusters harbour in their centre a population of BHs whose amount scales linearly with the cluster bound mass.The older the cluster is, the smaller the peak of the BH mass spectrum and the average BH mass.

Figure 1 .
Figure 1.Initial density, calculated at the half-mass radius, as a function of number of stars for several grids of direct N -body (blue points) and Monte Carlo (red boxes) simulations.The Dragon-II cluster database is represented by the green star.

Figure 2 .
Figure 2. Time evolution of the mass (top panel) and half-mass radius (bottom) of Dragon-II clusters (black lines) compared to observed massive clusters in the Milky Way (MW, blue stars), the two Magellanic Clouds (green squares and red points), the Andromeda galaxy (M31, green squares), the Henize 2-10 starburst galaxy (He2-10, grey diamonds), and the M83 galaxy (light green triangles).

Figure 3 .
Figure 3. Half-mass radius calculated at the end of each simulation, normalised to the initial value, as a function of the total simulated time normalised to the cluster initial relaxation time.The dashed line is a power-law fit of the data, while each point represents one simulation, with the colour map representing the initial cluster mass.

Figure 4 .
Figure 4. Binary fraction as a function of time normalised to the relaxation time, taking into account all binaries.

Figure 5 .
Figure 5.Initial values of total mass (x-axis) and period (y-axis) for all primordial binaries in one of the models with R HM = 1.75 pc, N = 120k, and f b = 0.2.The colour map marks the mean value of the mass ratio inside each pixel.

Figure 6 .
Figure6.Upper panels: Distribution of the semi-major axis of ejected binaries normalised to the hard-binary separation measured at the moment of ejection for binaries with a mass m bin < 2 M (blue steps) or heavier (orange steps).Lower panels: ejection velocities for the same classes of binaries.The plots refer to two simulations with (R HM , N ) = (1.78pc, 120k) (left-hand panels) and (0.8 pc, 300k) (right-hand panels).

Figure 7 .
Figure 7. Top: fraction of objects in binaries for different mass bins and different times.Bottom: number of binaries with mass in a given bin at different times.We show a simulation with f b = 0.2, N = 120k stars, and R HM = 1.75 pc.

Figure 8 .
Figure 8. Number of binaries with at least one WD (upper panel), a NS (central panel), or a BH (lower panel) as a function of time for all Dragon-II clusters.Different colours correspond to different values of the initial half-mass radius.There are two lines for each colour, corresponding to two realisations of the same cluster.

Figure 10 .
Figure10.Upper panel: mass (left) and mass ratio (right) distribution of ejected binaries containing either two compact objects, i.e.BHs, NSs, or WDs.The mass distribution refers to the primary (black straight steps) and companion (black dashed steps) in the case of double compact object binaries, and to the compact object (red straight steps) and the star (red dashed steps) in the case of binaries with one compact object.Lower panels: semi-major axis (left) and eccentricity (right) distribution of ejected binaries containing two (black straight steps) or one (red dashed steps) compact objects. name

Figure 12 .
Figure 12.Orbital period as a function of the ratio between the MS and BH masses for all BH-MS binaries formed inside the clusters.The colour coding marks the binary eccentricity.

Figure 13 .
Figure 13.Masses of the MS star and BH in primordial (grey squares) and dynamical (blue triangles) BH-MS binaries.The red and green stars represent observed binaries as in Figure 11.The plot refers to both in-cluster and ejected binaries.

Figure 14 .
Figure14.Fraction of mass in BHs within the cluster half-mass (straight lines) and core radius (dashed lines).Different colours correspond to different simulations.

Figure 15 .
Figure 15.Left-hand panel: number of BHs (straight lines) and BBHs (dashed lines) inside the core radius normalised to the initial cluster mass (BH formation efficiency) as a function of time normalised to the initial cluster relaxation time.The large jumps follow the jumps observed in the core calculation.The colour-map identifies the cluster mass.Right-hand panel: BH formation efficiency for single (points) and binary BHs (diamonds) calculated at the end of each simulation as a function of the ratio between the simulated time and the cluster relaxation time.
Figure 19.Mass distribution of BHs in a cluster model with N = 120k stars, half-mass radius R HM = 1.75 pc, and initial binary fraction f b = 0.2.The black straight line and shaded area correspond to the mass distribution of primary BHs in merging binaries inferred from observed BBH mergers during the first, second, and third observing runs of the LIGO-Virgo-Kagra collaboration The LIGO Scientific Collaboration et al. (2021b).

Figure 20 .
Figure20.Average BH mass for BHs within the cluster half-mass radius.Note that the spikes can be due to heavy BHs orbiting on the edge of the half-mass radius, or on eccentric trajectories that occasionally enter the half-mass radius.

Figure 21 .
Figure 21.Average BH mass inside the half-mass radius as a function of the time normalised to the cluster initial relaxation time.

Table 1 .
Col. 1-4: initial number of stars, cluster mass, half-mass radius, and primordial binary fraction.Col. 5: number of independent realisations.Col. 6-7: initial relaxation and segregation time.Col. 8: simulated time.Col. 9-10: number of 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.

Table 4 .
Orbital properties of the observed BH-MS binaries.For the sources labelled with an * only lower limits to the BH mass are available.