Evidence for Hierarchical Black Hole Mergers in the Second LIGO–Virgo Gravitational Wave Catalog

The second LIGO–Virgo Gravitational Wave Transient Catalog (GWTC-2) has significantly expanded our set of gravitational wave (GW) observations (Abbott et al. 2021a). It contains a total of 46 binary black hole (BBH) candidates, excluding GW190814 (Abbott et al. 2020a) whose source could also be a neutron star–black hole binary, whereas the previous catalog only contained 10 BBHs (Abbott et al. 2019). Multiple astrophysical formation channels have been suggested to explain the population of BBHs, and each of these have uncertainties in their underlying physics (e.g., Kruckow et al. 2016; Rodriguez et al. 2016; Abbott et al. 2016a; Klencki et al. 2018; Sasaki et al. 2018; Kumamoto et al. 2020; Tagawa et al. 2020; Tang et al. 2020; Di Carlo et al. 2020b; Zevin et al. 2020b). GW observations can constrain the relative contribution of formation channels and their uncertain physics, and as the catalog grows these constraints become more precise (Stevenson et al. 2017; Talbot & Thrane 2017; Vitale et al. 2017; Zevin et al. 2017; Barrett et al. 2018; Fishbach et al. 2018; Zevin et al. 2021).

Among the GWTC-2 systems there are high-mass BBHs that have components with masses of ≳45 M_⊙ (Abbott et al. 2021a), the most massive being the source of GW190521 (Abbott et al. 2020b). Black holes of ∼45–135M_⊙ are not typically expected to form via standard stellar evolution as the pair-instability process either limits the maximum mass of the progenitor star's core or completely disrupts the star entirely (Fryer et al. 2001; Heger & Woosley 2002; Belczynski et al. 2016; Spera & Mapelli 2017; Farmer et al. 2019; Stevenson et al. 2019; Farmer et al. 2020; Woosley & Heger 2021). Potential (non-mutually exclusive) astrophysical formation mechanisms for black holes in this mass gap include hierarchical mergers, where the remnant of a previous merger becomes part of a new binary (Miller & Hamilton 2002; Antonini & Rasio 2016; Gerosa & Berti 2017; Rodriguez et al. 2019; Yang et al. 2019; Anagnostou et al. 2020; Fragione & Silk 2020; Mapelli et al. 2020; Fragione et al. 2020a; Banerjee 2021); stellar mergers, which may result in a larger hydrogen envelope around a core below the pair-instability threshold (Spera et al. 2018; Kremer et al. 2020; Di Carlo et al. 2020a; González et al. 2021); formation of black holes from Population III stars that are able to retain their hydrogen envelopes (Farrell et al. 2021; Kinugawa et al. 2021; Vink et al. 2021), formation via stellar triples in the field (Vigna-Gómez et al. 2021); growth via accretion in an active galactic nucleus (AGN) disk (McKernan et al. 2012; Michaely & Perets 2020; Secunda et al. 2020; Tagawa et al. 2020), or growth via rapid gas accretion in dense primordial clusters (Roupas & Kazanas 2019).

Hierarchical mergers in globular clusters were considered as an origin for GW190521 (Abbott et al. 2020c) with inconclusive results. However, hints of eccentricity in follow-up analyses of GW190521 add weight to this explanation (Gayathri et al. 2020; Romero-Shaw et al. 2020a). To confidently identify hierarchical mergers, it is important to study events in the context of a population model that fits the mass distribution (and any mass cut-offs) for the first-generation (1G) black holes not formed through mergers (Doctor et al. 2020; Sedda et al. 2020; Kimball et al. 2020a; Tiwari & Fairhurst 2021).

We apply the population inference framework from Kimball et al. (2020b) to analyze the BBHs from GWTC-2. This framework assumes a phenomenological population model based on simulations of metal-poor globular clusters from Rodriguez et al. (2019). Considering a fiducial set of globular cluster masses, we simultaneously infer the properties of the 1G+1G BBH population—whose remnants are second-generation (2G) black holes—and the relative merger rates of hierarchical mergers. The expanded catalog enables the population parameters, including the mass distribution, to be more precisely determined (Abbott et al. 2021b). We find that several of the BBHs are likely to be the results of hierarchical mergers: the leading candidates are GW190519_153544 (GW190519) and GW190521.

