Evidence for hierarchical black hole mergers in the second LIGO--Virgo gravitational-wave catalog

We study the population properties of merging binary black holes in the second LIGO--Virgo Gravitational-Wave Transient Catalog assuming they were all formed dynamically in gravitationally bound clusters. Using a phenomenological population model, we infer the mass and spin distribution of first-generation black holes, while self-consistently accounting for hierarchical mergers. Considering a range of cluster masses, we see compelling evidence for hierarchical mergers in clusters with escape velocities $\gtrsim 100~\mathrm{km\,s^{-1}}$. For our most probable cluster mass, we find that the catalog contains at least one second-generation merger with $99\%$ credibility. We find that the hierarchical model is preferred over an alternative model with no hierarchical mergers (Bayes factor $\mathcal{B}>1400$) and that GW190521 is favored to contain two second-generation black holes with odds $\mathcal{O}>700$, and GW190519, GW190602, GW190620, and GW190706 are mixed-generation binaries with $\mathcal{O}>10$. However, our results depend strongly on the cluster escape velocity, with more modest evidence for hierarchical mergers when the escape velocity is $\lesssim 100~\mathrm{km\,s^{-1}}$. Assuming that all binary black holes are formed dynamically in globular clusters with escape velocities on the order of tens of $\mathrm{km\,s^{-1}}$, GW190519 and GW190521 are favored to include a second-generation black hole with odds $\mathcal{O}>1$. In this case, we find that $99\%$ of black holes from the inferred total population have masses that are less than $49\,M_{\odot}$, and that this constraint is robust to our choice of prior on the maximum black hole mass.

Among the GWTC-2 systems there are high-mass BBHs that have components with masses of 45M (Abbott et al. 2021), 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;Banerjee 2020;Fragione & Silk 2020;Mapelli et al. 2020;Fragione et al. 2020a); 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 cutoffs) for the first-generation (1G) black holes not formed through mergers (Doctor et al. 2020;Sedda et al. 2020;Tiwari & Fairhurst 2021;Kimball et al. 2020a).
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. 2020d). 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 Sec. 2 we review the key components of our population inference framework; the results of this are given in Sec. 3, with additional description of the population hyperparameters in Appendix A, and we discuss our findings in Sec. 4.

METHODS
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 Fig. 1 and Fig. 2. The 1G+1G model is nearly identical to population models used in Abbott et al. (2020d): 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 Appendix A.
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 5 M and r c = 1 pc, corresponding to a central escape velocity of ∼ 65 km s −1 We assume that the relative merger rates scale as R 1G+2G /R 1G+1G ∝ F ret , R 2G+2G /R 1G+1G ∝ F 2 ret , 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. 2020d). 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. (2021). 1 For the new GWTC-2 events we use results calculated with the IMRPhenomD and IMRPhenomPv2 waveforms (Khan et al. 2016;Hannam et al. 2014). We generate posterior samples for population hyperparameters Λ using the nested sampler dynesty (Speagle 2020) using the GWPopulation framework (Talbot et al. 2019), which takes ad-  Abbott et al. (2020d), we infer population hyperparameters for our mass and spin models (Fig. 6, Fig. 7, and Fig. 8 in Appendix A). For both Model TruncGauss and Model ZeroSubPop, we find that including the 1G+2G and 2G+2G population components is preferred, finding Bayes factors of 5 and 7, respectively, in favor of including versus excluding the hierarchical components.
In our inferred 1G mass distribution, the mean of the Gaussian mass component is well constrained to µ m = 32.0 +8.5 −6.8 M and µ m = 31.5 +23.0 −9.0 M for Model TruncGauss and Model ZeroSubPop, respectively. For both models we recover our prior on the maximum mass cut-off m max . In Fig. 1 and Fig. 2, we plot the posterior predictive distributions for the 1G+1G, 1G+2G, and 2G+2G populations. Using Model TruncGauss (Model ZeroSubPop), we find that 99% of 1G+1G black holes are less than 47M (47M ), and 99% of black holes in the total population are less than 49M (48M ), consistent with the results of Kimball et al. (2020b). These upper limits are lower than those found for the Power Law + Peak model in Abbott et al. (2020d), but that model does not include a high-mass hierarchical component, and requires a flatter power law to fit the heavier black holes in GWTC-2. Relaxing the prior on the maximum mass cut-off to a uniform prior out to 100M (Fig. 9, Fig. 10 and Fig. 11), we do not obtain stringent constraints on m max , but find that it peaks around ∼ 80M . In this case, we find that 99% of 1G+1G black holes are less than 49M (49M ), and 99% of black holes in the total population are less than 51M (50M ).
Using Model TruncGauss and Model ZeroSub-Pop, 90% of 1G+1G black holes have spins less than 0.65 and 0.50, respectively. With Model ZeroSub-Pop, the fraction λ 0 of BBHs originating from the zero spin channel is constrained to be less than 0.12 at the 99% credible level.

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 ZeroSub-Pop, which results in lower typical recoil velocities and hence higher relative hierarchical merger rates. In Fig. 3, we plot the posteriors for these rates, as well as for the fraction λ 0 of 1G+1G black holes in Model ZeroSub-Pop with zero spin. Using Model TruncGauss, we infer median relative rates R 1G+2G /R 1G+1G = 5.8 × 10 −3 and R 2G+2G /R 1G+1G = 1.7×10 −5 with 99% upper limits of 0.12 and 7.6 × 10 −3 , respectively. Taking into account selection effects, the detected population would have median relative rates of R det 1G+2G /R det 1G+1G = 0.01 and R 2G+2G /R 1G+1G = 7.0×10 −5 with 99% upper limits of 0.23 and 0.03, respectively. Using Model Ze-roSubPop, these rates become 5.3 × 10 −3 (0.01) and 1.4 × 10 −5 (5.7 × 10 −5 ), with 99% upper limits of 0.04 (0.08) and 9.8 × 10 −4 (4.0 × 10 −3 ), 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).

