Issue 
A&A
Volume 665, September 2022



Article Number  A20  
Number of page(s)  13  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202142331  
Published online  02 September 2022 
Binary black hole mergers from young massive clusters in the pairinstability supernova mass gap
^{1}
HelmholtzInstituts für Strahlen und Kernphysik, Nussallee 1416, 53115 Bonn, Germany
^{2}
ArgelanderInstitut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
email: sambaran@astro.unibonn.de
Received:
29
September
2021
Accepted:
14
June
2022
Context. The recent discovery of the binary black hole (BBH) merger event GW190521, between two black holes (BHs) of ≈100 M_{samp}, in addition to other massive BBH merger events involving BHs within the pairinstability supernova (PSN) mass gap have sparked widespread debate on the origin of such extreme gravitationalwave (GW) events. GW190521 simultaneously triggers two critical questions: how BHs can appear within the ‘forbidden’ PSN gap and, if they do, how they get to participate in generalrelativistic (GR) mergers.
Aims. In this study, I investigate whether dynamical interactions in young massive clusters (YMCs) serve as a viable scenario for assembling PSNgap BBH mergers.
Methods. To that end, I explore a grid of 40 new evolutionary models of a representative YMC of initial mass and size M_{cl} = 7.5 × 10^{4} M_{samp} (N ≈ 1.28 × 10^{5}) and r_{h} = 2 pc, respectively. The model grid ranges over metallicity 0.0002 ≤ Z ≤ 0.02 and comprises initial cluster configurations of King central concentration parameters W_{0} = 7 and 9. In each model, all BH progenitor stars are initially in primordial binaries following observationally motivated distributions. All cluster models are evolved with the direct, relativistic Nbody code NBODY7, incorporating uptodate remnant formation, BH natal spin, and GR merger recoil schemes.
Results. Binary black hole mergers from these model cluster computations agree well with the masses and effective spin parameters, χ_{eff}, of the events from the latest gravitationalwave transient catalogue (GWTC). In particular, GW190521like, that is to say ≈200 M_{samp}, low χ_{eff} events are produced via a dynamical merger among BHs derived from starstar merger products. GW190403_051519like, that is PSNgap, highly asymmetric, high χ_{eff} events result from mergers involving BHs that are spun up via matter accretion or a binary interaction. The resulting presentday, differential intrinsic merger rate density, within the PSN gap, accommodates that from GWTC well.
Conclusions. This study demonstrates that, subject to model uncertainties, the tandem of massive binary evolution and dynamical interactions in ≲100 Myrold, low metallicity YMCs in the Universe can plausibly produce GR mergers involving PSNgap BHs and in rates consistent with that from uptodate GW observations. Such clusters can produce extreme events similar to GW190521 and GW190403_051519. The upper limit of the models’ GW190521type event rate is within the corresponding LIGOVirgoKAGRA (LVK)estimated rate limits, although the typical model rate lies below LVK’s lower limit. The present YMC models yield a merger rate density of 0−3.8 × 10^{−2} yr^{−1} Gpc^{−3} for GW190521type events. They produce GW190403_051519like events at a rate within 0−1.6 × 10^{−1} yr^{−1} Gpc^{−3} and their total BBHmerger yield within the PSN gap is 0−8.4 × 10^{−1} yr^{−1} Gpc^{−3}.
Key words: stars: black holes / stars: massive / stars: kinematics and dynamics / supernovae: general / methods: numerical / gravitational waves
© S. Banerjee 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Within just a few years after the first detection (Abbott et al. 2016) of a gravitational wave (hereafter GW) signal from a merging binary black hole (hereafter BBH), we are already approaching a ‘golden era’ of GW and multimessenger astronomy (Branchesi 2016; Mapelli 2018; Mészáros et al. 2019; Mandel & Broekgaarden 2022). Recently, the LIGOVirgoKAGRA collaboration (hereafter LVK; Aasi et al. 2015; Acernese et al. 2015; KAGRA Collaboration 2020) has published, in their GW transient catalogue (hereafter GWTC)^{1}, about 90 candidates of general relativistic (hereafter GR) compact binary merger events. The uptodate GWTC includes all event candidates from LVK’s first, second (‘O1’, ‘O2’; Abbott et al. 2019), and third (‘O3’; Abbott et al. 2021a; The LIGO Scientific Collaboration 2021c) observing runs, including those in their Deep Extended Catalogue (The LIGO Scientific Collaboration 2021a).
It is important to note that GWTC contains several events that have triggered widespread debate regarding the events’ origin or formation channel. Perhaps the most explored event since its announcement has been GW190521 (Abbott et al. 2020): a merger of two black holes (hereafter BHs) of estimated masses and and effective aligned spin parameter (Abbott et al. 2021a; the limits correspond to 90% credibility). The χ_{eff} (Ajith et al. 2011) is a measure of the spinorbit alignment of a merging binary and is defined as
Here, the GRinspiralling masses M_{1}, M_{2}, with mass ratio q ≡ M_{2}/M_{1}, have Kerr vectors a_{1}, a_{2} that project with angles θ_{1}, θ_{2} on the orbital angular momentum vector just before the merger, respectively. The Kerr vector or Kerr parameter (also referred to as the dimensionless spin vector or dimensionless spin parameter), a, is defined as (Kerr 1963)
where S_{BH} is the total angular momentum vector of a Kerr BH of (nonspinning) mass M_{BH}.
Another notable PSNgap event is GW190403_051519 (hereafter GW190403), which is lighter than GW190521 but – in contrast – is highly massasymmetric () and spinorbit aligned (). GW190426_190642 has masses and χ_{eff} similar to GW190521, although this event is of a much higher false alarm rate (hereafter FAR), 4.1 yr^{−1}, as opposed to FAR = 2 × 10^{−4} yr^{−1} for GW190521. The current GWTC includes several additional candidates of BBH mergers with the primary being within the PSN gap.
The excitement is natural since both components of this BBH merger lie well within the ‘forbidden’ pairinstability supernova (hereafter PSN) mass gap between 45 M_{samp} ≲ M_{rem} ≲ 120 M_{samp} (Belczynski et al. 2016; Woosley 2017). Starting to evolve from their zero age main sequence (hereafter ZAMS), massive single stars are not expected to produce any compact remnant over this remnantmass range and, instead, explode completely due to PSN occurring in their late evolutionary stages (Langer et al. 2007; Woosley 2017; Mapelli et al. 2020). The lower limit of the PSN gap occurs as a result of pulsation pairinstability supernovae (hereafter PPSN). A PPSN episodically sheds the hydrogen envelope of a parent star until a helium core of ≈45 M_{samp} is retained (Woosley 2017; Woosley et al. 2020), which then directly collapses into a BH. The upper limit of the gap results when the progenitor star becomes massive enough so that PPSN and PSN quenches (see, e.g. Ziegler & Freese 2021) and the evolved star can directly collapse into an intermediate mass BH (hereafter IMBH). Readers can refer to Fig. 2 of Banerjee et al. (2020) for the ‘standard’ ZAMS massfinal mass relation, exhibiting the PSN gap at different metallicities.
Therefore, the observation of BHs within the PSN gap, at least at first glance, is a signature of an outlier BH population (Baibhav et al. 2021; O’Brien et al. 2021) and creates interest in more exotic scenarios of formation of BHs and mergers involving them. Several scenarios have already been invoked to explain GW190521type events: primordial BHs (e.g. Clesse & GarcíaBellido 2017; Carr & Silk 2018; De Luca et al. 2019), BHs derived from Population III stars (e.g. Tanikawa et al. 2021a,b; Ziegler & Freese 2021), gas accretion onto BHs in dense protoclusters (e.g. Roupas & Kazanas 2019), GR coalescences in field hierarchical systems (e.g. Fragione et al. 2020; VignaGómez et al. 2021; Hamers et al. 2021), and stellar collisions and dynamical interactions inside dense stellar clusters.
In this study, a young massive cluster (hereafter YMC) origin of PSNgap BBH mergers is investigated. Dynamical interactions among stellarremnant BHs inside dense star clusters is an intriguing scenario for generating PSNgap BBH mergers as well as other types of massive BBH mergers (Perna et al. 2019; Arca Sedda et al. 2020; Baibhav et al. 2020; Banerjee 2021a; Rizzuto et al. 2021, 2022; Gerosa & Fishbach 2021; Mapelli et al. 2021). This is because the scenario naturally allows for the prolonged postprocessing of BHs and BBHs that formed as a result of massive binary evolution. Dynamical interactions also enable kinds of starstar and starremnant mergers that would not be possible with the evolution of isolated massive binaries alone (Di Carlo et al. 2020a; González et al. 2021; Rastello et al. 2021). Such interactions would potentially lead to ‘forbidden’ BHs and mergers involving them depending on a cluster’s properties and its stellar and multiplicity content. Since the most massive members in a star cluster undergo dynamical processing the earliest, the most massive BBH mergers are expected to occur early in the cluster’s evolution (Banerjee 2017). As shown in earlier works (Banerjee et al. 2010; Banerjee 2017; Rodriguez et al. 2018), dynamical BBH mergers in a YMC can occur inside a cluster due to hierarchicalsystem interactions (e.g. Kozai–Lidov oscillation, chaotic triple interaction; Kozai 1962; Katz et al. 2011; Lithwick & Naoz 2011) or close flyby interactions, within the centrally segregated BH subsystem (Banerjee et al. 2010; Morscher et al. 2015). A fraction of the BBHs (and other types of compact binaries) ejected from the cluster, either dynamically or due to natal kick, can also merge within the Hubble time (Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000).
It should, however, be borne in mind that it is possible that the ‘standard’ PSN mass gap (see above) is simply an artefact of the ‘standard’ theoretical stellar evolution models. The PPSN mass cap can potentially be much higher or the PSN mass gap can even be nonexistent. This is supported by both stellar evolutionary models (Belczynski et al. 2020a; Farmer et al. 2020; Woosley & Heger 2021; Vink et al. 2021; Costa et al. 2021) and the observed population of merging BHs (The LIGO Scientific Collaboration 2021b; Edelman et al. 2021). In that case, GW190521like events can simply be explained by the standard isolated binary evolution scenario (Belczynski 2020). The present study is based on the existence of the standard PSN gap and involves a stellarremnant model that exhibits such a gap.
This paper is organized as follows. Section 2 describes a new set of direct Nbody computations of model YMCs. Section 3 discusses the outcomes, focussing on PSNgap mergers. Section 3.1 describes the inferred rate of GW190521like events and Sect. 3.2 covers that of GW190403like mergers. Section 4 summarizes the findings and discusses the implications of various model assumptions and the limitations.
2. Computations
Close dynamical interactions, starstar mergers, and starremnant mergers are generally favoured in massive star clusters with profiles of high central concentration. This is due to the resulting high central stellar density and efficient mass segregation. Therefore, massive BH formation and their mergers can be expected to preferably occur in such systems. Since massive, PSNgap BH formation and their involvement in GR mergers is the focus of this study, massive, concentrated model clusters are specifically considered. In this work, model star clusters with an initial mass of M_{cl} = 7.5 × 10^{4} M_{samp} (initial number of stars N ≈ 1.28 × 10^{5}) are considered to be representatives of YMCs. Such a cluster mass is comparable to those of most massive Galactic and localgroup YMCs. It is also representative of moderate mass ‘super star clusters’ (Portegies Zwart et al. 2010). The stellar content of these star cluster models is motivated by the cluster models described in Banerjee (2021a, hereafter Ba21). The initial density and kinematic profiles of the clusters follow the King (1966) model with a halfmass radius of r_{h} = 2 pc and a King dimensionless potential of W_{0} = 7 and 9. Such size is typical for YMCs (Portegies Zwart et al. 2010; Krumholz et al. 2019).
Previous results in Ba21 and Banerjee (2021b), where the models use the Plummer (1911) initial conditions, that is to say with the central concentration comparable to W_{0} ≈ 5−6 (Heggie & Hut 2003), suggest that massive clusters (M_{cl} ≥ 5 × 10^{4} M_{samp}) are capable of producing and merging PSNgap BHs. An objective of the present study is to see how somewhat more centrally peaked initial models fair in this, which is why W_{0} = 7 and 9 are chosen here (a less concentrated model would not be interesting in this regard due to the reasons given above). The King profile offers the freedom to alter the central concentration without altering the cluster’s overall length scale (as measured by r_{h}).
Another objective is to see how parsecscale clusters fair as opposed to highly concentrated, subparsecscale clusters of r_{h} comparable to widths of molecularcloud filaments, the latter conditions having been recently studied by several authors (Di Carlo et al. 2020a; Rastello et al. 2021; Rizzuto et al. 2022). (Near) monolithic, (near) gasfree, parsecscale YMCs are something that we do observe and can measure with relative unambiguity in various surveys (Mathieu 2000; Portegies Zwart et al. 2010; Kuhn et al. 2013; Krumholz et al. 2019). The gasfree nature of such clusters implies that any potential violentrelaxation phase and its impact on the cluster (Banerjee & Kroupa 2018) can be ignored and the cluster can be assumed to evolve solely via secular dynamical evolution, as is the case with the present model clusters.
As in Ba21, the initial models are composed of ZAMS stars of masses 0.08 M_{samp} ≤ m_{*} ≤ 150.0 M_{samp} that are distributed according to the canonical initial mass function (hereafter IMF; Kroupa 2001), and they have an overall (see below) primordialbinary fraction of f_{bin} = 5%. However, as in Ba21, the initial binary fraction of the Otype stars (m_{*} ≥ 16.0 M_{samp}), which are initially paired only among themselves, is considered to be f_{Obin}(0) = 100%, which is consistent with the observed high binary fraction among Ostars in young clusters and associations (e.g. Sana & Evans 2011; Sana et al. 2013; Moe & Di Stefano 2017). The Ostar binaries are considered to initially follow the observed orbitalperiod distribution of Sana & Evans (2011) and a uniform massratio distribution. The initial orbital periods of the nonOstar primordial binaries follow the period distribution of Duquennoy & Mayor (1991) and their massratio distribution is also uniform. The initial eccentricity of the Ostar binaries follows the Sana & Evans (2011) eccentricity distribution and that for the rest of the binaries obeys the thermal eccentricity distribution (Spitzer 1987). As explained in Banerjee (2018a), such a scheme for including primordial binaries provides a reasonable compromise between the economy of computing and consistencies with observations.
The model clusters are evolved with the starbystar, direct Nbody evolution code NBODY7 (Aarseth 2012), which has been updated in several aspects as detailed in Banerjee et al. (2020, hereafter Ba20) and Ba21. These updates mainly enable uptodate stellar wind mass loss and remnant formation in the code and also implement numerical relativity (hereafter NR)based GR merger recoil and spin recycling of merged BHs. The latter allows for onthefly and consistent treatment of BBH merger products. The primary ‘engine’ for stellar and binary evolution in NBODY7 is BSE (Hurley et al. 2000, 2002) and that for postNewtonian (hereafter PN) evolution of compact binaries and higher order systems is ARCHAIN (Mikkola & Tanikawa 1999; Mikkola & Merritt 2008). In the present computed models, the ‘F12rapid+B16PPSN/PSN’ (see Ba20) remnantmass prescription is applied. For single star evolution, this remnant prescription allows for the formation of stellar remnants, of mass M_{rem}, maintaining a neutron star (NS)BH mass gap between 2 M_{samp} ≲ M_{rem} ≲ 5 M_{samp} (Fryer et al. 2012) and a PSN mass gap between 45 M_{samp} ≲ M_{rem} ≲ 120 M_{samp} (Belczynski et al. 2016; Woosley 2017; see also Fig. 2 of Ba20).
In the present computations, a zero Kerr parameter is assigned to all BHs derived from single stars or from members of nonmasstransferring or noninteracting binaries, as per the BHformation models of Fuller & Ma (2019). This scheme is hereafter referred to as the FM19 BH natal spin model. If, after formation, a BH undergoes matter accretion due to a BHstar merger or mass transfer in a binary or if a BH is born in a tidally interacting (or symbiotic) binary, its Kerr parameter is set to the maximally spinning value of a = 1. In the event of a BHstar merger (the formation of a BH Thorne–Zytkow object, Thorne & Zytkow 1975 or BHTZO), the f_{TZ} = 0.95 fraction of the merging star’s mass is assumed to be accreted onto the BH. Also, in starstar mergers, the f_{mrg} = 0.2 fraction of the secondary’s mass (≤10% of the total merging stellar mass) is assumed to be lost in the merger process (Gaburov et al. 2008; Glebbeek et al. 2009). In the present runs, the NR treatment is updated to the more recent GWrecoil formulae of Lousto et al. (2012) and finalspin formulae of Hofmann et al. (2016).
For each W_{0} of the initial King profile, five metallicities are taken, namely, Z = 0.0002, 0.001, 0.005, 0.01, and 0.02. Four random model realizations of each (W_{0}, Z) pair, which are subjected to a solarneighbourhoodlike external field^{2}, are evolved for 300 Myr. These 40 newly computed models are listed in Table A.1.
3. Results
The bottom panels of Fig. 1 show the primary mass, M_{1}, versus secondary mass, M_{2} (M_{1} ≥ M_{2}), of the compactbinary mergers from the computed models. All these events are incluster or ejected BBH mergers, except for one that is an NSBH merger. The BBHmerger masses are, overall, well consistent with those from GWTC (black, filled squares with error bars; events with FAR ≤1 yr^{−1} and > 1 yr^{−1} are distinguished by different shades). The NSBH merger is also consistent with the LVKdetected NSBH merger events (Abbott et al. 2021b), as seen near the bottom left corners of the bottom panels in Fig. 1. For comparison, analogous plots are obtained from the computed models of Ba21 (Fig. 1; top panels), whose mergers show similarly good agreement with the LVK events.
Fig. 1. Comparison of BBH merger masses from computed model star clusters with those from GWTC. Top panels: primary mass, M_{1}, versus secondary mass, M_{2} (M_{1} ≥ M_{2}), of all BBH mergers from the computed model star clusters of Ba21 (filled circles). These data points are colourcoded according to the merger delay time, t_{mrg} (left panel), and the parent model cluster’s metallicity, Z (right panel). The GWTC events’ data points (black, filled squares) and their corresponding 90% credible intervals (horizontal and vertical error bars) are indicated where events with FAR ≤1 yr^{−1} and > 1 yr^{−1} are distinguished by different shades. The shaded regions over 45 M_{samp} − 120 M_{samp} on each panel represent the PSN mass gap. The shadings have been placed in the background so that they do not influence the colours of the data points. Bottom panels: same as the top panels, except that the mergers from the present computed clusters (Table A.1) with King concentration parameters W_{0} = 7.0 (filled triangles) and 9.0 (filled circles) are shown. 
The presently computed cluster models are more efficient in producing PSNgap BBH mergers: the 40 models produce a similar number of mergers with M_{1} and/or M_{2} within the gap as the 65 models of Ba21. This is due to the more centrally peaked initial profiles of the present models (Sect. 2) than those of Ba21 (despite similar r_{h} for both model sets). However, in both sets, mergers involving a PSNgap BH happen only in models with Z ≤ 0.005 (Fig. 1). In the present models, BHs in the PSN gap appear due to (i) starstar mergers and (ii) BHTZO accretion. Depending on the stellarevolutionary age of the merging stars, a starstar merger can result in an overmassive Henvelope of the merged star, with its Hecore mass being below the PPSN–PSN threshold. With sufficient Henvelope mass, such a merged star would evolve into a directcollapse BH within the PSNgap (e.g. Spera et al. 2019; Banerjee et al. 2020; Di Carlo et al. 2020b). Such a remnant mass range is ‘forbidden’ for single stars evolving directly from ZAMS. The high, 95% BHTZO accretion, as adopted here, would also push BHs into the PSN gap as a result of sufficiently massive BHstar mergers.
In the Ba21 models, a third channel has also produced PSNgap BBH mergers in the model runs with FM19 (vanishing) BH natal spins, namely, second generation (hereafter 2G) mergers^{3}. In those models, most 2G BBH mergers have occurred with delay times of t_{mrg} > 300 Myr. In the present models, although a few 2G BHs remain in the clusters until the end of the runs at 300 Myr evolutionary time (with the majority of the 2G BHs being ejected from their host clusters at their formation due to GW recoil kick; see Sect. 2), none of them could get involved in another merger by this time, as can be expected based on the Ba21 models.
The most massive, GW190521like merger obtained from the present model set (Fig. 1; bottom panels) is a 1G, incluster merger between two PSNgap BHs, both of which are derived from starstar merger products. The merger happened in one of the Z = 0.0002, W_{0} = 9 models. Figure 2 plots the mass distributions of stellar remnants retained within the cluster at different evolutionary times for all the 40 computed models. To improve counts over the mass bins and also facilitate comprehension, distributions – at a given cluster age – from all four models with a given W_{0} and Z are combined in one histogram. The figure demonstrates the overall trend that the maximum mass of the BHs retained in a cluster at a given age increases with decreasing Z and increasing W_{0}, as can be expected.
Fig. 2. Mass distributions of stellar remnants remaining in the computed model clusters (Table A.1) at 5 Myr, 10 Myr, 150 Myr, and 300 Myr clusterevolutionary times (legend). The distributions, at a given time, from all (four) models for a particular W_{0} and Z are combined in one histogram. The panels in the left (right) column correspond to the models with W_{0} = 7 (W_{0} = 9) with Z as indicated in each panel’s title. 
Figure 3 plots the effective (aligned) spin parameter, χ_{eff}, of the GR mergers from the present models against their chirp mass, M_{chirp} (top panels), and primary mass (bottom panels). A merging system’s chirp mass is defined as
Fig. 3. Effective aligned spin parameter, χ_{eff}, versus chirp mass, M_{chirp} (top panels), and primary mass, M_{1} (bottom panels), of all BBH mergers from the computed model clusters of Table A.1. For each merger, χ_{eff} values for a number of random orientations of the merging BHs’ spins are plotted (Sect. 3). The legends are the same as in Fig. 1, except that a black edge has been applied to the filled triangle symbol for improved visibility. 
For every GR merger from the computed models, a range of potential values of χ_{eff} is plotted in Fig. 3 by assigning independent, random values to θ_{1} and θ_{2} in Eq. (1). If the components of the merging binary are derived from stellar progenitors that were members of different primordial binaries or if they were single stars initially, that is to say if the merging binary is dynamically assembled, then the merging members are uncorrelated. Accordingly, (θ_{1}, θ_{2}) are considered to be isotropically oriented and values over the range 0 ≤ (θ_{1}, θ_{2})≤2π are assigned to them. On the other hand, if the merging components are derived from stellar progenitors that were members of the same primordial binary, then they would, plausibly, be at least partially correlated. Therefore, in this case, values over the range 0 ≤ (θ_{1}, θ_{2})≤π/2 are assigned. For every merger from the models, the corresponding values of M_{1}, M_{2}, a_{1}, and a_{2} are taken directly from the model^{4}.
For comparison, event data from GWTC are also plotted in Fig. 3. The LVK data points do exhibit an overall bias towards positive χ_{eff}, at least when the most probable values (squares) are considered. This is indicative of a contribution from additional channels which tend to produce spinaligned mergers, such as in isolated binary evolution. Mixing with additional channels will be addressed in a future work (see Banerjee 2022 for more details). Nevertheless, considering the 90% credible limits (error bars), the computed data points are consistent with the LVK data. In particular, the computed data well encompass the most massive observed PSNgap mergers (of M_{1} ≥ 55 M_{samp}) with nearly vanishing, mildly aligned, and highly aligned χ_{eff} values.
Figure 4 replots the model χ_{eff} as mean values with error bars, against M_{1} (upper panels) and M_{2} (lower panels), and it highlights the events GW190521 and GW190403. Figure 4 shows that the present model clusters produce only one GW190521like and five GW190403like events, in the sense that the model mergers’ M_{1}, M_{2}, and χ_{eff} all lie within or enter the 90% credible intervals (the LVKdata error bars) of the respective parameters of these events of interest.
Fig. 4. Effective aligned spin parameter, χ_{eff}, versus primary mass, M_{1} (top panels), and secondary mass, M_{2} (bottom panels), of all BBH mergers from the computed model clusters of Table A.1. For each model merger, χ_{eff} values for a number of random orientations of the merging BHs’ spins are obtained (Sect. 3), which are shown as a data point (mean value) with a vertical error bar (range). The legends of the data points are the same as in Fig. 1. The error bars (in a dashed line and solid line for mergers from the W_{0} = 7 and W_{0} = 9 models, respectively) share the same colour code as the corresponding data points. The LVK events GW190521 and GW190403 are highlighted with two empty squares, which are coloured differently for improved legibility and they do not correspond to the colour code. 
3.1. Rate of GW190521like events
The present computed cluster models yield BBH mergers that nicely resemble notable PSNgap events discovered by LVK, in terms of both merging masses (M_{1}, M_{2}) and spinorbit alignment (χ_{eff}). This is evident from Figs. 1, 3, and 4. The GW190521like merger event is an outcome of an incluster, dynamically assembled and triggered merger between two BHs of ≈110 M_{samp} and 57 M_{samp}, in one of the Z = 0.0002, W_{0} = 9 models. Both of these BHs are formed due to evolution of starstar merger products with an overmassive Henvelope (see above): both of these stellar mergers happened between components of massive primordial binaries. Since the BHs are derived directly from the (merged) stellar progenitors, they have zero natal spins as per the adopted BH natal spin model (Sect. 2). After forming as single BHs, they eventually get paired up via exchange interactions and merge due to a dynamical triple interaction (Banerjee 2018b), resulting in an event with masses and χ_{eff} well within the 90%credible limits of those of GW190521 and GW190426_190642. The detailed timeline of this event, from the formation of the merging BHs up to their merger, is depicted in Fig. 5 based on data directly from the Nbody simulation and is described further in the figure’s caption.
Fig. 5. Timeline of the GW190521like merger event that occurred inside one of the computed cluster models with W_{0} = 9, Z = 0.0002. In both panels the evolutions for the three cluster members, which are most relevant for the event, are shown from the start of the computation until the merger. Top panel: mass evolution of these objects: the trace of the merger primary and secondary mass is indicated by a circle and a square, respectively, and that of the third ‘perturber’ is represented by a triangle. The symbols are colourcoded according to the stellarevolutionary stage (main sequence, MS, or beyond main sequence, evolved, or BH). If an object is a member of a binary, then its symbol is filled; if it is not, then its symbol is open. Bottom panel: traces the radial position of these objects in units of the cluster’s instantaneous core radius r_{c}, which r_{c} was determined from the cluster’s stellar distribution by applying the Casertano & Hut (1985) method. The time evolution of r_{c}, 10% Lagrangian radius, and halfmass radius of all model clusters (Table A.1) are shown in Fig. A.1. The same symbol shapes and colour coding as in the top panel are used in the bottom panel: for better visibility of overlapping symbols, the symbol filling was not applied in the bottom panel. Both the primary and the secondary BHs are born from single stars, thereby having vanishing natal spins (Sect. 2). Both the progenitor single stars are merged primordial binaries: the BH can be more massive than the plotted ZAMS progenitor depending on the mass of the latter’s companion, the age of the primordial binary’s merger, and the mass loss during the merger (Sect. 2). Although they have formed well separated within the cluster, the primary and secondary BHs segregate to the cluster centre and pair up dynamically from ≈50 Myr. The third BH becomes bound to this BBH dynamically only a few million years before the merger (lower panel). The BBH merger happens due to this triple interaction (Samsing & RamirezRuiz 2017; Banerjee 2018b, 2021a). 
The merger events from the present computed models were utilized to evaluate the differential merger rate density within the PSN gap. This was done by applying the same cluster population synthesis approach as described in Banerjee (2021b). The methodology is detailed in this reference and, therefore, it is not repeated here. In the present application, a computed M_{cl} = 7.5 × 10^{4} M_{samp} cluster is considered as the representative YMC. Accordingly, when applying Eq. (5) of Banerjee (2021b), the YMC birth mass range is considered to be [M_{cl, low}, M_{cl, high}]=[5.7 × 10^{4} M_{samp}, 10^{5} M_{samp}]. This gives an average cluster mass of ⟨M_{cl}⟩ = 7.5 × 10^{4} M_{samp} for a powerlaw cluster mass function of index α = −2.
Population syntheses were performed with the merger outcomes from four different model sets with [W_{0} = 7 : Z = 0.0002, 0.001, 0.005, 0.01, 0.02] and four different model sets with [W_{0} = 9 : Z = 0.0002, 0.001, 0.005, 0.01, 0.02] (Table A.1). For each cluster model set, three independent Model Universe sample cluster populations, each with a size of N_{samp} = 5 × 10^{5}, were constructed to obtain three values for the presentday event count, N_{mrg}. These, in turn, yield six different differential merger rate density profiles (a reference and a pessimistic rate for each N_{mrg}; see the abovementioned study) per model set. As in the abovementioned study, the ‘moderateZ’ metallicityredshift dependence of Chruslinska & Nelemans (2019) was applied and the detector visibility horizon was taken at z_{max} = 1.0.
Figure 6 shows the resulting presentday, differential intrinsic merger rate density profiles, dℛ/dM_{1}, with respect to merger primary mass, within the PSN gap. Shown are the Model Universe average rates (blackfilled squares) and the upper and lower limits (blue error bars) over the different model sets at each M_{1} bin. As can be seen, the presentday, PSNgap intrinsic merger rate density profiles from the model clusters well accommodate for those that were estimated from GWTC (Abbott et al. 2021c; orange dots). However, the lower limit of the Model Universe dℛ/dM_{1} is zero for all the mass bins in Fig. 6. The upper limit of the Model Universe presentday merger rate density of GW190521like events is also well within the corresponding LVK limits (Abbott et al. 2022; hashed rectangle), but the model average rate is deficient by several factors.
Fig. 6. Presentday, differential intrinsic BBH merger rate density (blackfilled squares with blue error bars), from the Model Universe (Sect. 3.1) composed of the computed model clusters of Table A.1, as a function of merger primary mass, M_{1}, within the PSN mass gap (here considered to be 45 M_{samp} ≤ M_{1} ≤ 120 M_{samp}). At each M_{1} bin, the average rate over the model sets (Sect. 3.1) is indicated with the blackfilled square and the corresponding vertical, blue error bar indicates the maximum and minimum rates for the bin. The orange dots are random draws (300 per bin) of the posteriors of BBH differential intrinsic merger rate density, over the same range of M_{1}, as obtained by LVK (The LIGO Scientific Collaboration 2021b, their power law + peak model). The hashed region spans over the 90% credible intervals of the intrinsic merger rate density of GW190521like events and GW190521’s primary mass, as estimated by LVK (Abbott et al. 2020, 2022). The horizontal, blackdashed line indicates the LVKestimated upper limit of the merger rate density of equal mass, aligned spin (χ_{eff} > 0.8) BBHs of total mass 200 M_{samp} (Abbott et al. 2022). The line spans from the lower primary mass limit of GW190403 to the upper primary mass limit of GW190521. The redfilled circles and the corresponding vertical, red error bars represent the average and the limits of the aligned spin fractions of the Model Universe differential BBH merger rate density over 55 M_{samp} ≲ M_{1} ≲ 75.0 M_{samp}, respectively. The LVKestimated total rates were divided by the bin width (4 M_{samp}), so as to represent them in the dℛ/dM_{1} − M_{1} plane. The lower limit of the Model Universe rate is zero in each bin. Due to the logarithmic scale along the vertical axis, rate values below 10^{−4} yr^{−1} Gpc^{−3} are not shown in this figure. 
3.2. Rate of GW190403like events
The present computed models also produce several BBH mergers of M_{1}, M_{2}, and χ_{eff} similar to GW190403, lying well within the event’s 90% confidence limits, as seen in Figs. 1, 3, and 4. In such model events, at least one of the merging BHs has a = 1 due to undergoing BHTZO accretion or mass accretion in a binary (Sect. 2), prior to participating in dynamical pairing and the merger. The high χ_{eff} is due to chance alignment in dynamical pairing. Figure 7 depicts the timeline of one such event that occurred in one of the W_{0} = 9, Z = 0.005 models.
Fig. 7. Timeline of a GW190403like merger event that occurred inside one of the computed cluster models with W_{0} = 9, Z = 0.005. The meanings of the symbols are the same as in Fig. 5. Here, early in the cluster evolution, the merger primary BH was born in a primordial binary and gained mass via BHTZO accretion from its stellar companion. The merger secondary BH was born in a symbiotic binary. This is how both of the merging BHs started maximally spinning (a_{1} = a_{2} = 1), given the present scheme of assigning BH spins (Sect. 2). As in the example of Fig. 5, the BBH merger happens due to the dynamical pairing of the merging BHs and, thereafter, a triple interaction involving a third perturber BH. In this case, both interactions happened only a few million years before the merger (lower panel). All other GW190403type mergers in the present models are also incluster dynamical mergers, involving at least one such ‘spunup’ BH, except for one which is an ejected merger (see Fig. 3). 
The redfilled circles with error bars in Fig. 6 indicate the Model Universe, spinaligned differential merger rate density of GW190403like events in the mass range 55 M_{samp} ≲ M_{1} ≲ 75 M_{samp} over which such events occur in the computed models. These rates were obtained by multiplying the full differential rates, over the same M_{1} range, with the pvalue for χ_{eff} ≥ 0.8 of an isotropic χ_{eff} distribution corresponding to a_{1} = a_{2} = 1 and q = 0.25 (GW190403like mass ratio). The Model Universe, spinaligned merger rate remains below the LVKestimated upper limit of aligned spin (χ_{eff} ≥ 0.8), equal mass, 200 M_{samp} IMBHIMBH merger rate (Abbott et al. 2022; blackdashed line).
The above, LVKreported IMBHIMBH merger rate density should only be considered as a reference upper limit in this study. Equal mass IMBHIMBH mergers of 200 M_{samp} and high χ_{eff} neither occur in the present models nor have they been detected by LVK yet. The upper merger rate limit is inferred by LVK based on artificially injected events (Abbott et al. 2022). As of now, LVK has not reported a merger rate that applies specifically to GW190403like events. The vast majority of the 2G BHs (such BHs can approach 100 M_{samp} and be of high Kerr parameter) leave the clusters right after the 1G1G BBH mergers that form them due to the associated GW recoil kick (Sect. 2). Although 80 M_{samp} − 120 M_{samp} starstarmergerproduct and massaccreted BHs remain in the clusters, they get dynamically ejected within ≲100 Myr (Fig. 2) unless they get engaged in a massive 1G merger (Sect. 3.1).
Table 1 lists the Model Universe limits of the presentday, intrinsic merger rate densities of the various event types. The lower limit of these total rates and also of the differential rates (see above) is zero since some of the cluster model sets used in the population synthesis (Sect. 3.1) do not produce any PSNgap BBH mergers. The settoset variation of the merger outcome also results in the rather large uncertainty in the Model Universe differential rate profile within the PSN gap (Fig. 6). The rate densities of GW190403 and GW190521type events in Table 1 incorporate only those fractions of Model Universe mergers whose M_{1}, M_{2}, and χ_{eff} all lie within the 90% credible intervals of the respective event parameters.
Intrinsic merger rate densities from the computed models.
4. Summary and concluding remarks
This study investigates dynamical interactions among stars and BHs in YMCs as a potential mechanism for forming BHs within the PSN mass gap and engaging them in GR mergers. To that end, evolutionary models of YMCs (Sect. 2, Table A.1) of a representative initial mass and size of 7.5 × 10^{4} M_{samp} and 2 pc, respectively, are explored. The model clusters also include an observationally motivated population of primordial binaries. The model evolutions were computed with an updated version of the direct Nbody code NBODY7 that incorporates uptodate schemes of stellar mass loss, stellar remnant formation, BH natal spin, and GRmerger recoil.
These computed models produce BBH mergers that are well consistent with the merger masses (M_{1}, M_{2}) of the GWTC events (Sect. 3; Fig. 1). This is also true for the evolutionary cluster models of Ba21, but the present models tend to produce BBH merges involving PSNgap BHs more efficiently, owing to their initial configurations of higher central concentration (Sect. 2, Table A.1). A high central concentration accelerates the central segregation of the most massive cluster members, thereby facilitating starstar and starBH mergers which, in turn, facilitate the formation of PSNgap BHs (Sect. 3 and references therein). The model mergers’ spinorbit alignments (χ_{eff}) are also consistent with the uptodate LVK events (Fig. 3). In particular, the present models do produce BBH mergers resembling GW190521 and GW190403, in all aspects (M_{1}, q, χ_{eff}; Fig. 4).
Here, the large, ≈100 M_{samp}, component masses and vanishing χ_{eff} of the GW190521like event arise due to the involvement of BHs that were directly derived from starstar merger products whose BHs have zero natal spins (Sects. 2 and 3.1). In contrast, the high χ_{eff} of the GW190403like events is due to the involvement of BHs that are spun up via matter accretion (via mass transfer in a binary or from a BHTZO) or due to them being a former member of a tidally interacting binary (Sects. 2 and 3.2). The high central concentration and efficient mass segregation in the model clusters also facilitate the high massasymmetry (see Ba21) in the GW190403type mergers. The present models produce mergers with mass ratios down to ≈0.25 (similar to Ba21). However, beyond M_{1} ≳ 30 M_{samp}, the present, pure dynamical mechanism does not reproduce the bias towards χ_{eff} > 0, as is apparent in the LVK data (Fig. 3). This suggests a contribution from additional merger channels, for example isolated evolution of field massive binaries (e.g. Belczynski et al. 2020b; Wong et al. 2021; Olejak & Belczynski 2021). Further GW events and improved spin measurement in GW observations will better suggest the channel(s) for BBH mergers of high χ_{eff} magnitudes such as for GW190403.
The key uncertainties in the cluster model ingredients are the adopted small starstar merger mass loss and large BHTZO accretion fractions (≤10% and 95%, respectively; see Sect. 2). These settings optimize massive BH formation in the computed cluster models. Some support for a ≲10% mass loss in a starstar merger comes from ‘stickystar’type hydrodynamic calculations (e.g. Lombardi et al. 2002; Gaburov et al. 2008). The nature of BHTZO and the resulting matter accretion onto the BH is, currently, mostly elusive. However, an analogy with the direct collapse BH formation scenario may provide a qualitative justification for the assumed high BHTZO accretion (Banerjee 2019). Another important uncertainty is the spinup of BHs: at present, all BHs undergoing potential scenarios for gaining angular momentum (mass transfer from a stellar companion, membership of a symbiotic binary, and BHTZO accretion) are assigned the maximal spin (a = 1), irrespective of the BH’s natal spin (Sect. 2). Angular momentum deposition onto a BH is, largely, poorly understood and a = 1 is likely an extreme value. A qualitative basis for this model setting is that the typical spin angular momentum of a star (e.g. Lang 1992; Hurley et al. 2000) or orbital angular momentum of a massivestellar binary (e.g. Tauris et al. 2017) is much higher than that of a maximally spinning stellar mass BH^{5}.
The GWmerger recoil of the 2G BHs, as incorporated in the present models based on recent NR results (Sect. 2), ejects most of them from these model clusters which have a moderate, v_{esc} ≈ 40 km s^{−1}, central escape speed. This is why 2G BBH mergers do not contribute to highly aligned, PSNgap BBH mergers (GW190403type events) or any other GR merger involving a highly spinning BH in the present YMC models. (The much larger number of and much longerterm evolutionary cluster models in Ba21 have produced a few 2G BBH mergers.) However, in galactic nuclear clusters, where v_{esc} ≳ 100 km s^{−1}, 2G BHs are more likely to engage in hierarchical BBH mergers despite their GW recoil (Mahapatra et al. 2021). Such 2G BHs can even be of ≈100 M_{samp} and low spin (Antonini et al. 2019; Belczynski & Banerjee 2020).
The present YMC models, nevertheless, well accommodate the LVK presentday, differential intrinsic BBH merger rate density within the PSN gap (Sect. 3.1, Fig. 6). The alignedspin (χ_{eff} > 0.8) fraction of GW190403type mergers, as obtained from these models, make up a presentday rate density that is consistent with the LVKestimated reference upper limit of the rate density of 200 M_{samp}, equalmass, alignedspin BBH mergers (Sect. 3.2, Fig. 6). These results suggest that the tandem of massive binary evolution and dynamical interactions in ≲100 Myr old YMCs can plausibly produce BBH mergers involving BHs within the PSN gap and in rates consistent with what is inferred from uptodate GW observations.
However, the merger rate density of GW190521like events from the present models (0 ≈ 3.8 × 10^{−2} yr^{−1} Gpc^{−3}) is small compared to the models’ overall PSNgap BBH merger rate density (Table 1, Fig. 6). Also, the model merger rate density of GW190521like events fail to reach the corresponding LVKestimated upper limit, as opposed to the comparison for the less massive PSNgap mergers, and its average value is several factors smaller than LVK’s lower limit (see Fig. 6). This means that the present YMC models are rather inefficient in producing GW190521type mergers (see also Fragione & Banerjee 2021), although the upper limit of such events’ yield is still consistent with and within range of the corresponding LVKestimated rate limits. Alternative formation channels (see Sect. 1 and references therein) may be more efficient in producing such massive BBH mergers.
The results presented here do not comprise a ‘full story’, but they represent the yield from YMCs. The modelled clusters represent the most massive young clusters that we observe, such as R136, Westerlund 1. Also, the clusters are evolved up to 300 Myr, when they are still in their young phase, dense, and contain a substantial BH population. They would, therefore, continue to produce latetime GR mergers as they evolve into open clusters, as seen in Ba21. In the near future, such a cluster model set with a range of mass, size, and longterm evolution will be investigated.
As seen in recent studies such as Webb & Sills (2021), the dynamical evolution of YMCs similar to the present models remains practically unaffected by large alterations to the galactocentric distance.
In this work, BHs directly derived from the nuclear evolution of stellar progenitors are called first generation or 1G BHs. If a 1G BH gains mass due to stellar matter accretion, it is still 1G. A 1G+1G BBH merger results in a 2G BH, a BBH merger involving at least one 2G BH yields a 3G BH, and so on. This convention is often followed in the literature.
Special outputs are arranged in the updated NBODY7 for recording M_{1}, M_{2}, a_{1}, and a_{2} of the incluster and ejected mergers. In Fig. 3, each merger from the computed models is plotted with 100 random orientations of θ_{1} and θ_{2}. For the convenience of plotting, a small value of a_{i} = 10^{−2} was assigned if a_{i} = 0.
Of course, a = 1 represents the extremal value of the Kerr parameter and this value is often considered to represent a highly or near maximally spinning BH. In reality, the exact value of unity will not be achieved since an increasing specific angular momentum of the BH would also enhance angular momentum draining via outflow of the accreting matter or GW radiation, for example. However, any value of a infinitesimally close to unity is feasible and in astrophysical scenarios a can practically become unity (see, e.g. Kesden 2008; Benson & Babul 2009). Since there is no physical bound on how close to unity a can get as a BH continues to accrete angular momentum from its surrounding or companion, a = 1 can be considered as a representative of a BH that is subjected to a reservoir of angular momentum. In practice, taking a = 1 or close to 1 causes a negligible difference to the present results (see, e.g. Arca Sedda et al. 2021).
Acknowledgments
The author (S.B.) thanks the reviewer for criticisms and suggestions that have improved the paper. S.B. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) through the individual research grant “The dynamics of stellarmass black holes in dense stellar systems and their role in gravitationalwave generation” (BA 4281/61; PI: S. Banerjee). S.B. acknowledges the generous support and efficient system maintenance of the computing teams at the AIfA and HISKP. This work and the illustrations presented therein are greatly benefited by the use of the Python packages NumPy, SciPy, and Matplotlib. This work has been benefited by discussions with Aleksandra Olejak, Francesco Rizzuto, Chris Belczynski, Giacomo Fragione, Kyle Kremer, Fabio Antonini, Mark Gieles, Thorsten Naab, Jeremiah Ostriker, and Rainer Spurzem. S.B. has solely performed, managed, and analysed all the Nbody computations presented in this work. S.B. has done all the coding necessary for this work, prepared the illustrations, and written the manuscript.
References
 Aarseth, S. J. 2012, MNRAS, 422, 841 [NASA ADS] [CrossRef] [Google Scholar]
 Aasi, J., Abbott, B. P., Abbott, R., et al. 2015, Class. Quant. Grav., 32, 074001 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 031040 [Google Scholar]
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, Phys. Rev. Lett., 125, 101102 [Google Scholar]
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, Phys. Rev. X, 11, 021053 [Google Scholar]
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, ApJ, 915, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2021c, ApJ, 913, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, R., Abbott, T. D., Acernese, F., et al. 2022, A&A, 659, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [Google Scholar]
 Ajith, P., Hannam, M., Husa, S., et al. 2011, Phys. Rev. Lett., 106, 241101 [NASA ADS] [CrossRef] [Google Scholar]
 Antonini, F., Gieles, M., & Gualandris, A. 2019, MNRAS, 486, 5008 [NASA ADS] [CrossRef] [Google Scholar]
 Arca Sedda, M., Mapelli, M., Spera, M., Benacquista, M., & Giacobbo, N. 2020, ApJ, 894, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Arca Sedda, M., Mapelli, M., Benacquista, M., & Spera, M. 2021, MNRAS, submitted [arXiv:2109.12119] [Google Scholar]
 Baibhav, V., Gerosa, D., Berti, E., et al. 2020, Phys. Rev. D, 102, 043002 [NASA ADS] [CrossRef] [Google Scholar]
 Baibhav, V., Berti, E., Gerosa, D., Mould, M., & Wong, K. W. K. 2021, Phys. Rev. D, 104, 084002 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S. 2017, MNRAS, 467, 524 [NASA ADS] [Google Scholar]
 Banerjee, S. 2018a, MNRAS, 473, 909 [CrossRef] [Google Scholar]
 Banerjee, S. 2018b, MNRAS, 481, 5123 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S. 2019, ArXiv eprints [arXiv:1912.06022] [Google Scholar]
 Banerjee, S. 2021a, MNRAS, 500, 3002 [Google Scholar]
 Banerjee, S. 2021b, MNRAS, 503, 3371 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S. 2022, Phys. Rev. D, 105, 023004 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., & Kroupa, P. 2018, in Formation of Very Young Massive Clusters and Implications for Globular Clusters, ed. S. Stahler, Astrophys. Space Sci. Lib., 424, 143 [NASA ADS] [Google Scholar]
 Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371 [Google Scholar]
 Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belczynski, K. 2020, ApJ, 905, L15 [Google Scholar]
 Belczynski, K., & Banerjee, S. 2020, A&A, 640, L20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020a, ApJ, 890, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Klencki, J., Fields, C. E., et al. 2020b, A&A, 636, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benson, A. J., & Babul, A. 2009, MNRAS, 397, 1302 [NASA ADS] [CrossRef] [Google Scholar]
 Branchesi, M. 2016, J. Phys. Conf. Ser., 718, 022004 [NASA ADS] [CrossRef] [Google Scholar]
 Carr, B., & Silk, J. 2018, MNRAS, 478, 3756 [NASA ADS] [CrossRef] [Google Scholar]
 Casertano, S., & Hut, P. 1985, ApJ, 298, 80 [Google Scholar]
 Chruslinska, M., & Nelemans, G. 2019, MNRAS, 488, 5300 [NASA ADS] [Google Scholar]
 Clesse, S., & GarcíaBellido, J. 2017, Phys. Dark Univ., 15, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Costa, G., Bressan, A., Mapelli, M., et al. 2021, MNRAS, 501, 4514 [NASA ADS] [CrossRef] [Google Scholar]
 De Luca, V., Desjacques, V., Franciolini, G., Malhotra, A., & Riotto, A. 2019, JCAP, 2019, 018 [CrossRef] [Google Scholar]
 Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020a, MNRAS, 498, 495 [NASA ADS] [CrossRef] [Google Scholar]
 Di Carlo, U. N., Mapelli, M., Bouffanais, Y., et al. 2020b, MNRAS, 497, 1043 [NASA ADS] [CrossRef] [Google Scholar]
 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
 Edelman, B., Doctor, Z., & Farr, B. 2021, ApJ, 913, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36 [CrossRef] [Google Scholar]
 Fragione, G., & Banerjee, S. 2021, ApJ, 913, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Fragione, G., Loeb, A., & Rasio, F. A. 2020, ApJ, 895, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
 Fuller, J., & Ma, L. 2019, ApJ, 881, L1 [Google Scholar]
 Gaburov, E., Lombardi, J. C., & Portegies Zwart, S. 2008, MNRAS, 383, L5 [NASA ADS] [Google Scholar]
 Gerosa, D., & Fishbach, M. 2021, Nat. Astron., 5, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Glebbeek, E., Gaburov, E., de Mink, S. E., Pols, O. R., & Portegies Zwart, S. F. 2009, A&A, 497, 255 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, E., Kremer, K., Chatterjee, S., et al. 2021, ApJ, 908, L29 [CrossRef] [Google Scholar]
 Hamers, A. S., Fragione, G., Neunteufel, P., & Kocsis, B. 2021, MNRAS, 506, 5345 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D., & Hut, P. 2003, The Gravitational MillionBody Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Hofmann, F., Barausse, E., & Rezzolla, L. 2016, ApJ, 825, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
 KAGRA Collaboration (Akutsu, T., et al.) 2020, Progr. Theor. Exp. Phys., 2021, 05A103 [CrossRef] [Google Scholar]
 Katz, B., Dong, S., & Malhotra, R. 2011, Phys. Rev. Lett., 107, 181101 [NASA ADS] [CrossRef] [Google Scholar]
 Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237 [Google Scholar]
 Kesden, M. 2008, Phys. Rev. D, 78, 084030 [NASA ADS] [CrossRef] [Google Scholar]
 King, I. R. 1966, AJ, 71, 64 [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., McKee, C. F., & BlandHawthorn, J. 2019, ARA&A, 57, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Kuhn, M. A., Povich, M. S., Luhman, K. L., et al. 2013, ApJS, 209, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Lang, K. R. 1992, Astrophysical Data I. Planets and Stars (Berlin, Heidelberg: SpringerVerlag) [CrossRef] [Google Scholar]
 Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lithwick, Y., & Naoz, S. 2011, ApJ, 742, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Lombardi, J. C., Jr., Warren, J. S., Rasio, F. A., Sills, A., & Warren, A. R. 2002, ApJ, 568, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Lousto, C. O., Zlochower, Y., Dotti, M., & Volonteri, M. 2012, Phys. Rev. D, 85, 084015 [NASA ADS] [CrossRef] [Google Scholar]
 Mahapatra, P., Gupta, A., Favata, M., Arun, K. G., & Sathyaprakash, B. S. 2021, ApJ, 918, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., & Broekgaarden, F. S. 2022, Liv. Rev. Relativ., 25, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M. 2018, J. Phys. Conf. Ser., 957, 012001 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Santoliquido, F., Bouffanais, Y., et al. 2021, Symmetry, 13, 1678 [NASA ADS] [CrossRef] [Google Scholar]
 Mathieu, R. D. 2000, in Stellar Clusters and Associations: Convection, Rotation, and Dynamos, eds. R. Pallavicini, G. Micela, & S. Sciortino, ASP Conf. Ser., 198, 517 [NASA ADS] [Google Scholar]
 Mészáros, P., Fox, D. B., Hanna, C., & Murase, K. 2019, Nat. Rev. Phys., 1, 585 [Google Scholar]
 Mikkola, S., & Merritt, D. 2008, AJ, 135, 2398 [Google Scholar]
 Mikkola, S., & Tanikawa, K. 1999, MNRAS, 310, 745 [NASA ADS] [CrossRef] [Google Scholar]
 Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [Google Scholar]
 Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9 [NASA ADS] [CrossRef] [Google Scholar]
 O’Brien, B., Szczepánczyk, M., Gayathri, V., et al. 2021, Phys. Rev. D, 104, 082003 [CrossRef] [Google Scholar]
 Olejak, A., & Belczynski, K. 2021, ApJ, 921, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Perna, R., Wang, Y.H., Farr, W. M., Leigh, N., & Cantiello, M. 2019, ApJ, 878, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
 Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17 [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 507, 3612 [CrossRef] [Google Scholar]
 Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2021, MNRAS, 501, 5257 [NASA ADS] [CrossRef] [Google Scholar]
 Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2022, MNRAS, 512, 884 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., AmaroSeoane, P., Chatterjee, S., & Rasio, F. A. 2018, Phys. Rev. Lett., 120, 151101 [NASA ADS] [CrossRef] [Google Scholar]
 Roupas, Z., & Kazanas, D. 2019, A&A, 632, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samsing, J., & RamirezRuiz, E. 2017, ApJ, 840, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., & Evans, C. J. 2011, in Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits, eds. C. Neiner, G. Wade, G. Meynet, & G. Peters, IAU Symp., 272, 474 [Google Scholar]
 Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sigurdsson, S., & Hernquist, L. 1993, Nature, 364, 423 [Google Scholar]
 Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889 [Google Scholar]
 Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton University Press) [Google Scholar]
 Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., & Kinugawa, T. 2021a, ApJ, 910, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Tanikawa, A., Kinugawa, T., Yoshida, T., Hijikawa, K., & Umeda, H. 2021b, MNRAS, 505, 2170 [CrossRef] [Google Scholar]
 Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170 [Google Scholar]
 The LIGO Scientific Collaboration (Abbott, R., et al.) 2021a, ArXiv eprints [arXiv:2108.01045] [Google Scholar]
 The LIGO Scientific Collaboration (Abbott, R., et al.) 2021b, ArXiv eprints [arXiv:2111.03634] [Google Scholar]
 The LIGO Scientific Collaboration (Abbott, R., et al.) 2021c, ArXiv eprints [arXiv:2111.03606] [Google Scholar]
 Thorne, K. S., & Zytkow, A. N. 1975, ApJ, 199, L19 [NASA ADS] [CrossRef] [Google Scholar]
 VignaGómez, A., Toonen, S., RamirezRuiz, E., et al. 2021, ApJ, 907, L19 [Google Scholar]
 Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, J. J., & Sills, A. 2021, MNRAS, 501, 1933 [Google Scholar]
 Wong, K. W. K., Breivik, K., Kremer, K., & Callister, T. 2021, Phys. Rev. D, 103, 083021 [Google Scholar]
 Woosley, S. E. 2017, ApJ, 836, 244 [Google Scholar]
 Woosley, S. E., & Heger, A. 2021, ApJ, 912, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E., Sukhbold, T., & Janka, H. T. 2020, ApJ, 896, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Ziegler, J., & Freese, K. 2021, Phys. Rev. D, 104, 043015 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Direct Nbody runs of model star clusters
Summary of the new direct Nbody evolutionary models of star clusters and their GRmerger yields in this work.
continued.
Fig. A.1. Time evolution of the core radius, r_{c}, 10% Lagrangian radius, r_{0.1}, and 50% Lagrangian radius (halfmass radius), r_{0.5}, of all the computed cluster models in Table A.1. The panels in the left (right) column correspond to the models with W_{0} = 7 (W_{0} = 9) with Z as indicated in each panel’s title. Furthermore, r_{c}(t), r_{0.1}(t), and r_{0.5}(t) (legend) are plotted for each model separately. 
All Tables
Summary of the new direct Nbody evolutionary models of star clusters and their GRmerger yields in this work.
All Figures
Fig. 1. Comparison of BBH merger masses from computed model star clusters with those from GWTC. Top panels: primary mass, M_{1}, versus secondary mass, M_{2} (M_{1} ≥ M_{2}), of all BBH mergers from the computed model star clusters of Ba21 (filled circles). These data points are colourcoded according to the merger delay time, t_{mrg} (left panel), and the parent model cluster’s metallicity, Z (right panel). The GWTC events’ data points (black, filled squares) and their corresponding 90% credible intervals (horizontal and vertical error bars) are indicated where events with FAR ≤1 yr^{−1} and > 1 yr^{−1} are distinguished by different shades. The shaded regions over 45 M_{samp} − 120 M_{samp} on each panel represent the PSN mass gap. The shadings have been placed in the background so that they do not influence the colours of the data points. Bottom panels: same as the top panels, except that the mergers from the present computed clusters (Table A.1) with King concentration parameters W_{0} = 7.0 (filled triangles) and 9.0 (filled circles) are shown. 

In the text 
Fig. 2. Mass distributions of stellar remnants remaining in the computed model clusters (Table A.1) at 5 Myr, 10 Myr, 150 Myr, and 300 Myr clusterevolutionary times (legend). The distributions, at a given time, from all (four) models for a particular W_{0} and Z are combined in one histogram. The panels in the left (right) column correspond to the models with W_{0} = 7 (W_{0} = 9) with Z as indicated in each panel’s title. 

In the text 
Fig. 3. Effective aligned spin parameter, χ_{eff}, versus chirp mass, M_{chirp} (top panels), and primary mass, M_{1} (bottom panels), of all BBH mergers from the computed model clusters of Table A.1. For each merger, χ_{eff} values for a number of random orientations of the merging BHs’ spins are plotted (Sect. 3). The legends are the same as in Fig. 1, except that a black edge has been applied to the filled triangle symbol for improved visibility. 

In the text 
Fig. 4. Effective aligned spin parameter, χ_{eff}, versus primary mass, M_{1} (top panels), and secondary mass, M_{2} (bottom panels), of all BBH mergers from the computed model clusters of Table A.1. For each model merger, χ_{eff} values for a number of random orientations of the merging BHs’ spins are obtained (Sect. 3), which are shown as a data point (mean value) with a vertical error bar (range). The legends of the data points are the same as in Fig. 1. The error bars (in a dashed line and solid line for mergers from the W_{0} = 7 and W_{0} = 9 models, respectively) share the same colour code as the corresponding data points. The LVK events GW190521 and GW190403 are highlighted with two empty squares, which are coloured differently for improved legibility and they do not correspond to the colour code. 

In the text 
Fig. 5. Timeline of the GW190521like merger event that occurred inside one of the computed cluster models with W_{0} = 9, Z = 0.0002. In both panels the evolutions for the three cluster members, which are most relevant for the event, are shown from the start of the computation until the merger. Top panel: mass evolution of these objects: the trace of the merger primary and secondary mass is indicated by a circle and a square, respectively, and that of the third ‘perturber’ is represented by a triangle. The symbols are colourcoded according to the stellarevolutionary stage (main sequence, MS, or beyond main sequence, evolved, or BH). If an object is a member of a binary, then its symbol is filled; if it is not, then its symbol is open. Bottom panel: traces the radial position of these objects in units of the cluster’s instantaneous core radius r_{c}, which r_{c} was determined from the cluster’s stellar distribution by applying the Casertano & Hut (1985) method. The time evolution of r_{c}, 10% Lagrangian radius, and halfmass radius of all model clusters (Table A.1) are shown in Fig. A.1. The same symbol shapes and colour coding as in the top panel are used in the bottom panel: for better visibility of overlapping symbols, the symbol filling was not applied in the bottom panel. Both the primary and the secondary BHs are born from single stars, thereby having vanishing natal spins (Sect. 2). Both the progenitor single stars are merged primordial binaries: the BH can be more massive than the plotted ZAMS progenitor depending on the mass of the latter’s companion, the age of the primordial binary’s merger, and the mass loss during the merger (Sect. 2). Although they have formed well separated within the cluster, the primary and secondary BHs segregate to the cluster centre and pair up dynamically from ≈50 Myr. The third BH becomes bound to this BBH dynamically only a few million years before the merger (lower panel). The BBH merger happens due to this triple interaction (Samsing & RamirezRuiz 2017; Banerjee 2018b, 2021a). 

In the text 
Fig. 6. Presentday, differential intrinsic BBH merger rate density (blackfilled squares with blue error bars), from the Model Universe (Sect. 3.1) composed of the computed model clusters of Table A.1, as a function of merger primary mass, M_{1}, within the PSN mass gap (here considered to be 45 M_{samp} ≤ M_{1} ≤ 120 M_{samp}). At each M_{1} bin, the average rate over the model sets (Sect. 3.1) is indicated with the blackfilled square and the corresponding vertical, blue error bar indicates the maximum and minimum rates for the bin. The orange dots are random draws (300 per bin) of the posteriors of BBH differential intrinsic merger rate density, over the same range of M_{1}, as obtained by LVK (The LIGO Scientific Collaboration 2021b, their power law + peak model). The hashed region spans over the 90% credible intervals of the intrinsic merger rate density of GW190521like events and GW190521’s primary mass, as estimated by LVK (Abbott et al. 2020, 2022). The horizontal, blackdashed line indicates the LVKestimated upper limit of the merger rate density of equal mass, aligned spin (χ_{eff} > 0.8) BBHs of total mass 200 M_{samp} (Abbott et al. 2022). The line spans from the lower primary mass limit of GW190403 to the upper primary mass limit of GW190521. The redfilled circles and the corresponding vertical, red error bars represent the average and the limits of the aligned spin fractions of the Model Universe differential BBH merger rate density over 55 M_{samp} ≲ M_{1} ≲ 75.0 M_{samp}, respectively. The LVKestimated total rates were divided by the bin width (4 M_{samp}), so as to represent them in the dℛ/dM_{1} − M_{1} plane. The lower limit of the Model Universe rate is zero in each bin. Due to the logarithmic scale along the vertical axis, rate values below 10^{−4} yr^{−1} Gpc^{−3} are not shown in this figure. 

In the text 
Fig. 7. Timeline of a GW190403like merger event that occurred inside one of the computed cluster models with W_{0} = 9, Z = 0.005. The meanings of the symbols are the same as in Fig. 5. Here, early in the cluster evolution, the merger primary BH was born in a primordial binary and gained mass via BHTZO accretion from its stellar companion. The merger secondary BH was born in a symbiotic binary. This is how both of the merging BHs started maximally spinning (a_{1} = a_{2} = 1), given the present scheme of assigning BH spins (Sect. 2). As in the example of Fig. 5, the BBH merger happens due to the dynamical pairing of the merging BHs and, thereafter, a triple interaction involving a third perturber BH. In this case, both interactions happened only a few million years before the merger (lower panel). All other GW190403type mergers in the present models are also incluster dynamical mergers, involving at least one such ‘spunup’ BH, except for one which is an ejected merger (see Fig. 3). 

In the text 
Fig. A.1. Time evolution of the core radius, r_{c}, 10% Lagrangian radius, r_{0.1}, and 50% Lagrangian radius (halfmass radius), r_{0.5}, of all the computed cluster models in Table A.1. The panels in the left (right) column correspond to the models with W_{0} = 7 (W_{0} = 9) with Z as indicated in each panel’s title. Furthermore, r_{c}(t), r_{0.1}(t), and r_{0.5}(t) (legend) are plotted for each model separately. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.