In Section 2 we review the key components of our population inference framework; the results of this are given in Section 3, with additional description of the population hyperparameters in the Appendix, and we discuss our findings in Section 4.

We perform Bayesian hierarchical inference to infer the the population properties of BBHs following Kimball et al. (2020b). We employ phenomenological models for the mass and spin distributions of 1G+1G, 1G+2G, and 2G+2G BBHs merging in a dense stellar environment; see Figures 1 and 2. The 1G+1G model is nearly identical to population models used in Abbott et al. (2021b): it is equivalent to the Power Law + Peak mass model (but omits the low-mass smoothing and adopts a Gaussian prior on the maximum mass cut-off) and is similar to the Default spin model. We consider two separate modifications to that spin model: one that adds a parameter that allows for a subpopulation of zero-spin BBHs, and one consisting of a truncated Gaussian with a broad prior on the mean that allows for distributions with sharp peaks at 0 or 1; we refer to these as Model ZeroSubPop and Model TruncGauss, respectively. The particulars of the mass and spin models are discussed further in the Appendix.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The 2G black holes are assumed to be roughly twice the mass of 1G black holes; the mass ratio distribution for 1G+2G binaries is peaked around q ∼ 1/2 while the 2G+2G distribution is similar to the 1G+1G model but with an increased preference for near equal-mass binaries to account for the more massive components in a strong encounter forming bound binaries (Sigurdsson & Phinney 1993; Heggie et al. 1996; Downing et al. 2011). The 1G+2G and 2G+2G spin models presume that 2G black holes have dimensionless spin χ ≈ 0.67 inherited from the orbital angular momentum of the progenitor binary (Pretorius 2005; Gonzalez et al. 2007; Buonanno et al. 2008). The population models are described as conditional priors π(θ∣Λ) where θ are the parameters of a single binary (e.g., mass and spin) while Λ refers to the population hyperparameters describing the shape of the mass and spin distributions (e.g., the power-law index of the primary black hole mass spectrum). Our goal is two-fold: estimate the population hyperparameters Λ and carry out model selection to evaluate the Bayesian odds that events in GWTC-2 are formed hierarchically.

The relative rates of 1G+1G, 1G+2G, and 2G+2G mergers depend upon the properties of their cluster environment as well as the masses and spins of the BBH population. GW recoil kicks may lead to remnants being ejected from a cluster; kick magnitudes are strongly dependent on progenitor spins with larger spins leading to larger kicks (Campanelli et al. 2007; Gonzalez et al. 2007; Bruegmann et al. 2008; Lousto & Zlochower 2011; Varma et al. 2019), as well as the mass ratio of the merging binary. We calculate the fraction of retained merger remnants F_ret given the population properties of the 1G black holes and assuming a cluster described by a Plummer potential (Plummer 1911) mass of M_c with Plummer radius r_c. For our default cluster, we assume that M_c = 5 × 10⁵ M_⊙ and r_c = 1 pc, corresponding to a central escape velocity of ∼65 km s⁻¹. We assume that the relative merger rates scale as R_1G+2G /R_1G+1G ∝ F_ret, ${R}_{2{\rm{G}}+2{\rm{G}}}/{R}_{1{\rm{G}}+1{\rm{G}}}\propto {F}_{\mathrm{ret}}^{2}$, with the constant of proportionality calibrated against globular cluster simulations (Rodriguez et al. 2019).

For our analysis, we consider the 44 BBH (excluding GW190814, for which the nature of the secondary component is unknown) used in the GWTC-2 population analysis (Abbott et al. 2021b). For GWTC-1 events, we use the same single-event posterior samples as Kimball et al. (2020b), for GW190412 we use samples from Zevin et al. (2020a), for GW190521 we use the preferred samples from Abbott et al. (2020c), and for the other GWTC-2 events we use the public samples from Abbott et al. (2021a). ¹³ For the new GWTC-2 events we use results calculated with the IMRPhenomD and IMRPhenomPv2 waveforms (Hannam et al. 2014; Khan et al. 2016). We generate posterior samples for population hyperparameters Λ using the nested sampler dynesty (Speagle 2020) using the GWPopulation framework (Talbot et al. 2019), which takes advantage of Bilby (Ashton et al. 2019; Romero-Shaw et al. 2020b).