Posterior odds for hierarchical origin
For each event in GWTC-2, we calculate the posterior odds 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 Ze-roSubPop in Fig. 4.
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 ori- Figure 1. Posterior predictive distributions for the primary mass m1 and mass ratio q. The solid, dashed, and dotted lines indicate the 1G+1G, 1G+2G, and 2G+2G distribution, respectively. In blue, we plot the distributions inferred when modeling the 1G spin distribution as a non-singular Beta distribution together with a delta function at zero. In orange, we plot the distributions inferred when using a truncated Gaussian.

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. Figure 2. Posterior predictive distributions for the component black hole spins. The solid, dashed and dotted lines indicate the 1G+1G, 1G+2G, and 2G+2G distribution, respectively. In blue, we plot the distributions inferred when modeling the 1G spin distribution as a non-singular Beta distribution together with a delta function at zero. In orange, we plot the distributions inferred when using a truncated Gaussian. Figure 3. Posteriors of the inferred branching ratios, and the fraction λ0 of 1G+1G black holes with zero spin for Model ZeroSubPop. The branching ratios give the relative 1G+2G versus 1G+1G and 2G+2G versus 1G+1G merger rates. We plot the results using Model TruncGauss and Model ZeroSubPop in orange and blue, respectively.
In Fig. 5, we show results when considering models with Plummer masses 10 4 -10 9 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 −1 ), 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 −1 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%.
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 −1 ), with the highest Bayes factors corresponding to models where the central escape velocities are ∼ 300 km s −1 . In Fig. 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 −1 , 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 −3 ) 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) Figure 4. Odds of events in GWTC-2 having 1G+2G versus 1G+1G origin, as a function of the inferred median primary black hole mass, mass ratio, and primary black hole spin. The results using Model TruncGauss and Model ZeroSubPop are plotted on the left and right, respectively. The gray vertical lines are draws from the corresponding inferred posterior over the maximum 1G black hole mass. and 0.18 (0.13). When V esc ∼ 300 km s −1 , 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.

CONCLUSIONS
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 pairinstability 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 50M 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 5 M and r c = 1 pc: 1. The 44 events in GWTC-2 are best modelled when allowing for hierarchical formation channels. For Model TruncGauss and Model ZeroSub-Pop, 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 which 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.
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 6 M and r c = 0.1 pc, which has the highest Bayes factor: 1. We overwhelmingly favor models including hierarchical formation channels. For Model Trunc-Gauss 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.
4. We find that GW190519, GW190602, GW190620, and GW190706 are most likely of 1G+2G origin for both Model TruncGauss and Model Ze-roSubPop, favored over 1G+1G origin at > 10:1 odds, while GW190517 is favored to be of 1G+2G origin above this threshold for Model ZeroSub-Pop.

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 −3 , 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. 2020d), and comparisons of observations with model predictions indicate that a mix of formation channels is probable (Wong et al. 2021;Zevin et al. 2021;Bouffanais 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 (Vitale et al. 2017;Nishizawa et al. 2016;Breivik et al. 2016;Stevenson et al. 2017;Talbot & Thrane 2017;Gondán & Kocsis 2019;Lower et al. 2018;Romero-Shaw et al. 2019;Abbott et al. 2020d;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. 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 λ 0 of the population to have spins drawn from a delta function at zero: 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 Fig. 6, we plot the parameters governing the mass and mass ratio distributions. When using the astrophysicallymotivated prior on m max (a Gaussian centered at 50M 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 < 6.8 M and m min < 7.0M 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 µ m = 32.0 +8.5 −6.8 M . With Model ZeroSubPop, where we infer support for lower relative hierarchical merger rates (Fig. 3), we find a long tail at high masses when m max is low, finding µ m = 31.5 +23.0 −9.0 M . With both Model TruncGauss and Model ZeroSubPop, we find a bimodality in the relative hierarchical merger rates as seen in Fig. 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 they are therefore 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.
In Fig. 7 and Fig. 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 λ 0 of black holes originating from the zero-spin subpopulation (plotted in Fig. 3) is constrained to be less than 0.06 (0.12) at the 90% (99%) credible level.
We plot the posteriors for the population hyperparameters governing the mass distributions inferred when we assume a flat prior on m max in Fig. 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 µ m = 32.1 +4.2 −6.1 M and µ m = 31.6 +4.5 −7.3 M , 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 .
In Fig. 10 and Fig. 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 λ 0 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 = 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.