Understanding how fast black holes spin by analysing data from the second gravitational-wave catalogue

The Advanced LIGO and Virgo detectors have now observed approximately 50 black-hole-binary mergers, from which we can begin to infer how rapidly astrophysical black holes spin. The LIGO-Virgo Collaboration (LVC) analysis of detections up to the end of the first half of the third observing run (O3a) appeared to uncover a distribution of spin magnitudes that peaks at $\sim$0.2. This is surprising: is there a black-hole formation mechanism that prefers a particular, non-zero spin magnitude, or could this be the cumulative effect of multiple formation processes? We perform an independent analysis of the most recent gravitational-wave catalogue, and find that (a) the support for the LVC spin-magnitude is tenuous; in particular, adding or removing just one signal from the catalogue can remove the statistical preference for this distribution, and (b) we find potential evidence for two spin sub-populations in the observed black holes; one with extremely low spins and one with larger spin magnitudes. We make the connection that these spin sub-populations could be correlated with the mass of the binary, with more massive binaries preferring larger spin magnitudes, and argue that this may provide evidence for hierarchical mergers in the second gravitational-wave catalogue.


INTRODUCTION
During the first and second gravitational-wave (GW) observing runs (O1 and O2) (Abbott et al. 2019c) of the Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2014) GW observatories, the LIGO scientific and Virgo collaborations announced eleven GW candidates; ten from merging black-hole binaries (BBHs) (Abbott et al. 2016c(Abbott et al. ,b, 2017a(Abbott et al. ,b,c, 2019c and one from a binary neutron star coalescence (Abbott et al. 2017d). Independent groups also reported additional GW candidates (Nitz et al. 2020;Venumadhav et al. 2020;Zackay et al. 2019b,a). By combining parameter estimates for these BBH observations, first attempts at deciphering the astrophysical distribution of black hole spins were conducted Farr et al. 2018;Tiwari et al. 2018;Abbott et al. 2019a;Fairhurst et al. 2020a;Biscoveanu et al. 2020;García-Bellido et al. 2021). The limited sample size meant that only weak constraints could be placed on the distribution of black hole spins, although even with only a handful of signals it was shown that high spin magnitudes were strongly disfavoured, and there was some evidence for spin mis-alignment in binaries.
Precise measurements of both the spin magnitudes and their mis-alignment with the binary's orbital angular momentum rely on identifying the presence (or lack) of the General Relativistic phenomenon of spin-induced orbital precession -the misalignment of the binary's orbital angular momentum and the spins of each compact object resulting in characteristic modulations to the GWs amplitude and phase (Apostolatos et al. 1994). A direct measurement of spin-induced orbital precession would then provide a unique insight into the astrophysical distribution of black hole spins (e.g. Gerosa et al. 2018;Rodriguez et al. 2016). During the first half of the third GW observing run (O3a), a further 39 GW candidates were announced (Abbott et al. 2020a). However, similar to those from O1 and O2 (Abbott et al. 2019b), most of these detections remained largely uninformative about the presence of precession (Abbott et al. 2020a). Abbott et al. (2021b) recently showed, through a hierarchical Bayesian analysis (Thrane & Talbot 2019;Mandel et al. 2019;Vitale 2020), that there is clear evidence that the population of known BBHs includes misaligned spins, despite no single event unambiguously exhibiting evidence for precession. By assuming a population model where the spin magnitude of each black hole is described by a beta function  and the orientation by a model allowing for both isotropic and aligned spins (Talbot & Thrane 2017), it was shown that the most likely spin distribution has a peak in the spin magnitude at ∼ 0.2, with preference for primarily aligned spins (although there is non-vanishing support for angles > 90 • indicating the presence of misaligned component spins). Roulet et al. (2021) and Galaudage et al. (2021) later challenged this point of view, showing that all binaries observed to date are consistent with two spin populations: one with negligible black hole spin and a second with spins preferentially aligned with the orbital angular momentum. Galaudage et al. (2021) estimated that 70 -90% of merging binaries contain black holes with negligible spin and the black holes in the second sub-population have spins ∼ 0.5 with orientations preferentially (but not exactly) aligned to the orbital angular momentum. Callister et al. (2021) also searched for correlations within the black hole spin distribution and found that the binaries mass ratio is correlated with the black hole spin at 98.7% credibility with more unequal mass binaries exhibiting systematically larger black hole spin.
In this paper we analyse the results from Abbott et al. (2021b) and link our findings to the conclusions found from other works. We opt to extend previous analyses (see e.g. Farr et al. 2017;Tiwari et al. 2018) by including the effects of spin-induced orbital precession and perform a detailed model selection analysis for two reasons: firstly we are able to provide odds ratios to show by how much one model is preferred to another, and secondly we are able to "open the black box" of the results presented in Abbott et al. (2021b) by showing how each GW candidate contributes to the final result.
Since there are still only a limited number of binary black hole mergers available, which makes inferring the true black hole spin distribution challenging, we select a series of spin models which aim to identify the core features of the spin distribution. We choose spin models with isotropic or preferentially aligned spins with several choices for the distribution of spin magnitudes and compare our results to the most likely distribution from Abbott et al. (2020a). We test the robustness of our results by repeating the analysis using two parameterisations of precession.
We show that if we use the same set of BBHs, we obtain the same conclusion as Abbott et al. (2021b) irrespective of which parameterisation of precession we use, i.e., the current population of binary black holes prefers a spin distribution model with mild preference for aligned spins and spin magnitudes peaking at ∼ 0.2. However, the preference for this distribution disappears if we include GW190814 in the analysis, and can be disfavoured by as much as 14:1 if GW190517 055101 is removed.
Secondly, we show that our results agree with those presented in Roulet et al. (2021) and Galaudage et al. (2021) and show that there is evidence that the black hole population is consistent with two spin subpopulations, which, when combined, gives an overall preference for the distribution in Abbott et al. (2021b). We highlight that the majority of events (∼ 80%) prefer a spin distribution with extremely low spin magnitudes while several show preference for larger spins. We demonstrate that these two spin sub-populations could be correlated with the mass of the system with low mass events preferring extremely low spin magnitudes and high mass events preferring larger spins. Of the models that we use, the heavier binaries prefer isotropic rather than aligned spin orientations. This paper is structured as follows: in Section 2 we briefly summarize the details behind a model selection analysis and document which GW candidates and spin distribution models are used in this study. In Section 3 we describe our results and provide an understanding into how fast black holes spin. We then conclude in Section 4.

METHOD
We use Bayesian model selection to calculate the odds ratio between nine different spin distributions using the publicly released posterior samples from GWTC-2, made available through the Gravitational Wave Open Science Center (GWOSC; Vallisneri et al. 2015;Abbott et al. 2021a). As in Abbott et al. (2021b), we only consider BBHs with false alarm rates (FARs) < 1yr −1 . This means we exclude 2 marginal events included in Abbott et al. (2020a): GW190719 215514 and GW190909 114149 (we consider GW190426 152155 to be a marginal neutron star black hole candidate). Unlike Abbott et al. (2021b) where GW190814 was excluded in their spin distribution analysis, we consider how our results change if GW190814 is included in the population. This is because GW190814 is most likely (71%) the result of a BBH merger (Abbott et al. 2020d). For each GW candidate considered, we randomly draw 10 4 samples from the 'PublicationSamples' dataset to ensure a consistent number of samples across the population. For all candidates this meant that we used posterior samples that were generated using waveform models that at least included precession (Hannam et al. 2014;Ossokine et al. 2020), and for the majority of candidates this meant that we used posterior samples that were generated using waveform models that also included subdominant multipole moments (Khan et al. 2020;Ossokine et al. 2020;Varma et al. 2019); see Table  VIII in Abbott et al. (2020a). Unless otherwise stated, the odds ratio is compared to the mean distribution displayed in Fig. 10 of Abbott et al. (2021b), denoted in this work as the "LVC" distribution, and we do not incorporate the calculated uncertainty on the inferred LVC spin distribution. Given that population studies of showing support across the Left: χp-χ eff and Right: ρp/ρχ eff parameter space for a selection of the different spin distributions used in our analysis -LVC refers to the mean inferred spin distribution from Abbott et al. (2021b), ELI, VLI and LI are all distributions with isotropic spin orientations with either extremely low, very low or low spin magnitudes respectively and LA is a low spin magnitude distribution with spins nearly aligned with the orbital angular momentum (maximum misalignment angles of 30 • ), see text for definitions. Each contour encloses 90% of the distribution. We show the preferred model from Abbott et al. (2021b) in green. Gaussian kernel density estimates are used to estimate the probability density. The ρp/ρ marginalized PDF for the LA spin distribution is truncated to allow for distributions of other models to be seen.
BBHs from O1, O2 and O3a disfavour highly spinning black holes Tiwari et al. 2018;Fairhurst et al. 2020a;Abbott et al. 2021b), we consider only spin distributions with either low (L) (consistent with Farr et al. 2017), very low (VL) (consistent with Abbott et al. 2019b) or extremely low (EL) spin magnitudes, where a = |cS/(Gm 2 )| is the spin magnitude of a black hole with mass m and spin angular momentum S. We consider two distributions for the tilt angles: nearly aligned (A) and isotropic (I). The nearly aligned distribution is triangular in cos θ, with a peak at 1 and taking values between 0.85 ≤ cos θ ≤ 1. The nearly aligned distribution resembles field binaries -binaries which are formed from isolated stellar progenitors and are expected to have spins distributed about the orbital angular momentum with some unknown misalignment angle (e.g. Kalogera 2000;Mandel & O'Shaughnessy 2010;Gerosa et al. 2018) -with maximum misalignment angles of 30 • . The isotropic distribution is uniform in cos θ between −1 and 1 and resembles dynamic binaries -binaries which are expected to have randomly orientated spins and formed when two black holes become gravitationally bound in dense stellar environments (e.g. Rodriguez et al. 2016).
As it is difficult to constrain the individual black hole spins at typical signal-to-noise ratios (SNRs) (Pürrer et al. 2016), we use a mass weighted effective spin χ eff to describe the average projection of spins parallel to the orbital angular momentum (Ajith et al. 2011), and we use two competing quantities to describe the projection of spins perpendicular to the orbital angular momentum: χ p (Schmidt et al. 2015) and ρ p (Fairhurst et al. 2020b). χ p is widely used for inferring the occurrence of precession in GW data (see e.g. Abbott et al. 2019cAbbott et al. , 2020a, although alternative metrics have also been proposed (e.g. Gerosa et al. 2021;Thomas et al. 2021), and ranges between 0 (no precession of the orbital plane) and 1 (maximal precesion). ρ p is the precession SNR and describes the contribution to the total SNR of the system that can be attributed to precession. It has been shown previously that ρ p is a useful quantity for inferring the presence of precession in population studies (Fairhurst et al. 2020a). ρ p is calculated by de-composing a GW into two non-precessing harmonics and isolating the SNR contained in the harmonic orthogonal to the dominant one. By deconstructing a precessing gravitational-wave in this form, the characteristic amplitude and phase modulations can be interpreted as the beating of these harmonics. If ρ p is small (ρ p 2), the amplitude of the second harmonic is insignificant and any beating of the harmonics is negligible. For this case, we would observe a GW that looks like the dominant non-precessing harmonic. We choose to perform two analyses, one describing where the black hole spins are described by χ eff and χ p and another with χ eff and ρ p . This is because it is notoriously difficult to accurately extract precession information from gravitational-wave signals, especially at relatively low SNRs (see e.g. Although ρ p is dependent on the GW detector network and its sensitivity, if the harmonics are close to orthogonal, ρ p can be scaled by the total SNR ρ to provide a detector invariant quantity ρ p /ρ (see Eq. 39 in Fairhurst et al. (2020b), noting that ρ is denoted as ρ 2harm ). We therefore choose to parameterise precession in this work by ρ p /ρ since it allows for results that are independent of the detector network and the chosen detector sensitivity. This implies that the ρ p 2 criterion used in previous works for quantifying the measurability of precession (see e.g. Abbott et al. 2020a,b,d;Fairhurst et al. 2020b,a;Green et al. 2021) becomes ρ p /ρ 2/ρ, which is bounded between 0 ≤ ρ p /ρ ≤ 1/ √ 2, where the upper bound implies equal power in both harmonics. Consequently, for systems with large ρ p /ρ, precession contributes significantly to the total SNR of the system. Figure 1 shows how a subset of these spin distributions vary across the χ p -χ eff and ρ p /ρ-χ eff parameter space. The aligned distributions can easily be distinguished from the isotropic distributions as χ eff > 0 and ρ p /ρ is small by definition.
Following the methodology described in Tiwari et al. (2018), we calculate the odds ratio between two spin distribution models λ 1 and λ 2 as, where p(λ|{d}) is the posterior distribution for the model λ given a set of BBH observations {d}, V pop (λ) is the sensitive volume for the model λ, θ j i is the jth posterior sample for observation i, j p(θ j i |λ)/π(θ j i ) is the sum over posterior samples re-weighted from the default prior universe used in the LVC analyses π(θ j i ) (LAL prior), to the universe assuming a given model λ p(θ j i |λ), p(λ) is the prior on the model and we restrict θ = (M, q, χ eff , [χ p , ρ p /ρ], ι). As with Tiwari et al. (2018), we assume all models are equally likely i.e., p(λ) = 1, ∀λ. In order to evaluate p(θ j i |λ), we generate a universe for each spin distribution model consisting of 10 7 randomly chosen binaries and compute a five dimensional Gaussian Kernel Density Estimate (KDE).
For LVC parameter estimation analyses, the prior on m 1 and m 2 is taken to be flat and spin vectors are assumed to be uniform in spin magnitude and isotropic on the sphere (see Appendix B.1 of Abbott et al. 2019c). When generating a universe for a given model λ, we use an identical mass distribution but vary the spin magnitude and orientation vectors (see Sec. 3.1 for an analysis that uses an astrophysical mass distribution). All other binary parameters are randomly drawn from the same distributions as used in Abbott et al. (2019c).
The sensitive volume V pop (λ) is essential for accounting for selection effects. It is estimated numerically by injecting GW signals drawn from model λ into GW strain data and searching for them assuming a given detection threshold (Tiwari 2018). Currently search pipelines employ non-precessing waveform approximants for matched filtering (Usman et al. 2016;Messick et al. 2017). This means that current techniques to estimate the sensitive volume omit precession (although see Gerosa et al. (2020a), which suggests an alternative method that includes precession). Since precessing signals will be recovered at lower probabilities than an equivalent precessing search pipeline, we can expect that the sensitive volume will be underestimated for systems where precession effects are observable (Calderón Bustillo et al. 2017). However, Fairhurst et al. (2020b) argue that for signals with low ρ p this effect is minimal. Given that for most models used in this paper ρ p /ρ is small, we approximate V pop (λ) by V pop (λ np ): the sensitive volume for the non-precessing equivalent λ.
V pop (λ np ) as a function of each spin distribution model shows the trends we expect: we find that as the spin magnitude increases and spin orientation becomes more aligned, the sensitive volume increases. This is expected since binaries with a larger aligned spin (larger χ eff ) can be observed at a greater distance (Ajith et al. 2011). Since the ELI spin distribution leads to a population with the smallest aligned spin, it has the lowest sensitive volume. The odds ratio calculation in Eq. 2 involves dividing by the model's sensitive volume. This means that for cases where the sum over re-weighted posterior samples ( j p(θ j i |λ)/π(θ j i )) is similar between two spin models, the model with the lower sensitive volume is preferred. Although for a single event this effect may be minuscule, for 44 observations the odds ratio may increase substantially, e.g., the smaller sensitive volume for ELI over LVC contributes a factor ∼ 10 to the Bayes factor calculation. Figure 2 shows odds ratios with respect to the LVC spin distribution for two different metrics for precession: χ p and ρ p /ρ. Models that are not shown have odds ratios < 10 −5 : 1. Our analysis shows that if we assume that all black holes originate from the same population, models with aligned spins are strongly disfavoured (> 10 7 : 1 in favour of isotropic spins) as well as "low" spin magnitudes (> 10 3 : 1 in favour of "very low" or "extremely low" magnitudes). Only spin distributions with isotropic spin orientations and very low, or extremely low spin magnitudes are broadly consistent with the observations. However, the details of which spin distribution best describes the population depends upon the detailed choice of which events are included in the analysis. Specifically, when GW190814 is excluded from our analysis, of the distributions considered, the LVC spin distribution is marginally preferred. Our analysis infers that when GW190814 is excluded, both VLI and ELI are disfavoured with odds ratios 2.8 : 1 and 10 1 10 0 10 1

RESULTS AND DISCUSSION
Odds ratio (ELI / LVC) Figure 3. Odds ratios for ELI against LVC for each binary black hole candidate considered in this analysis (Abbott et al. 2020a). An orange line shows an odds ratio of 1 meaning that neither model is preferred. An odds ratio greater than 1 shows preference for ELI over LVC. In both cases odds ratios are calculated using two different paramerisations of precession: χp (grey) as used in Abbott et al. (2020a) and ρp/ρ (red). Odds ratios are calculated using the posterior samples released as part of GWTC-2 (Abbott et al. 2020a;Vallisneri et al. 2015;Abbott et al. 2021a).
3.8 : 1 for the χ p analysis and 5.9 : 1 and 5.5 : 1 for the ρ p /ρ analyses respectively. This equates to around a 1σ preference for LVC over ELI and 1.4σ for VLI. When GW190814 is included in the population, the extremely low isotropic spin distribution and LVC are equally preferred by the data, with odds ratios ∼ 1 (1.2 : 1 and 0.9 : 1 for the χ p and ρ p /ρ analyses respectively) while the very low isotropic model is only slightly dis-favoured. This is the first main result from our study. The preference for the LVC, extremely low or very low isotropic spin distribution depends sensitively on which signals are included in the analysis. All of these results assume a single population; we will consider multiple subpopulations in Section 3.2. We now consider these results in more detail.
In Figure 3 we show how the odds ratio of ELI vs LVC changes as a function of GW candidate. For the majority of events (32/44) the extremely low spin magnitude distribution is preferred over LVC with 31/44 events having odds ratio between 1 and 3. GW190814 exhibits the strongest preference for extremely low spin magnitudes with an odds ratio 5 : 1. On the other hand, there are a handful of events which exhibit a strong preference for larger spins, most notably GW151226, GW190412 and GW190517 055101. It is instructive to examine these events in more detail.
GW190517 055101 has the largest χ eff observed so far with χ eff = 0.52 +0.19 −0.19 . This leads to the strongest preference for the LVC distribution among the observed events. As can be seen in Figure 1 there is significantly greater support for large, positive χ eff in the LVC distribution than ELI: the majority of the LVC distribution supports χ eff > 0 (70% compared to 50% for ELI) and has a longer tail up to larger χ eff (∼ 0.32 compared to ∼ 0.1 for ELI). If this single event is removed from the analysis, ELI is preferred over LVC with odds ratios 14 : 1 and 11 : 1 when GW190814 is included and 3 : 1 and 2 : 1 when GW190814 is excluded from the analysis for the χ p and ρ p /ρ analyses respectively. Meanwhile, GW190412 supports χ eff > 0.15 at 90% confidence and is consistent with a mildly precessing system, χ p = 0.31 +0.19 −0.16 . This event contributes a factor of 5 to the odds ratio in favour of the LVC distribution. GW190412's spin has been discussed at length in previous work (see e.g. Mandel & Fragos 2020;Zevin et al. 2020). Similarly GW151226 has χ eff = 0.18 +0.20 −0.12 and again contributes significantly to the preference for larger spins. Indeed, all candidates that support χ eff > 0 at more than 90% probability (GW151226, GW170729, GW190412, GW190517 055101, GW190519 153544, GW190620 030421, GW170706 222641, GW190720 000836, GW190728 064510, GW190828 063405, GW190930 133541) show a preference for the LVC distribution.
The only event with a large odds ratio in favour of the ELI distribution is GW190814. Uniquely among the events considered, due to the large SNR and the unambiguous identification of higher harmonics, for this binary both the χ eff and χ p were constrained to be close to zero (Abbott et al. 2020d). As a result it is no surprise that ELI is preferred. In this particular region of parameter space ELI has ∼ 7× more support than LVC. Since comparing support in a given region of the parameter space is effectively computing a simplified version of Eq. 2, we expect this calculation to be indicative of the odds ratio. It is therefore a good sanity check for our results.
Next we consider the robustness of our results. We see from Figure 2 that we obtain the same conclusions when repeating the analysis using two different metrics of precession: χ p and ρ p /ρ. As can be seen in Figure 3, the largest difference between these analyses is for GW190521 (Abbott et al. 2020c,e) where the χ p and ρ p /ρ analyses prefer LVC over ELI by 1.2 : 1 and 1.4 : 1 respectively. GW190521 is consistent with a merger of two black holes with masses 85 +21 −14 M and 66 +17 −18 M , with effective spins χ p = 0.68 +0.26 −0.44 and χ eff = 0.08 +0.27 −0.36 . Owing to the large total mass, GW190521 is very short in duration, with only 4 cycles (2 orbits) within the sensitive frequency band of the GW observatories. In contrast to the χ p measurement, the inferred ρ p /ρ demonstrates a lack of measurable precession, ρ p /ρ = 0.16 +0.33 −0.13 . This is because the short signal implies almost degenerate non-precessing harmonics (with overlap 0.97 +0.01 −0.03 ) which leads to a near-zero SNR orthogonal to the dominant harmonic. This difference is explored in detail in . Consequently, the χ p analysis infers that GW190521 is consistent with a much larger spin magnitude distribution than the ρ p /ρ analysis.
Next we comment briefly on the difference between the χ p and ρ p /ρ analyses on a population level. From Figure 2 we see that the difference in analyses becomes larger for larger spin magnitudes with a 0.4σ, 0.7σ and 1.7σ difference between the χ p and ρ p /ρ analyses for the ELI, VLI and LI spin distributions respectively. This is expected since χ p and ρ p are more likely to give alternative descriptions regarding the presence of precession in a gravitational-wave signal at larger rather than lower spin magnitudes, see e.g. GW190521.

Reweighting to an astrophysical mass distribution
Up until now, all results have used a flat in m 1 and m 2 (with the condition that m 1 > m 2 ) mass distri-   (Abbott et al. 2020a;Vallisneri et al. 2015;Abbott et al. 2021a). Light grey and red traces show the posterior distributions for events which prefer ELI over LVC and LVC over ELI respectively (see Figure 3). Solid black and red curves shows the average of the light grey and red traces respectively. The orange curves show the default χ eff and χp priors used in the LVC analyses.
bution. Here, we investigate how the conclusions vary when the posterior samples are reweighted to an astrophysical mass distribution. We select a mass distribution where the probability of drawing the primary mass Odds ratios for LA (black) and LI (red) against ELI for each binary black hole candidate considered in this analysis. An orange line shows an odds ratio of 1 meaning that neither model is preferred. An odds ratio less than 1 shows preference for ELI. In both cases the quoted odds ratios are an average of the χp and ρp/ρ analyses. Candidates below the purple line have primary mass m1 > 50 M . follows a simple power law p(m 1 ) ∝ m −α 1 and the probability of drawing the secondary mass is uniform between 5M and m 1 Abbott et al. 2016a). For each α ∈ [1.5, 2.0, 2.35, 3.0] we repeat the analysis above and identify which spin distribution is preferred.
When reweighting the posterior samples to an astrophysical mass distribution, we expect to see a preference for lower effective spins. As explained in detail in Tiwari et al. (2018), this is because close to equal mass ratio binaries are significantly more likely, which, owing to known degeneracy between the mass ratio and aligned-spin (Cutler & Flanagan 1994;Poisson & Will 1995;Baird et al. 2013), leads to a preference for lower effective spins.
As expected, we find that when reweighting to an astrophysical mass prior, the ELI and VLI spin distributions describe the data well with odds ratios approaching unity. When GW190814 is included in the analysis, both ELI and VLI are preferred with odds ratios ∼ 4 : 1 and ∼ 2 : 1 respectively. We also draw similar conclusions to Tiwari et al. (2018) and find that the result is only mildly dependent on the chosen value of α. This demonstrates that even when reweighting to an astrophysically motivated mass prior, there is still no strong evidence to prefer one spin model over another (when GW190814 is excluded from the analysis), and that the unknown mass distribution of black holes does not cause a significant effect on the inferred spin distribution.

Structure in the preferred spin distribution
In Figure 4 we plot the posterior distributions for χ eff , χ p and ρ p /ρ for all events used in this analysis. On av-erage the χ eff distribution for the events that prefer ELI over LVC (see Figure 3) is strongly peaked at zero with width comparable to ELI (see Figure 1). Meanwhile, the average χ eff distribution for the events that prefer LVC over ELI peaks at ∼ 0.2 with little support for χ eff ≤ 0.
On average there is no information from precession for binaries that prefer ELI over LVC or binaries that prefer LVC over ELI with the average χ p posterior resembling the prior and near zero ρ p /ρ. This hints at possible subpopulations in the preferred spin distribution: one with extremely low spins (EL) and one with larger spins (L). This is the same conclusion found by Roulet et al. (2021) and Galaudage et al. (2021). We therefore investigate these possible spin sub-populations by calculating odds ratios between models with extremely low spin (ELI) and models with larger spins both nearly aligned with the orbital angular momentum (LA) and spins isotropically distributed (LI).
We show in Figure 5 that most events in our population prefer a distribution with EL spins (∼ 80%) while several prefer larger spin magnitudes. Of those events that prefer larger spins, the nearly aligned distribution (LA) is preferred to the isotropic (LI) with an odds ratio of 9 : 1. This result is consistent with the conclusion from Galaudage et al. (2021) where it is inferred that a) 70 -90% of merging black hole binaries contain black holes with negligible spins and b) the high spin sub-population has spins preferentially aligned with the orbital angular momentum. However, the preferentially aligned spin conclusion is primarily driven by GW190517 055101 with the odds ratio reducing from 9 : 1 to 2 : 1 when GW190517 055101 is excluded from the population.
Interestingly, we see a positive correlation between the binaries primary mass and the preferred black hole spin distribution with most low mass binaries preferring distributions with EL spins and most high mass binaries preferring larger spin magnitudes. This suggests that the sub-populations found in Roulet et al. (2021) and Galaudage et al. (2021) could be correlated with the primary mass of the binary. In fact, we find that most binaries with primary mass less than 50M tend to prefer EL spin magnitudes (33/37) while most binaries with primary mass greater than 50M tend to prefer larger spins (5/8), see Figure 6. This result is in agreement with the conclusions found by Tiwari & Fairhurst (2021), which hinted at a possible correlation between the aligned spin magnitude and the chirp mass of the binary since all events in this high mass sub-population have chirp mass greater than 32 M -the point at which the aligned spin magnitude starts to increase (see Figure 2 of Tiwari & Fairhurst (2021)). This positive correlation may suggest evidence for hierarchical mergers in GWTC-2 (Kimball et al. 2020(Kimball et al. , 2021Tiwari & Fairhurst 2021;Gerosa et al. 2020b;Abbott et al. 2020e;Fishbach et al. 2017;Doctor et al. 2019), where the remnant of a previous "first generation" binary becomes part of a new one (e.g. Antonini & Rasio 2016); although see Jaraba & Garcia-Bellido (2021) for an alternative mechanism which allows for heavier black holes to have larger spins owing to close hyperbolic encounters spinning up black holes in dense clusters. This is because hierarchical mergers are expected to have a) larger black hole mass and b) larger spins since the remnant of a first generation binary is expected to have mass nearly equal to the sum of its components and spin a ≈ 0.7 (inherited from the orbital angular momentum of the previous binary) ( Buonanno et al. 2008). Similarly it is expected that merging black hole binaries with black hole mass m 50M can only be formed through hierarchical mergers since pairinstability supernova theory (see e.g. Woosley 2017;Belczynski et al. 2016;Gerosa & Berti 2017;Abbott et al. 2020c,e) prohibits black holes forming from direct stellar collapse with masses within the range ∼ 50 − 120M .
From Figure 6 we see that even when GW190814 is excluded from the analysis, low mass binaries prefer distributions with EL spins and high mass binaries prefer distribution with larger spin magnitudes. In fact we calculate that for low mass binaries, EL spin magnitudes are preferred to low spin magnitudes by ∼ 10 3 : 1 while for high mass binaries, low spin magnitudes are preferred to EL spin magnitudes by ∼ 20 : 1. Although for high mass binaries there is significantly larger support for aligned spins than for low mass binaries, the high mass sub-population prefers more isotropic spin orientations. This is primarily driven by GW190701 203306 and GW190929 012149 since they both have support for negative χ eff , χ eff = −0.07 +0.23 −0.29 and χ eff = 0.01 +0.34 −0.33 respectively, a region of parameter space which is not permitted by the aligned spin models used in this analysis (see Figure 1). Since the high mass sub-population prefers isotropic spin distributions this suggests that if originated from hierarchical mergers, they are likely to have formed in dense stellar clusters, where the spins are predicted to be isotropic, rather than in the accretion disks surrounding Active Galactic Nuclei where the spins are predicted to be aligned with the orbital angular momentum (Wang et al. 2021). Abbott et al. (2021b) also investigated whether there is evidence for a mass dependence in the BH spin distribution through a hierarchical Bayesian analysis of the population of known BBHs. Similar to the work presented here, Abbott et al. (2021b) also found a preference for higher spin magnitudes in higher mass events, although weaker than what we find in this work (see e.g. Figure 13 of Abbott et al. (2021b)). However, since the uncertainty on their measurement was broad, a mass dependence could not be confidently claimed.
We note that this potential mass dependence could arise from systematic errors since higher mass systems have far fewer cycles in the sensitive frequency band of the LIGO-Virgo detectors, making it significantly harder to accurately infer the black hole spin (see e.g. Abbott et al. 2020e). For example, Haster et al. (2016) recently found that binaries with a larger total mass tend to infer a larger positive black hole spin. We propose that in order to test this possible correlation, this analysis should be repeated on a simulated population in which all systems, high-and low-mass, have very small spins. If the correlation between mass and spin is not found, our conclusion is not an systematic artifact. We leave this for future work. Callister et al. (2021) recently identified an anticorrelation between the binaries mass ratio and black hole spin where more unequal mass binaries exhibit a systematically larger χ eff at 98.7% confidence. Our analysis draws similar conclusions, and finds evidence for a correlation between the binaries mass ratio and black hole spin. We find that more unequal mass binaries prefer larger spin magnitudes than more equal mass binaries. For instance, for binaries with component masses more unequal than GW190527 092055 (based on their median value), 8/15 prefer EL spins while 7/15 prefer larger spin magnitudes 1 This results in a preference for low spin magnitudes over EL spins by 2 : 1 (8 : 1 if GW190814 is excluded from the analysis). For binaries with more equal masses than GW190527 092055, EL spins are preferred for all binaries except for GW190517 055101. For this subset of binaries the preference for EL spins is significant -10 3 : 1. Of those unequal mass binaries with a preference for larger spin magnitudes, Figure 5 shows that several have primary mass m 1 > 50M . It is possible that the apparent anti-correlation between the binaries mass ratio and black hole spin is a consequence of the observed mass dependence described above, although unlikely on theoretical grounds since hierarchical mergers are predicted to generate χ eff distributions centered at 0 (see e.g. Figure 7 of Rodriguez et al. (2019)). We note that there is also the possibility that our conclusion is a model dependent effect because there exists an inherent correlation between the χ eff and mass ratio of the binary for our models. For example, because our models assume uncorrelated spins, it is harder to produce χ eff = 0 if the mass ratio is equal than if it is extreme. This is a different effect from the known mass ratio and aligned-spin degeneracy in the gravitational wave likelihood (Cutler & Flanagan 1994;Poisson & Will 1995;Baird et al. 2013) and we leave an investigation into this effect for future work.