3.1. Inferred Populations

3.2. Relative Merger Rates

GW recoil kick velocities generally increase with the spin magnitudes of merging black holes. Figure 2 illustrates that while the spin distribution inferred using Model TruncGauss does not explicitly include a subpopulation at χ = 0, a larger portion of the population has spins less than ∼0.1 than for Model ZeroSubPop, which results in lower typical recoil velocities and hence higher relative hierarchical merger rates. In Figure 3, we plot the posteriors for these rates, as well as for the fraction λ₀ of 1G+1G black holes in Model ZeroSubPop with zero spin. Using Model TruncGauss, we infer median relative rates ${{ \mathcal R }}_{1{\rm{G}}+2{\rm{G}}}/{{ \mathcal R }}_{1{\rm{G}}+1{\rm{G}}}=5.8\times {10}^{-3}$ and ${{ \mathcal R }}_{2{\rm{G}}+2{\rm{G}}}/{{ \mathcal R }}_{1{\rm{G}}+1{\rm{G}}}\,=1.7\times {10}^{-5}$ with 99% upper limits of 0.12 and 7.6 × 10⁻³, respectively. Taking into account selection effects, the detected population would have median relative rates of ${{ \mathcal R }}_{1{\rm{G}}+2{\rm{G}}\,}^{\det }/{{ \mathcal R }}_{1{\rm{G}}+1{\rm{G}}\,}^{\det }=0.01$ and ${{ \mathcal R }}_{2{\rm{G}}+2{\rm{G}}}/{{ \mathcal R }}_{1{\rm{G}}+1{\rm{G}}}=7.0\times {10}^{-5}$ with 99% upper limits of 0.23 and 0.03, respectively. Using Model ZeroSubPop, these rates become 5.3 × 10⁻³ (0.01) and 1.4 × 10⁻⁵ (5.7 × 10⁻⁵), with 99% upper limits of 0.04 (0.08) and 9.8 × 10⁻⁴ (4.0 × 10⁻³), respectively, for the astrophysical (detected) population. The median inferred relative rates are roughly twice those found using the same model in Kimball et al. (2020b), though consistent with the upper limits reported there. The results for both models are consistent with the results of Monte Carlo modeling of black hole populations in globular clusters: Rodriguez et al. (2019) found that the ≈14% of merging BBHs from the underling population in their models contain 2G black holes in the extreme case where all 1G black holes have zero spin (this fraction drops to ≲1% when they increase 1G black hole spins to χ = 0.5).

Zoom In Zoom Out Reset image size

3.3. Posterior Odds for Hierarchical Origin

For each event in GWTC-2, we calculate the posterior odds ${ \mathcal O }$ in favor of hierarchical versus 1G+1G origin. We plot the odds in favor of 2G+2G versus 1G+1G origin assuming Model TruncGauss and Model ZeroSubPop in Figure 4.

Zoom In Zoom Out Reset image size

Assuming Model TruncGauss, we find that across all 44 BBHs in GWTC-2 the probability that at least one binary contains a 2G black hole is 99%. GW190519 and GW190521 are most likely of 1G+2G origin, favored over a 1G+1G origin with 1.2:1 and 2:1 odds, respectively. We also favor a 2G+2G origin over 1G+1G for GW190521 with odds of 1.2:1. We find roughly even odds for GW190602_175927 (GW190602) and GW190706_222641 (GW190706) being of 1G+2G origin. As in Kimball et al. (2020b), we find that GW170729 is most likely of 1G+1G origin, though at slightly higher odds of 1:10 of being of 1G+2G origin rather than 1G+1G.

Using Model ZeroSubPop, which finds lower relative hierarchical merger rates, odds decrease across all events. We find that the probability that at least one binary in GWTC-2 contains a 2G black hole is 96%. GW190519 is marginally favored to have a 1G+2G origin with 1.1:1 odds over a 1G+1G origin. Meanwhile, GW190521 has roughly even odds of being 1G+2G versus 1G+1G at 1.0:1 odds, and a 2G+2G origin is disfavored to a 1G+1G origin at 1:4 odds.

3.4. Varying Cluster Parameters

