Stellar-mass black holes in young massive and open stellar clusters – IV. Updated stellar-evolutionary and black hole spin models and comparisons with the LIGO-Virgo O1/O2 merger-event data

ABSTRACT

I present a set of long-term, direct, relativistic many-body computations of model dense stellar clusters with up-to-date stellar-evolutionary, supernova (SN), and remnant natal-kick models, including pair instability and pulsation pair instability supernova (PSN and PPSN), using an updated version of |${\rm{\small NBODY7}}$|  N-body simulation program. The N-body model also includes stellar evolution-based natal spins of black holes (BHs) and treatments of binary black hole (BBH) mergers based on numerical relativity. These, for the first time in a direct N-body simulation, allow for second-generation BBH mergers. The set of 65 evolutionary models have initial masses |$10^4{\!-\!}10^5\, \mathrm{M}_{\odot }$|⁠, sizes 1–3 pc, metallicity 0.0001–0.02, with the massive stars in primordial binaries and they represent young massive clusters (YMC) and moderately massive open clusters (OC). Such models produce dynamically paired BBH mergers that agree well with the observed masses, mass ratios, effective spin parameters, and final spins of the LVC O1/O2 merger events, provided BHs are born with low or no spin but spin-up after undergoing a BBH merger or matter accretion on to it. In particular, the distinctly higher mass, effective spin parameter, and final spin of GW170729 merger event is naturally reproduced, as also the mass asymmetry of the O3 event GW190412. The computed models produce intermediate-mass, |$\sim 100\, \mathrm{M}_{\odot }$| BBH mergers with primary mass within the ‘PSN gap’ and also yield mergers involving remnants in the ‘mass gap’. They also suggest that YMCs and OCs produce persistent, Local-Universe GW sources detectable by LISA. Such clusters are also capable of producing eccentric LIGO-Virgo mergers.

gravitational waves, methods: numerical, stars: black holes, stars: kinematics and dynamics, globular clusters: general, open clusters and associations: general

1 INTRODUCTION