CONCLUSION
In this work we performed an independent analysis and recomputed the black hole spin distribution using data from the second gravitational-wave catalog. We demonstrated that the surprising spin magnitude distribution obtained from Abbott et al. (2021b) is unlikely to be robust since the inclusion of GW190814 or exclusion of GW190517 055101 changes the preferred spin distribution to one with extremely low spins. We then demonstrated that our results are consistent with those from Roulet et al. (2021); Galaudage et al. (2021) and established that there is potential evidence for two spin sub-populations in the observed black holes -one with extremely low spins and one with larger spin magnitudes. We then made the argument that these spin sub-populations could be correlated with the primary mass of the binary, where we see an increase in spin magnitude for systems with higher masses, and argued 1 GW190814, GW190929 012149, GW190828 065509, GW190513 205428, GW190924 021846, GW190512 180714, GW151012, GW190930 133541 prefer EL spins while GW190412, GW151226, GW190720 000836, GW190706 222641, GW190519 153544, GW190620 030421, GW170729 prefer larger spin magnitudes. that this may provide evidence for hierarchical mergers in GWTC-2. Unlike recent works where hierarchical Bayesian inference has been used to infer the spin distribution of black holes, we chose to perform a detailed model selection analysis. We suggest that since there are still only a limited number of binary black hole observations, a much deeper understanding of the inferred black hole spin distribution can be achieved with this far simpler approach where it is clear how each GW candidate contributes to the final result. With the fourth gravitational-wave observing run anticipated to provide a plethora of additional binary black hole observations, with hopefully many more discoveries at high mass, we may soon be able to scrutinize the potential mass-spin correlation as well as deciphering the underlying black hole spin distribution.
Plots were prepared with the Matplotlib (Hunter 2007) and PESummary (Hoy & Raymond 2021) and both NumPy (Travis E 2006) and Scipy (McKinney 2010) were used in the analysis.