Our default Plummer model is chosen as representative of a typical globular cluster environment such as those in the vicinity of the Milky Way today, where central escape velocities are on the order of tens of kilometers per second (Baumgardt & Hilker 2018). However, globular clusters in the Milky Way may have been up to a few times more massive at formation than at present (Webb & Leigh 2015). Furthermore, hierarchical mergers may occur in a wide range of dynamical environments with significantly different escape velocities, including AGN disks and nuclear star clusters. Although our phenomenological models are tuned to the results of simulations of typical present-day globular clusters (Rodriguez et al. 2019), we can get an illustrative idea of how results scale with the mass and compactness of the assumed dynamical environment by varying the parameters of our simple Plummer model. We do not expect all BBHs to come from a single type of cluster, but our results let us explore a range of different average cluster sizes.

In Figure 5, we show results when considering models with Plummer masses 10⁴–10⁹ M_⊙ and radii 0.01–1 pc; both these parameters vary cluster escape velocities and thus the retention rate of hierarchical mergers. At low escape velocities (∼10–50 km s⁻¹), almost no 1G+1G merger products are retained and the relative 1G+2G and 2G+2G rates are negligible. In this case, the inferred posterior on the maximum black hole mass m_max shifts away from the astrophysical prior toward higher masses in order to accommodate massive GWTC-2 events as 1G+1G BBHs. As we move toward models with higher central escape velocities, the fraction of retained 1G+1G merger products and the relative 1G+2G and 2G+2G rates rapidly increase, and the odds in favor of GWTC-2 events being of hierarchical origin grow. It is expected that the rate of hierarchical mergers strongly depends on the escape velocity of their dynamical environments (Holley-Bockelmann et al. 2008; Moody & Sigurdsson 2009; Antonini et al. 2019; Gerosa & Berti 2019; Rodriguez et al. 2020; Sedda et al. 2020; Fragione et al. 2020b; Mapelli et al. 2021). We find that even a modest increase in the central escape velocity to ∼90 km s⁻¹ leads to the conclusion that GW190521 is favored to be a 2G+2G versus 1G+1G merger at >10:1 odds, and that the probability that at least one of the BBHs in GWTC-2 contains a 2G black hole is >99.99%.

Zoom In Zoom Out Reset image size

We find that for all of our assumed cluster models, the events in GWTC-2 are better fit including the hierarchical channels than when excluding those channels (equivalent to setting V_esc = 0 km s⁻¹), with the highest Bayes factors corresponding to models where the central escape velocities are ∼300 km s⁻¹. In Figure 5, we show Bayes factors in favor of our hierarchical model versus a model with only 1G+1G BBHs; taking the ratio of these Bayes factors gives cluster-wise Bayes factors comparing how well the data are supported by different cluster models. For clusters with escape velocities of ∼300 km s⁻¹, which have the highest Bayes factors, we find that 99% of 1G+1G black holes are below 40M_⊙ (40M_⊙) and that 99% of all black holes are below 67M_⊙ (66M_⊙) using Model TruncGauss (Model ZeroSubPop). We infer median relative rates for Model TruncGauss (Model ZeroSubPop) of 1G+2G and 2G+2G versus 1G+1G mergers of 0.15 (0.14) and 0.01 (9.2 × 10⁻³) respectively, with 99% upper limits of 0.29 (0.25) and 0.04 (0.03). When accounting for selection effects, we infer median relative rates for the detected population with Model TruncGauss (Model ZeroSubPop) of 1G+2G and 2G+2G versus 1G+1G mergers of 0.3 (0.26) and 0.05 (0.04) respectively, with 99% upper limits of 0.57 (0.48) and 0.18 (0.13). When V_esc ∼ 300 km s⁻¹, we find that GW190521 is most likely of 2G+2G origin, with 1200:1 and 700:1 odds in favor of being 2G+2G versus 1G+1G using Model TruncGauss and Model ZeroSubPop, respectively, with both spin models favoring 2G+2G over 1G+2G origin at ∼3.5:1. For both models, we find that GW190602, GW190620_030421 (GW190620), GW190706, and GW190519 are most likely of 1G+2G origin, favored over 1G+1G origin at >10:1 odds. Using Model ZeroSubPop, we also find that GW190517_055101 (GW190517) is favored to be of 1G+2G over 1G+1G origin at >10:1 odds.

