Open Access
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/0004-6361/202142331
Published online 02 September 2022

© S. Banerjee 2022

Licence Creative CommonsOpen 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 Subscribe-to-Open 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 multi-messenger astronomy (Branchesi 2016; Mapelli 2018; Mészáros et al. 2019; Mandel & Broekgaarden 2022). Recently, the LIGO-Virgo-KAGRA 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 up-to-date 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 M 1 = 95 . 3 18.9 + 28.7 M samp $ M_1=95.3_{-18.9}^{+28.7}\,M_{\mathrm{samp}} $ and M 2 = 69 . 0 23.1 + 22.7 M samp $ M_2=69.0_{-23.1}^{+22.7}\,M_{\mathrm{samp}} $ and effective aligned spin parameter χ eff = 0 . 03 0.39 + 0.32 $ \chi_{\mathrm{eff}}=0.03_{-0.39}^{+0.32} $ (Abbott et al. 2021a; the limits correspond to 90% credibility). The χeff (Ajith et al. 2011) is a measure of the spin-orbit alignment of a merging binary and is defined as

χ eff M 1 a 1 cos θ 1 + M 2 a 2 cos θ 2 M 1 + M 2 = a 1 cos θ 1 + q a 2 cos θ 2 1 + q · $$ \begin{aligned} \chi _{\rm eff} \equiv \frac{M_1a_1\cos \theta _1 + M_2a_2\cos \theta _2}{M_1 + M_2} = \frac{a_1\cos \theta _1 + qa_2\cos \theta _2}{1 + q}\cdot \end{aligned} $$(1)

Here, the GR-inspiralling masses M1, M2, with mass ratio q ≡ M2/M1, have Kerr vectors a1, a2 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)

a = c S BH G M BH 2 , $$ \begin{aligned} \boldsymbol{a} = \frac{c\boldsymbol{S}_{\rm BH}}{GM_{\rm BH}^2}, \end{aligned} $$(2)

where SBH is the total angular momentum vector of a Kerr BH of (non-spinning) mass MBH.