The LIGO-Virgo collaboration (hereafter LVC) has so far published 10 binary black hole (hereafter BBH) and one binary neutron star (hereafter BNS) merger events, in their first gravitational-wave transient catalogue (Abbott et al. 2019a, GWTC-1), through ground-based, interferometric detection of gravitational waves (Abbott et al. 2016a, hereafter GW) during their first (O1) and second (O2) observing runs. In their recently concluded third observing run (O3; https://gracedb.ligo.org/superevents/public/O3/), 56 additional candidates of compact-binary mergers are detected, a few of those being of BNS and even neutron star-black hole (hereafter NSBH) mergers. Additionally, a few runs are labelled as ‘mass gap’ in the sense that one or both of the merging members lie in the potential gap between the masses of neutron stars (hereafter NS) and stellar-remnant black holes (hereafter BHs) that certain supernova (hereafter SN) models (e.g. Fryer et al. 2012) predict. The handful of events from O1/O2 already suggest that detection of GW transients from compact binary merger events not only is interesting by its own right (Abbott et al. 2016a, b) but also has the potential to provide unprecedented information regarding masses, spins, and their boundaries (Abbott et al. 2019b), of stellar-remnant BHs and NSs. Such information would provide strongest constraints on the formation mechanisms of compact stellar remnants and of the environment in which their parent stars form and evolve (Belczynski et al. 2020; Olejak et al. 2020).

It is as well of wide interest and diverse implications (e.g. Abadie et al. 2010; Mandel & Farmer 2017) to consider how and under which conditions NSs and BHs would pair up in tight-enough binaries so that they can spiral in by emitting GW radiation and merge within the Hubble time. Recent numerical studies based on analytical (Hénon 1975), direct N-body integration (Aarseth 2003), and Monte Carlo approach (Hénon 1971; Joshi, Rasio & Portegies Zwart 2000; Hypki & Giersz 2013) show that the retention of BHs in dense stellar clusters of wide mass range, beginning from low-/medium-mass young and open clusters (e.g. Banerjee, Baumgardt & Kroupa 2010; Ziosi et al. 2014; Mapelli 2016; Banerjee 2017, 2018a, b; Park et al. 2017; Di Carlo et al. 2019; Kumamoto, Fujii & Tanikawa 2019; Rastello et al. 2019) through globular clusters (e.g. Breen & Heggie 2013; Morscher et al. 2013; Sippel & Hurley 2013; Arca-Sedda 2016; Hurley et al. 2016; Rodriguez, Chatterjee & Rasio 2016; Wang et al. 2016; Askar et al. 2017; Chatterjee, Rodriguez & Rasio 2017a; Chatterjee et al. 2017b; Fragione & Kocsis 2018; Rodriguez et al. 2018; Antonini & Gieles 2020; Kremer et al. 2020) to galactic nuclear clusters (e.g. Antonini & Rasio 2016; Hoang et al. 2018, 2019; Antonini, Gieles & Gualandris 2019; Arca-Sedda & Capuzzo-Dolcetta 2019; Arca Sedda 2020), comprise environments where BHs can pair up through close dynamical interactions, which, furthermore, lead to general-relativistic (hereafter GR) coalescences of these BBHs. Being much more massive than the rest of the stars, the BHs, which remain gravitationally bound to a cluster after their birth, spatially segregate and remain highly concentrated in the cluster’s innermost (and densest) region (e.g. Banerjee et al. 2010; Morscher et al. 2015) due to dynamical friction (Chandrasekhar 1943; Spitzer 1987) from the stellar background. This is essentially an early core collapse of the cluster leading to its post-core-collapse behaviour (Hénon 1975; Spitzer 1987; Heggie & Hut 2003), i.e. energy generation in the ‘collapsed’ BH core leading to an overall expansion of the cluster with time (Breen & Heggie 2013; Antonini & Gieles 2020). Inside this core, BHs undergo close binary–single and binary–binary encounters giving rise to compact subsystems (triples, quadruples, or even higher order multiples) whose resonant evolution can lead to GR inspiral and merger of their innermost binaries (Leigh & Geller 2013; Samsing, MacLeod & Ramirez-Ruiz 2014; Geller & Leigh 2015; Banerjee 2018b; Samsing 2018; Zevin et al. 2019), through the binaries’ eccentricity pumping. The breakup of such subsystems or simply close, flyby encounters may also lead to a sufficient boost in eccentricity of a BBH such that it merges either promptly, in between two close encounters (e.g. Kremer et al. 2019) or within a Hubble time if it gets ejected from the cluster as a result of the interaction (e.g. Rodriguez et al. 2015; Park et al. 2017; Kumamoto et al. 2019). Note that such GR mergers can also happen in hierarchical systems, containing NSs and BHs, in a galactic field that derive from field massive-stellar multiplets (e.g. Toonen, Hamers & Portegies Zwart 2016; Antonini, Toonen & Hamers 2017; Fragione & Loeb 2019; Fragione, Loeb & Rasio 2020b).

Alternatively, compact-binary mergers can happen in a galactic field either via common-envelope (hereafter CE) evolution of massive-stellar binaries with hydrogen-rich envelope (Dominik et al. 2012; Belczynski et al. 2016a; Stevenson et al. 2017; Giacobbo, Mapelli & Spera 2018; Baibhav et al. 2019) or in close, tidally interacting,‘over-contact’ binaries composed of chemically homogeneous members (De Mink et al. 2009; De Mink & Mandel 2016; Marchant et al. 2016). In such studies, BBH inspiral detection rate of ∼10 to ∼1000 yr⁻¹ has been estimated for the full-sensitivity LIGO, i.e. from similar to up to two orders of magnitudes higher detection rate than what is estimated for dynamically formed BBHs (Banerjee et al. 2010; Rodriguez et al. 2016; Askar et al. 2017; Banerjee 2017).

Recently, studies by Banerjee (2017, hereafter Paper I), Banerjee (2018a, Paper II), and Banerjee (2018b, Paper III) have shown the importance of triple-/higher order-dynamical interactions, including post-Newtonian (hereafter PN; Blanchet 2014) terms, in triggering in-cluster¹ GR BBH coalescences, particularly inside relatively low-velocity-dispersion systems such as open clusters and lower mass globular clusters (hereafter GC). These studies have utilized the direct (star-by-star) N-body integration code nbody7 (Aarseth 1999; Aarseth 2012; Nitadori & Aarseth 2012) that couples the fourth-order Hermite orbit integration technique with advanced subsystem identification and their regularization (Aarseth 2003), PN treatment for binaries, triples, and higher order subsystems containing NS or BH by adopting the |$\rm{\small ARCHAIN}$| algorithm (Mikkola & Tanikawa 1999; Mikkola & Merritt 2008), and semi-analytical, recipe-based population synthesis of single and binary stars by adopting the |$\rm{\small BSE}$| program (Hurley, Pols & Tout 2000; Hurley, Tout & Pols 2002). Note that nbody7 uses |$\rm{\small ARCHAIN}$| to treat triples and higher order subsystems even if they do not contain NS or BH. For close star–star binaries or those containing a white dwarf (hereafter WD) or close hyperbolic passages, KS regularization (Kustaanheimo & Stiefel 1965) is applied along with energy loss due to tidal interaction and gravitational radiation (PN-2.5 term; Peters 1964).

In the nbody7 computations presented in Paper I, II, and III, a somewhat old SN remnant mass scheme, as in Belczynski et al. (2008), is applied along with stellar wind mass-loss recipes of Belczynski et al. (2010), as implemented in the currently official, public version of nbody7 (in its |${\rm {\small BSE}}$| sector). In particular, the remnant scheme did not take into account the BH mass ceiling and the ‘upper mass gap’ due to pulsation pair-instability supernova and pair-instability supernova respectively (hereafter PPSN/PSN; Langer et al. 2007; Woosley 2017; Mapelli et al. 2020). Also, the possibility of having a ‘lower mass gap’ between NS and BH masses due to ‘rapid’ core-collapse SN (Fryer et al. 2012) was not incorporated. The schemes adopted for assigning natal kicks of NSs and BHs, that are crucial for their retention inside clusters, were rather basic. As GW observations continue to provide mass measurements of an increasing number BHs and NSs (Abbott et al. 2019a), signatures of such theoretically predicted features are being widely investigated and discussed (e.g. Fishbach & Holz 2017; Fishbach, Holz & Farr 2017; Giacobbo et al. 2018; Abbott et al. 2019a; Chatziioannou et al. 2019; Rodriguez et al. 2019; Spera et al. 2019; Fishbach, Farr & Holz 2020; Kimball, Berry & Kalogera 2020; Olejak et al. 2020). Full-fledged, long-term, direct, relativistic N-body (or many-body) computations of dense stellar clusters with up-to-date stellar-evolutionary and remnant-formation and retention recipes is therefore of urgent need and highly desired.

In this work, such a set of 65 computations is presented for the first time using an updated version of nbody7 that explicitly incorporates up-to-date stellar remnant formation models. It additionally includes schemes for assigning spins of stellar-remnant (first-generation²) BHs based on detailed stellar-evolutionary models, run-time tracking of the BHs’ spins and assigning final spins and GW recoil velocities to in-cluster BBH mergers based on numerical-relativity (hereafter NR) results. Note that one or more of these aspects have been incorporated in recent Monte Carlo (that use cmc or mocca codes; Joshi et al. 2000; Fregeau & Rasio 2007; Hypki & Giersz 2013; Giersz et al. 2013; Morscher et al. 2015) and direct N-body (that use nbody+ +gpu code; Spurzem et al. 2008; Wang et al. 2015) studies (Morawski et al. 2018; Di Carlo et al. 2019; Rodriguez et al. 2019; Kremer et al. 2020). The new BH-spin aspects of the updated nbody7 now automatically allows the formation and tracking of second-generation BHs and BBH mergers involving them as in latest cmc-based Monte Carlo studies (Rodriguez et al. 2018, 2019). This work introduces the new nbody7 computations, focusing on the compact-remnant merger outcomes from them and comparing them with LVC O1/O2 merger-event data.

In Section 2, elements of the updated nbody7 are described. Section 3 describes the cluster model characteristics and the set of direct N-body computations. Section 4 describes the compact-remnant merger outcomes from these new N-body simulations and compares them with the O1/O2 events. Section 5 summarizes the results and indicates future prospects.

2 THE UPDATED nbody7

Both the standalone bse code (Hurley et al. 2000, 2002) and its version that is integrated with nbody7 have recently been updated in parallel with up-to-date stellar wind mass-loss and remnant-formation recipes. The updated versions are thoroughly tested to agree well with the remnant outcomes of |$\rm{\small STARTRACK}$| population-synthesis code (Belczynski et al. 2008; Belczynski et al. 2016a). Furthermore, material fallback during a core-collapse SN is considered explicitly in determining natal kicks of NSs and BHs. The reader is advised to consult Banerjee et al. (2020, hereafter Ba20) where all these new implementations and |${\rm {\small BSE}}$|-|$\rm{\small STARTRACK}$| comparisons are elaborated. Below, these updates to nbody7 are summarized and further new ingredients are described. At present, the updated nbody7 is private (which will be publicized at a later occasion) but the parallel standalone bse is available publicly with Ba20.

2.1 Stellar wind

The stellar wind mass-loss follows the recipe of Belczynski et al. (2010, hereafter B10) as described in section 2.1 of Ba20. The main difference with respect to the currently public version of standalone |$\rm{\small BSE}$| is the application of Vink, de Koter & Lamers (2001) wind model, that exhibits a line-driven bi-stability jump feature at |$\approx 25\, 000$| K surface temperature, for massive, hot stars (of surface temperature |$\gt 12\, 500$| K). Other new ingredients are the use of the luminous-blue-variable (LBV) wind of Humphreys & Davidson (1994) and, for Helium stars, the Wolf–Rayet (WR) wind of Hamann & Koesterke (1998), Vink & de Koter (2005). For lower mass, colder stars, the original Hurley et al. (2000) wind model is maintained. The wind mass-loss depends on metallicity, Z, of the star through the use of the bi-stability jump, WR, and Hurley et al. (2000) wind models. Note that all these wind ingredients are in principle implemented in the current public nbody7. However, the implementation is revised so that the sequencing of the various wind elements are parallel to that of B10, as elaborated in Ba20.

2.2 Stellar remnant formation

The remnants (NSs and BHs) are formed according to the rapid- or delayed-SN mass and fallback models of Fryer et al. (2012, hereafter F12) which are available as options. The options for the existing older remnant models are also retained. Additionally, there is now the option of having a BH mass ceiling due to PPSN and the upper mass gap due to PSN, according to the prescriptions of Belczynski et al. (2016b, hereafter B16). These new implementations and the corresponding initial mass-final mass relations are detailed in Ba20. Optionally, PPSN-derived BHs can also be formed with masses according to the ‘moderate’ and ‘weak’ PPSN models of Leung, Nomoto & Blinnikov (2019), following their implementations in Belczynski et al. (2020, hereafter B20). For obtaining the gravitational mass of the remnant from its baryonic mass, a neutrino mass-loss of 10 per cent is assumed for BH formation and the mass-loss is according to Lattimer & Yahil (1989), Timmes, Woosley & Weaver (1996) for NS formation. NS formation through electron capture supernova (ECS; Podsiadlowski et al. 2004) is incorporated, as default in bse (which scheme is analogous to the scheme in Belczynski et al. 2008, producing the characteristic |$m_{\rm ECS,NS}=1.26\, \mathrm{M}_{\odot }$| ECS-NSs). In this work, F12-rapid and F12-delayed remnant mass models, with B16-PPSN/PSN (indicated as, e.g. F12-rapid+B16-PPSN/PSN remnant model), are applied in most of the model computations, except for a few where the weak PPSN is applied.

Note that with B16-PPSN/PSN, the (baryonic) mass of the pre-collapse Helium star in a PPSN is taken to be of |$45\, \mathrm{M}_{\odot }$| (B16; Woosley 2017) which collapses directly to a |$40.5\, \mathrm{M}_{\odot }$| BH, taking into account the 10 per cent neutrino mass-loss. |$45\, \mathrm{M}_{\odot }$| is widely considered as the theoretical upper mass limit of BHs due to PPSN and such a mass limit is also supported by BH masses measured in LVC O1/O2 BBH mergers (Abbott et al. 2019a; Chatziioannou et al. 2019; Kimball et al. 2020). Hence, a |$\gt 45\, \mathrm{M}_{\odot }$| BH is generally considered as a BH in the PSN (upper) mass gap. However, if one applies the weak PPSN model instead of the B16 model (in B20, B16-PPSN model is referred to as ‘strong’ PPSN), PPSN-derived BHs can reach up to |$\approx 50\, \mathrm{M}_{\odot }$| (assuming 10 per cent neutrino mass-loss) for very low metallicities, i.e. BH mass can lie in the ‘classical’ PSN mass gap. This is demonstrated in Fig. 1 (left-hand panel) where remnant masses are plotted against zero-age main-sequence (hereafter ZAMS) masses directly from nbody7 outputs (metallicity Z = 0.0002 is assumed). The characteristic lower mass gap between NSs and BHs, between |${\approx }2.5$| and |$5.0\, \mathrm{M}_{\odot }$|⁠, of F12-rapid remnant scheme is also indicated. See Ba20 for further examples of initial–final relations for various remnant mass schemes and metallicities and their comparisons with |${\rm{\small STARTRACK}}$|⁠.

Figure 1.

Left-hand panel: Examples of ZAMS mass-remnant mass relations as obtained from updated nbody7 used in this work. The outcomes for F12-rapid+B16-PPSN/PSN (filled, black squares) and F12-rapid+weak-PPSN/PSN (empty, brown squares) remnant-mass schemes (Section 2.2) are shown for the metallicity Z = 0.0002. For F12-rapid+B16-PPSN/PSN remnant model, the comparison with the corresponding relation from |$\rm{\small STARTRACK}$| (solid, magenta line) is demonstrated. In these N-body models (initially of total cluster mass |$M_{cl}(0)=5.0\times 10^4\, \mathrm{M}_{\odot }$| and half-mass radius r_h(0) = 2.0 pc), all stars are initially single whose ZAMS masses range from 0.08 to |$150.0\, \mathrm{M}_{\odot }$| and which are distributed according to the standard IMF. Right-hand panel: ZAMS mass-remnant mass relations with updated nbody7 from a model cluster of the same M_cl(0) and r_h(0) (Z = 0.0001; F12-rapid+B16-PPSN/PSN; standard IMF over |$0.08{\!-\!}150.0\, \mathrm{M}_{\odot }$|⁠), where all stars with |$M_{\rm ZAMS}\ge 16\, \mathrm{M}_{\odot }$| are in primordial binaries (see Section 3; filled, black squares). Models involving such a massive primordial-binary population are indicated by ‘+MB’ in the legends and axis labels. For progenitor stars that have undergone a star–star merger before the remnant formation, the ZAMS mass of the primary (the more massive of the members participating in the star–star merger, at the time of the merger) is plotted along the abscissa. In all star–star mergers, f_mrg = 0.5 of the secondary’s mass is assumed to be lost in the merger process (see Section 2.6). The corresponding single-star nbody7 outcomes (empty, orange squares) and |$\rm{\small STARTRACK}$| outcomes (solid, magenta line) are shown for comparison. Also shown for comparison is the corresponding outcome (filled, red circles) when the initial population of massive binaries from the N-body model is evolved individually using the updated standalone bse of Ba20. Likewise the N-body model with primordial binaries, the primary mass is plotted along the abscissa if the outcome of the binary evolution is a single remnant (f_mrg = 0.5 is assumed).

Even with B16-PPSN/PSN, the remnant BH mass can enter the PSN mass gap due to star–star mergers at low metallicities producing stars with overmassive Hydrogen envelope, when, as observations suggest, the BH-progenitor stars are in tight massive primordial binaries instead of being single. This possibility is discussed in detail and demonstrated in Ba20; see also Spera et al. (2019). Fig. 1 (right-hand panel) demonstrates this from nbody7 outputs (filled, black squares) for a model star cluster with massive primordial binaries, where, with F12-rapid+B16-PPSN/PSN remnant model, BHs up to |$\approx 60\, \mathrm{M}_{\odot }$| form at Z = 0.0001. Also shown in Fig. 1 (right-hand panel) is the initial–final outcome when the same primordial-binary population from the model cluster is evolved as isolated binaries using the standalone |$\rm{\small BSE}$| of Ba20 that is updated in parallel (filled, red circles). Although the overall pattern of the initial–final points is similar for the isolated and the in-cluster binary evolutions, individual points do differ. This difference can be attributed to the stochastic, weak and strong perturbations that the binaries receive during their evolutions, altering their orbital parameters, when they are cluster members. However, differences can also arise due to the way |$\rm{\small BSE}$|’s binary evolution engine is integrated with the direct N-body evolution engine in nbody7. (The binary-evolution physics and its nbody7 integration remains unaltered in the updated standalone |${\rm {\small BSE}}$| and nbody7/bse.) A detailed investigation on the differences between dynamically influenced and isolated evolution of a massive binary population is beyond the present scope; see Di Carlo et al. (2019) in this context.

As detailed in Ba20 (see their section 2.3.1), for numerical convergence with the updated wind and remnant recipes, both standalone bse and nbody7 are run with the bse time-step parameters |${hboxpts1}=0.001$|⁠, |${\tt pts2}=0.01$| (along with |${\tt pts3}=0.02$| for standalone bse). These choices of the time-step parameters practically do not slow down nbody7 or isolated-binary runs.

2.3 Remnant natal kick

The remnant natal kick is based on observed kick distribution of NSs (for single NSs in the Galactic field, the one-dimensional kick velocity dispersion is v_{kick, NS} ≈ 165 km s⁻¹; Hobbs et al. 2005) that is reduced linearly by the fraction of material fallback in the SN. This ‘momentum conserving’ formulation, as given by equation (1) of Ba20, is identical to the natal kick treatment in other widely used population synthesis programs such as |$\rm{\small STARTRACK}$| (Belczynski et al. 2008), tres (Toonen et al. 2016), and |${\tt MOBSE}$| (Giacobbo et al. 2018). The amount and fraction of the SN material fallback are provided by the chosen remnant-mass scheme (see Section 2.2). Fig. 2 shows the fallback fraction, f_fb, as a function of carbon–oxygen core mass, M_CO, and remnant mass, M_rem, for F12-rapid+B16-PPSN/PSN remnant mass model (Section 2.2), at different metallicities.

Figure 2.

Supernova fallback fraction (alternatively, fallback factor), f_fb, as a function of the progenitor star’s carbon–oxygen core mass, M_CO (left-hand panel), and remnant mass, M_rem (right-hand panel). The data points are obtained directly from nbody7 computations in this work, of model clusters composed of only single stars initially whose ZAMS masses range from 0.08 to |$150.0\, \mathrm{M}_{\odot }$| and which are distributed according to the standard IMF. The f_fb−M_CO and f_fb−M_BH dependencies are shown at the four metallicities Z = 0.0002, 0.001, 0.01, and 0.02 (legend) for the remnant-mass model F12-rapid+B16-PPSN/PSN (Section 2.2). f_fb = 1.0 corresponds to direct-collapse BHs.

Apart from this standard momentum-conserving recipe, two of its variants can be opted for, namely, models for ‘convection-asymmetry-driven’ kick (Scheck et al. 2004, 2008; Fryer & Young 2007) and ‘collapse-asymmetry-driven’ kick (Burrows & Hayes 1996; Fryer 2004; Meakin & Arnett 2006, 2007). These alternative recipes are given by equations (2) and (3) of Ba20. A recipe for ‘neutrino-driven’ (Fuller et al. 2003; Fryer & Kusenko 2006) kick is also available as an option (equation 4 of Ba20). The effects of these various natal kick models on the retention of NSs and BHs in clusters, right after their birth, is demonstrated and discussed in detail in Ba20.

To summarize, taking into account the slow down of natal kicks due to fallback, the momentum-conserving and convection-asymmetry-driven kicks lead to similar BH and NS retention in all types of clusters. The collapse-asymmetry-driven kick, including fallback reduction, would retain most of the BHs in young massive clusters (hereafter YMC; Portegies Zwart, McMillan & Gieles 2010), open clusters (hereafter OC), and globular clusters (hereafter GC) which have central escape speeds, v_esc, ranging from a few km s⁻¹ (for low-mass OCs), through ∼10 km s⁻¹ (for YMCs, massive OCs, and GCs; Portegies Zwart et al. 2010), up to 100 km s⁻¹ (for massive GCs; Georgiev et al. 2009; Baumgardt & Hilker 2018). The relatively slower NSs in the collapse-asymmetry-driven kick model would lead to a higher NS retention in GCs (e.g. Goswami, Kiel & Rasio 2014; Kremer et al. 2020) but still similar retention as in the momentum-conserving and convection-asymmetry-driven cases, in YMCs and OCs. Note that in the above three kick models, the majority of the NSs that retain after birth are products of ECS, assuming that ECS produces low kicks ∼ few km s⁻¹ (Podsiadlowski et al. 2004; Gessner & Janka 2018). BHs of |$\gtrsim 10\, \mathrm{M}_{\odot }$| would also retain in clusters since they receive small or zero natal kicks as, for them, f_fb ≈ 1 (see Fig. 2). f_fb = 1 implies a failed SN, leading to the formation of a BH via direct collapse.

Finally, since neutrino-driven kick is independent of fallback, natal kicks of all core-collapse-SN produced NSs and BHs would be much higher and they would escape from YMCs, OCs, and GCs. ECS-produced NSs would be the only retaining remnants in the case of neutrino-driven natal kick mechanism. See section 3.2.4, figs 9 and 11 of Ba20 for further discussions. The retention fractions would be moderately affected if the BH and NS progenitor stars are in primordial binaries as demonstrated in Section 4 and table 1 of Ba20. Note that irrespective of the natal kick mechanism, a large fraction of NSs and BHs will always retain in nuclear star clusters (hereafter NSCs) owing to their very high v_esc between 300 and 500 km s⁻¹ (Schödel et al. 2014; Georgiev et al. 2009, 2016).

Based on these considerations, all computations in this work uses either momentum-conserving or collapse-asymmetry-driven natal kick prescription. The kick-velocity components of core-collapse-SN NSs (without fallback) have a Gaussian distribution with the Hobbs et al. (2005) dispersion of v_{kick, NS} = 265 km s⁻¹. The counterpart for ECS-NSs is assigned a much lower value of v_{kick, ECS} = 3 km s⁻¹ (Gessner & Janka 2018). In this work, no natal kick is applied to WDs; ∼ few km s⁻¹ natal kicks of WDs would not significantly influence the outcomes of the present models except for the least massive ones (see Section 3; Table C1), where WD kicks would have led to a somewhat higher mass-loss from the cluster over its long-term evolution (Fellhauer et al. 2003).

2.4 Natal spins of black holes

The O-type parent stars of BHs typically possess significant spin angular momentum, as observations of young, massive stars suggest (Ramírez-Agudelo et al. 2013, 2015). However, this doesn’t necessarily translate into high spins of their remnant BHs. The spin angular momentum of a BH is often expressed through its dimensionless counterpart, namely, the Kerr vector/parameter (dimensionless spin vector/parameter; Kerr 1963) |$\boldsymbol{a}$|⁠, which is defined as

$$\begin{eqnarray*} \boldsymbol{a} = \frac{c\boldsymbol{S}_{\rm BH}}{GM_{\rm BH}^2}. \end{eqnarray*}$$

(1)

Here, |$\boldsymbol{S}_{\rm BH}$| is the total angular momentum vector of a Kerr BH of (non-spinning) mass M_BH.

In the absence of any angular momentum supply on to the evolved parent star (e.g. due to mass accretion or tidal interaction from a binary companion), the remnant BH can potentially have a low spin if the angular momentum of the inner stellar core is carried away along with the stellar wind. Therefore, the magnitude of |$\boldsymbol{a}$| (a) depends on (i) how efficient is the transport of angular momentum from the core of the pre-SN star to its envelope and (ii) how high is the wind mass-loss rate. In this work, the BH spin estimates of B20 are adopted for BHs at their birth i.e. for those BHs that are non-recycled, having not undergoing any mass accretion or GR coalescence after their formation (see Sections 2.5 and 2.6 for the treatment of spins of second/higher generation and mass-accreted BHs). B20 estimate BH spins from detailed evolutionary models of fast-rotating single stars.

Fig. 3 (top row) shows the outcomes, as a function of metallicity Z, of the present adoption of the BH natal spin model of B20 that is based on rotating stellar models using the Geneva stellar-evolution code (Eggenberger et al. 2008; Ekström et al. 2012; hereafter ‘Geneva’ BH-spin model). Since the Geneva code does not include magnetic field, core-to-envelope angular momentum transport is purely convective and therefore inefficient leading to most of the BHs forming with a high spin, a = 0.85. Only the strong wind of the most massive stars are effective in eliminating the core’s angular momentum so that the most massive BHs form with 0.25 ≥ a ≥ 0.0 as see in Fig. 3 (top row; cf. fig. 1 of B20). Note that the M_CO–a relation (and hence the M_BH–a relation) depends on Z: here, the M_CO–a function is taken piecewise over Z-ranges in the same way as in Morawski et al. (2018, see their table 1).

Figure 3.

Magnitude of dimensionless spin parameter, a, of stellar-remnant BHs at birth (i.e. of BHs that have not undergone any mass accretion or GR coalescence after their formation) as a function of the progenitor star’s carbon–oxygen core mass, M_CO (left-hand column), and the BH mass, M_BH (right-hand column). The data points are obtained directly from nbody7 computations in this work. Top panels: The N-body models corresponding to these panels employ the ‘Geneva model’ of Belczynski et al. (2020, B20) for BH spin (Section 2.4) and comprise only single stars initially, whose ZAMS masses range from 0.08 to |$150.0\, \mathrm{M}_{\odot }$| and which are distributed according to the standard IMF. These models use the F12-rapid+B16-PPSN/PSN remnant mass prescription (Section 2.2). The models are for four metallicities, Z = 0.0002, 0.001, 0.01, and 0.02 as indicated in the legends. Middle panels: The N-body models corresponding to these panels employ the ‘mesa model’ of B20 for BH spin (Section 2.4). The other model characteristics are the same as those in the top panels except that the ‘weak’ PPSN mass prescription (Section 2.2) of Leung et al. (2019) is utilized (resulting in the non-monotonic behaviour w.r.t. M_BH which, here, extends up to |$\approx 50\, \mathrm{M}_{\odot }$| as opposed to the models in the top panels where M_BH is capped at |$40.5\, \mathrm{M}_{\odot }$| due to the use of B16-PPSN/PSN). Bottom panels: These M_CO−a and M_BH−a relations, employing the Geneva and mesa BH-spin prescriptions (see legend), are from model computations with Z = 0.0001 where all stars with |$M_{\rm ZAMS}\ge 16\, \mathrm{M}_{\odot }$| are in primordial binaries (see Section 3). The other model characteristics are the same as those in the top panels. Here, as a result of star–star mergers occurring in the massive binaries, M_BH exceeds the widely accepted PPSN upper limit of |$45\, \mathrm{M}_{\odot }$|⁠. Such ‘mass-gap BH’s possess high (low) spins when Geneva (mesa) models are applied.

Fig. 3 (middle row) shows the M_CO−a and M_BH−a relations in nbody7 for the BH natal spin model of B20 that is based on rotating stellar models using the |$\rm{\small MESA}$| stellar-evolution code (Paxton et al. 2011, 2015; hereafter ‘mesa’ BH-spin model). Unlike the Geneva code, |$\rm{\small MESA}$| includes magnetic field that makes the outwards angular momentum transport from the core much more efficient by forming a Tayler–Spruit magnetic dynamo (Spruit 2002; Fuller, Piro & Jermyn 2019). This causes BHs of all masses to form with a small residual spin, a ∼ 0.1, as seen in Fig. 3 (middle row; cf. fig. 2 of Ba20). Likewise the Geneva BH-spin model, the piecewise Z-dependence of the mesa  M_CO−a relation is divided over the same Z-ranges as in Morawski et al. (2018).

When a primordial massive-binary population as in the current models (see Section 3) is present, star–star mergers can form BHs in the PSN mass gap (see Section 2.7; Ba20). In conjunction with the Geneva (mesa) BH-spin model, such BHs will also form with high (low) spins as demonstrated in Fig. 3 (bottom row).

Note that according to the recent work by Fuller et al. (2019), the Tayler–Spruit magnetic dynamo can essentially extract all of the angular momentum of the proto-remnant core, leading to nearly non-spinning BHs. To consider this possibility, the option of assigning a = 0 irrespective of M_BH is kept in nbody7 (hereafter Fuller BH-spin model).

Note that, for the time being, the above natal spin models are applied only to BHs. The NS spins are still treated as defaulted in nbody7/bse.³

2.5 General relativistic merger recoil and final spin

The final merged BH out of a BBH merger would receive a recoil kick, |$\boldsymbol{{v}_{\rm k}}$|⁠, due to asymmetric radiation of GW. The magnitude and direction (w.r.t. the BBH orbit) of the recoil depends on the BBH’s mass ratio and on the magnitudes and directions w.r.t. the orbital angular momentum, |$\boldsymbol{{L_{\rm orb}}}$|⁠, of the component BHs’ |$\boldsymbol{a}$|s (Pretorius 2005; Campanelli et al. 2007; Hughes 2009). If the merging BHs’ spins are zero, |$\boldsymbol{{v}_{\rm k}}$| will be aligned along the line joining the BHs just before the merger. Its magnitude, v_k, is zero (small) for equal-mass (extreme-mass-ratio) components and maximizes to ≈170 km s⁻¹ at the mass ratio of ≈1/2.9 (Baker et al. 2007, 2008).

The GW emission is particularly asymmetric if the merging BHs are spinning. In that case, depending on the spins’ magnitudes and orientations, |$\boldsymbol{{v}_{\rm k}}$| will as well have an in-orbital-plane component perpendicular to the mass axis and a component perpendicular to the orbital plane. For (near-) maximally spinning BHs, this off-plane component typically dominates and can well exceed 500 km s⁻¹; for certain configurations, it can reach ≈3000 km s⁻¹ (Campanelli et al. 2007; Baker et al. 2008; van Meter et al. 2010; Lousto & Zlochower 2013).

For typical configurations, the orbital angular momentum is the primary contributor to the final (dimensionless) spin, |$\boldsymbol{{a}_{\rm f}}$|⁠, of the merged BH (Rezzolla et al. 2008). If the merging BHs’ spins are zero, then the only source of angular momentum in the BBH system is |$\boldsymbol{{L_{\rm orb}}}$|⁠. Accordingly, |$\boldsymbol{{a}_{\rm f}}$|⁠, will be aligned with |$\boldsymbol{{L_{\rm orb}}}$|⁠, and, for equal-mass components, will have the magnitude a_f ≈ 0.7 (Pretorius 2005). With finite, misaligned spins of the merging BHs, |$\boldsymbol{{a}_{\rm f}}$| will be misaligned relative to |$\boldsymbol{{L_{\rm orb}}}$| and its magnitude will be augmented (suppressed) w.r.t. the non-spinning-merger value if the spins are pro-aligned (anti-aligned). Since for most configurations, |$\boldsymbol{{L_{\rm orb}}}$| dominates the BBH system’s angular momentum budget, the misalignment is, typically, <10°.

As a first attempt to incorporate on-the-fly NR-based treatments of a BBH merger that takes place while the BBH is bound to the cluster, nbody7 utilizes the NR-based fitting formulae of van Meter et al. (2010) for evaluating the components of |$\boldsymbol{{v}_{\rm k}}$|⁠. These formulae incorporate cases where the BHs’ spins are inclined w.r.t. |$\boldsymbol{{L_{\rm orb}}}$| and hence would undergo spin-orbit precession during the in-spiral and merger phases. These formulae agree with NR outcomes within |$5{{\ \rm per\ cent}}$|⁠. As for |$\boldsymbol{{a}_{\rm f}}$|⁠, the NR-based fitting functions of Rezzolla et al. (2008) are applied. When the members of the BBH do not derive from the same primordial binary (i.e. the BBH is dynamically assembled), the orientations of the BHs’ spins are chosen randomly and isotropically (their dimensionless magnitudes being according to the chosen BH-spin model; see Section 2.4) to evaluate |$\boldsymbol{{v}_{\rm k}}$| and |$\boldsymbol{{a}_{\rm f}}$|⁠. If both BH members retain the membership of the same parent primordial binary, then, considering the parent stars’ close birth locations and potential interactions thereafter, the BHs’ random relative spin orientations are restricted within ≤90°. The same procedure is followed for computing |$\boldsymbol{{v}_{\rm k}}$| and |$\boldsymbol{{a}_{\rm f}}$| of BBH mergers that take place within a Hubble time after getting ejected from the cluster, using a standalone version of the above-mentioned|$\boldsymbol{{v}_{\rm k}}=0$| NR formulae.

For all BNS and NSBH mergers, is assigned and |$\boldsymbol{{a}_{\rm f}}$| is evaluated in the same way as for BBH mergers, as a preliminary treatment. In both cases, the merger product is assumed to be a BH of mass equal to the total binary mass.

For all in-cluster and ejected GR mergers, the ‘effective spin parameter’ χ_eff, which is a measure of the spin-orbit alignment of the merging system (Ajith et al. 2011; Abbott et al. 2017), is evaluated according to the expression

$$\begin{eqnarray*} \chi _{\rm eff}= \frac{{ M}_1{ a}_1\cos \theta _1 + { M}_2{ a}_2\cos \theta _2}{{ M}_1+{ M}_2}. \end{eqnarray*}$$

(2)

Here, M₁, M₂ are, respectively, the masses of the merging members with dimensionless spins |$\boldsymbol{{a}_1}$|⁠, |$\boldsymbol{{a}_2}$| that project with angles θ₁, θ₂ on |$\boldsymbol{{L_{\rm orb}}}$|⁠.

After assigning |$\boldsymbol{{v}_{\rm k}}$|⁠, the escape algorithm of nbody7 (Aarseth 2003) determines whether the merged BH should remain bound in the cluster or escape, in the usual way. If the BH remains bound, it is treated alike other cluster members. The a_f of the BH is also remembered and recycled if the BH engages itself in further GR mergers. How |$\boldsymbol{{v}_{\rm k}}$| is assigned and spins of first and higher generation BHs are tracked are detailed in Appendix  A.

2.6 Star–star and star–remnant mergers

In a dynamically active dense stellar system, mergers among stellar entities can take place due to either binary evolution or dynamical interactions (in-orbit collisions, due to binary–single and binary–binary encounters boosting the binaries’ eccentricities, and hyperbolic collisions). In both standalone and nbody7/bse, a complete mixing is assumed in composing the merged star (Hurley et al. 2002). In this work, all star–star mergers are assigned a mass-loss, as indicated in theoretical studies (e.g. de Mink et al. 2013). The amount of mass-loss is taken to be equal to a fixed fraction, f_mrg, of the secondary’s mass (the less massive of the merging stars, at the time of the merger). The merged product is a rejuvenated star of type and age according to the merger-product-type scheme of Hurley et al. (2002). The total mass of the merged product is equal to the total stellar mass just before the merger minus the merger mass-loss.

Furthermore, in a BH-star merger (forming a BH-Thorne-Zytkow-object; hereafter BH-TZO), a fixed fraction, f_TZ, of the star’s mass is assumed to be accreted on to the BH, increasing the BH’s final mass. The post-stellar-accretion BH is assumed to be maximally spinning, i.e. has a = 1 irrespective of the BH’s natal spin, due to the associated accretion of angular momentum. Observations of high-mass X-ray binaries indeed suggest near-maximal spins of accreting BHs (Miller & Miller 2015). (This treatment is different from that of a BBH merger product whose spin is assigned based on NR; see Section 2.5.) As shown in Banerjee (2019), a significant BH-TZO accretion can lead to BHs well within the PSN mass gap, at all metallicities. Finally, if a BH forms during a tidal-interaction (symbiotic), mass-transfer, or CE episode with a stellar binary companion, then also a = 1 is assigned to it, as supported by observations and theoretical modelling (e.g. Qin et al. 2019).

No matter accretion is assumed on to the NS, in an NS-star merger (forming an NS-TZO; Thorne & Zytkow 1975), whose final outcome is just the NS (the stellar material of the TZO is assumed to be completely irradiated away by the NS).

See Appendices  A and  B for the details of the implementations of these aspects in updated nbody7.

2.7 Black holes in the ‘upper’ and ‘lower’ mass gaps

Combining with the discussions from Sections 2.2–2.6, a BH can be in the PSN (upper) mass gap, in the classical sense, due to (i) weak PPSN, (ii) a star–star merger, (iii) a BBH merger, (iv) a BH-TZO accretion. In a dynamically active environment of a dense stellar cluster, all the four processes can, in principle, take place and can even combine. Among these, a PSN-gap BH from channel (i) or (ii) is a first-generation BH, that from channel (iii) is a second-generation BH, and that from channel (iv) is a first (second) generation BH if the accreting BH is of first (second) generation.⁴ A PSN-gap BH out of channel (i) or (ii) will have its spin as per the adopted BH natal spin model (Section 2.4) and that from channel (iv) will necessarily be maximally spinning, as assumed here (Section 2.6). For most BBH mergers, channel (iii) will also produce highly spinning BHs (a ≥ 0.5). However, depending on the spin-orbit configuration of the merging BHs, the second-generation (PSN-gap) BH can also be of moderate or low spin (a < 0.5), as discussed in Belczynski & Banerjee (2020).

If the PSN-gap BH is retained in the cluster, it can, in principle, participate in (further) GR mergers and be detectable by ground- and space-based GW detectors such as LIGO-Virgo and LISA (laser interferometer space antenna Amaro-Seoane et al. 2017). Among the above four channels, channel (iii) will typically eject the PSN-gap BH from any cluster when one or both of the merging, first-generation BHs are highly spinning (but see Belczynski & Banerjee 2020). When they are equal (close) in mass and are (nearly) non-spinning (Fuller BH-spin model in Section 2.4), they would retain in most OCs and more massive systems. For small BH spins (mesa BH-spin model in Section 2.4) or for unequal-mass non-spinning BHs, the recoiled BH would typically retain in NSCs (Miller & Lauburg 2009; Antonini et al. 2019) but escape from OCs, YMCs, and GCs. These inferences are based on the GR merger recoils as discussed in Section 2.5. The other channels will retain the PSN-gap BH in any cluster (all PPSN BHs and core-collapse-SN BHs of |$\gtrsim 10\, \mathrm{M}_{\odot }$| are direct collapse BHs, in the present remnant schemes; see Section 2.3).

Note that if a seed BH grows in mass via extreme-mass-ratio inspirals (hereafter EMRI) with many small BHs, as proposed in, e.g. Miller & Hamilton (2002), Hughes & Blandford (2003), then a cluster could sustain many generations of BBH mergers, due to the small GR recoils in EMRIs. That way, BHs can enter the PSN gap with low spins. In the present computed model clusters (Section 3), BHs derive from a continuous stellar initial mass function (hereafter IMF) and comprise a (nearly) continuous mass distribution. In that case, seed BH growth will be ineffective. (The present models also do not undergo an early runaway stellar merger episode as in, e.g. Portegies Zwart & van den Heuvel 2007; Fujii & Portegies Zwart 2013 that would produce a seed BH.)

As for the lower mass gap between NSs and BHs (Section 2.2), BHs can lie in the gap only in the case of F12-delayed remnant mass scheme that does not produce this gap (Section 2.2). This is due to the current assumptions regarding NS-star mergers (Section 2.6). Note that such low-mass BHs are formed with small SN fallback (Section 2.2), so that they receive large natal kicks in the momentum-conserving kick scenario (Section 2.3). Therefore, to retain such low-mass BHs in OCs, YMCs, and GCs, it is also necessary to adopt the collapse-asymmetry-driven kick model (Section 2.3) instead.

With F12-rapid remnant model, for which this lower gap is inherent (Section 2.2; Fig. 1), BNS mergers can produce BHs in the gap that would then retain in the clusters (Section 2.5) and potentially participate in further GR mergers. However, in-cluster BNS mergers are rare in the present computed models (see Section 3, Table C1).

3 DIRECT, POST-NEWTONIAN MANY-BODY COMPUTATIONS WITH UPDATED STELLAR-EVOLUTIONARY AND BLACK HOLE SPIN MODELS

Using the updated nbody7 as described in Section 2, a set of 65 long-term N-body computations of model stellar clusters are performed with varied remnant-mass, BH natal spin, metallicity, and initial cluster parameters. The initial cluster parameters and model choices are given in Table C1. These runs are performed over about an year whilst the code developments in nbody7 which is why, at present, the evolutionary model grid in Table C1 is somewhat heterogeneous. However, the present set is still composed of representative examples of the various remnant and natal spin schemes and their observable outcomes (see Section 4).

The initial model clusters have Plummer (1911) profiles with masses |$1.0\times 10^4\, \mathrm{M}_{\odot }\le M_{cl}(0) \le 1.0\times 10^5\, \mathrm{M}_{\odot }$|⁠, half-mass radii 1.0 pc ≤ r_h(0) ≤ 3.0 pc, overall (see below) primordial-binary fractions 0.0 ≤ f_bin(0) ≤ 0.1, and have metallicities 0.0001 ≤ Z ≤ 0.02. All the computed models are taken to be initially virialized (Spitzer 1987) and unsegregated and they are subjected to a solar-neighbourhood-like external galactic field. All initial models are composed of ZAMS stars. All models are evolved until 11.0 Gyr unless the cluster is dissolved earlier.⁵

Note that the values of f_bin(0) > 0.00 quoted in Table C1 represent the overall initial binary fraction in the cluster for the entire stellar-mass distribution, the latter being taken to be the standard IMF (Kroupa 2001) and in the ZAMS mass range |$0.08\!-\!150.0\, \mathrm{M}_{\odot }$|⁠. However, as in Papers II and III, the initial binary fraction of the O-type stars (ZAMS mass ≥M_crit; |${M}_{\rm crit}=16.0\, \mathrm{M}_{\odot }$|⁠), which are paired among themselves, is taken to be f_Obin(0) = 100 per cent, to be consistent with the observed high binary fraction among O-stars in young clusters and associations (e.g. Sana & Evans 2011; Sana et al. 2013). The initial binary fraction among the non-O-type binaries would still be ≈f_bin(0) since for the adopted standard stellar IMF, the O-stars comprise only a small fraction of the total stellar population. The O-star binaries are taken to initially follow the observed orbital-period distribution of Sana & Evans (2011) and a uniform mass-ratio distribution. The orbital periods of the non-O-star primordial binaries are taken to follow the period distribution of Duquennoy & Mayor (1991) and their mass-ratio distribution is also taken to be uniform. The initial eccentricities of the O-star binaries follow the Sana & Evans (2011) eccentricity distribution and those for the rest of the binaries are drawn from the thermal eccentricity distribution (Spitzer 1987). As explained in Paper II, such a scheme of including primordial binaries provides a reasonable compromise between the economy of computing and consistencies with observations.

When f_bin(0) = 0.00 is quoted in Table C1, an initial model with all single ZAMS star (standard IMF) is implied. Although having an observationally motivated primordial-binary population is more realistic, large-N runs, as in here, with primordial binaries is significantly more time consuming (and tedious). Therefore, to obtain reasonable statistics, initially single-star models are also performed. In a few primordial-binary runs, |${M}_{\rm crit}=5.0\, \mathrm{M}_{\odot }$| is taken instead of |$16.0\, \mathrm{M}_{\odot }$| as test cases; it is mentioned in the table’s footnote when so. Unless stated otherwise in Table C1, f_mrg = f_TZ = 0.5 (Section 2.6) is taken. In several runs, a reduced f_mrg and an enhanced f_TZ are assumed to facilitate BHs enter the PSN mass gap (Section 2.7) and explore its consequences on observable GR mergers (Section 4).

The star cluster models computed here represent young massive and open clusters that continue to form and dissolve throughout a gas-rich galaxy, as in, e.g. the Milky Way and Local Group galaxies. Note that the 1–3 pc initial size of the present models is consistent with the typical sizes observed in gas-free YMCs in the Milky Way and neighbouring galaxies (Portegies Zwart et al. 2010; Banerjee & Kroupa 2017). It is therefore implicit in the present models that these clusters have survived their assembling and violent-relaxation phases and have expanded to parsec scale sizes from sub-parsec sizes (due to, e.g. gas expulsion), as observed in newly formed, gas-embedded and partially embedded clusters and associations; see Banerjee & Kroupa (2018) and references therein. The long-term evolution explored here would essentially wash out any imprint (Heggie & Hut 2003) of the complex birth history of the clusters. The present models bridge between lower mass embedded-cluster-type initial conditions explored in works such as Kumamoto et al. (2019), Di Carlo et al. (2019) and GC-progenitor systems as in Askar, Arca Sedda & Giersz (2018), Kremer et al. (2020).

4 RESULTS

In the present paper, focus is given on how the GR-merger outcomes from the present N-body computations (Section 3; Table C1) agree with the measured parameters of those observed in LVC O1/O2 (Abbott et al. 2019a) and how they comply with the O3 candidates (https://gracedb.ligo.org/superevents/public/O3/).

4.1 Merger masses

Fig. 4 plots the primary and secondary masses, M₁, M₂ (M₁ ≥ M₂), of all the GR mergers from the computed models in Table C1. The vast majority these mergers are in-cluster BBH mergers as is the case with somewhat older prescriptions and treatments, for similar types of clusters, in Papers I, II, and III. See the discussions and analyses in these papers.

Figure 4.

Primary mass, M₁, versus secondary mass, M₂ (M₁ ≥ M₂), of GR compact-binary mergers from all the computed models in Table C1 (filled circles and triangles). They are colour-coded according to their delay times (colour bar). The circles represent mergers inside a model cluster while being bound to it and the triangles represent those happening in compact-binary systems after they get ejected from the cluster. The filled, black squares with error bars (90 per cent credible intervals) represent the observed merger-event data from LVC O1/O2 (GWTC-1, Abbott et al. 2019a). The ‘lower mass gap’ over |$3.0\, \mathrm{M}_{\odot }\le {M}_1,{M}_2\le 5.0\, \mathrm{M}_{\odot }$| is indicated with the grey shades.

Fig. 4 demonstrates that the trend of the M₁−M₂ combination, their ranges, and their scatter, in the mergers from the computed models, agree well with those of the GR merger events from LVC O1/O2 (taking into account the latter’s 90 per cent confidence intervals). For |${M}_1\gtrsim 30\, \mathrm{M}_{\odot }$|⁠, all BBH mergers happen at short delay times (the time interval between the beginning of the cluster evolution and the occurrence of the merger), t_mrg < 1.0 Gyr, except for a few ejected mergers (see also Fig. 9). In this work, an ejected double-compact binary is labelled as a ‘merger’ if its GR coalescence time, τ_insp, is less than the Hubble time, t_Hubble = 13.7 Gyr. The adopted τ_insp for this purpose is given by (Peters 1964)

$$\begin{eqnarray*} \tau _{\rm insp}\approx \frac{5}{64}\frac{c^5a_{\rm ej}^4\left(1-e_{\rm ej}^2\right)^{7/2}}{G^3 {M}_1{M}_2({M}_1+{M}_2)} \left(1+\frac{73}{24}e_{\rm ej}^2+\frac{37}{96}e_{\rm ej}^4\right)^{-1}, \end{eqnarray*}$$

(3)

a_ej and e_ej being the binary’s semimajor axis and eccentricity, respectively, at the time of its ejection.

A handful of mergers are ejected BNS mergers (left lower corner of Fig. 4). All these BNS mergers are ‘isolated’ mergers in the sense that the merging NS members have maintained their original primordial-binary membership. These binaries got ejected as eccentric BNSs due to either the later-born NS member’s natal kick or a dynamical encounter. The ‘dynamical heating’ effect of the BH core (see Section 1) prevents mass segregation of lighter members including NSs, inhibiting strong dynamical interactions and dynamical pairing among themselves until all BHs are nearly depleted (Paper II; Fragione, Pavlík & Banerjee 2018; Ye et al. 2019). The present set of models do not yield dynamically paired BNS mergers or NSBH mergers (but see Rastello et al. 2020). In contrast, as in Papers I–III and other studies involving intermediate-mass and massive stellar clusters (see Section 1 and references therein), most of the BBH mergers are dynamically assembled and take place inside the clusters (like the BNS mergers, there are only a few BBH mergers where the primordial binary membership is maintained until the merger).

4.1.1 Mergers in the ‘upper’ and ‘lower’ mass gaps

Fig. 4 shows that in the present models, BBH mergers take place with the primary being in the PSN mass gap (see Section 2.2), i.e. of |${M}_1\gt 45\, \mathrm{M}_{\odot }$|⁠. All these mergers happen inside clusters and at delay times t_mrg < 1.0 Gyr (see also Fig. 9). Section 2.7 discusses the various ways in which a BH can appear in the PSN gap. Among the four most prominent PSN-gap mergers in Fig. 4 (⁠|${M}_1\gtrsim 55\, \mathrm{M}_{\odot }$|⁠), only the one with |${M}_1\approx 80\, \mathrm{M}_{\odot }$| is an in-cluster second-generation merger; this particular sequential BBH merger event is demonstrated in Example 1 of Appendix  A (outcome from model 49 of Table C1). The rest are mass-gained BHs through significant BH-TZO accretion (f_TZ ≥ 0.5; see Table C1), which is the most common form of participation in PSN-gap BBH mergers, in the present models. The |${M}_1\approx 100\, \mathrm{M}_{\odot }$| BBH-merger primary has originated due to a f_TZ = 0.95 BH-TZO accretion in a low Z model (model 61 of Table C1). All the PSN-gap mergers took place in models with GC-like metallicities, of Z ≤ 0.001, where the BHs’ birth mass and mass gain are the largest.

Note that M₁−M₂s of the BBH mergers from the computed models show a linear trend up to |${M}_1\approx 40\, \mathrm{M}_{\odot }$|⁠, with mass ratios typically M₂/M₁ > 0.5 (but see Section 4.1.2; see also Fig. 9) and agreeing well with the O1/O2 merger-event data. This preferred pairing of similar masses as in the O1/O2 events (Fishbach & Holz 2020) comes out naturally of the present computations and is a consequence of dynamical pairing. As already explained in Paper II, despite pairing of unrelated members (e.g. members not belonging to the same primordial binary or that have formed well separated from each other), commonly termed as ‘random pairing’, the most massive members segregate to the innermost part of the host cluster and, hence, interact preferentially among each other. In fact, this selection effect operates at all levels: among the segregated members, the most massive ones would engage themselves in close encounters most frequently due to their stronger gravitational focusing (Spitzer 1987) and such encounters (e.g. three-body, binary–single, or binary–binary encounters) would most likely result in pairing up the most massive participants since such an outcome would be energetically favourable (Heggie & Hut 2003). Since the BH mass distribution up to |$\approx 40\, \mathrm{M}_{\odot }$| is continuous (see Fig. 1; figs 8 and 12 of Ba20), the pairing is automatically biased towards similar masses, typically resulting in 0.5 < M₂/M₁ < 1.0.

However, this trend deviates sharply for larger M₁, in the PSN gap, as seen in Fig. 4. (Since most of the computed models assume B16-PPSN/PSN model, the PSN gap actually begins from |$\approx 40\, \mathrm{M}_{\odot }$|⁠, in the present model set; see Section 2.2; for the sake of discussions, PSN gap will still refer to |$\gt 45\, \mathrm{M}_{\odot }$|⁠.) Within the PSN gap, the mass ratio tends to decline from 0.5 with M₁. This is, again, a consequence of random, dynamical pairing of the merging BBHs and of the fact that mass-gained PSN-gap BHs are, in most models, distinctively the most massive members. The computed data points, nevertheless, well encompass the most massive and marginally PSN-gap (but see Chatziioannou et al. 2019; Kimball et al. 2020) BBH merger event GW170729 (the most massive O1/O2 M₁−M₂ pair in Fig. 4).

For the typical ∼10 km s⁻¹ escape speeds from the present models, the only way a second-generation BH retains in a cluster is via mergers of two similarly massive, first-generation BHs, that did not undergo any previous matter interaction (see Section 2.6), with the Fuller BH-spin model (i.e. BHs having practically zero spin at birth; see Section 2.4). In all other cases, one or both of the merging BHs will have spins of a > 0.1 and the merger recoil speed will eject the merged BH from the cluster (see Sections 2.5 and 2.7).⁶ Even if the second-generation BH retains in a cluster, it does not necessarily participate in further (in-cluster or ejected) mergers and may get ejected via dynamical encounters, as seen in the present models that employ the Fuller BH spin. Nevertheless, such models do tend to produce a larger number of in-cluster BBH mergers; see Table C1. The probability of a retained second-generation BH to actually participate in further mergers would increase in higher escape-velocity systems like massive GCs, as seen in, e.g. Rodriguez et al. (2018, 2019). Also, owing to their ∼100 km s⁻¹ escape speed, the chances of retaining second and even higher generation BHs and of their participation in further mergers increases in NSCs, as explored in, e.g. Miller & Lauburg (2009), Antonini et al. (2019), and Arca Sedda et al. (2020).

Fig. 4 also demonstrates BBH mergers in the lower NS–BH mass gap (the shaded strips; see Sections 2.2 and 2.7), happening both inside clusters and after ejection. As explained in Section 2.7, under the present assumptions and adopted models, only the combination of F12-delayed remnant mass model (Section 2.2) with collapse-asymmetry-driven natal kick (Section 2.3) would retain such lower mass gap BHs in the clusters, after their birth. The lower mass gap mergers in Fig. 4 are from such models (see Table C1), exclusively. Note that these mergers happen at long delay times of 1 Gyr ≲ t_mrg ≲ 10 Gyr. This is expected since such low-mass BHs have to wait until their more massive counterparts are dynamically depleted from the cluster, before they can effectively participate in dynamical encounters, as also found in previous studies (e.g. Morscher et al. 2015; Chatterjee et al. 2017b). Note, further, that the clusters with F12-rapid remnant model, that have actually produced such lower mass gap mergers, are of higher Z (0.01 or 0.02), which produce and retain the least number of massive, |$\gtrsim 10\, \mathrm{M}_{\odot }$| BHs. Interestingly, LVC has indeed detected a few of such (lower) ‘mass gap’ mergers in their O3. Note that although in all model clusters the BH population depletes with time due to dynamical ejections (plus BBH mergers), BHs continue to be present in them until late evolutionary times, as demonstrated in Fig. C1. This happens due to the fact that the energetic dynamical interactions, that form BBHs and eject BHs, also lead to heating and expansion of the clusters which effect, in turn, moderates the BHs’ dynamical activities (see e.g. Breen & Heggie 2013; Morscher et al. 2015; Banerjee 2017).

4.1.2 Mass-asymmetric mergers: on GW190412

In the light of the above discussions, it would be worth highlighting the BBH merger event GW190412 of the LVC O3, published very recently by the collaboration (The LIGO Scientific Collaboration & the Virgo Collaboration 2020). With respect to all the LVC GW merger events so far, the most striking aspect of GW190412 is its highly asymmetric mass ratio of |${M}_2/{M}_1=0.28^{+0.13}_{-0.07}$|⁠. Although relatively rare, similarly low-mass ratio BBH mergers, with the primary below the PSN gap, do take place in the present computed models. In Fig. 4 (see also Fig. 9), they are the few lower mass ratio outliers of the overall O1/O2/computed-merger trend. Apart from these, the BBH mergers with M₁ within the PSN gap are also of M₂/M₁ ≲ 0.5, as discussed in Section 4.1.1. Overall, 3 per cent of all GR mergers from the present models are highly asymmetric in mass, in the sense of having M₂/M₁ < 0.5, with the most asymmetric one being of M₂/M₁ ≈ 0.3 (primary below the PSN gap). All but one of such mergers are in-cluster mergers; since the majority of the mergers are in-cluster (see above), this is also reflected in their mass-asymmetric subset (the asymmetric merger fraction is 4 per cent w.r.t. the in-cluster mergers).

Despite the ‘dynamical selection’ effect discussed in Section 4.1.1, asymmetric outliers are statistically possible in random pairing. The most concentrated ones among the present model set produce the asymmetric mergers; those with |$M_{cl}(0)\ge 7.5\times 10^5\, \mathrm{M}_{\odot }$| and with |$M_{cl}(0)=5\times 10^4\, \mathrm{M}_{\odot }$|⁠, r_h(0) = 1.0 pc (Table C1). Higher stellar density allows central mass segregation of BHs down to lower mass BHs (Section 1 and references therein) promoting close interactions and pairings among more dissimilar masses. If the GR inspiral (Section 4.4) begins early enough in the cluster evolution that many exchange interactions do not equalize the in-cluster BBH component masses (e.g. Antonini & Gieles 2020), highly mass-asymmetric mergers are possible, as the present computations demonstrate. Indeed, all the mass-asymmetric BBH mergers, in the present computations, take place within t_mrg < 1 Gyr (see Figs 4 and 9; the ejected one is ejected from the parent cluster at ≈500 Myr cluster-evolutionary time).

This is further supported by the very recent work of Di Carlo et al. (2020), who have indeed obtained a larger fraction of such asymmetric, GW190412-like BBH mergers in their computations of a large number of cluster models with r_h(0) 0.2–1.5 pc and M_cl(0) |$10^3\!-\!3\times 10^4\, \mathrm{M}_{\odot }$|⁠. Such models have shorter two-body relaxation and hence BH mass-segregation times (Spitzer 1987; Banerjee et al. 2010; Breen & Heggie 2013; Antonini & Gieles 2020) and, hence, are expected to produce more of such asymmetric BBH mergers. The shorter mass-segregation times in their models also induce runaway stellar mergers, enabling the formation and participation in mergers of PSN-gap BHs which are also, typically, asymmetric (Di Carlo et al. 2019; see also Section 4.1.1).

Therefore, YMCs and OCs serve as sites that naturally enable mass-asymmetric BBH mergers. Apart from star clusters, other astrophysical environments that enable close interactions and mergers among unequal-mass BHs are active galactic nucleus (AGN) gas discs (e.g. Secunda et al. 2019) and field triple and higher order systems (e.g. Antonini et al. 2017; Silsbee & Tremaine 2017; Fragione & Kocsis 2019).

4.2 Merger final spins and effective spin parameters: on GW170729

Figs 5 and 6 show the final spins, a_fs, and the effective spin parameters, χ_eff (Section 2.5), of all the mergers from the computed models in Table C1. In these figures, the outcomes from the models employing high BH natal spins, i.e. Geneva BH-spin scheme, and low BH natal spins, i.e. mesa and Fuller BH-spin schemes, are distinguished (left-hand and right-hand panels, respectively).

Figure 5.

Total mass, M_tot, versus final spin, a_f, of GR mergers from those computed models in Table C1 that employ the high BH natal spin, i.e. Geneva BH-spin model (left-hand panel) and that employ the low BH natal spin, i.e. mesa and Fuller BH-spin models (right-hand panel). The legend for the points is the same as in Fig. 4. The a_f values for both the in-cluster and ejected mergers are computed as described in Section 2.5.

Figure 6.

Total mass, M_tot, versus effective spin parameter, χ_eff, of GR mergers from those computed models in Table C1 that employ the high BH natal spin, i.e. Geneva BH-spin model (left-hand panel) and that employ the low BH natal spin, i.e. mesa and Fuller BH-spin models (right-hand panel). The legend for the points is the same as in Fig. 4. The χ_eff values for both the in-cluster and ejected mergers are computed as described in Section 2.5.

As one expects, high BH natal spin would result in significant scatter of a_f and χ_eff values around their equal-mass, non-spinning merger values, a_f ≈ 0.7 and χ_eff = 0.0, due to the random mass pairing and spin orientations of the merging members, in dynamical mergers (Section 2.5). This is demonstrated in Figs 5 and 6 (their left-hand panels). The notable tapering of the scatter with increasing M_tot is due to the fact that even in Geneva BH-spin model, depending on the BH-progenitor star’s metallicity, the most massive BHs still receive low natal spins (Section 2.4; top panel of Fig. 3). This leads to the convergence to the non-spinning a_f and χ_eff values with increasing M_tot. Beyond |$M_{\rm tot}\approx 70\, \mathrm{M}_{\odot }$|⁠, mass-gained BHs become the abundant participants of BBH mergers (Section 4.1.1), which BHs have high or maximal spins (Sections 2.6 and 2.7). Therefore, the scatter of a_f and χ_eff resumes for |$M_{\rm tot}\gtrsim 70\, \mathrm{M}_{\odot }$| (left-hand panels of Figs 5 and 6). Overall, the scatter in a_f and χ_eff, for the computed mergers with Geneva BH-spin model, well exceeds that in the O1/O2 merger events (left-hand panels of Figs 5–7).

Figure 7.

Distribution of effective spin parameter, χ_eff (shaded histogram), of GR mergers from those computed models in Table C1 that employ the high BH natal spin, i.e. Geneva BH-spin model (left-hand panel) and that employ the low BH natal spin, i.e. mesa and Fuller BH-spin models (right-hand panel). The χ_eff values for in-cluster and ejected mergers are computed as described in Section 2.5. The filled, black squares with error bars (90 per cent credible intervals) represent the observed χ_eff values in LVC O1/O2 merger events (GWTC-1, Abbott et al. 2019a). For clarity, the LVC O1/O2 data points are displaced along the Y-axis in the increasing order of their date of discovery.

On the other hand, BBH mergers from the models with low BH natal spin (mesa and Fuller BH-spin model; Section 2.4; Fig. 3) yield χ_eff and a_f values that agree well with the observed values of these quantities from O1/O2 events and as well with their overall trend and scatter (except for GW170729, see below, the O1/O2 values are similar to the non-spinning values). This is demonstrated in Figs 5–7 (their right-hand panels). None the less, as discussed above, highly spinning, mass-gained BHs cause the model mergers to fluctuate between spin-orbit-aligned and spin-orbit-anti-aligned values of a_f and χ_eff for |$M_{\rm tot}\gtrsim 70\, \mathrm{M}_{\odot }$|⁠, as seen in Figs 5 and 6 (their right-hand panels). In particular, the distinctively high (positive) value of a_f (χ_eff) for GW170729 could be reproduced (taking into account their 90 per cent confidence intervals).

For |$M_{\rm tot}\lesssim 70\, \mathrm{M}_{\odot }$|⁠, a few mergers still show large scatter in mesa and Fuller BH-spin models. These are either lower mass, mass-gained BHs or those mergers where the primordial-binary membership is maintained and therefore a partial spin-orbit-alignment is imposed (Section 2.5). These include the ejected BNS mergers (Section 4.1).

Overall, the mergers from the computed models can explain well the masses, final spins, and spin-orbit alignments of the BNS and BBH mergers observed by LVC in their O1/O2, including the distinctively large masses, pro spin-orbit alignment, and high final spin of GW170729, provided the natal spins of the BHs are taken to be a ≲ 0.1 (mesa or Fuller BH natal spin model). The results also imply that beyond |$M_{\rm tot}\approx 70\, \mathrm{M}_{\odot }$|⁠, both high and low final spins and (correspondingly) pro- and anti-aligned spin-orbit configurations can generally be expected in BBH mergers despite low natal spins of BHs. This is also occasionally possible for |$M_{\rm tot}\lesssim 70\, \mathrm{M}_{\odot }$|⁠. Example 3 (from model 61 of Table C1) of Appendix  A shows a specific example of a GW170729-like merger between a mass-gained (through BH-TZO accretion) and an intact BH.

In the low BH natal spin models of the current set, very massive BBH mergers, up to |$M_{\rm tot}\approx 140\, \mathrm{M}_{\odot }$|⁠, are obtained that lead to (final) BHs of intermediate-mass and moderate to near-maximal spins (a_f ≳ 0.5; see right-hand panel of Fig. 5). These mergers happen in low-Z models through either second-generation BBH merger (when Fuller BH-spin is applied) or significant BH-TZO accretion (f_TZ ≳ 0.7). The current set with Geneva BH natal spins assume smaller f_TZ (Table C1) and second-generation mergers in them is unlikely (see Section 4.1.1), which is why such massive BBH mergers are missing in them (left-hand panel of Fig. 5).

4.3 General-relativistic coalescences inside clusters driven by compact subsystems

As seen in Table C1, the majority of the GR mergers (most of which are dynamically driven BBH mergers) take place inside clusters. While, overall, the total number of mergers increases with increasing M_cl(0) and decreasing r_h(0), the growth is mainly in in-cluster mergers, over the cluster mass and size ranges explored here, namely, |$10^4\, \mathrm{M}_{\odot }\le M_{cl}(0)\le 10^4\, \mathrm{M}_{\odot }$|⁠, 1.0 pc ≤ r_h(0) ≤ 3.0 pc. As explained and demonstrated in Paper III (see sections 1 and 3.1 therein), most of these in-cluster mergers are driven by the short time-scale, few-body dynamics of dynamically formed compact subsystems that are made of one or more of the retained stellar remnants, typically of BHs. As shown therein (and found in the present models as well), the in-cluster triples or higher order subsystems are strongly perturbed so that it is, typically, the chaotic part of the few-body evolution (e.g. Antonini et al. 2016; Samsing & D’Orazio 2018) of the subsystem rather than its secular evolution (e.g. Katz, Dong & Malhotra 2011; Lithwick & Naoz 2011) that drives the eccentricity boost of its innermost binary, leading to the latter’s GR inspiral.

Fig. 8 shows the hierarchy of the (innermost) triple configuration, |${\mathcal {R}}$| (X-axis), against their respective Kozai–Lidov (hereafter KL) time-scale, T_KL (Y-axis), of the in-cluster subsystems hosting a GR merger, from all the computed models. The triple hierarchy is measured by the ratio of the triple’s outer periastron to inner apoastron, |${\mathcal {R}}\equiv a_0(1-e_0)/a_i(1+e_i)$|⁠, a₀ and e₀ being the triple’s outer semimajor axis and eccentricity, respectively, and a_i and e_i being the inner semimajor axis and eccentricity, respectively. The triple’s secular KL time period, T_KL, is given by (Kiseleva, Eggleton & Mikkola 1998)

$$\begin{eqnarray*} T_{\rm KL}=\frac{2P_0^2}{3\pi P_i}\left(1-e_0^2\right)^{3/2}\frac{{M}_1+{M}_2+m_0}{m_0}, \end{eqnarray*}$$

(4)

where P_i and P₀ are the Keplerian periods of the triple’s inner and outer orbits, respectively, and M₁, M₂, and m₀ are the masses of the inner binary’s components and the outer member, respectively.

Figure 8.

The configurations of the subsystems hosting a GR binary merger occurring inside the computed model clusters in Table C1 (filled circles). If the merger happens to the innermost binary of a triple (higher order) subsystem, then the abscissa is the ratio, |${\mathcal {R}}$|⁠, of the outer periastron to the inner apoastron of the triple (innermost triple), which measures the extent of hierarchy of this triple. The ordinate represents the Kozai–Lidov time period (Myr), T_KL, of the triple (innermost triple). The colour coding (colour bar) represents the total membership, n_p, of the merger-hosting subsystem. If the merger happens in an in-cluster binary (n_p = 2; filled, black circles), then the binary’s apoastron to periastron ratio is plotted along the abscissa and its orbital period (Myr), T_orb, is plotted along the ordinate. These are the configurations of the host subsystems at the start of the binary’s final GR in-spiral. The grey, vertical line indicates the canonical stability limit, |${\mathcal {R}}=3$|⁠, for hierarchical triple systems (Mardling & Aarseth 1999).

Fig. 8 is essentially similar to fig. 1 of Paper III but with a much larger number of data points from the present, much larger set of computed models. The data points are colour-coded with the total membership, n_p (n_p = 2 is a binary), of the subsystem containing the GR in-spiralling binary. Note that n_p simply represents the total number of members that are bound to the merging binary at the instant the binary is labelled COALESCENCE during the PN integration of its host subsystem via |$\rm{\small ARCHAIN}$| (Mikkola & Tanikawa 1999; Mikkola & Merritt 2008), based on the GR merger criteria in |$\rm{\small ARCHAIN}$| (Aarseth 2012); see Examples 1, 2, and 3 of Appendix  A. These criteria boil down to the arrival of an (a_i, e_i) combination of the innermost PN binary such that it can safely be considered to be spiralling in up to its merger, unperturbed.⁷ For each in-cluster merger, n_p is determined by searching for Chain members and Chain perturbers close to the COALESCENCE, as detailed in Paper III (see section 2.3.1 therein). At other times, n_p can be different or comprise different members since subsystems are, typically, strongly perturbed (see above, Paper III).

Note that in Fig. 8, n_p = 2 implies that the binary has merged without being a part of a higher order subsystem: due to a recent strong dynamical interaction, the binary has shrunk and/or become sufficiently eccentric that it has merged before meeting the next encounter.⁸ In such cases, the apoastron to periastron ratio of the binary at its COALESCENCE stage, (1 + e_i)/(1 − e_i), and the corresponding orbital period, T_orb = P_i, are plotted along the X-axis and Y-axis, respectively.

Fig. 8 suggests that binary GR mergers inside clusters take place in compact binaries that are, at the time of the beginning of the final inspiral, either a part of a dynamically formed, higher order subsystem or without any host subsystem, in comparable probabilities (merger inside a subsystem being of moderately higher probability). When part of a subsystem, a merger happens more often in triples than in higher order systems: in the present model set, mergers inside quadruples and those inside n_p > 4 systems are comparable in number. The figure also suggests that the majority of the (innermost) triple systems hosting a merger have hierarchy ratio |${\mathcal {R}}\lesssim 10$|⁠, implying that they are typically of unstable to hierarchical configuration (Mardling & Aarseth 1999; Mardling 2008). For further insight, it would be worth studying the secular evolutionary properties of these systems (as in, e.g. Fragione et al. 2020a) or integrate them, after extracting them from the clusters, through PN few-body integration (as in, e.g. Antonini, Murray & Mikkola 2014; Antonini et al. 2016; Kimpson et al. 2016) which will be taken up in a forthcoming study.

4.4 General-relativistic inspirals

Fig. 9 shows the merger delay times, t_mrg (X-axis), against M_tot (left-hand panel) and M₂/M₁ (right-hand panel) for all the mergers from the computed models in Table C1. The computed merger data points well encompass the 90 per cent confidence limits of M_tots and M₂/M₁s of the O1/O2 merger events (horizontal lines), except for a few |$M_{\rm tot}\gtrsim 80\, \mathrm{M}_{\odot }$| mergers with t_mrg < 1.0 Gyr as consistent with Fig. 4 (Section 4.1). The colour coding represents the inspiral time, τ_insp, corresponding to a merger as estimated from equation (3) (a_ej, e_ej replaced by a_i, e_i, see Section 4.3, for in-cluster mergers). All in-cluster mergers (filled circles), consistently, have τ_insp  |$\ll$| 1 Myr (see colour coding), which are therefore much shorter than ∼Myr strong interaction time-scale among BHs in such model clusters with ≲10 km s⁻¹ velocity dispersion (Spitzer 1987; Bacon, Sigurdsson & Davies 1996).

Figure 9.

Total mass, M_tot (left-hand panel), and secondary-to-primary mass ratio, M₂/M₁ (right-hand panel), versus merger delay time, t_mrg, of GR compact-binary mergers from all the computed models in Table C1 (filled circles and triangles). They are colour-coded according to their GR in-spiral times, τ_insp (colour bar). The circles represent mergers inside a model cluster while being bound to it and the triangles represent those happening in compact-binary systems after they get ejected from the cluster. The 90 per cent credible intervals of the total masses and mass ratios of the LVC O1/O2 merger events are shown as horizontal lines.

Fig. 10 shows the unperturbed orbital evolution up to coalescence (Section 4.3), starting from the a_i, e_i (a_ej, e_ej) for in-cluster (ejected) mergers at a cluster-evolutionary time/delay time t_insp. The orbits are obtained by integrating Peters (1964) orbit-averaged orbital shrinkage and circularization equations due to quadrupole GW radiation (PN-2.5 term) and are plotted in the log₁₀(e) − log₁₀(f_GWp) plane, f_GWp being the peak-power GW frequency. f_GWp is given by (Wen 2003)

$$\begin{eqnarray*} f_{\rm GWp}=\frac{\sqrt{G({M}_1+{M}_2)}}{\pi } \frac{(1+e_t)^{1.1954}}{\left[a_t\left(1-e_t^2\right)\right]^{1.5}}, \end{eqnarray*}$$

(5)

a_t, e_t being the instantaneous orbital parameters of the in-spiralling binary. The curves are colour-coded with the time since t_insp.

Figure 10.

Final (unperturbed) orbital inspiral curves of all the GR binary mergers (both, in-cluster and ejected) from the computed models in Table C1 in the log₁₀(e) − log₁₀(f_GWp) plane, e being the binary’s eccentricity and f_GWp being the peak GW frequency (see the text). The gold and the blue shades represent the characteristic detection frequency bands of the LISA and the LIGO-Virgo instruments. The curves are obtained by integrating the Peters (1964) orbit-averaged semimajor axis and eccentricity decay equations. They are colour-coded according to the time since the beginning of the final inspiral (colour bar). When e ≲ 0.7 (when f_GWp corresponds to the ≈10th or a lower harmonic), GW from a binary is detectable by LISA. LIGO-Virgo will detect an eccentricity in the inspiral waveform if e ≳ 0.1. These limits in e are indicated by the grey, horizontal lines.

Fig. 10 (cf. fig. 7 of Paper II) clearly shows that the vast majority of the inspirals begin with f_GWp in the LISA’s characteristic frequency band or at an even lower frequency (especially, the ejected mergers). Then, after a time of |$\tau _{\rm insp}^\prime$|⁠,⁹ the merger happens after traversing through the decihertz and LIGO-Virgo’s frequency bands. As in Fig. 9, |$\tau _{\rm insp}^\prime \ll 1$| Myr for the in-cluster mergers and |$1{\rm ~Gyr}\lesssim \tau _{\rm insp}^\prime \lesssim 14{\rm ~Gyr}$| for the ejected mergers. Note that |$\tau _{\rm insp}^\prime$| reaches down to seconds for the few in-cluster inspirals that begin right at the LIGO-Virgo band; see Fig. 10. These type of mergers begin with eccentricity e_i  |$\gg$| 0.1, so that they can be detected as highly eccentric mergers in the LIGO-Virgo band. All other in-cluster and ejected mergers enter the LIGO-Virgo band (f_GWp > 10 Hz) with e < <0.1 so that they would appear as circular inspirals by LIGO-Virgo (Huerta & Brown 2013). In the present model set, such eccentric LIGO-Virgo mergers comprise ≈5 per cent of the total number of mergers (≈6 per cent of all in-cluster mergers), which fraction is comparable to the estimate by Samsing (2018) and somewhat lower than that of Rodriguez et al. (2018). These studies involve typical GC models, which are an order of magnitude more massive than the present YMC and OC-type models, and would therefore contain tighter compact binaries, promoting eccentric mergers (Samsing et al. 2014).

Some in-cluster inspirals begin as highly eccentric in the decihertz frequency range (see Fig. 10) and would nearly circularize when they reach the LIGO-Virgo band. These inspirals serve as the ‘missing link’ dynamical mergers (Sesana 2016; Chen & Amaro-Seoane 2017) and would be of interest for the future, space-based decihertz (0.1–10 Hz) GW antenna such as DECIGO (Kawamura et al. 2008). The rest of the in-cluster mergers begin with f_GWp spread over the LISA band. Although they also initiate with e_i > 0.9 and therefore would be undetectable by LISA, they circularize within the LISA band to become of e ≲ 0.7, when they are detectable by LISA (Nishizawa et al. 2016, 2017; Chen & Amaro-Seoane 2017). Note that the typical time-scale of eccentricity decline from e_i, in the LISA band, is |$\tau _{\rm insp}^\prime \sim 0.1$| Myr (colour coding in Fig. 10). This means that, depending on the parent cluster’s formation epoch, light traveltime from the cluster’s distance, and instrument sensitivity, such binaries either are already in a state of being observable by LISA as persistent or semipersistent sources over its mission lifetime or will remain inaudible by LISA. It is not possible that such an in-spiralling binary stellar remnant ‘becomes’ visible by LISA over its mission lifetime (typically, 5 yr; Amaro-Seoane et al. 2017).

This is as opposed to all in-cluster and ejected mergers that either initiate in or enter DECIGO (LIGO-Virgo) band which have |$\tau _{\rm insp}^\prime \sim$| year to minute (⁠|$\tau _{\rm insp}^\prime \lesssim$| minute). Irrespective of inspiral time-scales, the numbers of persistent, semipersistent, and transient sources from YMCs and OCs in the Local Universe at the current cosmic epoch, in various instrument frequency bands, depend on a convolution of star formation history, cosmology, and instrument sensitivity. Such estimates will be taken up in a forthcoming study.

Fig. 11 (left-hand panel) shows the orbital inspiral curves of Fig. 10 in the log₁₀(h_c)−log₁₀(f_GWp) plane. In Fig. 11 (left-hand panel), the X-axis is the detector-frame peak GW frequency and h_c (Y-axis) is the GW characteristic strain on the detector at that frequency, assuming a comoving distance of D = 150.0 Mpc (redshift z ≈ 0.03) of the in-spiraling binary (see the text). h_c is given by (see Kremer et al. 2019 and references therein for a derivation)

$$\begin{eqnarray*} h_c^2 = \frac{2}{3\pi ^{4/3}}\frac{G^{5/3}}{c^3}\frac{M_{\rm chirp}^{5/3}}{D^2}\frac{1}{f_{\rm GWp}^{1/3}} \left(\frac{2}{n_p}\right)^{2/3}\frac{g(n_p,e)}{F(e)}, \end{eqnarray*}$$

(6)

where, f_GWp is the source-frame peak GW frequency (as given by equation 5 and plotted along the X-axis of Fig. 10), M_chirp ≡ (M₁M₂)^3/5/(M₁ + M₂)^1/5 is the source-frame chirp mass, n_p is the harmonic corresponding to f_GWp, g(n, e), with n = n_p, is the relative GW power function as in Peters & Mathews (1963), and F(e) is the eccentricity correction function as in Peters (1964) with (1 − e²)^−7/2 factored in. The source-frame values from the computed inspiral orbits in Fig. 10 are directly used in equation (6). In Fig. 11 (left-hand panel), the corresponding detector-frame, redshifted frequency, f_GWp/(1 + z), is plotted along the X-axis; since this quantity is very close to f_GWp for the chosen D, a new symbol is not used in Fig. 11 (left-hand panel). Note that n_p decreases with the orbital evolution such that n_p ≲ 10 for e ≲ 0.7 and n_p = 2 for e = 0. The design sensitivity (noise floor) curves of Advanced LIGO (LIGO document number LIGO-T1800044-v5), its proposed A+upgrade (LIGO document number LIGO-T1800042-v5), and LISA (Amaro-Seoane et al. 2017; Robson, Cornish & Liu 2019) are shown in the same plane.

Figure 11.

Left-hand panel: Final (unperturbed) orbital inspiral curves of all the GR binary mergers (both, in-cluster and ejected) from the computed models in Table C1, in the log₁₀(h_c)–log₁₀(f_GWp) plane. Here, the redshifted, detector-frame peak GW frequency is plotted along the X-axis and h_c, along the Y-axis, is the GW characteristic strain at that frequency, assuming a comoving distance of D = 150.0 Mpc (redshift z ≈ 0.03) of the in-spiraling binary (see the text). The same analytic GR orbital evolution as in Fig. 10 are utilized to obtain these curves and they bear the same colour coding. The design sensitivity curves of Advanced LIGO, its proposed A+upgrade, and LISA are shown in the same plane (legend). Right-hand panel: Same plot as in the left-hand panel, except that h_c (Y-axis) is multiplied with a correction factor to take into account the GW frequency evolution over a LISA mission lifetime of T_LISA = 5 yr.

Fig. 11 (right-hand panel) repeats the plot of the left-hand panel with

$$\begin{eqnarray*} h_c \times \min \left(\sqrt{\frac{\mathrm{ d}f_{\rm GWp}}{\mathrm{ d}t}\frac{T_{\rm LISA}}{f_{\rm GWp}}}, {\rm ~~}1 \right) \end{eqnarray*}$$

(7)

along the Y-axis. The multiplicative factor takes into account the lesser detectability in the LISA band the slower the evolution of f_GWp is (Sesana et al. 2005; Willems et al. 2007; Kremer et al. 2019) over a LISA mission lifetime of T_LISA = 5 yr. The factor diminishes the characteristic strain over the LISA and decihertz frequency bands; that over the LIGO-Virgo band remains unaffected. Figs 10 and 11 suggest that BBH mergers from YMCs and open clusters within the ∼100 Mpc Local Universe would be detectable by LISA (LIGO-Virgo) with S/N ratio ∼ few to 10s (∼10s–100s).

4.5 Caveats

There are several caveats in the current model computations that deserve improvements in the near future. An important improvement would be to implement a more consistent scheme for recycling BH spins due to matter accretion on to it. At present, such BHs are simply assigned a = 1 (see Section 2.6 and Appendix  A). Future models would also incorporate GW recoil kick and final spins in BBH mergers based on more recent NR studies such as those of Lousto & Zlochower (2013), Hofmann, Barausse & Rezzolla (2016). Although these improvements are unlikely to alter the main aspects of the present computations and their outcomes, they would, nevertheless, enrich nbody7 with even more updated and consistent prescriptions and treatments. The ongoing grid of N-body simulations, the current of which is presented here (Table C1), is somewhat heterogeneous, especially, in r_h(0), BH-spin model, and natal-kick models. They, nevertheless, encompass all the ingredients and options of the current model. In the near future, the set of computed models will be expanded, also taking into account external galactic field in dwarf-galaxy-like environments. The plethora of merger events from the recently concluded O3 of LVC will suggest which ingredients/options to focus on or whether further physical aspects need to be included.

5 SUMMARY AND OUTLOOK

In this section, the work presented in the previous sections are summarized and future prospects are indicated.

• This work presents a set of 65 long-term, direct, post-Newtonian many-body nbody7 computations of model dense stellar clusters with up-to-date stellar wind (B10 wind; Section 2.1), SN remnant-mass (F12-rapid and delayed SN; Section 2.2; Fig. 1), and SN natal-kick (momentum-conserving and collapse-asymmetry-driven, including slow down due to SN material fallback; Section 2.3; Fig. 2) prescriptions, including PPSN and PSN (B16 and weak PPSN/PSN prescriptions; Sections 2.2 and 2.7; Fig. 1). They also incorporate models of natal spins of stellar-remnant (i.e. first-generation) BHs (Section 2.4; Fig. 3), runtime tracking of NR-based GW recoils and final spins of BBH mergers (Section 2.5 and Appendix  A), and preliminary treatments of star–star mergers, star–remnant mergers, and spin-up of BHs due to matter accretion (Section 2.6, Appendices  A and  B). These ingredients and treatments allow BBH mergers inside clusters involving second-generation and mass-gained BHs. The model clusters, initially, have Plummer density profiles with masses |$1.0\times 10^4\, \mathrm{M}_{\odot }\le M_{cl}(0) \le 1.0\times 10^5\, \mathrm{M}_{\odot }$|⁠, half-mass radii 1.0 pc ≤ r_h(0) ≤ 3.0 pc, and ZAMS stars with masses ranging over |$0.08\!-\!150.0\, \mathrm{M}_{\odot }$| that are distributed according to the standard IMF (Section 3, Table C1). They range in metallicity over 0.0001 ≤ Z ≤ 0.02, initially unsegregated, are subjected to an external galactic field, and about half of them have their O-type stars paired among themselves with an observationally motivated distribution of primordial binaries (Section 3, Table C1). The structure and stellar content of such star cluster evolutionary models are consistent with those observed in YMCs and moderately massive OCs, that bridge low-mass OCs and GCs.

• These computed models produce GR mergers of BBHs that, primarily, take place while being bound to the cluster (i.e. are in-cluster mergers) and either being a part of a triple or higher order subsystem or by itself following a close encounter (Section 4.3; Fig. 8). The vast majority of the BBH mergers, in these models, take place following dynamical pairing, i.e. pairing among BHs whose parent stars were not members of the same primordial binary or of any binary (Section 4.1). The models also produce a few BNSs (that maintain the primordial-binary membership) that merge after getting ejected from their parent clusters, likewise for a few ejected BBH mergers (Section 4.1).

• These mergers, collectively, agree well with the observed masses, mass ratios, effective spin parameters, and final spins of the LVC O1/O2 merger events and also with the overall trends and 90 per cent confidence limits of these quantities in O1/O2, provided first-generation BHs are born with low or no spin (mesa or Fuller BH-spin model; Section 2.4) but spin-up after undergoing a (first-generation) BBH merger (Section 2.5 and Appendix  A) or matter accretion on to it (Section 2.6 and Appendices  A and  B); see Sections 4.1 and 4.2; Figs 4–7 and 9. In particular, the distinctly higher mass, effective spin parameter, and final spin of GW170729 merger event is naturally reproduced.

• The computed models also produce massive, |$M_{\rm tot}\sim 100\, \mathrm{M}_{\odot }$| BBH mergers with primaries within the ‘PSN gap’ (i.e. with |${M}_1\gt 45\, \mathrm{M}_{\odot }$|⁠; Section 4.1.1; Fig. 4). Irrespective of the BHs’ natal spins, such PSN-gap mergers would show spin signatures, i.e. typically have pro- or anti-aligned effective spin parameters leading to high or moderate final spins (Figs 5 and 6). If the SN natal kick is driven by collapse asymmetry (Section 2.3), the cluster models with the F12-delayed remnant scheme (Section 2.2) and metallicity Z ≥ 0.01 also yield GR mergers involving remnants (BHs) with masses within the (lower, NS–BH) ‘mass gap’ (Section 4.1.1; Fig. 4). Candidates of such ‘mass-gap’ mergers have been detected during the O3 of LVC. Furthermore, mass-asymmetric BBH mergers with distinctly low-mass ratio, similar to that of the LVC O3 event GW190412, is also naturally obtained in the present models (Section 4.1.2).

• The computed models also produce BBH inspirals that have eccentricities >0.1 in the LIGO-Virgo frequency band so that their eccentricity is detectable (Section 4.4; Fig. 10). In the present computations, ≈5 per cent of all GR mergers are such eccentric BBH mergers, which fraction is comparable to that estimated in recent studies of more massive GC models.

• These computations produce persistent and semipersistent GW sources detectable by LISA from within the ∼100 Mpc Local Universe (Section 4.4; Fig. 11).

• The computed evolutionary model set therefore shows that with state-of-the-art stellar-evolutionary and remnant-formation ingredients, YMCs and medium-mass OCs are able to produce dynamical GR mergers that are well consistent with the observed characteristics of the O1/O2 events. Like GCs and NSCs, such clusters are also capable of producing second-generation BBH mergers, PSN-gap mergers, mass-gap mergers, and LIGO-Virgo-detectable eccentric mergers. Such clusters also assemble highly mass-asymmetric mergers. The mergers from such clusters are distributed over a wide range of delay times (Figs 4 and 9) and they take place in clusters with a wide range of metallicity. However, this is compensated by the fact that such clusters form over a wide range of cosmic epochs and environments, so that it is, in principle, possible for all such events, from YMCs and OCs, to occur at the current epoch and be observable (Section 4.4; see also Kumamoto, Fujii & Tanikawa 2020).

The current model set and their physical ingredients would provide a variety of opportunities for comparing models with observations and make estimates and predictions. This study has focused on the GR-merger aspects but, of course, many other aspects are yet to be studied. A Monte Carlo approach is being designed to estimate the number of persistent GW sources and the detection rate of transient GW events from YMCs and OCs, based on the present (or a somewhat more expanded) model set, taking into account star formation history, cosmology, and instrument sensitivities (Banerjee 2020). Another important task is to explicitly compare the size, structure, and kinematics of the current model clusters, under various remnant-mass and natal-kick scenarios, with those in individual, well-observed YMCs and OCs, which will be taken up in the near future (see also Mackey et al. 2008; Paper I, II). Rooms for improvements on the current model ingredients are discussed in Section 4.5.

ACKNOWLEDGEMENTS

SB is thankful to the anonymous referee for constructive comments and useful suggestions that have helped to improve the manuscript. SB acknowledges the support from the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) through the individual research grant ‘The dynamics of stellar-mass black holes in dense stellar systems and their role in gravitational-wave generation’ (BA 4281/6-1; PI: S. Banerjee). SB acknowledges relevant discussions with Sverre Aarseth, Chris Belczynski, Chris Fryer, Mirek Giersz, Rainer Spurzem, Peter Berczik, Bhusan Kayastha, Sukanta Bose, Parameswaran Ajith, Deirdre Shoemaker, Pablo Laguna, Achamavedu Gopakumar, Sourav Chatterjee, Harald Pfeiffer, and Philipp Podsiadlowski. SB acknowledges the generous support and efficient system maintenance of the computing team at the AIfA and HISKP. SB acknowledges support from the Silk Road Project while visiting the National Astronomical Observatories of China (NAOC), Beijing, that has facilitated discussions and independent testings with nbody6+ +gpu (thanks to Peter Berczik). A few of the N-body computations presented here have been partially performed on the KEPLER cluster of the Astronomisches Rechen-Institut (ARI), Heidelberg. SB has solely performed, managed, and analysed all the N-body computations presented in this work. SB has done all the coding necessary for this work and has prepared the manuscript, except for a few output manoeuvres which are borrowed from NBODY6|${{+\,+}}$|GPU (thanks to Rainer Spurzem).

DATA AVAILABILITY

The simulation data underlying this article will be shared on reasonable request to the corresponding author. Some of the data are meant for further research and publications and may not be shared immediately.

Footnotes

1

Like in Papers I–III, ‘in-cluster coalescence’ in this work will refer to any GR coalescence taking place inside a cluster, while being gravitationally bound to the cluster.

2

In this work, BHs that are direct descendants of stars will be referred to as ‘first generation’, their merger product with other first-generation BHs as ‘second generation’, and so on.

3

The new implementations described in Sections 2.1–2.4 are also adopted in and independently tested to work correctly in a version of NBODY6|${{+\,+}}$|GPU (in preparation).

4

The outcome of a BH-TZO accretion can also be referred to as a ‘recycled’ BH but, in this work, we simply stick to referring to BH generations.

5

For a few models, the run had to be concluded prematurely as it could not be recovered after a crash or could be run only with impractically small regular time-steps, after arriving at its core-collapsed state (Spitzer 1987). This is, of course, not the case with all core-collapsed models.

6

An exception is the application of the Geneva BH-spin model with Z = 0.001 which would also make the most massive BHs non-spinning; see Fig. 3 (top panel; see also Morawski et al. 2018).

7

Until this condition is reached, |$\rm{\small ARCHAIN}$| continues to integrate the subsystem, applying PN treatment to the innermost binary. This seldom results in very long integration time of the subsystem. All |$\rm{\small ARCHAIN}$| PN sub-integration continues until the subsystem is resolved via either a GR merger of the innermost binary or the dynamical disintegration of the subsystem. Note that although spins are assigned to BHs at their formation and after experiencing matter interaction (Sections 2.4 and 2.6), BH spins are not used in the |$\rm{\small ARCHAIN}$| PN integration (both PN members are treated as non-spinning) and are used only for determining the GW recoil kick and the final spin (Section 2.5) upon COALESCENCE declaration. Since the vast majority of the final inspirals begin when the binary is of sub-Hz peak GW frequency (Fig. 10), spin PN terms would not have a significant effect on most of the sub-integrations. Although the speed of light can be ‘chosen’ in |$\rm{\small ARCHAIN}$|⁠, the physical speed of light has been used in all computations presented here.

8

In some works such as Rodriguez et al. (2018) and Kremer et al. (2019), it is only this type of mergers that have been labelled ‘in-cluster’.

9

For highly eccentric binaries, |$\tau _{\rm insp}^\prime$|⁠, as obtained by integrating Peters (1964) orbital-decay equations, can be up to twice of τ_insp, as obtained from equation (3).

APPENDIX A: IMPLEMENTATION OF BH SPIN AND BBH MERGER RECOIL IN UPDATED nbody7

A preliminary arrangement has been made to assign spins to BHs, at their birth from a stellar progenitor, based on stellar-evolutionary models, as described in Section 2.4. The Geneva, mesa, and Fuller BH spin models are implemented in subroutine Block/kick.f where the remnant natal kick (Section 2.3) is also evaluated and assigned (see Ba20). An input integer parameter bhflag toggles between the various BH-spin models (bhflag=2/3/4 for Geneva/mesa/Fuller model). The formulae representing the Geneva and mesa models are functions of carbon–oxygen core mass, M_CO, of the progenitor star (see B20) which value (along with fallback fraction, fallback mass, and ECS indicator; see Ba20) is imported from Block/hrdiag.f (the bse/nbody7 subroutine that assigns remnant mass; see Ba20) via a common block.

After computing the magnitude of the dimensionless spin parameter, a, the BH spin magnitude, S_BH, is evaluated according to equation (1). S_BH is then scaled to ‘N-body unit’, |$\hat{S}_{\rm BH}$|⁠, by dividing the former by the spin scaling parameter SPNFAC. |$\hat{S}_{\rm BH}$| and a are then copied to a newly added common block

COMMON /SPIN2/ SPN(NMAX), ASPN(NMAX) placed in the program-wide header file common6.h for easy accessibility from other subroutines. For maintaining the association of the spin values with the BH, the values are stored in the above arrays against the BH’s NAME, the latter being an integer that uniquely identifies a member. This is why the existing SPIN(NMAX) array is not used which is continued to be utilized independently and as before by bse-based routines that treat Roche lobe overflow, CE, and tidal circularization (except, in Block/kick.f, |$\hat{S}_{\rm BH}$| is also copied to SPIN(NMAX) against the newly formed remnant’s array index, once for all). At present, a fully consistent treatment of the recycling the BH’s natal spin, in the event of mass transfer on to the BH or a merger of the BH with a stellar member, is unavailable. Note that, at the start of a run, SPIN2 common block is initialized, once for all for all stars, in subroutine Block/instar.f. The newly formed BHs’ spins are then overwritten on to it from Block/kick.f, as described above.

In the event of a BH-star merger, a = 1 is simply set for the merger-product BH, as discussed in Section 2.6, which is done in the procedures in nbody7 for assembling a new, merged body. Specifically, the merged object (a BH) assumes the NAME of the merging stellar member, whose spin would typically correspond to a > 1 (due to the way stellar spins are assigned in instar.f following bse prescriptions; Hurley et al. 2002), making the BH maximally spinning when the BH’s spin parameter is later searched for against its NAME in ASPN array of SPIN2 common block. (In the NR subroutines described below, a = 1 is set if a > 1.) For the relatively rare case of BH formation during an interacting-binary phase (a symbiotic, mass transfer, or CE phase), a = 1 is set for the BH by a ‘serendipitous bug’ that prevents the newly formed BH, while being treated in subroutine Block/roche.f, to be accessed by Block/kick.f and hence prevents the reassignment of the BH’s spin in the latter routine (see above). Such a BH will also be found maximally spinning due to its progenitor star’s a > 1. (The BH, thereby, also receives zero natal kick but such a BH is typically in the direct collapse regime, i.e. it would get zero kick, anyway, even if it would have been formally treated in Block/kick.f; see Section 2.3.) It is checked in the computations presented in this work (Table C1) that the handful of BHs with a history of matter interaction with stars inevitably show up with a = 1, if they are encountered in the NR routines (see below).

During execution of an in-cluster GR merger in archain/chain.f, the NAMEs, masses, and stellar types of the merging members are collected and stored in a private (as opposed to code-wide) common block

COMMON/EXTRA3/ NMOBJ1,NMOBJ2,BOBJ1,BOBJ2,KOBJ1,KOBJ2.

This common block is shared with the ‘chain termination’ routines ARint/chterm.f and ARint/chterm2.f. It is in the latter two routines where the GR merger recoil kick is evaluated for BBH mergers, along with the final dimensionless spin parameter, and assigned to the merged BH (the computed recoil kick is vector-added to the instantaneous velocity of the BBH’s centre of mass and the final spin magnitudes are stored in the SPIN2 common block against the NAME of the merged BH).

This is done through calling a newly added routine GWREC3 which, in turn, calls the newly added routine GWKICK. Routine GWKICK serves as the ‘NR engine’ where the NR-based GR recoil (as in van Meter et al. 2010) and final spin (as in Rezzolla et al. 2008) formulae are implemented. Before calling GWKICK, GWREC3 makes use of EXTRA3 common block to get the names and masses of the merging bodies. The dimensionless spin magnitudes of the merging members are then obtained, against the names, from SPIN2 common block (see above). The spins are then assigned random orientations as described in Section 2.5. The masses, dimensionless spins, and spin orientations serve as arguments to GWKICK routine. GWREC3 also prints detailed information of the merging members and the merged object.

GWREC3 is also called from Block/brake4.f which routine treats GR inspiral and merger when these processes are not treated via |$\rm{\small ARCHAIN}$|⁠. In nbody7, brake4.f is never used for treating NS- or BH-containing binaries; it is used only whilst running in the ‘nbody6 mode’, that uses the classical or KS |${\tt CHAIN}$| (Mikkola & Aarseth 1993) instead of |$\rm{\small ARCHAIN}$|⁠.

Note that at the moment, only BH spins are assigned and tracked in the ways described above. The treatment naturally enables second-generation BHs to potentially be present and undergo second-generation GR mergers inside a cluster. In the following, two such examples of in-cluster, second-generation BBH mergers are presented from a computed model with |$M_{cl}(0)=5.0\times 10^4\, \mathrm{M}_{\odot }$|⁠, r_h(0) = 1.0 pc, f_bin(0) = 0, Z = 0.001, delayed+B16-PPSN/PSN remnant model, and Fuller BH spin model (model 49 of Table C1). Here, relevant text outputs from NBODY7|$\rm{\small ARCHAIN}$| are shown during BBH inspirals. The carryover of the final mass, NR-based final spin, and the identity of the first-generation BBH merger product (i.e. of the newly formed second-generation BH) to the second-generation BBH merger is highlighted with * symbols. The delay times (N-body unit; Heggie & Mathieu 1986) of the merger events are highlighted with square brackets.

In Example 1, the second-generation BBH merger happens well within the PSN mass gap, the second-generation BH involved being of |$78.3\, \mathrm{M}_{\odot }$|⁠. The final outcome of this BBH merger sequence is a BH of |$118.8\, \mathrm{M}_{\odot }$| and dimensionless spin 0.82 which crosses the cluster’s tidal radius in a few dynamical times (N-body times) and permanently escapes the system, due to the large recoil (>300 km s⁻¹) it receives during the second merger. This is an example of a massive, second-generation BBH merger early in the cluster’s evolution (first- and second-generation mergers at 522.2 and 1152.4 N-body times or 51.2 and 112.9 Myr). In this work, the ∼ few per cent mass-loss due to quadrupole GW energy radiation, which has a complex dependence on mass ratio and spin orientations (e.g. Sperhake 2015), is ignored and the mass of the merged BH is simply the sum of the merging BH masses.

In Example 2, such BBH merger sequence happens at late evolutionary times: the first-generation merger at t = 43476.8 (4260.7 Myr) and the second-generation merger at t = 103684.3 (10161.1 Myr). The second-generation merger happened within a compact triple system (as verified by the runtime data of in-cluster compact subsystems; see Banerjee 2018b), with an ECS-derived NS of mass |$1.26\, \mathrm{M}_{\odot }$| as the outer member. After the coalescence, an NSBH binary is formed (the NEW KSREG statement) which escapes the cluster after a few dynamical times, due to the GW recoil kick (the BINARY  escape statement).

Example 3 (from model 61 of Table C1) is an example for a merger between a mass-gained BH (a = 1) and a zero-spin BH (Fuller BH natal spin model; Section 2.4), leading to a GW170729-like final BH.

Example 1:

––––––––––––––––––––––––––––––––––––––––––––

INSPIRAL    T NP IPN E A TZ TKOZ DW  [522.247]   0   3  0.99971  3.12E−05  5.98E−01  1.00E+04  2.95E−03

REVERSE INFALL    IBH JBH I* MI MJ TZ    2   1  14  14  7.56E−04  8.10E−04  5.98E−01

COALESCENCE    T IC IPN NAM E EB PM TZ     522.2467591   0   3    *20286*  28908  0.99970502 −9.80E−  03  9.21E−09  5.98E−01

SWALLOWED STAR/BH    NAMC K* M1 M2  28908  14  37.78  *78.28*

INFALL CHECK:   T CG     522.25 −2.1E−22 −6.4E−23 −1.9E−22 −8.8E−16 −1.7E−16 −5.0E−16

CHAIN CHECK    ENERGY ENER0 DE ECH EnGR  −9.8047E−03  0.0000E+00 −9.8047E−03 −9.8033E−03 −1.4395E−06

END CHAIN    T # N NBH ECC SEMI RX ECH G   522.2468       90   1   1  0.99971  3.12E−05  5.81E−  06  0.00E+00  0.00E+00

BH1,2: m s theta phi a k

   3.78E+01   0.00E+00    271.977    240.010   0.000000   14

   4.05E+01   0.00E+00    132.642     24.660   0.000000   14

Merged BH: vprp1 vprp2 vpar xeff afin sfin thfin

     18.435255        0.000000        0.000000    0.000000    *0.686259*    1.21E+05       0.000

COALESCENCE KICK    VF ECDOT VCM VESC      1.24    0.000303   11.0     *13.6*

TERMINATE ARC    NNB SI SR BCM   133  2.98E−08  2.98E−08  1.57E−03

..............

..............

..............

INSPIRAL    T NP IPN E A TZ TKOZ DW [1152.357]   0   2  0.80386  1.61E−07  1.00E+00  1.00E+04  1.45E−03

REVERSE INFALL    IBH JBH I* MI MJ TZ    2   1  14  14  8.10E−04  1.57E−03  1.00E+00

COALESCENCE    T IC IPN NAM E EB PM TZ    1152.3570814   0   3    *20286*  32673  0.80381349   −3.95E+00  3.15E−08  1.00E+00

SWALLOWED STAR/BH    NAMC K* M1 M2  32673  14  40.50 *118.78*

INFALL CHECK:   T CG    1121.16  2.8E−24 −7.7E−25 −4.7E−24 −6.0E−15 −3.9E−17 −7.4E−16

CHAIN CHECK    ENERGY ENER0 DE ECH EnGR  −3.9466E+00  0.0000E+00 −3.9466E+00 −2.2870E−02 −3.9237E+00

END CHAIN    T # N NBH ECC SEMI RX ECH G  1152.3571   336794   1   1  0.80381  1.61E−07  1.98E−  07  0.00E+00  0.00E+00

BH1,2: m s theta phi a k

   4.05E+01    0.00E+00     176.869      40.068    0.000000    14

  *7.83E+01*   1.21E+05     337.076      47.201   *0.686259*   14

Merged BH: vprp1 vprp2 vpar xeff afin sfin thfin

     14.801988      −90.962949      361.319037    0.416548    0.818478    3.33E+05       8.155

COALESCENCE KICK    VF ECDOT VCM VESC     60.98    0.972884     6.0   *368.6*

TERMINATE ARC    NNB SI SR BCM   117  6.10E−05  1.22E−04  2.38E−03

..............

..............

..............

ESCAPE    N  =  84392  1265   0  0.7437  −12.665655   1.04  81.38  0.062   0.23  0.00001 84230 *20286*

––––––––––––––––––––––––––––––––––––––––––––

Example 2:

––––––––––––––––––––––––––––––––––––––––––––

INSPIRAL    T NP IPN E A TZ TKOZ DW  [43476.799] 0   2  0.99790  2.89E−06  1.00E+00  1.00E+04  1.56E−03

COALESCENCE    T IC IPN NAM E EB PM TZ   43476.7996768   0   3   *16743*   34163  0.99789827 −1.28E−  02  6.06E−09  1.00E+00

SWALLOWED STAR/BH    NAMC K* M1 M2  34163  14  13.78  *27.20*

INFALL CHECK:   T CG   43448.87 −2.6E−23 −2.3E−23 −3.9E−23 −3.5E−16 −7.7E−18 −1.5E−17

CHAIN CHECK    ENERGY ENER0 DE ECH EnGR  −1.2835E−02  0.0000E+00 −1.2835E−02 −1.7925E−04 −1.2655E−02

END CHAIN    T # N NBH ECC SEMI RX ECH G 43476.7997    66559   1   1  0.99790  2.88E−06  4.04E−  06  0.00E+00  0.00E+00

BH1,2: m s theta phi a k

   1.34E+01    0.00E+00     246.902      37.005    0.000000    14

   1.38E+01    0.00E+00     223.879       7.597    0.000000    14

Merged BH: vprp1 vprp2 vpar xeff afin sfin thfin

      6.809817       −0.000000        0.000000   −0.000000    *0.686827*    1.47E+04       0.000

COALESCENCE KICK    VF ECDOT VCM VESC      0.19   −0.000061     6.2     *1.2*

TERMINATE ARC    NNB SI SR BCM   255  2.44E−04  2.44E−04  5.44E−04

..............

..............

..............

ARC SWITCH    T TB−T S1 R DW G *************  4.37E−11  5.82E−11  2.83E−08  1.10E−03  3.85E−07

EINSTEIN SHIFT    # IPN IX E A DW       1   2  14  0.29411  2.82E−08  1.11E−03

RELATIVISTIC    ECC AX PM RZ TZ TPOM DW   0.2941  2.82E−08  1.99E−08  8.02E−12  1.61E+00  1.09E−  05  1.11E−03

REVERSE INFALL    IBH JBH I* MI MJ TZ    2   1  14  14  2.72E−04  5.44E−04  1.61E+00

COALESCENCE    T IC IPN NAM E EB PM TZ  [103684.3414642] 0   2   *16743*   63694  0.29410903   −2.63E+00  1.99E−08  1.61E+00

SWALLOWED STAR/BH    NAMC K* M1 M2  63694  14  13.60 *40.80*

INFALL CHECK:   T CG  103684.34  3.1E−26  2.1E−25 −1.1E−24  3.5E−15 −5.2E−16 −1.9E−16

CHAIN CHECK    ENERGY ENER0 DE ECH EnGR  −2.6260E+00  0.0000E+00 −2.6260E+00 −2.6260E+00  0.0000E+00

END CHAIN    T # N NBH ECC SEMI RX ECH G **********        1   1   1  0.29411  2.82E−08  3.63E−  08  0.00E+00  0.00E+00

BH1,2: m s theta phi a k

   1.36E+01    0.00E+00      35.571     330.895    0.000000    14

  *2.72E+01*   1.47E+04     111.447     287.338   *0.686827*   14

Merged BH: vprp1 vprp2 vpar xeff afin sfin thfin

     200.227384       35.752667     −374.437814   −0.167434    0.609463    2.93E+04      27.792

COALESCENCE KICK    VF ECDOT VCM VESC     27.77    0.415896    14.8   *411.4*

TERMINATE ARC    NNB SI SR BCM   160  1.86E−09  1.86E−09  8.16E−04

NEW KSREG   TIME[NB] 1.0368434146E+05 NM1,2,S=    *16743*   (63069)   102405 KW1,2,S =   14  13   0

IPAIR        1 DTAU  3.83E−03 M1,2[NB]  8.16E−04  2.52E−05 R12[NB]  3.89E−06

e,a,eb[NB]=   4.8387081E−01  2.6195E−06 −3.93E−03 P[d] =   3.30E+01 H −1.61E+02

GAMMA  0.00E+00 STEP(ICM)  2.33E−10 NPERT    0 NB(ICM)  159 M1,2[*] *4.08E+01* (1.26E+00)

RAD1,2,S[*]  5.69E−05  1.40E−05  2.23E+02 RI,VI[NB] =   1.91E−01  3.16E+01

..............

..............

..............

BINARY ESCAPE    KS =     1  NM =    *16743*  (63069) K* = 14 13  0 −1  M = *40.80* (1.26) EB =  −0.0039  R*/PM =  0.000  V/<V> =  28.89  E =  0.48387081  EI =   0.42009  P  =  3.3E+01

––––––––––––––––––––––––––––––––––––––––––––

Example 3:

––––––––––––––––––––––––––––––––––––––––––––

WATCH    # IPN E EN EGR A NAM K* M GP ES R     1190   2   0.9998418  –0.0043159  −0.0001085  4.2151E−  05      6313      6355  14  14    6.84E−04    5.32E−04    0.00E+00   −9.67E−05    7.96E−05

INSPIRAL    T NP IPN E A TZ TKOZ DW   1114.497   0   2  0.99984  4.22E−05  1.00E+00  1.00E+04  2.37E−03

REVERSE INFALL    IBH JBH I* MI MJ TZ    2   1  14  14  5.32E−04  6.84E−04  1.00E+00

COALESCENCE    T IC IPN NAM E EB PM TZ    1114.4967764   0   3     6313     6355  0.99984181 −4.32E−  03  6.67E−09  1.00E+00

SWALLOWED STAR/BH    NAMC K* M1 M2   6355  14  39.91  91.19

INFALL CHECK:   T CG    1114.47  1.2E−22  2.2E−21  1.2E−21  7.3E−18  9.2E−18  1.8E−17

CHAIN CHECK    ENERGY ENER0 DE ECH EnGR  −4.3161E−03  0.0000E+00 −4.3161E−03 −4.2074E−03 −1.0874E−04

END CHAIN    T # N NBH ECC SEMI RX ECH G  1114.4968     1193   1   1  0.99984  4.21E−05  8.22E  −05  0.00E+00  0.00E+00

BH1,2: m s theta phi a k

  *3.99E+01*   0.00E+00     209.119      73.663   *0.000000*   14

  *5.13E+01*   8.72E+04      31.900     120.558   *1.000000*   14

Merged BH: vprp1 vprp2 vpar xeff afin sfin thfin

    −113.815020     −125.059397      548.802952    0.477386   *0.849598*   2.04E+05      11.342

COALESCENCE KICK    VF ECDOT VCM VESC    283.34    1.617051     2.0   575.1

TERMINATE ARC    NNB SI SR BCM   194  7.63E−06  6.10E−05  1.22E−03

––––––––––––––––––––––––––––––––––––––––––––

APPENDIX B: MERGER MASS-LOSS, BH-TZO ACCRETION, AND, bse user inputs in updated nbody7

As discussed in Section 2.6, in star–star mergers within a cluster, a constant fraction, f_mrg, of the instantaneous secondary mass is eliminated. Also, in a BH-star merger (initially forming a BH-TZO object), a constant fraction, f_TZ, of the stellar mass is assumed to be accreted on to the BH.

The star–star merger mass-loss is implemented in the nbody7 routines Block/mix.f and Block/coal.f. BH-TZO accretion is also implemented in the same routines. The values of f_mrg and f_TZ are supplied to these routines via common blocks and are read from input at the start of the run (see below).

Apart from nbody7’s standard parameter and option inputs, all the bse input parameters as described in Ba20 (see also Hurley et al. 2002), along with f_mrg and f_TZ, are read from a special input file called input_bse. All the common blocks associated with these input parameters are placed in the code-wide header file common6.h for supplying their input-read values to the relevant bse routines. Care is taken that the values of all these parameters, as well as the arrays in SPIN2 common block (Appendix  A), are correctly saved in the ‘common dump’ file for resuming the runs.

APPENDIX C: N-BODY RUNS WITH UPDATED nbody7

Table C1 lists the long-term direct N-body computations performed in this work with updated nbody7 (Section 2). These runs are performed on server computers containing NVIDIAFermi-, Kepler-, or Turing-series GPUs and 4- or 16-thread CPUs. These calculations, altogether, took about 1 yr to complete including the times taken for continued code developments, their testing, the management of the runs, essential data analyses of the runs, and maintenance of the computer servers.

Figure C1.

Number of BHs, N_{BH, bound}, bound to the model clusters of Table C1 as a function of the models’ evolutionary time, t. The left-hand and the right-hand panel shows N_{BH, bound}(t) for the models with momentum-conserving and collapse-asymmetry-driven natal kicks, respectively. The lines are colour-coded according to the models’ initial mass, M_cl(0) (colour bar). The initial growth of BH population, due to the retention of BHs in the clusters at birth (Section 2.3), scales, overall, with M_cl(0). All evolutionary models exhibit long-term retention of BHs despite the decay of their population with time due to dynamical ejections. Since clusters with collapse-asymmetry-driven natal kick retain larger numbers of BHs at birth (Section 2.3), they retain larger numbers of BHs at the end of the computation (right-hand panel).

Table C1.

Summary of model evolutionary calculations. The columns from left to right give the model cluster’s (a) ID number, (b) initial mass, M_cl(0), (c) initial half-mass radius, r_h(0), (d) metallicity, Z, (e) initial fraction of primordial binaries, f_bin(0), (f) model evolutionary time, T_evol, (g) remnant-mass and PPSN/PSN model (Section 2.2), (h) remnant natal kick model (Section 2.3), (i) BH natal spin model (Section 2.4), (j) number of GR mergers within the cluster, N_{mrg, in}, (k) number of GR mergers after getting ejected from the cluster, N_{mrg, out}. f_bin(0) = 0.0 implies that the model initially contains only single stars. When f_bin(0) > 0, the quoted value represents the overall initial binary fraction, with the binary fraction for stars with ZAMS mass ≥M_crit being initially f_Obin = 1.0 as described in Section 3. Unless otherwise indicated, |${M}_{\rm crit}=16.0\, \mathrm{M}_{\odot }$|⁠. Initially, all stellar members in a cluster are ZAMS stars whose masses are distributed over |$0.08\!-\!150.0\, \mathrm{M}_{\odot }$| according to the standard IMF. All initial models follow a Plummer profile and are unsegregated. All model clusters are subjected to a solar-neighbourhood-like external galactic field. See the footnotes associated with this table for further details and specifications.

No.	Mcl(0) (M⊙)	rh(0) (pc)	Z	fbin(0)	Tevol (Gyr)	Remnant model	SN kick	BH spin	Nmrg, in	Nmrg, out
1	1.0 × 104	1.0	0.001	0.10a	8.5	Rapid+B16	mom. cons.b	Geneva	0	0
2	1.0 × 104	1.0	0.001	0.10	5.2	Rapid+B16	col. asym.	Geneva	0	1
3	2.0 × 104	2.0	0.001	0.10	11.0	Rapid+B16	mom. cons.	Geneva	1	0
4	2.0 × 104	2.0	0.001	0.10	8.7	Rapid+B16	col. asym.	Geneva	1	0
5	2.0 × 104	2.0	0.01	0.10	4.4	Rapid+B16	col. asym.	Geneva	0	0
6	2.0 × 104	2.0	0.01	0.10	5.8	Rapid+B16	mom. cons.	Geneva	1	0
7	2.0 × 104	2.0	0.02	0.10	4.4	Rapid+B16	mom. cons.	Geneva	0	0
8	3.0 × 104	1.0	0.0002	0.00	11.0	Rapid+B16	mom. cons.	Geneva	3	0
9	3.0 × 104	1.0	0.01	0.00	7.2	Rapid+B16	mom. cons.	Geneva	2	0
10	3.0 × 104	1.0	0.01	0.00	11.0	Rapid+B16	col. asym.	Geneva	1	1
11	3.0 × 104	1.0	0.01	0.00	10.9	Delayed+B16	col. asym.	mesa	4	0
12	3.0 × 104	1.0	0.02	0.00	11.0	Delayed+B16	col. asym.	mesa	0	0
13	3.0 × 104	1.0	0.02	0.00	8.2	Rapid+B16	mom. cons.	Geneva	0	0
14	3.0 × 104	1.0	0.02	0.00	7.0	Rapid+B16	col. asym.	Geneva	1	0
15	3.0 × 104	2.0	0.0002	0.00	11.0	Rapid+B16	mom. cons.	Geneva	1	0
16	3.0 × 104	2.0	0.01	0.00	9.7	Rapid+B16	mom. cons.	Geneva	2	0
17	3.0 × 104	2.0	0.01	0.00	11.0	Rapid+B16	col. asym.	Geneva	2	0
18	3.0 × 104	2.0	0.02	0.00	11.0	Rapid+B16	col. asym.	Geneva	0	0
19	3.0 × 104	2.0	0.001	0.10	11.0	Rapid+B16	mom. cons.	Geneva	2	0
20	3.0 × 104	2.0	0.001	0.10	11.0	Rapid+B16	col. asym.	Geneva	0	0
21	3.0 × 104	2.0	0.01	0.10	11.0	Rapid+B16	mom. cons.	Geneva	0	0
22	3.0 × 104	2.0	0.01	0.10	11.0	Rapid+B16	col. asym.	Geneva	1	0
23	3.0 × 104	2.0	0.02	0.10	3.0	Rapid+B16	mom. cons.	Geneva	0	0
24	3.0 × 104	2.0	0.001	0.10c	11.0	Rapid+weak	mom. cons.	mesa	1	2
25	3.0 × 104	2.0	0.001	0.10d	9.3	Rapid+weake	mom. cons.	mesa	0	1
26	3.0 × 104	3.0	0.0002	0.00	6.6	Rapid+B16	mom. cons.	Geneva	0	0
27	3.0 × 104	3.0	0.01	0.00	11.0	Rapid+B16	mom. cons.	Geneva	3	0
28	3.0 × 104	3.0	0.01	0.00	11.0	Rapid+B16	col. asym.	Geneva	0	0
29	3.0 × 104	3.0	0.02	0.00	10.1	Rapid+B16	mom. cons.	Geneva	2	0
30	3.0 × 104	3.0	0.02	0.00	11.0	Rapid+B16	col. asym.	Geneva	0	0
31	5.0 × 104	2.0	0.0002	0.00	11.0	Rapid+B16	mom. cons.	Geneva	0	1
32	5.0 × 104	2.0	0.0002	0.00	11.0	Rapid+weak	mom. cons.	mesa	1	0
33	5.0 × 104	2.0	0.001	0.00	11.0	Rapid+B16	mom. cons.	mesa	2	0
34	5.0 × 104	2.0	0.001	0.00	11.0	Rapid+weak	mom. cons.	mesa	0	0
35	5.0 × 104	2.0	0.001	0.00	11.0	Rapid+B16	col. asym.	mesa	2	0
36	5.0 × 104	2.0	0.005	0.00	11.0	Rapid+weak	mom. cons.	mesa	1	0
37	5.0 × 104	2.0	0.005	0.00	11.0	Delayed+B16	col. asym.	mesa	4	0
38	5.0 × 104	2.0	0.01	0.00	11.0	Rapid+B16	mom. cons.	Geneva	1	0
39	5.0 × 104	2.0	0.01	0.00	11.0	Rapid+B16	col. asym.	Geneva	3	0
40	5.0 × 104	2.0	0.02	0.00	9.9	Rapid+B16	mom. cons.	Geneva	0	0
41	5.0 × 104	2.0	0.02	0.00	11.0	Delayed+B16	col. asym.	mesa	3	0
42	5.0 × 104	2.0	0.0001	0.05	11.0	Rapid+B16	mom. cons.	Geneva	1	1
43	5.0 × 104	2.0	0.001	0.05	10.0	Rapid+B16	mom. cons.	Geneva	4	2
44	5.0 × 104	2.0	0.001	0.05f	11.0	Rapid+weak	mom. cons.	mesa	1	2
45	5.0 × 104	2.0	0.0001	0.05g	11.0	Rapid+B16	mom. cons.	mesa	2	2
46	5.0 × 104	2.0	0.01	0.05h	11.0	Delayed+B16	col. asym.	mesa	2	0
47	5.0 × 104	1.0	0.001	0.00	11.0	Rapid+B16	mom. cons.	Geneva	1	2
48	5.0 × 104	1.0	0.001	0.00	11.0	Delayed+B16	col. asym.	mesa	3	0
49	5.0 × 104	1.0	0.001	0.00	11.0	Delayed+B16	mom. cons.	Fuller	5	0
50	5.0 × 104	1.0	0.01	0.00	11.0	Delayed+B16	col. asym.	mesa	7	0
51	5.0 × 104	1.0	0.02	0.00	11.0	Delayed+B16	col. asym.	mesa	4	4
52	7.5 × 104	2.0	0.001	0.00	11.0	Rapid+B16	mom. cons.	Geneva	2	0
53	7.5 × 104	2.0	0.001	0.00	11.0	Rapid+weak	mom. cons.	mesa	3	0
54	7.5 × 104	2.0	0.001	0.00	11.0	Rapid+B16	col. asym.	Geneva	3	0
55	7.5 × 104	2.0	0.005	0.00	11.0	Delayed+B16	col. asym.	mesa	4	0
56	7.5 × 104	2.0	0.01	0.00	11.0	Rapid+B16	mom. cons.	Geneva	6	0
57	7.5 × 104	2.0	0.02	0.00	11.0	Rapid+B16	mom. cons.	Geneva	6	0
58	7.5 × 104	2.0	0.02	0.00	11.0	Delayed+B16	col. asym.	mesa	3	1
59	7.5 × 104	2.0	0.0001	0.05i	11.0	Rapid+B16	mom. cons.	mesa	1	2
60	7.5 × 104	2.0	0.001	0.05j	11.0	Rapid+B16	mom. cons.	mesa	1	4
61	7.5 × 104	2.0	0.001	0.05k	9.8	Rapid+B16	mom. cons.	Fuller	5	1
62	7.5 × 104	2.0	0.01	0.05l	11.0	Rapid+B16	mon. cons.	Fuller	1	1
63	7.5 × 104	2.0	0.02	0.05m	11.0	Delayed+B16	col. asym.	mesa	2	0
64	1.0 × 105	2.0	0.001	0.00	11.0	Delayed+B16	col. asym.	Geneva	5	1
65	1.0 × 105	1.5	0.001	0.05n	1.0o	Rapid+B16	mom. cons.	Fuller	4	0

^aThe binary fraction is defined as f_bin = 2N_bin/N, N_bin being the total number of binaries and N being the total number of members.

^bmom. cons. = momentum conserving natal kick model, col. asym. = collapse-asymmetry-driven natal kick model.

^c|${M}_{\rm crit}=5.0\, \mathrm{M}_{\odot }$|⁠.

^d|${M}_{\rm crit}=5.0\, \mathrm{M}_{\odot }$|⁠.

^e30 per cent of the full B10 wind is applied.

^ff_TZ = 0.70, f_mrg = 0.3.

^gf_mrg = 0.2.

^hf_TZ = 0.90, f_mrg = 0.2.

ⁱf_mrg = 0.2.

^jf_TZ = 0.95, f_mrg = 0.2.

^kf_TZ = 0.95, f_mrg = 0.2.

^lf_TZ = 0.95, f_mrg = 0.2.

^mf_mrg = 0.2.

ⁿf_TZ = 0.95, f_mrg = 0.2.

^oOngoing run.