GW observations have demonstrated that black holes merge to form more massive black holes (Abbott et al. 2016b). If these merger products form new binaries, they may again merge as a detectable GW source. It is necessary to consider this hierarchical merger channel when using catalogs of GW sources to make inferences about the physics of black hole formation. For example, inference of the location of the lower edge of the pair-instability mass gap, which could potentially constrain nuclear reaction rates (Farmer et al. 2020) or beyond Standard Model physics (Croon et al. 2020; Straight et al. 2020; Baxter et al. 2021), using detections of black holes in the ≳50 M_⊙ regime would be contaminated by the presence of 2G black holes. In order to distinguish between 1G and 2G black holes, we must account simultaneously for the shapes of 1G and 2G populations and the relative rate of hierarchical mergers. Here, we apply the analysis of Kimball et al. (2020b) to 44 BBHs in GWTC-2, and self-consistently infer a black hole population that accounts for 1G+1G, 1G+2G, and 2G+2G binary mergers, as well as the relative branching ratios between them, in order to identify candidate hierarchical mergers in the current catalog of GW sources.

We find the following, assuming our nominal globular cluster model with M_c = 5 × 10⁵ M_⊙ and r_c = 1 pc.

• 1.  
  The 44 events in GWTC-2 are best modeled when allowing for hierarchical formation channels. For Model TruncGauss and Model ZeroSubPop, we find Bayes factors of 5 and 7, respectively, in favor of including hierarchical components.
• 2.  
  At least one BBH in GWTC-2 contains a 2G black hole with 99% and 96% probability using Model TruncGauss and Model ZeroSubPop, respectively.
• 3.  
  The two binaries that are most likely to contain a 2G black hole are GW190519 and GW190521, with 1.2:1 and 2.0:1 odds, respectively, of being 1G+2G versus 1G+1G assuming Model TruncGauss. Using Model ZeroSubPop, we find that both events have approximately equal odds of being 1G+2G and 1G+1G.
• 4.  
  The relative rates of hierarchical mergers are dependent on how the 1G+1G spin is modeled. Using Model TruncGauss, the median relative merger rates of 1G+2G and 2G+2G to 1G+1G mergers are inferred to be 5.8 × 10⁻³ and 1.7 × 10⁻⁵, respectively, with 99% upper limits of 0.12 and 7.6 × 10⁻³. While using Model ZeroSubPop, the relative rates drop to 5.3 × 10⁻³ and 1.4 × 10⁻⁵, with 99% upper limits of 0.04 and 9.8 × 10⁻⁴, respectively.
• 5.  
  Using Model TruncGauss (Model ZeroSubPop), we find that 99% of 1G+1G black holes are below 40M_⊙ (40M_⊙) and that 99% of all black holes are below 67M_⊙ (66M_⊙).

Since we do not believe that all BBHs come from a single type of cluster, we also consider a range of other typical cluster sizes, demonstrating that results depend upon the assumed escape velocity. For a cluster model with M_c = 10⁶ M_⊙ and r_c = 0.1 pc, which has the highest Bayes factor:

• 1.  
  We overwhelmingly favor models including hierarchical formation channels. For Model TruncGauss and Model ZeroSubPop, we find Bayes factors of 1400:1 and 25000:1, respectively, in favor of including hierarchical components.
• 2.  
  At least one BBH in GWTC-2 contains a 2G black hole with probability >99.99% for both Model TruncGauss and Model ZeroSubPop.
• 3.  
  GW190521 is most likely of 2G+2G origin, with 1200:1 and 700:1 odds in favor of being 2G+2G versus 1G+1G using Model TruncGauss and Model ZeroSubPop, with both models favoring 2G+2G over 1G+2G origin at ∼3.5:1.
• 4.  
  We find that GW190519, GW190602, GW190620, and GW190706 are most likely of 1G+2G origin for both Model TruncGauss and Model ZeroSubPop, favored over 1G+1G origin at >10:1 odds, while GW190517 is favored to be of 1G+2G origin above this threshold for Model ZeroSubPop.
• 5.  
  Using Model TruncGauss, the median relative merger rates of 1G+2G and 2G+2G to 1G+1G mergers are inferred to be 0.15 and 0.01, respectively, with 99% upper limits of 0.29 and 0.04. While using Model ZeroSubPop, the relative rates drop slightly to 0.14 and 9.2 × 10⁻³, with 99% upper limits of 0.25 and 0.03, respectively.