Another notable PSN-gap event is GW190403_051519 (hereafter GW190403), which is lighter than GW190521 but – in contrast – is highly mass-asymmetric ( M 1 = 88 . 0 32.9 + 28.2 M samp , M 2 = 22 . 1 9.0 + 23.8 M samp $ M_1=88.0_{-32.9}^{+28.2}\,M_{\mathrm{samp}}, M_2=22.1_{-9.0}^{+23.8}\,M_{\mathrm{samp}} $) and spin-orbit aligned ( χ eff = 0 . 70 0.27 + 0.15 $ \chi_{\mathrm{eff}}=0.70_{-0.27}^{+0.15} $). 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’ pair-instability supernova (hereafter PSN) mass gap between 45 Msamp ≲ Mrem ≲ 120 Msamp (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 remnant-mass 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 pair-instability supernovae (hereafter PPSN). A PPSN episodically sheds the hydrogen envelope of a parent star until a helium core of ≈45 Msamp 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 mass-final 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 GW190521-type events: primordial BHs (e.g. Clesse & García-Bellido 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 proto-clusters (e.g. Roupas & Kazanas 2019), GR coalescences in field hierarchical systems (e.g. Fragione et al. 2020; Vigna-Gó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 PSN-gap BBH mergers is investigated. Dynamical interactions among stellar-remnant BHs inside dense star clusters is an intriguing scenario for generating PSN-gap 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 post-processing of BHs and BBHs that formed as a result of massive binary evolution. Dynamical interactions also enable kinds of star-star and star-remnant 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 hierarchical-system 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, GW190521-like 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 stellar-remnant model that exhibits such a gap.

This paper is organized as follows. Section 2 describes a new set of direct N-body computations of model YMCs. Section 3 discusses the outcomes, focussing on PSN-gap mergers. Section 3.1 describes the inferred rate of GW190521-like events and Sect. 3.2 covers that of GW190403-like mergers. Section 4 summarizes the findings and discusses the implications of various model assumptions and the limitations.

2. Computations

Close dynamical interactions, star-star mergers, and star-remnant 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, PSN-gap 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 Mcl = 7.5 × 104Msamp (initial number of stars N ≈ 1.28 × 105) are considered to be representatives of YMCs. Such a cluster mass is comparable to those of most massive Galactic and local-group 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 half-mass radius of rh = 2 pc and a King dimensionless potential of W0 = 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 W0 ≈ 5−6 (Heggie & Hut 2003), suggest that massive clusters (Mcl ≥ 5 × 104Msamp) are capable of producing and merging PSN-gap BHs. An objective of the present study is to see how somewhat more centrally peaked initial models fair in this, which is why W0 = 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 rh).

Another objective is to see how parsec-scale clusters fair as opposed to highly concentrated, sub-parsec-scale clusters of rh comparable to widths of molecular-cloud 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) gas-free, parsec-scale 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 gas-free nature of such clusters implies that any potential violent-relaxation 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 Msamp ≤ m* ≤ 150.0 Msamp that are distributed according to the canonical initial mass function (hereafter IMF; Kroupa 2001), and they have an overall (see below) primordial-binary fraction of fbin = 5%. However, as in Ba21, the initial binary fraction of the O-type stars (m* ≥ 16.0 Msamp), which are initially paired only among themselves, is considered to be fObin(0) = 100%, which is consistent with the observed high binary fraction among O-stars in young clusters and associations (e.g. Sana & Evans 2011; Sana et al. 2013; Moe & Di Stefano 2017). The O-star binaries are considered to initially follow the observed orbital-period distribution of Sana & Evans (2011) and a uniform mass-ratio distribution. The initial orbital periods of the non-O-star primordial binaries follow the period distribution of Duquennoy & Mayor (1991) and their mass-ratio distribution is also uniform. The initial eccentricity of the O-star 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 star-by-star, direct N-body 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 up-to-date 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 on-the-fly 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 post-Newtonian (hereafter PN) evolution of compact binaries and higher order systems is ARCHAIN (Mikkola & Tanikawa 1999; Mikkola & Merritt 2008). In the present computed models, the ‘F12-rapid+B16-PPSN/PSN’ (see Ba20) remnant-mass prescription is applied. For single star evolution, this remnant prescription allows for the formation of stellar remnants, of mass Mrem, maintaining a neutron star (NS)-BH mass gap between 2 Msamp ≲ Mrem ≲ 5 Msamp (Fryer et al. 2012) and a PSN mass gap between 45 Msamp ≲ Mrem ≲ 120 Msamp (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 non-mass-transferring or non-interacting binaries, as per the BH-formation 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 BH-star 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 BH-star merger (the formation of a BH Thorne–Zytkow object, Thorne & Zytkow 1975 or BH-TZO), the fTZ = 0.95 fraction of the merging star’s mass is assumed to be accreted onto the BH. Also, in star-star mergers, the fmrg = 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 GW-recoil formulae of Lousto et al. (2012) and final-spin formulae of Hofmann et al. (2016).

For each W0 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 (W0, Z) pair, which are subjected to a solar-neighbourhood-like external field2, 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, M1, versus secondary mass, M2 (M1 ≥ M2), of the compact-binary mergers from the computed models. All these events are in-cluster or ejected BBH mergers, except for one that is an NS-BH merger. The BBH-merger 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 NS-BH merger is also consistent with the LVK-detected NS-BH 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.

thumbnail Fig. 1.

Comparison of BBH merger masses from computed model star clusters with those from GWTC. Top panels: primary mass, M1, versus secondary mass, M2 (M1 ≥ M2), of all BBH mergers from the computed model star clusters of Ba21 (filled circles). These data points are colour-coded according to the merger delay time, tmrg (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 Msamp − 120 Msamp 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 W0 = 7.0 (filled triangles) and 9.0 (filled circles) are shown.

The presently computed cluster models are more efficient in producing PSN-gap BBH mergers: the 40 models produce a similar number of mergers with M1 and/or M2 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 rh for both model sets). However, in both sets, mergers involving a PSN-gap BH happen only in models with Z ≤ 0.005 (Fig. 1). In the present models, BHs in the PSN gap appear due to (i) star-star mergers and (ii) BH-TZO accretion. Depending on the stellar-evolutionary age of the merging stars, a star-star merger can result in an over-massive H-envelope of the merged star, with its He-core mass being below the PPSN–PSN threshold. With sufficient H-envelope mass, such a merged star would evolve into a direct-collapse BH within the PSN-gap (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% BH-TZO accretion, as adopted here, would also push BHs into the PSN gap as a result of sufficiently massive BH-star mergers.

In the Ba21 models, a third channel has also produced PSN-gap BBH mergers in the model runs with FM19 (vanishing) BH natal spins, namely, second generation (hereafter 2G) mergers3. In those models, most 2G BBH mergers have occurred with delay times of tmrg > 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, GW190521-like merger obtained from the present model set (Fig. 1; bottom panels) is a 1G, in-cluster merger between two PSN-gap BHs, both of which are derived from star-star merger products. The merger happened in one of the Z = 0.0002, W0 = 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 W0 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 W0, as can be expected.

thumbnail 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 cluster-evolutionary times (legend). The distributions, at a given time, from all (four) models for a particular W0 and Z are combined in one histogram. The panels in the left (right) column correspond to the models with W0 = 7 (W0 = 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, Mchirp (top panels), and primary mass (bottom panels). A merging system’s chirp mass is defined as

M chirp ( M 1 M 2 ) 3 / 5 ( M 1 + M 2 ) 1 / 5 · $$ \begin{aligned} M_{\rm chirp} \equiv \frac{(M_1M_2)^{3/5}}{(M_1+M_2)^{1/5}}\cdot \end{aligned} $$(3)

thumbnail Fig. 3.

Effective aligned spin parameter, χeff, versus chirp mass, Mchirp (top panels), and primary mass, M1 (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 M1, M2, a1, and a2 are taken directly from the model4.

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 spin-aligned 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 PSN-gap mergers (of M1 ≥ 55 Msamp) with nearly vanishing, mildly aligned, and highly aligned χeff values.

Figure 4 re-plots the model χeff as mean values with error bars, against M1 (upper panels) and M2 (lower panels), and it highlights the events GW190521 and GW190403. Figure 4 shows that the present model clusters produce only one GW190521-like and five GW190403-like events, in the sense that the model mergers’ M1, M2, and χeff all lie within or enter the 90% credible intervals (the LVK-data error bars) of the respective parameters of these events of interest.

thumbnail Fig. 4.

Effective aligned spin parameter, χeff, versus primary mass, M1 (top panels), and secondary mass, M2 (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 W0 = 7 and W0 = 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 GW190521-like events

The present computed cluster models yield BBH mergers that nicely resemble notable PSN-gap events discovered by LVK, in terms of both merging masses (M1, M2) and spin-orbit alignment (χeff). This is evident from Figs. 1, 3, and 4. The GW190521-like merger event is an outcome of an in-cluster, dynamically assembled and triggered merger between two BHs of ≈110 Msamp and 57 Msamp, in one of the Z = 0.0002, W0 = 9 models. Both of these BHs are formed due to evolution of star-star merger products with an over-massive H-envelope (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 N-body simulation and is described further in the figure’s caption.

thumbnail Fig. 5.

Timeline of the GW190521-like merger event that occurred inside one of the computed cluster models with W0 = 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 colour-coded according to the stellar-evolutionary 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 rc, which rc was determined from the cluster’s stellar distribution by applying the Casertano & Hut (1985) method. The time evolution of rc, 10% Lagrangian radius, and half-mass 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 & Ramirez-Ruiz 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 Mcl = 7.5 × 104Msamp cluster is considered as the representative YMC. Accordingly, when applying Eq. (5) of Banerjee (2021b), the YMC birth mass range is considered to be [Mcl, low, Mcl, high]=[5.7 × 104Msamp, 105Msamp]. This gives an average cluster mass of ⟨Mcl⟩ = 7.5 × 104Msamp for a power-law cluster mass function of index α = −2.

Population syntheses were performed with the merger outcomes from four different model sets with [W0 = 7 : Z = 0.0002, 0.001, 0.005, 0.01, 0.02] and four different model sets with [W0 = 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 Nsamp = 5 × 105, were constructed to obtain three values for the present-day event count, Nmrg. These, in turn, yield six different differential merger rate density profiles (a reference and a pessimistic rate for each Nmrg; see the above-mentioned study) per model set. As in the above-mentioned study, the ‘moderate-Z’ metallicity-redshift dependence of Chruslinska & Nelemans (2019) was applied and the detector visibility horizon was taken at zmax = 1.0.

Figure 6 shows the resulting present-day, differential intrinsic merger rate density profiles, dℛ/dM1, with respect to merger primary mass, within the PSN gap. Shown are the Model Universe average rates (black-filled squares) and the upper and lower limits (blue error bars) over the different model sets at each M1 bin. As can be seen, the present-day, PSN-gap 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ℛ/dM1 is zero for all the mass bins in Fig. 6. The upper limit of the Model Universe present-day merger rate density of GW190521-like 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.

thumbnail Fig. 6.

Present-day, differential intrinsic BBH merger rate density (black-filled 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, M1, within the PSN mass gap (here considered to be 45 Msamp ≤ M1 ≤ 120 Msamp). At each M1 bin, the average rate over the model sets (Sect. 3.1) is indicated with the black-filled 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 M1, 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 GW190521-like events and GW190521’s primary mass, as estimated by LVK (Abbott et al. 2020, 2022). The horizontal, black-dashed line indicates the LVK-estimated upper limit of the merger rate density of equal mass, aligned spin (χeff > 0.8) BBHs of total mass 200 Msamp (Abbott et al. 2022). The line spans from the lower primary mass limit of GW190403 to the upper primary mass limit of GW190521. The red-filled 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 Msamp ≲ M1 ≲ 75.0 Msamp, respectively. The LVK-estimated total rates were divided by the bin width (4 Msamp), so as to represent them in the dℛ/dM1 − M1 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 M 1 $ M_\odot^{-1} $ are not shown in this figure.

3.2. Rate of GW190403-like events

The present computed models also produce several BBH mergers of M1, M2, 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 BH-TZO 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 W0 = 9, Z = 0.005 models.

thumbnail Fig. 7.

Timeline of a GW190403-like merger event that occurred inside one of the computed cluster models with W0 = 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 BH-TZO 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 (a1 = a2 = 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 GW190403-type mergers in the present models are also in-cluster dynamical mergers, involving at least one such ‘spun-up’ BH, except for one which is an ejected merger (see Fig. 3).

The red-filled circles with error bars in Fig. 6 indicate the Model Universe, spin-aligned differential merger rate density of GW190403-like events in the mass range 55 Msamp ≲ M1 ≲ 75 Msamp over which such events occur in the computed models. These rates were obtained by multiplying the full differential rates, over the same M1 range, with the p-value for χeff ≥ 0.8 of an isotropic χeff distribution corresponding to a1 = a2 = 1 and q = 0.25 (GW190403-like mass ratio). The Model Universe, spin-aligned merger rate remains below the LVK-estimated upper limit of aligned spin (χeff ≥ 0.8), equal mass, 200 Msamp IMBH-IMBH merger rate (Abbott et al. 2022; black-dashed line).

The above, LVK-reported IMBH-IMBH merger rate density should only be considered as a reference upper limit in this study. Equal mass IMBH-IMBH mergers of 200 Msamp 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 GW190403-like events. The vast majority of the 2G BHs (such BHs can approach 100 Msamp and be of high Kerr parameter) leave the clusters right after the 1G-1G BBH mergers that form them due to the associated GW recoil kick (Sect. 2). Although 80 Msamp − 120 Msamp star-star-merger-product and mass-accreted 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 present-day, 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 PSN-gap BBH mergers. The set-to-set 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 GW190521-type events in Table 1 incorporate only those fractions of Model Universe mergers whose M1, M2, and χeff all lie within the 90% credible intervals of the respective event parameters.

Table 1.

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 × 104Msamp 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 N-body code NBODY7 that incorporates up-to-date schemes of stellar mass loss, stellar remnant formation, BH natal spin, and GR-merger recoil.

These computed models produce BBH mergers that are well consistent with the merger masses (M1, M2) 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 PSN-gap 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 star-star and star-BH mergers which, in turn, facilitate the formation of PSN-gap BHs (Sect. 3 and references therein). The model mergers’ spin-orbit alignments (χeff) are also consistent with the up-to-date LVK events (Fig. 3). In particular, the present models do produce BBH mergers resembling GW190521 and GW190403, in all aspects (M1, q, χeff; Fig. 4).

Here, the large, ≈100 Msamp, component masses and vanishing χeff of the GW190521-like event arise due to the involvement of BHs that were directly derived from star-star merger products whose BHs have zero natal spins (Sects. 2 and 3.1). In contrast, the high χeff of the GW190403-like events is due to the involvement of BHs that are spun up via matter accretion (via mass transfer in a binary or from a BH-TZO) 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 mass-asymmetry (see Ba21) in the GW190403-type mergers. The present models produce mergers with mass ratios down to ≈0.25 (similar to Ba21). However, beyond M1 ≳ 30 Msamp, 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 star-star merger mass loss and large BH-TZO 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 star-star merger comes from ‘sticky-star’-type hydrodynamic calculations (e.g. Lombardi et al. 2002; Gaburov et al. 2008). The nature of BH-TZO 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 BH-TZO accretion (Banerjee 2019). Another important uncertainty is the spin-up 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 BH-TZO 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 massive-stellar binary (e.g. Tauris et al. 2017) is much higher than that of a maximally spinning stellar mass BH5.

The GW-merger 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, vesc ≈ 40 km s−1, central escape speed. This is why 2G BBH mergers do not contribute to highly aligned, PSN-gap BBH mergers (GW190403-type events) or any other GR merger involving a highly spinning BH in the present YMC models. (The much larger number of and much longer-term evolutionary cluster models in Ba21 have produced a few 2G BBH mergers.) However, in galactic nuclear clusters, where vesc ≳ 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 Msamp and low spin (Antonini et al. 2019; Belczynski & Banerjee 2020).

The present YMC models, nevertheless, well accommodate the LVK present-day, differential intrinsic BBH merger rate density within the PSN gap (Sect. 3.1, Fig. 6). The aligned-spin (χeff > 0.8) fraction of GW190403-type mergers, as obtained from these models, make up a present-day rate density that is consistent with the LVK-estimated reference upper limit of the rate density of 200 Msamp, equal-mass, aligned-spin 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 up-to-date GW observations.

However, the merger rate density of GW190521-like events from the present models (0 ≈ 3.8 × 10−2 yr−1 Gpc−3) is small compared to the models’ overall PSN-gap BBH merger rate density (Table 1, Fig. 6). Also, the model merger rate density of GW190521-like events fail to reach the corresponding LVK-estimated upper limit, as opposed to the comparison for the less massive PSN-gap 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 GW190521-type mergers (see also Fragione & Banerjee 2021), although the upper limit of such events’ yield is still consistent with and within range of the corresponding LVK-estimated 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 late-time 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 long-term evolution will be investigated.


2

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.

3

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.

4

Special outputs are arranged in the updated NBODY7 for recording M1, M2, a1, and a2 of the in-cluster 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 ai = 10−2 was assigned if ai = 0.

5

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 stellar-mass black holes in dense stellar systems and their role in gravitational-wave generation” (BA 4281/6-1; 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 N-body computations presented in this work. S.B. has done all the coding necessary for this work, prepared the illustrations, and written the manuscript.

References

  1. Aarseth, S. J. 2012, MNRAS, 422, 841 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aasi, J., Abbott, B. P., Abbott, R., et al. 2015, Class. Quant. Grav., 32, 074001 [Google Scholar]
  3. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [Google Scholar]
  4. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 031040 [Google Scholar]
  5. Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, Phys. Rev. Lett., 125, 101102 [Google Scholar]
  6. Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, Phys. Rev. X, 11, 021053 [Google Scholar]
  7. Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, ApJ, 915, L5 [NASA ADS] [CrossRef] [Google Scholar]
  8. Abbott, R., Abbott, T. D., Abraham, S., et al. 2021c, ApJ, 913, L7 [NASA ADS] [CrossRef] [Google Scholar]
  9. Abbott, R., Abbott, T. D., Acernese, F., et al. 2022, A&A, 659, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [Google Scholar]
  11. Ajith, P., Hannam, M., Husa, S., et al. 2011, Phys. Rev. Lett., 106, 241101 [NASA ADS] [CrossRef] [Google Scholar]
  12. Antonini, F., Gieles, M., & Gualandris, A. 2019, MNRAS, 486, 5008 [NASA ADS] [CrossRef] [Google Scholar]
  13. Arca Sedda, M., Mapelli, M., Spera, M., Benacquista, M., & Giacobbo, N. 2020, ApJ, 894, 133 [NASA ADS] [CrossRef] [Google Scholar]
  14. Arca Sedda, M., Mapelli, M., Benacquista, M., & Spera, M. 2021, MNRAS, submitted [arXiv:2109.12119] [Google Scholar]
  15. Baibhav, V., Gerosa, D., Berti, E., et al. 2020, Phys. Rev. D, 102, 043002 [NASA ADS] [CrossRef] [Google Scholar]
  16. Baibhav, V., Berti, E., Gerosa, D., Mould, M., & Wong, K. W. K. 2021, Phys. Rev. D, 104, 084002 [NASA ADS] [CrossRef] [Google Scholar]
  17. Banerjee, S. 2017, MNRAS, 467, 524 [NASA ADS] [Google Scholar]
  18. Banerjee, S. 2018a, MNRAS, 473, 909 [CrossRef] [Google Scholar]
  19. Banerjee, S. 2018b, MNRAS, 481, 5123 [NASA ADS] [CrossRef] [Google Scholar]
  20. Banerjee, S. 2019, ArXiv e-prints [arXiv:1912.06022] [Google Scholar]
  21. Banerjee, S. 2021a, MNRAS, 500, 3002 [Google Scholar]
  22. Banerjee, S. 2021b, MNRAS, 503, 3371 [NASA ADS] [CrossRef] [Google Scholar]
  23. Banerjee, S. 2022, Phys. Rev. D, 105, 023004 [NASA ADS] [CrossRef] [Google Scholar]
  24. 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]
  25. Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371 [Google Scholar]
  26. Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Belczynski, K. 2020, ApJ, 905, L15 [Google Scholar]
  28. Belczynski, K., & Banerjee, S. 2020, A&A, 640, L20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020a, ApJ, 890, 113 [NASA ADS] [CrossRef] [Google Scholar]
  31. Belczynski, K., Klencki, J., Fields, C. E., et al. 2020b, A&A, 636, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Benson, A. J., & Babul, A. 2009, MNRAS, 397, 1302 [NASA ADS] [CrossRef] [Google Scholar]
  33. Branchesi, M. 2016, J. Phys. Conf. Ser., 718, 022004 [NASA ADS] [CrossRef] [Google Scholar]
  34. Carr, B., & Silk, J. 2018, MNRAS, 478, 3756 [NASA ADS] [CrossRef] [Google Scholar]
  35. Casertano, S., & Hut, P. 1985, ApJ, 298, 80 [Google Scholar]
  36. Chruslinska, M., & Nelemans, G. 2019, MNRAS, 488, 5300 [NASA ADS] [Google Scholar]
  37. Clesse, S., & García-Bellido, J. 2017, Phys. Dark Univ., 15, 142 [NASA ADS] [CrossRef] [Google Scholar]
  38. Costa, G., Bressan, A., Mapelli, M., et al. 2021, MNRAS, 501, 4514 [NASA ADS] [CrossRef] [Google Scholar]
  39. De Luca, V., Desjacques, V., Franciolini, G., Malhotra, A., & Riotto, A. 2019, JCAP, 2019, 018 [CrossRef] [Google Scholar]
  40. Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020a, MNRAS, 498, 495 [NASA ADS] [CrossRef] [Google Scholar]
  41. Di Carlo, U. N., Mapelli, M., Bouffanais, Y., et al. 2020b, MNRAS, 497, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  42. Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
  43. Edelman, B., Doctor, Z., & Farr, B. 2021, ApJ, 913, L23 [NASA ADS] [CrossRef] [Google Scholar]
  44. Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36 [CrossRef] [Google Scholar]
  45. Fragione, G., & Banerjee, S. 2021, ApJ, 913, L29 [NASA ADS] [CrossRef] [Google Scholar]
  46. Fragione, G., Loeb, A., & Rasio, F. A. 2020, ApJ, 895, L15 [NASA ADS] [CrossRef] [Google Scholar]
  47. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
  48. Fuller, J., & Ma, L. 2019, ApJ, 881, L1 [Google Scholar]
  49. Gaburov, E., Lombardi, J. C., & Portegies Zwart, S. 2008, MNRAS, 383, L5 [NASA ADS] [Google Scholar]
  50. Gerosa, D., & Fishbach, M. 2021, Nat. Astron., 5, 749 [NASA ADS] [CrossRef] [Google Scholar]
  51. 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]
  52. González, E., Kremer, K., Chatterjee, S., et al. 2021, ApJ, 908, L29 [CrossRef] [Google Scholar]
  53. Hamers, A. S., Fragione, G., Neunteufel, P., & Kocsis, B. 2021, MNRAS, 506, 5345 [NASA ADS] [CrossRef] [Google Scholar]
  54. Heggie, D., & Hut, P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
  55. Hofmann, F., Barausse, E., & Rezzolla, L. 2016, ApJ, 825, L19 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [Google Scholar]
  57. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  58. KAGRA Collaboration (Akutsu, T., et al.) 2020, Progr. Theor. Exp. Phys., 2021, 05A103 [CrossRef] [Google Scholar]
  59. Katz, B., Dong, S., & Malhotra, R. 2011, Phys. Rev. Lett., 107, 181101 [NASA ADS] [CrossRef] [Google Scholar]
  60. Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237 [Google Scholar]
  61. Kesden, M. 2008, Phys. Rev. D, 78, 084030 [NASA ADS] [CrossRef] [Google Scholar]
  62. King, I. R. 1966, AJ, 71, 64 [Google Scholar]
  63. Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
  64. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  65. Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA&A, 57, 227 [NASA ADS] [CrossRef] [Google Scholar]
  66. Kuhn, M. A., Povich, M. S., Luhman, K. L., et al. 2013, ApJS, 209, 29 [NASA ADS] [CrossRef] [Google Scholar]
  67. Lang, K. R. 1992, Astrophysical Data I. Planets and Stars (Berlin, Heidelberg: Springer-Verlag) [CrossRef] [Google Scholar]
  68. Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Lithwick, Y., & Naoz, S. 2011, ApJ, 742, 94 [NASA ADS] [CrossRef] [Google Scholar]
  70. Lombardi, J. C., Jr., Warren, J. S., Rasio, F. A., Sills, A., & Warren, A. R. 2002, ApJ, 568, 939 [NASA ADS] [CrossRef] [Google Scholar]
  71. Lousto, C. O., Zlochower, Y., Dotti, M., & Volonteri, M. 2012, Phys. Rev. D, 85, 084015 [NASA ADS] [CrossRef] [Google Scholar]
  72. Mahapatra, P., Gupta, A., Favata, M., Arun, K. G., & Sathyaprakash, B. S. 2021, ApJ, 918, L31 [NASA ADS] [CrossRef] [Google Scholar]
  73. Mandel, I., & Broekgaarden, F. S. 2022, Liv. Rev. Relativ., 25, 1 [NASA ADS] [CrossRef] [Google Scholar]
  74. Mapelli, M. 2018, J. Phys. Conf. Ser., 957, 012001 [NASA ADS] [CrossRef] [Google Scholar]
  75. Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76 [NASA ADS] [CrossRef] [Google Scholar]
  76. Mapelli, M., Santoliquido, F., Bouffanais, Y., et al. 2021, Symmetry, 13, 1678 [NASA ADS] [CrossRef] [Google Scholar]
  77. 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]
  78. Mészáros, P., Fox, D. B., Hanna, C., & Murase, K. 2019, Nat. Rev. Phys., 1, 585 [Google Scholar]
  79. Mikkola, S., & Merritt, D. 2008, AJ, 135, 2398 [Google Scholar]
  80. Mikkola, S., & Tanikawa, K. 1999, MNRAS, 310, 745 [NASA ADS] [CrossRef] [Google Scholar]
  81. Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [Google Scholar]
  82. Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9 [NASA ADS] [CrossRef] [Google Scholar]
  83. O’Brien, B., Szczepánczyk, M., Gayathri, V., et al. 2021, Phys. Rev. D, 104, 082003 [CrossRef] [Google Scholar]
  84. Olejak, A., & Belczynski, K. 2021, ApJ, 921, L2 [NASA ADS] [CrossRef] [Google Scholar]
  85. Perna, R., Wang, Y.-H., Farr, W. M., Leigh, N., & Cantiello, M. 2019, ApJ, 878, L1 [NASA ADS] [CrossRef] [Google Scholar]
  86. Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
  87. Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17 [Google Scholar]
  88. Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
  89. Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 507, 3612 [CrossRef] [Google Scholar]
  90. Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2021, MNRAS, 501, 5257 [NASA ADS] [CrossRef] [Google Scholar]
  91. Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2022, MNRAS, 512, 884 [NASA ADS] [CrossRef] [Google Scholar]
  92. Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., & Rasio, F. A. 2018, Phys. Rev. Lett., 120, 151101 [NASA ADS] [CrossRef] [Google Scholar]
  93. Roupas, Z., & Kazanas, D. 2019, A&A, 632, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Samsing, J., & Ramirez-Ruiz, E. 2017, ApJ, 840, L14 [NASA ADS] [CrossRef] [Google Scholar]
  95. 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]
  96. Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Sigurdsson, S., & Hernquist, L. 1993, Nature, 364, 423 [Google Scholar]
  98. Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889 [Google Scholar]
  99. Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton University Press) [Google Scholar]
  100. Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., & Kinugawa, T. 2021a, ApJ, 910, 30 [NASA ADS] [CrossRef] [Google Scholar]
  101. Tanikawa, A., Kinugawa, T., Yoshida, T., Hijikawa, K., & Umeda, H. 2021b, MNRAS, 505, 2170 [CrossRef] [Google Scholar]
  102. Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170 [Google Scholar]
  103. The LIGO Scientific Collaboration (Abbott, R., et al.) 2021a, ArXiv e-prints [arXiv:2108.01045] [Google Scholar]
  104. The LIGO Scientific Collaboration (Abbott, R., et al.) 2021b, ArXiv e-prints [arXiv:2111.03634] [Google Scholar]
  105. The LIGO Scientific Collaboration (Abbott, R., et al.) 2021c, ArXiv e-prints [arXiv:2111.03606] [Google Scholar]
  106. Thorne, K. S., & Zytkow, A. N. 1975, ApJ, 199, L19 [NASA ADS] [CrossRef] [Google Scholar]
  107. Vigna-Gómez, A., Toonen, S., Ramirez-Ruiz, E., et al. 2021, ApJ, 907, L19 [Google Scholar]
  108. Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146 [NASA ADS] [CrossRef] [Google Scholar]
  109. Webb, J. J., & Sills, A. 2021, MNRAS, 501, 1933 [Google Scholar]
  110. Wong, K. W. K., Breivik, K., Kremer, K., & Callister, T. 2021, Phys. Rev. D, 103, 083021 [Google Scholar]
  111. Woosley, S. E. 2017, ApJ, 836, 244 [Google Scholar]
  112. Woosley, S. E., & Heger, A. 2021, ApJ, 912, L31 [NASA ADS] [CrossRef] [Google Scholar]
  113. Woosley, S. E., Sukhbold, T., & Janka, H. T. 2020, ApJ, 896, 56 [NASA ADS] [CrossRef] [Google Scholar]
  114. Ziegler, J., & Freese, K. 2021, Phys. Rev. D, 104, 043015 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Direct N-body runs of model star clusters

Table A.1.

Summary of the new direct N-body evolutionary models of star clusters and their GR-merger yields in this work.

Table A.1.

continued.

thumbnail Fig. A.1.

Time evolution of the core radius, rc, 10% Lagrangian radius, r0.1, and 50% Lagrangian radius (half-mass radius), r0.5, of all the computed cluster models in Table A.1. The panels in the left (right) column correspond to the models with W0 = 7 (W0 = 9) with Z as indicated in each panel’s title. Furthermore, rc(t), r0.1(t), and r0.5(t) (legend) are plotted for each model separately.

All Tables

Table 1.

Intrinsic merger rate densities from the computed models.

Table A.1.

Summary of the new direct N-body evolutionary models of star clusters and their GR-merger yields in this work.

All Figures

thumbnail Fig. 1.

Comparison of BBH merger masses from computed model star clusters with those from GWTC. Top panels: primary mass, M1, versus secondary mass, M2 (M1 ≥ M2), of all BBH mergers from the computed model star clusters of Ba21 (filled circles). These data points are colour-coded according to the merger delay time, tmrg (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 Msamp − 120 Msamp 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 W0 = 7.0 (filled triangles) and 9.0 (filled circles) are shown.

In the text
thumbnail 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 cluster-evolutionary times (legend). The distributions, at a given time, from all (four) models for a particular W0 and Z are combined in one histogram. The panels in the left (right) column correspond to the models with W0 = 7 (W0 = 9) with Z as indicated in each panel’s title.

In the text
thumbnail Fig. 3.

Effective aligned spin parameter, χeff, versus chirp mass, Mchirp (top panels), and primary mass, M1 (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
thumbnail Fig. 4.

Effective aligned spin parameter, χeff, versus primary mass, M1 (top panels), and secondary mass, M2 (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 W0 = 7 and W0 = 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
thumbnail Fig. 5.

Timeline of the GW190521-like merger event that occurred inside one of the computed cluster models with W0 = 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 colour-coded according to the stellar-evolutionary 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 rc, which rc was determined from the cluster’s stellar distribution by applying the Casertano & Hut (1985) method. The time evolution of rc, 10% Lagrangian radius, and half-mass 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 & Ramirez-Ruiz 2017; Banerjee 2018b, 2021a).

In the text
thumbnail Fig. 6.

Present-day, differential intrinsic BBH merger rate density (black-filled 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, M1, within the PSN mass gap (here considered to be 45 Msamp ≤ M1 ≤ 120 Msamp). At each M1 bin, the average rate over the model sets (Sect. 3.1) is indicated with the black-filled 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 M1, 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 GW190521-like events and GW190521’s primary mass, as estimated by LVK (Abbott et al. 2020, 2022). The horizontal, black-dashed line indicates the LVK-estimated upper limit of the merger rate density of equal mass, aligned spin (χeff > 0.8) BBHs of total mass 200 Msamp (Abbott et al. 2022). The line spans from the lower primary mass limit of GW190403 to the upper primary mass limit of GW190521. The red-filled 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 Msamp ≲ M1 ≲ 75.0 Msamp, respectively. The LVK-estimated total rates were divided by the bin width (4 Msamp), so as to represent them in the dℛ/dM1 − M1 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 M 1 $ M_\odot^{-1} $ are not shown in this figure.

In the text
thumbnail Fig. 7.

Timeline of a GW190403-like merger event that occurred inside one of the computed cluster models with W0 = 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 BH-TZO 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 (a1 = a2 = 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 GW190403-type mergers in the present models are also in-cluster dynamical mergers, involving at least one such ‘spun-up’ BH, except for one which is an ejected merger (see Fig. 3).

In the text
thumbnail Fig. A.1.

Time evolution of the core radius, rc, 10% Lagrangian radius, r0.1, and 50% Lagrangian radius (half-mass radius), r0.5, of all the computed cluster models in Table A.1. The panels in the left (right) column correspond to the models with W0 = 7 (W0 = 9) with Z as indicated in each panel’s title. Furthermore, rc(t), r0.1(t), and r0.5(t) (legend) are plotted for each model separately.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.