Our analysis indicates that there are plausible hierarchical merger candidates in GWTC-2, meriting further study.

There are a number of possible extensions to this analysis. Most importantly, we have assumed that all merging binaries are formed dynamically in clusters with a specific mass and density. While illustrative, this is unrealistic as (i) the observed BBH population may come from a mixture of formation channels including isolated field evolution, and (ii) dynamically formed binaries may occur in a wide range cluster types ranging from young open clusters to nuclear star clusters. An excess of events with aligned spin in GWTC-2 suggests that at least some binaries are assembled in the field (Abbott et al. 2021b), and comparisons of observations with model predictions indicate that a mix of formation channels is probable (Bouffanais et al. 2021; Wong et al. 2021; Zevin et al. 2021). The potential for multiple formation channels could be accounted for by including an additional mixture model for dynamically formed binaries versus those formed in isolation (Kimball et al. 2020b). Previous analyses have suggested using the distribution of spin orientations or eccentricities to measure the fraction of binaries formed dynamically (Breivik et al. 2016; Nishizawa et al. 2016; Stevenson et al. 2017; Talbot & Thrane 2017; Vitale et al. 2017; Lower et al. 2018; Gondán & Kocsis 2019; Romero-Shaw et al. 2019; Abbott et al. 2021b; Zevin et al. 2021). Relaxing the assumption that all dynamically formed binaries form in identical environments requires a model for the distribution of globular cluster properties and other dense environments, e.g., AGNs. It is possible that future GW observations will allow us to directly probe the distribution of cluster masses if we obtain sufficient observations to reconstruct the population of host environments. We leave incorporating these extensions to future work.

The authors thank Kyle Kremer, Carl Rodriguez, Mario Spera, and Zoheyr Doctor for their expert advice in constructing this study, and Isobel Romero-Shaw for comments on a draft manuscript. This research has made use of data obtained from the Gravitational Wave Open Science Center (www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the US National Science Foundation (NSF). Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes. This work is supported by the NSF Grant PHY-1607709 and through the Australian Research Council (ARC) Centre of Excellence CE170100004. C.K. is supported supported by the National Science Foundation under grant DGE-1450006. C.P.L.B. is supported by the CIERA Board of Visitors Research Professorship. M.Z. is supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51474.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. E.T. is supported through ARC Future Fellowship FT150100281 and CE170100004. T.D. acknowledges support from the María de Maeztu Unit of Excellence MDM-2016-0692, by Xunta de Galicia under project ED431C 2017/07, by Consellería de Educacíon, Universidade e Formacíon Profesional as Centros de Investigacíon do Sistema universitario de Galicia (ED431G 2019/05), and by FEDER. This research was supported in part through the computational resources from the Grail computing cluster at Northwestern University—funded through NSF PHY-1726951—and staff contributions provided for the Quest high performance computing facility at Northwestern University, which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by NSF Grants PHY-0757058 and PHY-0823459. This document has been assigned LIGO document number LIGO-P2000466.

Here, we include the full sets of inferred population hyperparameter posteriors for Model TruncGauss and Model ZeroSubPop. The 1G+1G primary mass distribution consists of two components. The first is a truncated power law with minimum mass m_min, maximum mass m_max, and power-law index α. The second is a Gaussian component with mean μ_m and standard deviation σ_m . The mixing parameter λ_m gives the fraction of BHs drawn from the Gaussian component. The mass ratio distribution is governed by a power law with index β_q . The 1G+1G spin distributions for Model TruncGauss are modeled as truncated Gaussians, with standard deviation σ_χ ∈ [0.1,10] and mean μ_χ ∈ [ − 3σ_χ , 1 + 3σ_χ ]. For Model ZeroSubPop, we use a Beta distribution with shape parameters α_χ > 1 and β_χ > 1 and allow a fraction λ₀ of the population to have spins drawn from a delta function at zero:

$$\begin{eqnarray}&&\begin{array}{l}\pi (\chi | {\lambda }_{0},{\alpha }_{\chi },{\beta }_{\chi },1{\rm{G}}+1{\rm{G}}\,)={\lambda }_{0}\delta (\chi )\\ \quad +\,(1-{\lambda }_{0}){\rm{B}}(\chi | {\alpha }_{\chi },{\beta }_{\chi }).\end{array}\end{eqnarray} \tag{ A1 }$$

This zero-spin subpopulation is inspired by simulations of massive stars with efficient angular momentum transfer, where black holes in effective isolation would be born with spins of ∼0.01 (Qin et al. 2018; Fuller & Ma 2019), and may also describe primordial black holes (De Luca et al. 2020). The 1G+2G and 2G+2G mass and spin distributions are obtained using the transfer functions defined in Kimball et al. (2020b).

In Figure 6, we plot the parameters governing the mass and mass ratio distributions. When using the astrophysically motivated prior on m_max (a Gaussian centered at 50 M_⊙ with standard deviation 10M_⊙), we mostly recover this prior; we do not yet have an informative enough catalog to measure this within our phenomenological model. As in Kimball et al. (2020b), we find that m_max is restricted at small values of the power-law index α where the mass distribution is flatter and more sensitive to the upper mass cut-off. We are able to place stronger constraints on the minimum black hole mass, finding ${m}_{\min }\lt 6.8$ M_⊙ and ${m}_{\min }\lt 7.0{M}_{\odot }$ at the 99% credible level using Model TruncGauss and Model ZeroSubPop, respectively. For Model TruncGauss, we find that the Gaussian component of the mass spectrum is well constrained to ${\mu }_{m}={32.0}_{-6.8}^{+8.5}{M}_{\odot }$. With Model ZeroSubPop where we infer support for lower relative hierarchical merger rates (Figure 3), we find a long tail at high masses when m_max is low, finding ${\mu }_{m}={31.5}_{-9.0}^{+23.0}{M}_{\odot }$. With both Model TruncGauss and Model ZeroSubPop we find a bimodality in the relative hierarchical merger rates, as seen in Figure 3. The peak at higher relative rates is associated with small values of σ_m . When we restrict the Gaussian component to peak sharply, it is unable to accommodate the highest-mass events as 1G+1G in its tail, and therefore they are fit as hierarchical mergers. If we omit the five highest-mass events in GWTC-2, this peak at higher relative hierarchical rates disappears. Overall, the inferred mass distributions are largely consistent between our two spin models.

Zoom In Zoom Out Reset image size

In Figures 7 and 8, we plot the parameters governing the component spin distributions for Model TruncGauss and Model ZeroSubPop, respectively. In Model TruncGauss, we find μ_χ < 0, and therefore preference for spin distributions that peak at 0. In Model ZeroSubPop, we prefer low values of α_χ , which increases support at low component spin, but find that β_χ is unconstrained. With Model ZeroSubPop, we also find that the fraction λ₀ of black holes originating from the zero-spin subpopulation (plotted in Figure 3) is constrained to be less than 0.06 (0.12) at the 90% (99%) credible level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We plot the posteriors for the population hyperparameters governing the mass distributions inferred when we assume a flat prior on m_max in Figure 9. Using the flat prior on m_max, we no longer constrain the maximum mass cut-off, and the posterior peaks at around ∼80M_⊙. For Model TruncGauss and Model ZeroSubPop, we constrain the mean of the Gaussian component to be ${\mu }_{m}={32.1}_{-6.1}^{+4.2}{M}_{\odot }$ and ${\mu }_{m}={31.6}_{-7.3}^{+4.5}{M}_{\odot }$, respectively. The preference for high m_max means that the high-mass tail on μ_m is no longer required to fit the more massive events in GWTC-2, even though the relative 1G+2G and 2G+2G versus 1G+1G merger rates (as well as the events-wise odds of hierarchical merger) drop by an order of magnitude under the flat prior on m_max.

Zoom In Zoom Out Reset image size

In Figures 10 and 11, we plot the posteriors of the population hyperparameters governing the component spin distributions inferred when we assume a flat prior on m_max. We find that the inferred spin distributions are consistent across choices of prior on m_max. For Model ZeroSubPop the fraction λ₀ of BBHs originating from the zero-spin subpopulation is constrained to be less than 0.04 (0.09) at the 90% (99%) credible level, and is still consistent with λ₀ = 0. We also find that α_χ and β_χ are less constrained; when we allow for high values of m_max, it is easier to explain higher mass systems as 1G+1G, and hence we do not require a low-spin population to enable high relative hierarchical merger rates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size