<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 15. The distributions of component spin magnitudes $\chi$ (left) and spin-orbit misalignment angles $\theta$ (right) among binary black hole mergers, inferred using the DEFAULT component spin model described further in Sect. B2a; e.g., both spin magnitudes are drawn from the same distribution. In each figure, solid black lines denote the median and central 90% credible bounds inferred on $p(\chi)$ and $p(\cos\theta)$ using GWTC-3. The light grey traces show individual draws from our posterior distribution on the DEFAULT model parameters, while the blue traces show our previously published results obtained using GWTC-2. As with GWTC-2, in GWTC-3 we conclude that the spin magnitude distribution peaks near $\chi_i \approx 0.2$, with a tail extending towards larger values. Meanwhile, we now more strongly favor isotropy, obtaining a broad $\cos\theta_i$ distribution that may peak at alignment ($\cos\theta_i = 1$) but that is otherwise largely uniform across all $\cos\theta$.</figcaption></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 16. Left panel: Inferred distribution of $\chi_{eff}$ for our latest full analysis in black. For comparison, the blue distribution and interval shows our inferences derived from GWTC2. Right panel: Corresponding result for $\chi_p$. While both panels in this figure are derived using the Gaussian spin model, we find similar conclusions with the other spin models used to analyze GWTC-2.</figcaption></figure>

at 99.7% credibility.

This interpretation was challenged in [194] and [195], which argued that no evidence for extreme spin misalignment exists if BBH spin models are expanded to allow the existence of a secondary subpopulation with vanishingly small spins. Other avenues of investigation are also in tension with the identification of extreme spin-orbit misalignment. When the $\chi_{eff}$ distribution is allowed to correlate with other BBH parameters, like the binary mass ratio (see Sec. VII B), evidence for negative $\chi_{eff}$ values diminishes [196]. Motivated by the concerns raised in [194] and [195], we repeat our inference of $\chi_{eff,min}$ but under an expanded model that allows for a narrow subpopulation of BBH events with extremely small effective

inspiral spins:

$$p(\chi_{eff}|\mu_{eff}, \sigma_{eff}, \chi_{eff,min}) = \zeta_{bulk}\mathcal{N}_{[\chi_{eff,min}, 1]}(\chi_{eff}|\mu_{eff}, \sigma_{eff}) + (1 - \zeta_{bulk})\mathcal{N}_{[-1, 1]}(\chi_{eff}|0, 0.01). \quad (18)$$

Here, $\zeta_{bulk}$ is the fraction of BBHs in the wide bulk population, truncated above $\chi_{eff,min}$, while $(1 - \zeta_{bulk})$ is the fraction of events residing in the vanishing spin subpopulation, which formally extends from -1 to 1. When repeating our inference of $\chi_{eff,min}$ under this expanded model, our data still prefer a negative $\chi_{eff,min}$ but with lower significance. As seen in Fig. 18, we now infer that $\chi_{eff,min} < 0$ at 92.5% credibility. This expanded model allows us to additionally investigate evidence for the existence
<figure><img src="image_20.png" /><figcaption>FIG. 17. Distribution of magnitudes of the most ($\chi_A$) and least ($\chi_B$) rapid component spin among BBHs in GWTC-3. Traces show individual draws from our posterior on the spin population under the DEFAULT model, while dark curves bound 90% credible bounds on $p(\chi_A)$ and $p(\chi_B)$.</figcaption></figure>

tence of a sub-population of BBHs with vanishingly small spins. GWTC-3 prefers but does not require such a sub-population to exist. We measure $\zeta_{\text{bulk}} = 0.54^{+0.36}_{-0.26}$, with $\zeta_{\text{bulk}} > 0.2$ at 99% credibility, but also find that our posterior remains consistent with $\zeta_{\text{bulk}} = 1$.

A. Spin distribution consistent as mass increases

Our previous analysis adopted the same spin distribution at all masses. The spins of low-mass binaries dominate the reconstructed spin distribution. However, the binaries with the most extreme values of spins have heavier masses: observations GW170729, GW190517, GW190519, GW190620, GW190706, GW190805, and GW191109, constitute 70% of the binaries with moderate to high spins. This preponderance of massive binaries with large spin suggests a one-size-fits-all approach might not fully capture how well we can predict black hole spins *given their mass*; conversely, this preponderance can also reflect the increased impact on our search sensitivity at the highest masses. Too, astrophysical formation scenarios often predict correlations between mass and spin, both from isolated and dynamical formation [173, 197]. Using the FM model for the aligned spin components we reconstruct the trend of $|s_z|$ versus mass. Figure 19 shows the aligned spin magnitude distribution versus binary chirp mass. At low masses, the aligned spin is consistent with (and well constrained to be close to) zero, (i.e., maximum aligned spin magnitude averaged over chirp masses $30M_{\odot}$ or less is 0.38 at 90% credibility). At heavier masses, the aligned spin is still consistent with zero, albeit with larger dispersion (i.e., maximum aligned spin magnitude averaged over chirp masses $30M_{\odot}$ or more is 0.5 at 90% credibility). This trend is qualitatively consistent with the relative proportion of

<figure><img src="image_21.png" /><figcaption>FIG. 18. Cumulative probabilities on the minimum truncation bound on the $\chi_{\text{eff}}$ distribution as inferred using GWTC-2 and GWTC-3. When modeling the effective inspiral spin distribution as a Gaussian truncated on $\chi_{\text{eff,min}} \le \chi_{\text{eff}} \le 1$, we inferred using GWTC-2 that $\chi_{\text{min,eff}} < 0$ at 99.1% credibility, and hence that the data support the existence of BBH mergers with negative effective inspiral spins. Using GWTC-3, this same analysis more strongly infers that $\chi_{\text{min,eff}} < 0$, now at 99.7% credibility. As discussed further below, evidence for negative effective inspiral spins is diminished under an expanded model that allows for a subset of BBHs to possess vanishing effective inspiral spins. When instead modeling the $\chi_{\text{eff}}$ distribution as a mixture between a broad Gaussian and a narrow Gaussian sub-population centered at $\chi_{\text{eff}} = 0$ (e.g., the second consistent with zero spin), we infer $\chi_{\text{min,eff}} < 0$ at 92.5% credibility.</figcaption></figure>

events versus chirp mass: very few observations have high chirp masses, providing relatively little leverage to constrain spins. At high chirp masses, the spin distribution is poorly constrained by only a handful of measurements, closer to our broad prior assumptions, in contrast to the better-constrained distribution at low mass. We have no evidence to support or refute a trend of aligned spin with chirp mass.

Figure 19 suggests aligned spin magnitude remains constrained to be close to zero independently of the most well-identified peaks in the mass distribution, contrary to what would be expected from hierarchical formation scenarios for these peaks [110, 197–201].

B. High spin correlates with asymmetric binaries

BBHs may exhibit an anti-correlation between their mass ratios and spins, such that binaries with $q \sim 1$ favor effective inspiral spin parameters near zero, while binaries with more unequal mass ratios exhibit preferentially positive $\chi_{\text{eff}}$ values [196]. To evaluate the degree to which $q$ and $\chi_{\text{eff}}$ are (or are not) correlated, following prior work [196] we adopt a Gaussian model for the $\chi_{\text{eff}}$ distribution with a mean and standard deviation that are
<figure><img src="image_22.png" /><figcaption>FIG. 19. The dependence of aligned spin magnitude on the chirp mass. The light/dark shaded regions are the aligned spin magnitude at a credibility 90%/50%. The distribution is consistent with small values for lower chirp mass binaries, however, the spin magnitude is less tightly constrained for chirp masses of 30 $M_{\odot}$ and higher.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIG. 21. Posteriors on the mass ratios and effective inspiral spins of BBHs in GWTC-3, reweighted to a population-informed prior allowing for a correlation between $q$ and $\chi_{\text{eff}}$. We infer that the mean of the BBH $\chi_{\text{eff}}$ distribution shifts towards larger values with decreasing mass ratios. Accordingly, reweighted events shift considerably, such that events with $q \sim 1$ contract about $\chi_{\text{eff}} \approx 0$ while events with $q < 1$ shift towards larger effective inspiral spins.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIG. 20. Posterior constraints on the mean (top) and standard deviation (bottom) of the $\chi_{\text{eff}}$ distribution as a function of mass ratio $q$. At 97.5% credibility, we find that the mean of the $\chi_{\text{eff}}$ shifts towards larger values for more unequal mass systems. The grey region in the lower panel shows the area artificially excluded by our prior on the parameters $\sigma_0$ and $\beta$; see Eq. (20).</figcaption></figure>

allowed to evolve with $q$:

$$p(\chi_{\text{eff}}|q) \propto \exp \left[ - \frac{(\chi_{\text{eff}} - \mu(q))^2}{2\sigma^2(q)} \right], \quad (19)$$

with

$$\mu(q) = \mu_0 + \alpha(q - 1) \quad (20a)$$

$$\log_{10} \sigma(q) = \log_{10} \sigma_0 + \beta(q - 1). \quad (20b)$$

The new hyperparameters $\alpha$ and $\beta$ measure the extent to which the location or width of the $\chi_{\text{eff}}$ distribution changes as a function of mass-ratio.

We repeat hierarchical inference of the BBH population, adopting the fiducial model for the primary mass and redshift distribution. At 97.5% credibility, we constrain $\alpha < 0$, indicating that more unequal-mass binaries preferentially possess larger, more positive $\chi_{\text{eff}}$. Figure 20 illustrates our constraints on the mean and standard deviation of the $\chi_{\text{eff}}$ distribution as a function of mass ratio. Each light trace represents a single sample from our hyperposterior, and the solid black lines denote the median values and central 90% bounds on $\mu(q)$ and $\sigma(q)$ at a given value of $q$. If we adopt these hierarchical results as a new, population-informed prior, Fig. 21 shows the resulting reweighted posteriors for the BBHs among GWTC-3. Each filled contour bounds the central 90% region for a given event in the $q-\chi_{\text{eff}}$ plane, while black points mark events' one-dimensional median $q$ and $\chi_{\text{eff}}$ measurements.
VIII. COMPARISON WITH OTHER GW
CATALOGS

In this paper, we have presented population inferences
based upon events identified by the LIGO Scientific,
Virgo and KAGRA Collaborations in data taken by the
Advanced LIGO and Advanced Virgo instruments during
their first three observing runs [1, 2]. We have imposed a
FAR threshold of < 0.25 yr⁻¹ across all analyses incorpo-
rating NS binaries and a lower threshold FAR < 1 yr⁻¹ for
BBH analyses. This excludes several events which pass
the threshold of p<sub>astro</sub> > 0.5 for inclusion in GWTC-3
. In addition, a number of analyses of the public GW
data from O1, O2 and O3a [24, 66–69, 71] have identi-
fied additional candidate binary merger events. In the
remainder of the paper, we have restricted the primary
analysis to events included in GWTC-3. The overriding
reason for this is that differences in the analysis meth-
ods prevent a detailed evaluation of search sensitivity, as
described in Section III, which is critical to interpreting
the population. In this Section, we investigate the con-
sistency of the remaining GWTC-3 events and additional
non-GWTC events with the population models inferred
in this paper.

For concreteness, when referring to results reported
by external groups, we include all events identified
as GWs in their catalogs. In O1, there is one ad-
ditional event, GW151216 identified in [68]. The
additional events from O2 are GW170304, GW170425,
and GW170403 which are identified in [69, 71], and
GW170121, GW170202, and GW170727 which were also
then independently found in open gravitational-wave
catalog two (2-OGC) [67]. In O3, we include 16 ad-
ditional events. These include GW190916.200658 and
GW190926.050336 which were originally identified in
open gravitational-wave catalog three (3-OGC) [24] and
independently identified in Gravitational Wave Tran-
sient Catalog 2.1 (GWTC-2.1) [3]; GW190403.051519,
GW190426.190642 and GW190514.065416 which are in-
cluded in GWTC-2.1 but have a FAR below our < 1 yr⁻¹
threshold; GW191113.071753, GW191126.115259,
GW191204.110529, GW191219.163120,
GW200208.222617, GW200210.092254,
GW200220.061928, GW200220.124850,
GW200306.093714, GW200308.173609,
GW200322.091133 which are included in GWTC-3
but again have a FAR below our < 1 yr⁻¹ threshold.

In Figure 22, we show the additional gravitational wave
events which were not included in the sample used in
this paper. The additional events are broadly consis-
tent with the population presented here although several
events lie at the boundaries of the identified population.
Specifically, two of the events have aligned spins that lie
outside the inferred population. These are GW151216
with a mean $\chi_{\text{eff}} = 0.82$ and GW170403 with a mean
$\chi_{\text{eff}} = -0.58$. The analysis in [71] used a prior which is
constant in $\chi_{\text{eff}}$ which is significantly different from the
uniform in spin-magnitude prior used in the GWTC pa-

pers. A re-analysis of GW151216 and GW170403 [202]
leads to inferred χ<sub>eff</sub> distributions which are more consis-
tent with the population inferred here. Specifically, this
gives χ<sub>eff</sub> = 0.5<sup>+0.2</sup><sub>-0.5</sub> for GW151216 and χ<sub>eff</sub> = -0.2<sup>+0.4</sup><sub>-0.3</sub>
for GW170403. In addition, the sub-threshold events
from GWTC-3 extend the distribution to both higher
masses and higher mass ratios. However, only low sig-
nificance events currently populate these regions. Addi-
tional observations in future runs will allow us to deter-
mine whether these low significance events are more likely
spurious, or were the first hints of a broader population
in the mass space.

With regard to events potentially containing NS,
GWTC-3 contains several candidates that do not sat-
isfy our FAR< 1/yr threshold but do have m<sub>2</sub> poten-
tially consistent with NS masses, namely GW191219 and
GW200210. Both events are inferred to have highly
asymmetric masses and could possibly be an indication
of additional NSBH sources, or asymmetric BBH similar
to GW190814. Further observations in future runs will
again allow us to investigate these interesting regions of
the binary parameter space in greater detail.

IX. ASTROPHYSICAL INTERPRETATION

A. Implications for binary black hole formation

1. *Mass distribution*

The statistical distribution of BH source properties
such as their mass, spin and redshift can be used to
probe the astrophysics of BH binary formation and evo-
lution [94, 110, 113, 173, 180, 184, 199, 201, 203–214].
The analysis performed in Sec. VI has identified struc-
tures in the mass distribution of BBHs that go beyond
a standard power-law model and can help to shed light
on formation processes. These features were previously
identified in [47], but we are now more confident that
they are statistically significant (see Sec. VI).

The underlying mass distribution of BBHs inferred in
this paper peaks at a primary mass ~10M⊙, withthe ma-
jority of BBHs having a primary BHs with a mass lower
than this value (e.g., see Fig. 11). Formation in glob-
ular clusters has been long recognized as an important
channel for merging BBHs [215–223]. In this scenario,
BBHs are assembled during three body dynamical inter-
actions in a low metallicity environment. The resulting
BH mass distribution is generally predicted to peak at
> 10M⊙. Three recent studies of globular cluster for-
mation find that the BBH merger rate is severely sup-
pressed where we observe a peak: one study [159] finds
that the BBH merger rate is severely suppressed below
about m ≃ 13M⊙ with a corresponding realistic merger
rate at this mass value of ~0.5 Gpc⁻³yr⁻¹M⊙⁻¹ (see their
Fig. 2); another recent study [224] finds similar results,
with a peak in their mass distribution at about m ≃ 15–
20M⊙ (see their Fig. 5); a third analysis [217] finds the
<figure><img src="image_25.png" /><figcaption>FIG. 22. The measured properties of the BBH candidates not included in the population study presented in this paper (shaded regions), compared to the inferred population from the PP model presented in Section VI A (black contours). These include both events which fall below our FAR threshold as well as events identified by other groups. The events are color coded based upon the search which first identified them: catalogs from O1 and O2 [68, 69, 71] in red, events from open gravitational-wave catalog three (3-OGC) (which incorporates events in O1–O3a) [24] in green, and from GWTC-3 with FAR> 1 yr⁻¹ threshold in blue.</figcaption></figure>

peak at $m \simeq 20M_{\odot}$. Taking these results at face value, the inferred high merger rate of sources with $\lesssim 10M_{\odot}$ may suggest that globular clusters contribute subdominantly to the detected population. Dynamical formation in young clusters is also disfavored to explain the whole BH population at $m \simeq 10M_{\odot}$ because lighter BHs are ejected by supernova kicks and do not participate to the dynamical evolution of the cluster [225–227].

Galactic nuclei can produce a BBH population with a much wider mass spectrum than both young and globular clusters [228–234]. Because of their high metallicities and escape velocities, nuclear star clusters can form and retain a significant number of lighter BHs, which can then pair and merge. BBH formation near an AGN disk can produce a significant population of BBH mergers with a wide mass spectrum [167, 197, 201, 235, 236]. In such scenarios, the observed low mass overdensities without counterparts in spin could be reflections of supernova physics; by contrast, in these hierarchical formation models no evident mechanism can impart them without a corresponding signature in spin. If the BBHs are formed near an AGN disk, this process might select heavier BHs, hardening the BBH mass function and driving the peak of the mass distribution towards values higher than observed [237].

Isolated binary evolution models often predict a peak near $m \simeq 10M_{\odot}$ [173, 238–241]. Recent population models find component masses of merging BBHs that peak at $8–10M_{\odot}$ and come from $\sim 20–30M_{\odot}$ progenitors [173, 175, 240]. The overall merger rate normalization of the peak remains, however, poorly constrained. Moreover, the peak of the mass distribution can shift significantly depending on the adopted supernova, natal kick, mass transfer, and wind prescriptions, and star formation history of the Universe [32, 166, 175, 227, 239, 242, 243].

The analysis in Section VI suggests two additional peaks in the mass distribution at $m \sim 17M_{\odot}$ and at

$35M_{\odot}$. The three most significant mass peaks are therefore separated by roughly a factor of two from each other [47]. Assuming these peaks exist, an explanation consistent with our constraints on BH spins is that they originate either from the initial BH mass function, or that they are produced by different populations formed by separate physical processes or formation channels.

The other feature of the inferred BH mass distribution that was shown in our analysis is the apparent lack of truncation at $m \sim 40M_{\odot}$, which confirms our results based on GWTC-1 and GWTC-2 [2, 25]. A mass gap between approximately $50_{-10}^{+20}M_{\odot}$ and $\sim 120M_{\odot}$ is predicted by stellar evolution models as the result of the pair instability process in the cores of massive stars [147–153]. However, due to our limited knowledge of the evolution of massive stars, the formation of BHs heavier than $\sim 40M_{\odot}$ from stellar collapse cannot be fully excluded [153, 173, 244–247]. The location of the mass gap is sensitive to the uncertain $^{12}$C($\alpha$, $\gamma$)$^{16}$O reaction rate, which governs the production of oxygen at the expense of carbon [152, 246, 248]. Moreover, BHs formed from progenitor stars at low metallicities ($Z/Z_{\odot} \lesssim 0.1$) might avoid all together the mass limit imposed by pair-instability [178, 244, 247]. The lack of a sharp truncation at high masses might indicate a dynamical process, such as the hierarchical merger of BHs [110, 154, 198–200, 224, 229, 234, 249–253] or stars [254–257] in dense clusters or in the gaseous disk surrounding a massive BH [213]. In a hierarchical scenario we would expect the more massive BHs to also have the larger spins [199, 200]. While we do observe such a mass-spin correlation above $m \sim 40M_{\odot}$ (Fig. 19), the binaries with a signature that $\chi_{\text{eff}}$ is not zero all prefer $\chi_{\text{eff}} > 0$ (see Table IV), while hierarchical formation in dynamical environments would lead to isotropically oriented spins. BHs can also increase their birth mass beyond the pair-instability mass gap through the efficient accretion of gas from a stel-
lar companion or from a surrounding gaseous disk [258–262]. Highly coherent accretion on one of the BHs could also explain the negative correlation between $\chi_{\text{eff}}$ and $q$ shown in Fig. 20, although accretion in gas-free scenarios should be highly super-Eddington in order to impart significant spin [258, 262]. Alternatively, primordial BHs can have masses above the pair-instability mass threshold, although this most likely requires efficient accretion before the reionization epoch in order not to violate current constraints [263].

## 2. Redshift distribution

In Section VI we showed that the BBH merger rate increases with redshift, as $(1+z)^\kappa$, with $\kappa \sim 3$. Although error bars are large, current data prefer a model in which the merger rate evolves steeply with redshift and at a rate that is consistent with the growth in star formation. For binary formation in the field, the predicted value of $\kappa$ is sensitive to the assumed efficiency of common envelope ejection: values between $\kappa = 0.2$ and 2.5 are all possible, although relatively small values $\kappa \sim 1$ are preferred [160, 166, 175, 227, 264]. Delay times in the field are also dependent upon stability of mass transfer, e.g., [265]. Similarly, $\kappa \lesssim 2$ is often found in models of BBH formation in open and young clusters [225, 266]. Dynamical formation in globular clusters predicts $\kappa \lesssim 2$ [164, 267], e.g., [159] find $\kappa = 1.6^{+0.4}_{-0.6}$, and show that the most important parameter affecting the value of $\kappa$ in the globular cluster scenario is the initial cluster half-mass density, $\rho_h$, while uncertainties in other model parameters (e.g., natal kicks, black hole masses, metallicity) have a small effect. Only models in which globular clusters are formed with a high half-mass density, $\rho_h > 10^5 M_\odot \text{pc}^{-3}$, lead to $\kappa \gtrsim 2$ [159]. While uncertainties are large, improved constraints on the merger rate evolution have the potential to unveil important information about the physics of massive binaries [96, 173, 239, 264] and the initial conditions of clusters across cosmic time [164, 208, 267].

## 3. Spin distribution

We observe evidence that the spin distribution both requires spin-orbit misalignment and also includes events with anti-aligned spins. BBHs with a large spin-orbit misalignment can be formed in dynamical environments such as globular, nuclear, and young star clusters, or active galactic nuclei [185, 198, 229, 230]. In these systems, two single BHs are paired together during a three body interaction and/or undergo a number of subsequent dynamical interactions before merging. Their spins have directions that are therefore uncorrelated with each other and with the orbital angular momentum of the binary, leading to an isotropic spin-orbit alignment [183, 185]. The evolution of BH spins in AGN disks depends on several uncertain factors, such as the importance of accre-

tion and dynamical encounters, the initial spin orientation, and the efficiency of migration [201, 268]. If radial migration of BHs is inefficient, the distribution of $\chi_{\text{eff}}$ skews toward higher values because scattering encounters that randomize spin directions become less frequent. On the other hand, efficient migration would imply more frequent dynamical encounters, producing a $\chi_{\text{eff}}$ distribution centered around zero [268]. However, the dispersion of the $\chi_{\text{eff}}$ distribution also increases characteristically with mass, as with other hierarchical formation scenarios [110, 197, 201].

Formation from field binaries is thought to produce components with preferentially aligned spins [181, 183, 185]. Such an alignment, however, is not certain. In fact, all population models of isolated binaries customarily start with the stellar progenitor spins initially perfectly aligned with the orbital angular momentum of the binary. This assumption is made due to simplicity and partly because tidal interactions are thought to quickly remove any spin-orbit misalignment prior to BH formation. However, the observational evidence of close massive binaries with highly inclined spin axes suggests that close massive binaries can form with misaligned spins and that tides might not in all cases be able to realign the spins [269–272]. Moreover, a large spin-orbit misalignment can be produced if a binary is the inner component of a triple system [267, 273–276], where the tertiary component can be either a star, another BH, or even a massive BH [231, 277–279]. In this scenario, the secular gravitational interaction of the binary with an external companion can randomize the orbital plane of the binary. The complex precessional dynamics of the BH spins in triple systems also changes the spin orientation and lead to a distribution of $\chi_{\text{eff}}$ peaked near zero, although with a marginal preference for aligned spins [273, 276]. A spin-orbit misalignment can also be produced in field binaries by a stable episode of mass-transfer prior to the formation of the BHs [280] or by asymmetric mass and neutrino emission during core-collapse [193, 281], although these processes are unlikely to produce a large misalignment for a significant fraction of the population [186, 212, 282]. We conclude that the presence of systems with misaligned spins is not in contradiction with a scenario in which the majority, if not all, BBHs form in the field of galaxies. On the other hand, the fact that the $\chi_{\text{eff}}$ distribution is not symmetric around zero, if confirmed, can be used to rule out a model in which all BBHs are formed through dynamical encounters in star clusters [194, 204].

Corroborating our previous conclusion based on GWTC-2, we find that the BH population is typically described by small spins. Predictions for BH spin magnitudes vary depending on the assumptions about stellar winds and their metallicity dependence, tides, and are particularly sensitive to the efficiency of angular momentum transport within the progenitor star [188, 189]. If the stellar core remains strongly coupled to the outer envelope during the stellar expansion off the main sequence, then a significant amount of spin can be carried from the
core to the envelope. In this case, a BH formed from stellar collapse may be born with nearly zero spins. For formation in isolated binaries, this implies that the first-born BH will essentially be a Schwarzschild black hole. The second-born BH can still form with significant spin as tidal interactions may realign and increase the spin of its stellar progenitor in between the two supernova explosions [283, 284]. If the binary undergoes chemically homogeneous evolution [285, 286], its components may both be tidally spun up to near break-up velocity, and keep this rotation rate throughout main sequence evolution, evolving into BHs with large and aligned spins. Black holes that form during the QCD phase transition in the early Universe will all have essentially zero natal spins [287, 288]. However, a significant spin can be attained through subsequent gas accretion [263]. Finally, if BBHs are formed or migrate within the accretion disk of a supermassive BH, they can accrete from the surrounding gaseous environment and spin up [268].

We observe neither evidence for nor against an increase in spin magnitude for systems with higher masses [213, 289, 290] and more unequal mass ratios [196]. Current stellar evolution models suggest that larger BH masses should correlate with smaller spins because larger BHs originate from more massive stars which undergo more extensive mass loss, carrying away most of the angular momentum and producing BHs with small spins [173, 181, 284]. A consequence of this should be either a decrease in spin magnitude (and $\chi_{eff}$) with mass above $\sim 20M_{\odot}$ or no correlation, where the predicted trend depends on the specific stellar evolution models adopted [173, 181]. Predictions remain uncertain and are strongly dependent on modeling assumptions about the angular momentum transport within the star, spin dissipation during the supernova and the treatment of binary interaction prior to BH formation. If future observations identify such a trend, then an increase in spin magnitude with mass might suggest a hierarchical formation scenario. However, as mentioned above, this scenario seems currently at odds with the fact that binaries with more unequal mass ratios and massive components exhibit preferentially positive $\chi_{eff}$.

## B. Implications for neutron stars

One result from gravitational wave (GW) observations is tension with the strong preference for $1.35M_{\odot}$ mass objects which has been recovered in galactic BNS [41]. Instead, conservatively assuming all objects below the maximum neutron star mass are neutron stars, our unmodeled analysis of the lowest-mass compact objects is consistent with a broad unimodal Gaussian, allowing for highly asymmetric binaries. Our analysis of all individual low-mass (assumed NS) objects suggests a wide NS mass distribution, without the bimodal structure seen in the galactic NS population. The GW-observed population of low-mass mergers is still small. If this tension

persists, however, several avenues exist to explain a discrepancy between Galactic and GW observations, including but not limited to additional formation channels for GW systems; strong observational selection effects, like those used to explain GW190425 [21, 291, 292] and the smaller body in GW200105; or the prospect that BHs form below the maximum NS mass.

Our conclusions about the compact object mass spectrum in general and the mass spectrum of NS in particular will have substantial impact on the understanding of the stellar explosions that generate such compact objects [32, 293, 294] and the binary interactions that carry these objects towards merger, assuming a stellar origin for low-mass binary mergers.

Our analyses show no evidence for or against the presence of a mass distribution feature closely corresponding to the maximum neutron star mass. Rather, the shape of the neutron star mass distribution, the existence of GW190814, and our results for the mass distribution for compact objects between $3M_{\odot}$ and $7M_{\odot}$ may instead suggest a continuous mass spectrum, albeit strongly suppressed above the masses of known NS.

Fortunately, the comparatively high prevalence of objects close to the maximum neutron star mass suggests that we will likely observe several objects near this region in the future, potentially providing several avenues to connect features in the NS mass distribution to fundamental nuclear physics. Our analysis of NS in merging binaries alone alone suggests the NS mass distribution extends to the maximum NS mass $M_{TOV}$ expected from the EOS.

Our analyses are also consistent with both symmetric ($q \simeq 0.8$) and significantly asymmetric ($q < 0.8$) binaries containing NS in BNS, and modestly ($q \in [0.5, 0.8]$) to highly ($q < 0.5$) asymmetric binaries in NSBH. Compared to equal-mass mergers [295, 296], modestly asymmetric NS mergers (with either NS or BH counterparts) are potentially strong candidates for multimesenger counterparts [297], since an asymmetric merger can eject more mass [298], produce a larger remnant disk [299, 300], and potentially produce significant associated gamma ray burst emission [301–303]. For BNS, our analyses are consistent with a significant fraction of highly asymmetric events. For NSBH, the discovery of GW200105 and GW200115 demonstrate the existence of asymmetric binaries containing NS with a range of mass ratios. Based on these events, our inferences about the low-mass compact object distribution suggests that EM-bright NSBH mergers could occur at a significant fraction of the overall NSBH rate. Generally, a broad mass ratio distribution suggests modestly more favorable prospects for electromagnetic follow-up observations. Conversely, a broad mass ratio distribution complicates simple efforts to interpret existing GW observations which were developed under the assumption that low-mass binary mergers are very frequently of comparable mass [304].

Finally, our analyses here leave GW190814 as an outlier both from BBH systems and from systems that con-
tain a likely NS. Neither component of this binary has exceptional masses; for example, the secondary component could easily be produced from conventional supernova engines [32]. However, based on the merger rates versus mass identified in our study, this system (and the larger sample of high-mass-ratio binaries available at a lower threshold) may require a different formation pathway [230, 305, 306].

## X. THE GW BACKGROUND FROM BINARY MERGERS

The observation of binaries with masses in the NSBH range allows us to provisionally complete a census of the different classes of compact binaries that contribute to an astrophysical gravitational-wave background, assuming our existing surveys are sensitive to all relevant sources (i.e., not accounting for frequent subsolar mass mergers). We have previously predicted the contributions of BBH and BNS mergers to the gravitational-wave background, based on the compact binary population observed in GWTC-2 [172]. In Fig. 23, we update this forecast with our latest knowledge of the BBH and BNS population and the newly measured rate of NSBH mergers.

The shaded bands on the left side of Fig. 23 shows estimates of and uncertainties on the dimensionless energy-density spectra

$$ \Omega(f) = \frac{1}{\rho_c} \frac{d\rho}{d\ln f} \qquad (21) $$

of gravitational waves radiated by each class of compact binary. In Eq. (21), $d\rho$ is the gravitational-wave energy density per logarithmic frequency interval $d\ln f$, while $\rho_c$ is the critical energy density required to close the Universe. We adopt the same model for the merger history of compact binaries used previously [172, 307], assuming that compact binary formation rate traces a metallicity-weighted star formation rate model [308–310] with a $p(t_d) \propto t_d^{-1}$ distribution of time delays $t_d$ between binary formation and merger. Time delays are restricted to 20 Myr $\le t_d \le 13.5$ Gyr for BNS and NSBH mergers and 50 Myr $\le t_d \le 13.5$ Gyr for BBHs [175, 311], with binary formation restricted to redshifts below $z_{max} = 10$. The birth rate of BBH progenitors is further weighted by the fraction of star formation at metallicities $Z < 0.1 Z_{\odot}$ [174, 312].

Within Fig. 23, the stochastic energy-density due to BBHs has been marginalized over our uncertainty on both the local merger rate and mass distribution of the BBH population, as measured using the PP mass model. At 25 Hz, we estimate the energy-density due to BBHs to be $\Omega_{BBH}(25\,Hz) = 5.0_{-1.8}^{+1.4} \times 10^{-10}$. To estimate the contribution due to BNS systems, we adopt the simple rate measurement presented in Section IV A under a fixed mass distribution, and correspondingly assume a uniform distribution of neutron star masses between 1

and $2.5 M_{\odot}$, giving $\Omega_{BNS}(25\,Hz) = 0.6_{-0.5}^{+1.7} \times 10^{-10}$. The contribution due to NSBH systems, meanwhile, is estimated using the BGP rate reported in Table II. For simplicity, we again assume a uniform distribution of neutron star masses between 1 and $2.5 M_{\odot}$ and a logarithmically uniform distribution of black hole masses between 5 and 50 $M_{\odot}$ among NSBH mergers. Under these assumptions, we find $\Omega_{NSBH}(25\,Hz) = 0.9_{-0.7}^{+2.2} \times 10^{-10}$.

The blue band in the right side of Fig. 23 denotes our estimate of the total gravitational-wave background due to the superposition of these three source classes; we expect a total energy-density of $\Omega(25\,Hz) = 6.9_{-2.1}^{+3.0} \times 10^{-10}$. For comparison, the solid black curve marks our present sensitivity to the gravitational-wave background [172, 313]. Although our estimate for the background amplitude lies well below current limits, it may be accessible with future detectors, such as the planned “A+” LIGO configuration.

## XI. CONCLUSIONS

The third LIGO–Virgo gravitational wave transient catalog (GWTC-3) [1] has increased our census of the population of compact mergers by nearly a factor of two, compared to our analysis of the first half of O3 [11]. We simultaneously employ all observations with FAR< 0.25 yr$^{-1}$ to infer the merger rate versus both component masses across the observed mass spectrum. For NS, we find a broad mass distribution, extending up to $2.0_{-0.3}^{+0.3} M_{\odot}$, in contrast to the unimodal mass distribution observed for Galactic BNS. We find the BBH mass distribution is nonuniform, with overdensities at BH masses of $10 M_{\odot}$ and $35 M_{\odot}$. These overdensities may reflect the astrophysics associated with generating coalescing binaries, potentially reflecting properties of stellar physics or astrophysical environments. These features may assist future applications of gravitational wave astronomy. As an example, these sharp features could be redshift-independent and, if so, used as standard candles for cosmology [314, 315]. We find the BH mass distribution exhibits an interval between $2.2 M_{\odot}$ and $6.1 M_{\odot}$ where merger rates are suppressed, which could be consistent with past X-ray observations suggesting a mass gap [28–31]. Our analysis lacks sufficient sensitivity to probe the structure of the mass distribution at the highest masses $m_1 > 70 M_{\odot}$ in detail; however, so far, we find no evidence for or against an upper mass gap.

We find that observed BH spins are typically small (half less than 0.25). We still conclude that at least some of these spins exhibit substantial spin-orbit misalignment. We corroborate a correlation between BBH effective aligned inspiral spins and mass ratio.

Using parametric models to infer the distribution of BBH merger rate with redshift, we find the BBH merger rate likely increases with redshift; we cannot yet assess more complex models where the shape or extent of the mass distribution changes with redshift.
<figure><img src="image_26.png" /><figcaption>FIG. 23. Forecast of astrophysical gravitational-wave background due to binary mergers following O3. (Left): The individual contributions expected from BNS, NSBH and BBH mergers. While uncertainties on the energy-density due to BNS and NSBH are due to Poisson uncertainties in their merger rates, our forecast for the stochastic background due to BBHs additionally includes systematic uncertainties associated with their imperfectly known mass distribution. (Right): Estimate of the total gravitational-wave background (blue), as well as our experimental current sensitivity (solid black) [172, 313]. For comparison, we additionally show the expected sensitivities of the LIGO-Virgo network at design sensitivity, as well as that of LIGO's anticipated “A+” configuraton.</figcaption></figure>

Analyses presented in our previous work [11] and in a companion paper [316] employ coarse-grained models for the BBH population, smoothing over some of the subtle features identified above. We find that these coarse-grained models draw similar conclusions on current data as our previous studies; see Sec. VI A. Applications that focus on large-scale features of the mass distribution (e.g., the stochastic background, as described in Sec. X) only require these coarse-grained results. Nonetheless, the mass distribution remains a critical source of systematic uncertainty in any merger rate integrated over any mass interval, particularly in mass intervals with few observations. We specifically find the BNS and NSBH merger rates exhibit considerable uncertainty in the mass distribution, with relative merger rate errors within (and between) models far in excess of the expected statistical The Poisson error associated with the count of these events. These systematics propagate directly into our most conservative estimates for their merger rates.

The next GW survey could have a BNS detection range increased by approximately 15–40% [317]. Even without allowing for increased merger rates at higher redshift, the next survey should identify roughly 3 times more events of each class then used in this study, including several new events from the BNS and BHNS category. We continuously revise our assessment of future observing prospects [317].

ACKNOWLEDGMENTS

This material is based upon work supported by NSF's
LIGO Laboratory which is a major facility fully funded
by the National Science Foundation. The authors also

gratefully acknowledge the support of the Science and
Technology Facilities Council (STFC) of the United
Kingdom, the Max-Planck-Society (MPS), and the State
of Niedersachsen/Germany for support of the construc-
tion of Advanced LIGO and construction and operation
of the GEO600 detector. Additional support for Ad-
vanced LIGO was provided by the Australian Research
Council. The authors gratefully acknowledge the Italian
Istituto Nazionale di Fisica Nucleare (INFN), the French
Centre National de la Recherche Scientifique (CNRS)
and the Netherlands Organization for Scientific Research
(NWO), for the construction and operation of the Virgo
detector and the creation and support of the EGO consor-
tium. The authors also gratefully acknowledge research
support from these agencies as well as by the Council
of Scientific and Industrial Research of India, the De-
partment of Science and Technology, India, the Science
& Engineering Research Board (SERB), India, the Min-
istry of Human Resource Development, India, the Span-
ish Agencia Estatal de Investigación (AEI), the Span-
ish Ministerio de Ciencia e Innovación and Ministerio de
Universidades, the Conselleria de Fons Europeus, Uni-
versitat i Cultura and the Direcció General de Política
Universitaria i Recerca del Govern de les Illes Balears,
the Conselleria d'Innovació, Universitats, Ciència i Soci-
etat Digital de la Generalitat Valenciana and the CERCA
Programme Generalitat de Catalunya, Spain, the Na-
tional Science Centre of Poland and the European Union
– European Regional Development Fund; Foundation for
Polish Science (FNP), the Swiss National Science Foun-
dation (SNSF), the Russian Foundation for Basic Re-
search, the Russian Science Foundation, the European
Commission, the European Social Funds (ESF), the Eu-
ropean Regional Development Funds (ERDF), the Royal
Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the French Lyon Institute of Origins (LIO), the Belgian Fonds de la Recherche Scientifique (FRS-FNRS), Actions de Recherche Concertées (ARC) and Fonds Wetenschappelijk Onderzoek - Vlaanderen (FWO), Belgium, the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, the Natural Science and Engineering Research Council Canada, Canadian Foundation for Innovation (CFI), the Brazilian Ministry of Science, Technology, and Innovations, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Leverhulme Trust, the Research Corporation, the Ministry of Science and Technology (MOST), Taiwan, the United States Department of Energy, and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, INFN and CNRS for provision of computational resources. This work was supported by MEXT, JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, JSPS Grant-in-Aid for Scientific Research on Innovative Areas 2905: JP17H06358, JP17H06361 and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grant-in-Aid for Scientific Research (S) 17H06133 and 20H05639, JSPS Grant-in-Aid for Transformative Research Areas (A) 20A203: JP20H05854, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, National Research Foundation (NRF) and Computing Infrastructure Project of KISTI-GSDC in Korea, Academia Sinica (AS), AS Grid Center (ASGC) and the Ministry of Science and Technology (MoST) in Taiwan under grants including AS-CDA-105-M06, Advanced Technology Center (ATC) of NAOJ, Mechanical Engineering Center of KEK. We would like to thank all of the essential workers who put their health at risk during the COVID-19 pandemic, without whom we would not have been able to complete this work.

## Appendix A: Sensitivity estimation

A key ingredient in Eqs. (3) and (4) is the detection fraction $\xi(\Lambda)$, which estimates the fraction of systems that we expect to successfully detect from some prior volume that extends past our detector's reach. The detection fraction quantifies selection biases, and so it is critical to accurately characterize. For a population described by parameters $\Lambda$, the detection fraction is

$$ \xi(\Lambda) = \int P_{\text{det}}(\theta)\pi(\theta|\Lambda)d\theta. \quad (A1) $$

Here, $P_{\text{det}}(\theta)$ is the detection probability: the probability that an event with parameters $\theta$ would be detected by a

particular search. The detection probability depends on the angular/sky position and orientation of the source binary, and crucially for our purposes, on the masses and redshift of a system, and, to a lesser degree, on the spins.

Given the non-ideal nature of the detector data, the variation in network sensitivity over time, and the complexity of both the signal waveforms and the search pipelines, an accurate estimate of $P_{\text{det}}(\theta)$ and $\xi(\Lambda)$ requires empirical methods, specifically the use of a large suite of simulated signals added to the data: injections. For analyses that focus on the BBH subpopulation in Section VI, we simulate compact binary signals from a reference BBH population and record which ones are successfully detected by the PYCBC, GSTLAL or MBTA search pipelines. We omit the cWB search from our volume estimate, since at present any detection of a binary merger was corroborated by a detection in the remaining pipelines. In addition, we also simulate compact binary signals from reference BNS, NSBH and IMBH populations. These injections include binaries with component masses in the range 1–600$M_{\odot}$, have spins that are isotropic in orientation and are uniform in comoving volume. Spins are drawn from a distribution that is uniform in the dimensionless spin magnitude up to a maximum of $\chi_{\text{max}} = 0.998$ for black holes and $\chi_{\text{max}} = 0.4$ for neutron stars. To control computational costs, the expected network signal to noise ratio (SNR) of each injection is computed using representative detector power spectral densities (PSDs) for O3. Injections with expected SNR below 6 are assumed not to be detected, and are thus removed from the set analyzed by the search pipelines. A thorough description of the injections and their underlying probability distribution is available in [318]. These injections are then combined into a single dataset as a mixture model [1, 318] in order to assess sensitivity across the entire parameter space and subpopulations. Our analyses in Sec. IV make use of these injections to estimate sensitivity.

Unlike previous synthetic simulation sets used in our population analysis following GWTC-2 [11], the injections used here model spins that are isotropically distributed in orientation and hence allow for orbital precession. Further, the maximum spin magnitude we assume for NS components, 0.4, is significantly larger than for previous injection sets used to estimate BNS merger rates [21]. That said, our injections have an effective $\chi_{\text{eff}}$ distribution that is narrow and centered at 0 while analyses using BNS populations with small NS spins inherently have $\chi_{\text{eff}} \approx 0$. Because the merger rate depends on spins primarily through the system's $\chi_{\text{eff}}$, the specific assumptions made about the spin distribution at low mass have modest impact on the inferred low-mass merger rate.

Following [94, 95, 105, 319], the point estimate for Eq. (A1) is calculated using a Monte Carlo integral over found injections:

$$ \hat{\xi}(\Lambda) = \frac{1}{N_{\text{inj}}} \sum_{j=1}^{N_{\text{found}}} \frac{\pi(\theta_j | \Lambda)}{p_{\text{draw}}(\theta_j)}, \quad (A2) $$
where $N_{inj}$ is the total number of injections, $N_{found}$ are the injections that are successfully detected, and $p_{draw}$ is the probability distribution from which the injections are drawn. When using this approach to estimate sensitivity, we marginalize over the uncertainty in $\xi(\Lambda)$ and ensure that the effective number of found injections remaining after population re-weighting is sufficiently high ($N_{eff} > 4N_{det}$) following [105]. We also compute (and some analyses like MS employ) semi-analytic approximations to the integrated network sensitivity $VT(\theta, \kappa) = \int dt dz dV_c/dz/(1+z) \langle P_{det}(\theta, z) \rangle (1+z)^\kappa$ for fiducial choices of $\kappa$, appropriate to characterize sensitivity to a population with a fixed redshift evolution.

For the O3 observing period, we characterize the found injections as those recovered with a FAR below the corresponding thresholds used in population analyses described in this paper (1 per year and 1 per 4 years) in either PYCBC, GSTLAL or MBTA. For the O1 and O2 observing periods, we supplement the O3 pipeline injections with mock injections drawn from the same distribution $p_{draw}$ above. For the mock injections, we calculate $P_{det}(m_1, m_2, z, \chi_{1,z}, \chi_{2,z})$ according to the semi-analytic approximation used in our analysis of GWTC-2 [25], based on a network signal-to-noise ratio threshold $\rho = 10$ and representative strain noise power spectral densities estimated from data recorded during the O1 and O2 observing runs. We combine O1, O2 and O3 injection sets ensuring a constant rate of injections across the total observing time [320].

## Appendix B: Population Model Details

In this section we provide details about the low-dimensional parameterized population models described

above in Section III. Each subsection includes a table with a summary of the parameters for that model and the prior distribution used for each parameter. The prior distributions are indicated using abbreviations: for example, $U(0, 1)$ translates to uniform on the interval $(0, 1)$, $LU(10^{-6}, 10^5)$ translates to log-uniform on the interval $10^{-6}, 10^5$, and $N(0, 1)$ translates to a Gaussian distribution with mean 0 and standard deviation of 1.

Using Monte Carlo summations over samples from each event's posterior distribution to approximate the integral in the likelihood (Eq. 4) results in statistical error in the likelihood estimates [105, 106]. In order to avoid including relics from to unconverged Monte Carlo integrals in the posterior distribution, we introduce a data-dependent constraint on the prior, determined by the number of effective samples used in the Monte Carlo integral. We define the effective number of samples as

$$N_{\text{eff}} = \frac{(\sum_i w_i)^2}{\sum_i w_i^2}, \qquad (\text{B1})$$

where $w_i$ is the weight for the $i$th event in the Monte Carlo integral.

For the TRUNCATED, POWER LAW + PEAK, POWER LAW + SPLINE, POWER LAW + DIP + BREAK, DEFAULT, GAUSSIAN, and POWER-LAW population models, we only assign nonzero likelihoods to points in parameter space with an effective sample size of at least the number of observed events in our event list. This is similar to the convergence constraints we enforce when computing sensitivity (see A) [105].

# 1. Details of mass population models

## a. TRUNCATED mass model

TRUNCATED mass model serves as the primary component for some of our mass models. The primary mass distribution for this model follows a power-law with spectral index $\alpha$, and with a sharp cut-off at the lower end $m_{min}$ and the upper end of the distribution $m_{max}$:

$$\pi(m_1|\alpha, m_{\min}, m_{\max}) \propto \begin{cases} m_1^{-\alpha} & m_{\min} < m_1 < m_{\max} \\ 0 & \text{otherwise,} \end{cases} \qquad (\text{B2})$$

Meanwhile, the mass ratio $q \equiv m_2/m_1$ follows a power-law distribution with spectral index $\beta_q$

$$\pi(q|\beta_q, m_{\min}, m_1) \propto \begin{cases} q^{\beta_q} & m_{\min} < m_2 < m_1 \\ 0 & \text{otherwise.} \end{cases} \qquad (\text{B3})$$

The parameters for this model are summarized in Table V. For this model, as well as further mass models where a prior on the total merger rate is not specified, the rate prior is proportional to $1/R$, or equivalently to $1/N$ in the notation of Eq. 3-4.
<table><caption>TABLE V. Summary of TRUNCATED model parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>α</td><td>Spectral index for the power-law of the primary mass distribution.</td><td>U(-4, 12)</td></tr><tr><td>β<sub>q</sub></td><td>Spectral index for the power-law of the mass ratio distribution.</td><td>U(-2, 7)</td></tr><tr><td>m<sub>min</sub></td><td>Minimum mass of the power-law component of the primary mass distribution.</td><td>U(2 M<sub>☉</sub>, 10 M<sub>☉</sub>)</td></tr><tr><td>m<sub>max</sub></td><td>Maximum mass of the power-law component of the primary mass distribution.</td><td>U(30 M<sub>☉</sub>, 100 M<sub>☉</sub>)</td></tr></tbody></table>

<table><caption>TABLE VI. Summary of POWER LAW + PEAK model parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>α</td><td>Spectral index for the power-law of the primary mass distribution.</td><td>U(-4, 12)</td></tr><tr><td>β<sub>q</sub></td><td>Spectral index for the power-law of the mass ratio distribution.</td><td>U(-2, 7)</td></tr><tr><td>m<sub>min</sub></td><td>Minimum mass of the power-law component of the primary mass distribution.</td><td>U(2 M<sub>☉</sub>, 10 M<sub>☉</sub>)</td></tr><tr><td>m<sub>max</sub></td><td>Maximum mass of the power-law component of the primary mass distribution.</td><td>U(30 M<sub>☉</sub>, 100 M<sub>☉</sub>)</td></tr><tr><td>λ<sub>peak</sub></td><td>Fraction of BBH systems in the Gaussian component.</td><td>U(0, 1)</td></tr><tr><td>μ<sub>m</sub></td><td>Mean of the Gaussian component in the primary mass distribution.</td><td>U(20 M<sub>☉</sub>, 50 M<sub>☉</sub>)</td></tr><tr><td>σ<sub>m</sub></td><td>Width of the Gaussian component in the primary mass distribution.</td><td>U(1 M<sub>☉</sub>, 10 M<sub>☉</sub>)</td></tr><tr><td>δ<sub>m</sub></td><td>Range of mass tapering at the lower end of the mass distribution.</td><td>U(0 M<sub>☉</sub>, 10 M<sub>☉</sub>)</td></tr></tbody></table>

b. POWER LAW + PEAK mass model

This is equivalent to Model C from [25]. The primary mass distribution is a truncated powerlaw, with the addition of tapering at the lower mass end of the distribution and a Gaussian component:

$$ \pi(m_1|\lambda_{\text{peak}}, \alpha, m_{\text{min}}, \delta_m, m_{\text{max}}, \mu_m, \sigma_m) = \left[ (1 - \lambda_{\text{peak}})\mathfrak{P}(m_1|\alpha, m_{\text{max}}) + \lambda_{\text{peak}}G(m_1|\mu_m, \sigma_m) \right] S(m_1|m_{\text{min}}, \delta_m). \quad (B4) $$

Here, $\mathfrak{P}(m_1|\alpha, m_{\text{max}})$ is a normalized power-law distribution with spectral index $-\alpha$ and high-mass cut-off $m_{\text{max}}$. Meanwhile, $G(m_1|\mu_m, \sigma_m)$ is a normalized Gaussian distribution with mean $\mu_m$ and width $\sigma_m$. The parameter $\lambda_{\text{peak}}$ is a mixing fraction determining the relative prevalence of mergers in $\mathfrak{P}$ and $G$. Finally, $S(m_1, m_{\text{min}}, \delta_m)$ is a smoothing function, which rises from 0 to 1 over the interval $(m_{\text{min}}, m_{\text{min}} + \delta_m)$:

$$ S(m | m_{\text{min}}, \delta_m) = \begin{cases} 0 & (m < m_{\text{min}}) \\ [f(m - m_{\text{min}}, \delta_m) + 1]^{-1} & (m_{\text{min}} \le m < m_{\text{min}} + \delta_m) \\ 1 & (m \ge m_{\text{min}} + \delta_m) \end{cases} \quad (B5) $$

with

$$ f(m', \delta_m) = \exp\left(\frac{\delta_m}{m'} + \frac{\delta_m}{m' - \delta_m}\right). \quad (B6) $$

The conditional mass ratio distribution in this model also includes the smoothing term:

$$ \pi(q | \beta, m_1, m_{\text{min}}, \delta_m) \propto q^{\beta q} S(qm_1 | m_{\text{min}}, \delta_m). \quad (B7) $$

The parameters for this model are summarized in Table VI.

c. POWER LAW + SPLINE mass model

The POWER LAW + SPLINE mass model explicitly applies a perturbation to a modified version of the fiducial POWER LAW + PEAK model that does not include the Gaussian peak [45]. Let $p(m_1|\alpha, m_{\text{min}}, m_{\text{max}}, \delta_m)$ be the modified POWER LAW + PEAK model without the Gaussian, then primary mass distribution for the POWER LAW + SPLINE model is given as:
<table><caption>TABLE VII. Summary of POWER LAW + SPLINE model parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>α</td><td>Spectral index for the power-law of the primary mass distribution.</td><td>U(-4, 12)</td></tr><tr><td>β<sub>q</sub></td><td>Spectral index for the power-law of the mass ratio distribution.</td><td>U(-2, 7)</td></tr><tr><td>m<sub>min</sub></td><td>Minimum mass of the power-law component of the primary mass distribution.</td><td>U(2 M<sub>☉</sub>, 10 M<sub>☉</sub>)</td></tr><tr><td>m<sub>max</sub></td><td>Maximum mass of the power-law component of the primary mass distribution.</td><td>U(30 M<sub>☉</sub>, 100 M<sub>☉</sub>)</td></tr><tr><td>δ<sub>m</sub></td><td>Range of mass tapering at the lower end of the mass distribution.</td><td>U(0 M<sub>☉</sub>, 10 M<sub>☉</sub>)</td></tr><tr><td>{f<sub>i</sub>}</td><td>Spline perturbation knot heights.</td><td>N(0, 1)</td></tr></tbody></table>

<table><caption>TABLE VIII. Summary of FLEXIBLE MIXTURES model parameters. All rates are in Gpc⁻³ yr⁻¹.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>w<sub>i</sub></td><td>Mixing weights.</td><td>Dirichlet(α), α<sub>1...N</sub> = 1/N</td></tr><tr><td>μ<sub>i</sub><sup>M</sup></td><td>Mean of the normal distribution modeling the chirp mass.</td><td>LU(5.2M<sub>☉</sub>, 65M<sub>☉</sub>)</td></tr><tr><td>σ<sub>i</sub><sup>M</sup></td><td>Scale of the normal distribution modeling the chirp mass.</td><td>U(0.02 μ<sub>i</sub><sup>M</sup>, 0.18 μ<sub>i</sub><sup>M</sup>)</td></tr><tr><td>μ<sub>i</sub><sup>sz</sup></td><td>Mean of the normal distribution modeling the aligned spin distribution.</td><td>U(-0.5, 0.5)</td></tr><tr><td>σ<sub>i</sub><sup>sz</sup></td><td>Scale of the normal distribution modeling the aligned spin distribution.</td><td>U(0.05, 0.6)</td></tr><tr><td>q<sub>i</sub><sup>min</sup></td><td>Minimum value of the mass ratio.</td><td>U(0.1, 0.95)</td></tr><tr><td>α<sub>i</sub><sup>q</sup></td><td>Slope of the power-law.</td><td>U(-7, 2)</td></tr><tr><td>R</td><td>Merger rate.</td><td>LU(1, 100)</td></tr></tbody></table>

$$p_{\text{PS}}(m_1|\alpha, m_{\text{min}}, m_{\text{max}}, \delta_m, \{f_i\}) = k p(m_1|\alpha, m_{\text{min}}, m_{\text{max}}, \delta_m) \exp(f(m_1|\{f_i\})). \quad (\text{B8})$$

Above, *k* is a normalization factor found by numerically integrating *p*<sub>PS</sub> over the range of allowed primary masses, *f*(*m*<sub>1</sub>|{*f*<sub>*i*</sub>) is the perturbation function we model with cubic splines, and {*f*<sub>*i*</sub>} are the heights of the *n* knots from which *f* is interpolated. The *n* knot locations are fixed, spaced linearly in log *m*<sub>1</sub> space from 2–100 M<sub>☉</sub>. We additionally restrict the perturbations to converge to the underlying distribution at the boundary nodes by fixing both *f*<sub>0</sub> and *f*<sub>*n*−1</sub> to be 0. We chose *n* = 20 to be the optimal number of knots for this analysis following the same procedure in [45], which adds a total of 18 additional parameters describing the perturbations to the underlying model. In addition to the primary mass, the conditional mass ratio distribution follows the same form as the POWER LAW + PEAK model defined in Eq. (B7). For each mass distribution inference with the POWER LAW + SPLINE model, we simultaneously fit the spin distribution with the DEFAULT model and the redshift evolution of the merger rate with the POWER LAW evolution model. The parameters and chosen prior distributions for the POWER LAW + SPLINE model are summarized in Table VII.

### d. FLEXIBLE MIXTURES model

The FLEXIBLE MIXTURES model, Vamana, predicts the population using a sum of weighted components. Each component is composed of a Gaussian, another Gaussian and a power-law to model the chirp mass, the aligned spins and the mass ratio respectively. The model is defined as

$$p(M, q, s_{1z}, s_{2z}|\lambda) = \sum_{i=1}^{N} w_i G(M|\mu_i^M, \sigma_i^M) G(s_{1z}|\mu_i^{sz}, \sigma_i^{sz}) \sigma_i^M G(s_{2z}|\mu_i^{sz}, \sigma_i^{sz}) P(q|\alpha_i^q, q_i^{\text{min}}, 1), \quad (\text{B9})$$

where *G* is the normal distribution and *P* is the truncated power-law. For the presented analysis we use *N* = 11 components. This choice maximises the marginal likelihood, however, the predicted population is robust for a wide range of *N*. For detailed description of this model see [49]. FLEXIBLE MIXTURES model uses a power-law to model the redshift evolution of the merger rate, as described in subsection B3. The merger rate has a uniform-in-log distributed prior; the prior distributions for parameters in Eq. B9 are summarized in Table VIII.
<table><caption>TABLE IX. Summary of Binned Gaussian Process model parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>μ</td><td>Mean log (Rate) in each bin.</td><td>N(0, 10)</td></tr><tr><td>σ</td><td>Amplitude of the covariance kernel.</td><td>N(0, 10)</td></tr><tr><td>log(l)</td><td>log (Length scale) of the covariance kernel.</td><td>N(-0.085, 0.93)</td></tr></tbody></table>

### e. Binned Gaussian Process model

The Binned Gaussian Process models the rate densities, $m_1 m_2 \frac{dR^i}{d\bar{m}_1 d\bar{m}_2} = n^i$, as a binned Gaussian Process where the index $i$ denotes a particular bin in the two-dimensional $\log m_1 - \log m_2$ parameter space [118, 119]. The bin edges in the analysis presented in the paper are located at [1, 2, 2.5, 3, 4, 5, 6.5, 8, 10, 15, 20, 30, 40, 50, 60, 70, 80, 100] M$_\odot$ with the assumption that $m_2 \le m_1$. The probabilistic model for the logarithm of the rate density in each bin is defined as

$$ \log n^i \sim N(\mu, \Sigma), \qquad (B10) $$

where $\mu$ is the mean of the Gaussian process and $\Sigma$ is the covariance matrix that correlates the bins. Each element of the covariance matrix $\Sigma$ is generated using a squared-exponential kernel $k(x, x')$ which is defined as

$$ k(x, x') = \sigma^2 \exp \left( -\frac{(x - x')^2}{2l^2} \right). \qquad (B11) $$

For the specific analysis here we take $x, x'$ to be the bin centers in $\log m$. The parameter $\sigma$ models the amplitude of the covariances while $l$ is a parameter that defines the length scales over which bins are correlated. The prior distribution chosen here for the length scale is a log-normal distribution with a mean that is the average between the minimum bin spacing

$$ \Delta_{\min} \equiv \min_{m_1, m_2} \Delta \log m \qquad (B12) $$

and the maximum bin spacing

$$ \Delta_{\max} \equiv \max_{m_1, m_2} \Delta \log m \qquad (B13) $$

with a standard deviation of $(\Delta_{\max} - \Delta_{\min})/\sqrt{4}$. This constrains (at “2-$\sigma$” in the prior) the correlation length for the GP to lie between “one bin” and “all the bins.” For our analyses presented in the paper, the mean and standard deviation are -0.085 and 0.93 respectively. The Binned Gaussian Process model assumes a redshift distribution such that the overall merger rate of compact binaries is uniform-in-comoving volume. The spin distributions for each component are isotropic in direction and uniform in the spin magnitude with a maximum spin of 0.998 for BHs and 0.4 for NSs; the prior distribution for the relevant parameters in Equations B10 and B11 is summarized in Table IX.

### f. Power Law + Dip + Break model

The Power Law + Dip + Break model explicitly searches for separation in masses between two subpopulations by employing a broken power law with a dip at the location of the power law break. As described in [108] and [109], the dip is modeled by a notch filter with depth $A$, which is fit along with the other model parameters in order to determine the existence and depth of a potential mass gap. No gap corresponds to $A = 0$, whereas $A = 1$ corresponds to precisely zero merger rate over some interval. Power Law + Dip + Break also employs a low-pass filter at high masses to allow for a tapering of the mass spectrum, which has the effect of a smooth second break to the power law.

The PDB model assumes a merger rate that is uniform in comoving volume. It also assumes a spin distribution with isotropically oriented component spins and uniform component spin magnitudes.
<table><caption>TABLE X. Summary of POWER LAW + DIP + BREAK model parameters. The first entries describe the mass distribution parameters, and the last two entries describe the spin distribution parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>α<sub>1</sub></td><td>Spectral index for the power-law of the mass distribution at low mass.</td><td>U(-8, 2)</td></tr><tr><td>α<sub>2</sub></td><td>Spectral index for the power-law of the mass distribution at high mass.</td><td>U(-3, 2)</td></tr><tr><td>A</td><td>Lower mass gap depth.</td><td>U(0, 1)</td></tr><tr><td>M<sub>low</sub><sup>gap</sup></td><td>Location of lower end of the mass gap.</td><td>U(1.4M<sub>☉</sub>, 3M<sub>☉</sub>)</td></tr><tr><td>M<sub>high</sub><sup>gap</sup></td><td>Location of upper end of the mass gap</td><td>U(3.4M<sub>☉</sub>, 9M<sub>☉</sub>)</td></tr><tr><td>η<sub>low</sub></td><td>Parameter controlling how the rate tapers at the low end of the mass gap</td><td>50</td></tr><tr><td>η<sub>high</sub></td><td>Parameter controlling how the rate tapers at the low end of the mass gap.</td><td>50</td></tr><tr><td>η</td><td>Parameter controlling tapering the truncated power law at high mass</td><td>U(-4, 12)</td></tr><tr><td>β</td><td>Spectral index for the power-law-in-mass-ratio pairing function.</td><td>U(-2, 7)</td></tr><tr><td>m<sub>min</sub></td><td>Minimum mass of the power-law component of the mass distribution.</td><td>U(1 M<sub>☉</sub>, 1.4 M<sub>☉</sub>)</td></tr><tr><td>m<sub>max</sub></td><td>Maximum mass of the power-law component of the mass distribution.</td><td>U(35 M<sub>☉</sub>, 100 M<sub>☉</sub>)</td></tr><tr><td>a<sub>max,NS</sub></td><td>Maximum allowed component spin for objects with mass &lt; 2.5 M<sub>☉</sub></td><td>0.4</td></tr><tr><td>a<sub>max,BH</sub></td><td>Maximum allowed component spin for objects with mass ≥ 2.5 M<sub>☉</sub></td><td>1</td></tr></tbody></table>

The joint mass distribution in this model has the form:

$$p(m_1, m_2) \propto p(m_1)p(m_2)(m_2/m_1)^\beta, \qquad (B14)$$

$$p(m) \propto p_{\text{pl}}(m)n(m)\ell(m), \qquad (B15)$$

$$n(m) = 1 - \frac{A}{(1 + (M_{\text{low}}^{\text{gap}}m)^{\eta_{\text{low}}})(1 + (M_{\text{high}}^{\text{gap}}m)^{\eta_{\text{high}}})}, \quad \text{and} \qquad (B16)$$

$$\ell(m) = \frac{1}{1 + (m/m_{\text{max}})^{\eta}}. \qquad (B17)$$

where $p_{\text{pl}}(m)$ is a broken power law with exponents $\alpha_1$ between $m_{\text{min}}$ and $M_{\text{low}}^{\text{gap}}$ and $\alpha_2$ between $M_{\text{low}}^{\text{gap}}$ and $m_{\text{max}}$. The parameters for this model are summarized in Table X.

### g. Neutron star mass models

The mass models adopted for the BNS and NSBH events in Sec. V assume a basic mass distribution that is common to all NSs, with random pairing into compact binaries. The basic mass distribution is taken to be either a power law or, inspired by the shape of the Galactic BNS mass distribution [39–41], a Gaussian. The BH mass distribution is fixed to be uniform between 3 and 60 $M_{\odot}$. The NS mass distribution analysis assumes definite source classifications for the events. Thus, the joint mass distribution takes the form

$$p(m_1, m_2) \propto \begin{cases} p(m_1)p(m_2) & \text{if BNS} \\ U(3M_{\odot}, 60M_{\odot})p(m_2) & \text{if NSBH,} \end{cases} \qquad (B18)$$

with $p(m)$ either a power law with exponent $\alpha$, minimum mass $m_{\text{min}}$ and maximum mass $m_{\text{max}}$, or a Gaussian with a peak of width $\sigma$ at $\mu$, plus sharp minimum and maximum mass cutoffs $m_{\text{min}}$, $m_{\text{max}}$. We call these models POWER and PEAK, respectively. Their hyper-parameters, and the choices for their prior distributions, are listed in Table XI. We additionally impose the constraint $m_{\text{min}} \le \mu \le m_{\text{max}}$ on the PEAK model. Besides the flat $m_{\text{max}}$ prior described in the table, for the analyses excluding GW190814 we use a prior proportional to the cumulative distribution function of $M_{\text{max,TOV}}$, i.e., $p(m_{\text{max}}) \propto \int_{m_{\text{max}}}^{\infty} dM_{\text{max,TOV}} p(M_{\text{max,TOV}})$. This enforces our expectation that the NS masses in the gravitational-wave population should not exceed $M_{\text{max,TOV}}$.
<table><caption>TABLE XI. Summary of POWER and PEAK NS mass model parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>α</td><td>Spectral index for the power-law in the POWER NS mass distribution.</td><td>U(-4, 12)</td></tr><tr><td>m<sub>min</sub></td><td>Minimum mass of the NS mass distribution.</td><td>U(1.0M<sub>☉</sub>, 1.5M<sub>☉</sub>)</td></tr><tr><td>m<sub>max</sub></td><td>Maximum mass of the NS mass distribution.</td><td>U(1.5M<sub>☉</sub>, 3.0M<sub>☉</sub>)</td></tr><tr><td>μ</td><td>Location of the Gaussian peak in the PEAK NS mass distribution.</td><td>U(1.0M<sub>☉</sub>, 3.0M<sub>☉</sub>)</td></tr><tr><td>σ</td><td>Width of the Gaussian peak in the PEAK NS mass distribution.</td><td>U(0.01M<sub>☉</sub>, 2.00M<sub>☉</sub>)</td></tr></tbody></table>

<table><caption>TABLE XII. Summary of DEFAULT spin parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>μ<sub>χ</sub></td><td>Mean of the Beta distribution of spin magnitudes.</td><td>U(0,1)</td></tr><tr><td>σ<sup>2</sup><sub>χ</sub></td><td>Variance of the Beta distribution of spin magnitudes.</td><td>U(0.005,0.25)</td></tr><tr><td>ζ</td><td>Mixing fraction of mergers from truncated Gaussian distribution.</td><td>U(0,1)</td></tr><tr><td>σ<sub>t</sub></td><td>Width of truncated Gaussian, determining typical spin misalignment.</td><td>U(0.1,4)</td></tr></tbody></table>

## 2. Details of spin population models

### a. DEFAULT spin model

This model was introduced in [25]. Following [113], the dimensionless spin magnitude distribution is taken to be a Beta distribution,

$$ \pi(\chi_{1,2}|\alpha_{\chi}, \beta_{\chi}) = \text{Beta}(\alpha_{\chi}, \beta_{\chi}), \qquad (B19) $$

where $\alpha_{\chi}$ and $\beta_{\chi}$ are the standard shape parameters that determine the distribution's mean and variance. The Beta distribution is convenient because it is bounded on (0,1). The distributions for $\chi_1$ and $\chi_2$ are assumed to be the same. Following [114], we define $z_i = \cos \theta_i$ as the cosine of the tilt angle between component spin and a binary's orbital angular momentum, and assume that $\mathbf{z}$ is distributed as a mixture of two populations:

$$ \pi(\mathbf{z}|\zeta, \sigma_t) = \zeta G_t(\mathbf{z}|\sigma_t) + (1-\zeta)\mathcal{I}(\mathbf{z}). \qquad (B20) $$

Here, $\mathcal{I}(z)$ is an isotropic distribution, while $G_t(z|\sigma_t)$ is a truncated two-dimensional Gaussian, peaking at $\mathbf{z} = 0$ (perfect alignment) with width $\sigma_t$. The mixing parameter $\zeta$ controls the relative fraction of mergers drawn from the isotropic distribution and Gaussian subpopulations. The isotropic subpopulation is intended to accommodate dynamically assembled binaries, while $G_t$ is a model for field mergers. The parameters for this model and their priors are summarized in Table XII. Additional constraints to the priors on $\mu_{\chi}$ and $\sigma_{\chi}^2$ are applied by setting $\alpha_{\chi}, \beta_{\chi} > 1$.

### b. GAUSSIAN spin model

In addition to the distribution of component spin magnitudes and tilts, we explore the distribution of the effective inspiral spin parameter $\chi_{\text{eff}}$ and the effective precession spin parameter $\chi_{\text{p}}$. In particular, we wish to measure the mean and variance of each parameter, and so model the joint distribution of $\chi_{\text{eff}}$ and $\chi_{\text{p}}$ as a bivariate Gaussian:

$$ \pi(\chi_{\text{eff}}, \chi_{\text{p}}|\mu_{\text{eff}}, \sigma_{\text{eff}}, \mu_{\text{p}}, \sigma_{\text{p}}, \rho) \propto G(\chi_{\text{eff}}, \chi_{\text{p}}|\mu, \Sigma). \qquad (B21) $$

The mean of this distribution is $\mu = (\mu_{\text{eff}}, \mu_{\text{p}})$, and its covariance matrix is

$$ \Sigma = \begin{pmatrix} \sigma_{\text{eff}}^2 & \rho \sigma_{\text{eff}} \sigma_{\text{p}} \\ \rho \sigma_{\text{eff}} \sigma_{\text{p}} & \sigma_{\text{p}}^2 \end{pmatrix}. \qquad (B22) $$

The population parameters governing this model and their corresponding priors are shown in Table XIII. Equation (B21) is truncated to the physically allowed range of each effective spin parameter, with $\chi_{\text{eff}} \in (-1, 1)$ and $\chi_{\text{p}} \in (0, 1)$. All results in the main text using the GAUSSIAN model are obtained while simultaneously fitting for the
<table><caption>TABLE XIII. Summary of GAUSSIAN spin parameters. The $\chi_{\text{eff,min}}$ and $\zeta$ parameters appear only in variants of the GAUSSIAN model, as discussed below.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>μ<sub>eff</sub></td><td>Mean of the χ<sub>eff</sub> distribution.</td><td>U(-1, 1)</td></tr><tr><td>σ<sub>eff</sub></td><td>Standard deviation of the χ<sub>eff</sub> distribution.</td><td>U(0.05,1)</td></tr><tr><td>μ<sub>p</sub></td><td>Mean of the χ<sub>p</sub> distribution.</td><td>U(0.05, 1)</td></tr><tr><td>σ<sub>p</sub></td><td>Standard deviation of the χ<sub>p</sub> distribution.</td><td>U(0.05, 1)</td></tr><tr><td>ρ</td><td>Degree of correlation between χ<sub>eff</sub> and χ<sub>p</sub>.</td><td>U(-0.75, 0.75)</td></tr><tr><td>χ<sub>eff,min</sub></td><td>Lower truncation bound on χ<sub>eff</sub>.</td><td>U(-1, μ<sub>eff</sub>)</td></tr><tr><td>ζ</td><td>Non-vanishing mixture fraction in Eq. (18).</td><td>U(0, 1)</td></tr></tbody></table>

BBH mass distribution, assuming the POWER LAW + PEAKmodel, and the evolving redshift distribution model in Appendix B3 below.

Two variants of this model are additionally discussed in Sect. VII. In the first, Eq. (B21) is modified such that the effective inspiral spin parameter is truncated not on the interval $(-1, 1)$, but on $(\chi_{\text{eff,min}}, 1)$, where $\chi_{\text{eff,min}}$ is another parameter to be inferred by the data. The second variant, inspired by [194] and [195] and defined in Eq. (18), alternatively treats the $\chi_{\text{eff}}$ distribution as a mixture between a bulk component with a variable mean and width and a narrow zero-spin component centered on $\chi_{\text{eff}} = 0$. In this second variant, we measure only the marginal $\chi_{\text{eff}}$ distribution, implicitly assuming that the remaining spin degree of freedom are distributed uniformly and isotropically. As $\chi_{\text{eff}}$ is the primary spin measurable, we do not expect this implicit prior to have a strong effect.

### 3. Redshift Evolution Model

The POWER-LAW redshift evolution model parameterizes the merger rate density per comoving volume and source time as [96]

$$ \mathcal{R}(z) = \mathcal{R}_0(1+z)^{\kappa}, \qquad (B23) $$

where $\mathcal{R}_0$ denotes the merger rate density at $z=0$. This implies that the redshift distribution is

$$ \frac{dN}{dz} = C \frac{dV_c}{dz} (1+z)^{\kappa-1}, \qquad (B24) $$

where $dV_c/dz$ is the differential comoving volume, and $C$ is related to $\mathcal{R}_0$ by

$$ \mathcal{R}_0 = C \frac{dV_c}{dz} \left[ \int_0^{z_{\text{max}}} \frac{dV_c}{dz} (1+z)^{\kappa-1} \right]^{-1}. \qquad (B25) $$

We adopt $z_{\text{max}} = 2.3$ as this is a conservative upper bound on the redshift at which we could detect BBH systems during O3, for both detection thresholds used in this work. We employ a uniform prior on $\kappa$ centered at $\kappa = 0$. We take a sufficiently wide prior so that the likelihood is entirely within the prior range, $\kappa \in (-6, 6)$.

### 4. Models with multiple independent components

#### a. MULTI SOURCE model

The MULTI SOURCE model, introduced in [321], extends the MULTI SPIN BBH model introduced in [11] to include additional subpopulations for BNS and NSBH systems. Each subpopulation (two for BBH, one for BNS, and one for NSBH) is assumed to have an independent rate parameter.

The BBH subpopulation is itself a mixture of two subpopulations, i. a power law mass distribution $m_1^{-\alpha} q^{\beta}$ truncated to a range $[m_{\text{min,BBH}}, m_{\text{max,BBH}}]$ which is inferred from the data, and ii. a Gaussian in $(m_1, m_2)$ with independent mean and standard deviation parameters $\mu_{m_1,\text{BBH}}$, $\mu_{m_2,\text{BBH}}$, $\sigma_{m_1,\text{BBH}}$, $\sigma_{m_2,\text{BBH}}$. Both subpopulations, and both binary components within them follow independent DEFAULT spin models, with $\zeta \equiv 1$.

Two more bivariate Gaussians in $m_1, m_2$ are used to model BNS and NSBH. The BH in NSBH follow a Gaussian mass distribution, with free parameters $\mu_{m,\text{NSBH}}$, $\sigma_{m,\text{NSBH}}$. As with BBH, these BH follow an independent DEFAULT
spin model with $\zeta \equiv 1$. All three types of NS (two in BNS and one in NSBH) follow the same Gaussian mass distribution, with free parameters $\mu_{m,NS}$, $\sigma_{m,NS}$, $m_{max,NS}$ (the minimum mass is assumed to be $1 M_{\odot}$). Each type of NS follows an independent DEFAULT spin model. To stay within astrophysically plausible spins, the magnitude distributions are scaled down to $\chi_{max} = 0.05$. Since NS spin tilts are not well measured, we set $\zeta \equiv 0$, assuming they are isotropic, which has the effect of not wasting any samples from parameter estimation.

In addition to any mass cutoffs mentioned above, all BHs component masses are assumed to lie on the range $[2, 100] M_{\odot}$, with those in NSBHs further restricted to $[2, 50] M_{\odot}$ due to our limited injections.

Priors for all parameters are given in Table XIV.

<table><caption>TABLE XIV. Summary of MULTI SOURCE model parameters. All rates are in Gpc$^{-3}$yr$^{-1}$, and all masses in $M_{\odot}$. Rate, mass, and spin hyperparameters are separated by horizontal lines.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>R<sub>BBH,pl</sub></td><td>Local merger rate for the BBH power-law subpopulation.</td><td>U(0, 1000)</td></tr><tr><td>R<sub>BBH,g</sub></td><td>Local merger rate for the BBH Gaussian subpopulation.</td><td>U(0, 1000)</td></tr><tr><td>R<sub>BNS</sub></td><td>Local merger rate for the BNS subpopulation.</td><td>U(0, 2000)</td></tr><tr><td>R<sub>NSBH</sub></td><td>Local merger rate for the NSBH subpopulation.</td><td>U(0, 500)</td></tr><tr><td>α</td><td>Primary mass spectral index for the BBH power-law subpopulation.</td><td>U(-4, 12)</td></tr><tr><td>β</td><td>Mass ratio spectral index for the BBH power-law subpopulation</td><td>U(-4, 10)</td></tr><tr><td>m<sub>min,BBH,pl</sub></td><td>Minimum mass of the BBH power-law subpopulation.</td><td>U(2, 10)</td></tr><tr><td>m<sub>max,BBH,pl</sub></td><td>Maximum mass of the BBH power-law subpopulation.</td><td>U(30, 100)</td></tr><tr><td>μ<sub>m1,BBH,g</sub> (μ<sub>m2,BBH,g</sub>)</td><td>Centroid of the primary (secondary) mass distribution for the BBH Gaussian subpopulation</td><td>U(20, 50)</td></tr><tr><td>σ<sub>m1,BBH,g</sub> (σ<sub>m2,BBH,g</sub>)</td><td>Width of the primary (secondary) mass distribution for the BBH Gaussian subpopulation</td><td>U(0.4, 20)</td></tr><tr><td>μ<sub>m,NSBH</sub></td><td>Centroid of the BH mass distribution for NSBH</td><td>U(3, 50)</td></tr><tr><td>σ<sub>m,NSBH</sub></td><td>Width of the BH mass distribution for NSBH</td><td>U(0.4, 20)</td></tr><tr><td>μ<sub>m,NS</sub></td><td>Centroid of the NS mass distribution</td><td>U(1, 3)</td></tr><tr><td>σ<sub>m,NS</sub></td><td>Width of the NS mass distribution</td><td>U(0.05, 3)</td></tr><tr><td>m<sub>max,NS</sub></td><td>Maximum mass of all NS.</td><td>U(2, 3)</td></tr><tr><td>μ<sub>χ1,BBH,pl</sub> (μ<sub>χ2,BBH,pl</sub>)</td><td>Mean of the Beta distribution of primary (secondary) spin magnitudes for the BBH Gaussian sub-population.</td><td>U(0, 1)</td></tr><tr><td>σ<sup>2</sup><sub>χ1,BBH,pl</sub> (σ<sup>2</sup><sub>χ2,BBH,pl</sub>)</td><td>Variance of the Beta distribution of primary (secondary) spin magnitudes for the BBH Gaussian sub-population.</td><td>U(0, 0.25)</td></tr><tr><td>σ<sub>t1,BBH,pl</sub> (σ<sub>t2,BBH,pl</sub>)</td><td>Width of truncated Gaussian, determining typical primary (secondary) spin misalignment for the BBH Gaussian sub-population.</td><td>U(0, 4)</td></tr><tr><td>μ<sub>χ1,BBH,g</sub> (μ<sub>χ2,BBH,g</sub>)</td><td>Mean of the Beta distribution of primary (secondary) spin magnitudes for the BBH Gaussian sub-population.</td><td>U(0, 1)</td></tr><tr><td>σ<sup>2</sup><sub>χ1,BBH,g</sub> (σ<sup>2</sup><sub>χ2,BBH,g</sub>)</td><td>Variance of the Beta distribution of primary (secondary) spin magnitudes for the BBH Gaussian sub-population.</td><td>U(0, 0.25)</td></tr><tr><td>σ<sub>t1,BBH,g</sub> (σ<sub>t2,BBH,g</sub>)</td><td>Width of truncated Gaussian, determining typical primary (secondary) spin misalignment for the BBH Gaussian sub-population.</td><td>U(0, 4)</td></tr><tr><td>μ<sub>χ,NSBH</sub></td><td>Mean of the Beta distribution of spin magnitudes for BH in the NSBH sub-population.</td><td>U(0, 1)</td></tr><tr><td>σ<sup>2</sup><sub>χ,NSBH</sub></td><td>Variance of the Beta distribution of spin magnitudes for BH in the NSBH sub-population.</td><td>U(0, 0.25)</td></tr><tr><td>σ<sub>t,NSBH</sub></td><td>Width of truncated Gaussian, determining typical primary (secondary) spin misalignment for BH in the NSBH sub-population.</td><td>U(0, 4)</td></tr><tr><td>μ<sub>χ,NSBH</sub></td><td>Mean of the Beta distribution of spin magnitudes for NS in the NSBH sub-population.</td><td>U(0, 0.05)</td></tr><tr><td>σ<sup>2</sup><sub>χ,NSBH</sub></td><td>Variance of the Beta distribution of spin magnitudes for NS in the NSBH sub-population.</td><td>U(0, 0.0125)</td></tr><tr><td>μ<sub>χ1,BNS</sub> (μ<sub>χ2,BNS</sub>)</td><td>Mean of the Beta distribution of primary (secondary) spin magnitudes in the BNS sub-population.</td><td>U(0, 0.05)</td></tr><tr><td>σ<sup>2</sup><sub>χ1,BNS</sub> (σ<sup>2</sup><sub>χ2,BNS</sub>)</td><td>Variance of the Beta distribution of primary (secondary) spin magnitudes in the BNS sub-population.</td><td>U(0, 0.0125)</td></tr></tbody></table>
# Appendix C: Validation studies

We employ several methods to validate our calculations, notably including comparing results from multiple independent analyses; reproducing previous work through O3a [11]; assessing the sensitivity of our results to threshold choices (changing from $1 \text{ yr}^{-1}$ to $0.2 \text{ yr}^{-1}$ for BBH; or from $0.25 \text{ yr}^{-1}$ to $1 \text{ yr}^{-1}$ for analyses containing NS); and performing posterior predictive checks as in our analysis of GWTC-2 [11]. Though these specific technical checks will not be described here, some of these checks can be reproduced with the data release associated with this paper.

Below, we describe additional validation studies we have performed to assess whether our results for merger rates are sensitive to the choice of threshold; waveform systematics; or updates to our sensitivity model.

## 1. Effects of the Spin Distribution on Merger Rates Across All Masses

In principle, the mass, spin, and redshift distributions of binaries should be fit simultaneously in order to avoid systematic biases in inferred distribution parameters caused by correlations between measurements of these intrinsic parameters [106, 322–328]. However, fixing one or more of these distributions to a realistic form typically introduces biases that have little impact on the parameters of interest. We therefore seek to determine if our choice to fix the spin distribution for the PDB and BGP models has introduced any significant biases in our inference of the mass distribution and overall merger rate.

We compare the PDB analysis presented in Section IV with an analysis that utilizes the same mass and redshift distribution but fits for the spin distribution rather than fixing it to one that assumes isotropic and uniformly distributed component spins. For this, we apply the DEFAULT [114] spin model described in Section III C 2 and Appendix B 2a. The resulting fit is compared to the fiducial analysis in Figure 24.

We find some differences between the fixed-spin and fit-spin analyses. Firstly, the hyperposterior for the fit-spin analysis is broader than that of the fit-spin analysis, presumably due to an increase in free parameters. Second, some hyperparameters exhibit a slight shift. The most notable shifts are in the rate and upper gap edge parameters. The shift in the rate is to be expected because the fit to the DEFAULT model favors lower spinning objects. Since the detectors are slightly less sensitive to low-spin objects, more support for those objects implies a higher astrophysical rate. Nonetheless, all hyperposterior differences are well within statistical uncertainty, so we conclude that both the fixed-spin and fit-spin cases are acceptable, and use the fixed-spin case for our fiducial results for simplicity.

In a preprint version of this paper, the rates reported

by the PDB model were higher than in this version of the paper. This was due to an incorrect approximation of detector sensitivity in the region of highly spinning low mass objects. This approximation has been removed in the current version, and the effect has been to lower the BNS and NSBH rates by nearly $1\sigma$, lowering the overall rate by a similar amount.

## 2. NS mass distribution including marginal events

If we loosen the FAR threshold to < $1 \text{ yr}^{-1}$ so as to include the marginal events GW190917 and GW190426, and repeat the analysis of Sec. V B, the inferred NS mass distribution is virtually unchanged. This can be seen in Fig. 25, which compares the posterior population distributions inferred with and without the marginal events. Traces from the posterior population distribution with respect to the original FAR threshold are also shown. This alternative analysis strongly suggests that substantial uncertainties in the merger rate versus mass dominate our error budget; the handful of observations made to date is not sufficient to overcome the strong impact of our highly uncertain model priors. Moreover, the masses of the NS secondaries in the marginal events are poorly constrained relative to those in GW170817, GW190425, GW200105 and GW200115, such that the FAR < $0.25 \text{ yr}^{-1}$ events continue to drive the inference.

## 3. Merger rates including subthreshold triggers

In the main text, our merger rates were calculated after adopting a fixed significance threshold to identify confident events, then fitting population model families to the recovered events’ posteriors. By design, such an approach depends on the threshold. Here we employ an alternative threshold-free method of rate estimation which lacks potential biases from an arbitrary choice of significance threshold [12].

We extend methods from GWTC-2.1 [3], also applied to the discovery of GW200105 and GW200115 [329], to estimate the event rate from the full set of triggers (including subthreshold triggers) from a specific binary merger search: here, GSTLAL [50–52]. In doing so, we allow for population distributions that fit our observations and account for still-considerable uncertainty in the mass distribution, rather than adopting a fixed population model with fixed model hyper-parameters. Compared to previous publications, the results presented in this section update the BBH merger rates presented in GWTC-2.1 by including O3b events [1]. We also update the NSBH rate quoted in [329] by incorporating all O3 triggers, rather than as previously truncating to the first 9 months of O3.

We use a multi-component mixture model [121] to construct the posterior of astrophysical counts of CBC events by assuming that foreground and background events are
<figure><img src="image_27.png" /><figcaption>FIG. 24. Corner plot of inferred PDB mass and rate hyperparameters under an analysis that fixes the spin distribution (blue) and simultaneously fits the spin distribution using the DEFAULT model (orange). The fit-spin hyperposterior is slightly shifted and widened when compared to the fixed-spin case, but all changes are within statistical uncertainties.</figcaption></figure>

independent Poisson processes. We then estimate the
space-time volume sensitivity of the pipeline using simu-
lated events which are re-scaled to an astrophysical popu-
lation model [319]. We then compute the rates as the ra-
tio of the counts to VT. In order to marginalize over pop-
ulation hyper-parameters we compute several VT's, each
corresponding to a population hyper-parameter sample
drawn from the inferred hyper-posterior for the astro-

physical population model. Finally, we integrate over the
count posterior obtained for each of these samples with
an appropriate weight, effectively marginalizing over the
<figure><img src="image_28.png" /><figcaption>FIG. 25. Inferred neutron star mass distribution with and without the marginal events GW190426 and GW190917. Top panel: Median and 90% confidence region of the inferred NS mass distribution for the POWER model, using the event list at a FAR threshold of 0.25 yr⁻<fnref n="1" /> (blue) and 1 yr⁻<fnref n="1" /> (orange). Traces from the posterior population distribution with respect to the stricter FAR threshold are plotted in grey. Bottom panel: Same as the top panel but for the PEAK model. The inclusion of the marginal events has a negligible impact on the inferred mass distribution.</figcaption></figure>

population hyper-parameters:

$$ 
\begin{aligned}
p(R|\vec{x}) &= \int p(R|\vec{\Lambda}, \vec{x})p(\Lambda|\{d\})d\vec{\Lambda} \\
&= \int VT p(N|\vec{x})p(VT|\vec{\Lambda})p(\vec{\Lambda}|\{d\})d(VT)d\vec{\Lambda} \\
&= \sum_{i,j} VT_{ij} \times p(N_{ij}|\vec{x}),
\end{aligned} 
\quad (C1) $$

where $\vec{x}$ is the complete set of triggers (including sub-threshold triggers) and $\{d\}$ is the set of data from gravitational-wave detections used in population model inference, as in IIIB. The astrophysical count posterior is given by $p(N|\vec{x})$, where $N = R \times VT$; we evaluate by sampling via $N_{ij} = R \times VT_{ij}$ where $VT_{ij}$ is the i'th $VT$ sample drawn from $p(VT|\vec{\Lambda}_j)$ for the j'th hyper-parameter sample $\vec{\Lambda}_j$ drawn from the inferred hyper-posterior $p(\vec{\Lambda}|\{d\})$.<fnref n="1" /> Following [121], we take the dis-

tribution $p(VT|\vec{\Lambda}_j)$ to be

$$ p(VT|\vec{\Lambda}_j) = \frac{1}{VT\sqrt{2\pi}\sigma^2} \exp \left[ -\frac{[\ln VT - \ln \langle VT \rangle(\vec{\Lambda}_j)]^2}{2\sigma^2} \right], \quad (C2) $$

where $\langle VT(\vec{\Lambda}_j) \rangle$ is calculated by re-weighting simulated sources to an astrophysical population with hyper-parameter $\vec{\Lambda}_j$, and $\sigma$ is the quadrature sum of a calibration error of 10% [91] and Monte-Carlo uncertainty.

Using hyper-parameter samples from the posterior inferred using the PP model with data through the end of O3, as in Section VI and imposing a Jeffreys prior $\propto N^{-1/2}$ on the astrophysical counts, we compute a BBH merger rate of 24.81–63.58 Gpc⁻³yr⁻<fnref n="1" />. A similar calculation for the BGP model, again with a Jeffreys prior $\propto N^{-1/2}$ imposed on the astrophysical counts, yields an NSBH merger rate of 14.57–187.96 Gpc⁻³yr⁻<fnref n="1" />, which is consistent with 11–140 Gpc⁻³yr⁻<fnref n="1" />, the joint inference for the NSBH merger rate presented in the main text. We also compute a BNS merger rate using a fixed population of BNS’s, distributed uniformly in component masses that lie with in 1 to 2.5 $M_\odot$. This uses the same multi-component mixture model [121] as described above, with the only difference being, that instead of marginalizing over population hyperparameters like with the BBH and NSBH merger rates, for BNS, we use a fixed population. Hence, it updates the BNS merger rate reported in GWTC-2.1 [3] by including all of O1 through O3 instead of truncating at O3a. We report a BNS merger rate of 28.76–462.23 Gpc⁻³yr⁻<fnref n="1" /> which is consistent with the GWTC 2.1 Rate as well as the other BNS rates quoted in this paper that were computed from only high-significance triggers.

## 4. Effect of Waveform Systematics on Population

All O3b BBH events analyzed in this paper have source properties inferred using two different waveform models: SEOBNRv4PHM [87] and IMRPhenomXPHM [88], both of which include effects of higher-order multipole moments and spin precession. The posterior distribution for each event is then checked for consistency between waveform models before use in our analyses.

The event GW200129 is the highest SNR event exhibiting notable inconsistencies between the source properties inferred with the two waveform models. The event analysis using IMRPhenomXPHM infers much more support for unequal masses and precessing spins relative to the analysis using SEOBNRv4PHM. See [1] for an extended discussion of these systematic differences.

To test if the inferred BBH spin population depends on the waveform model chosen for this event, we repeat our O3 population inference using the PP model for three different choices of waveform model for GW200129: IMR-PhenomXPHM, SEOBNRv4PHM, and a mix of the two.

<footnote n="1">¹ The extra $VT$ factor in Eq. (C1) arises from the Jacobian $dN/dR$.</footnote>
<figure><img src="image_29.png" /><figcaption>FIG. 26. Inferred differential merger rate as a function of the cosine of the tilt angle ($t_i$), where $i$ indexes the body of the binary. We demonstrate that the differences in the posterior distribution for GW200129's spin parameters have a minimal effect on the inferred spin tilt population. The population is inferred using posterior distributions for GW200129 using the IMRPhenomXPHM waveform model (orange), SEOB-NRv4PHM model (blue), and a mixture of both (green). Dashed lines are 90% credible intervals.</figcaption></figure>

As shown in Fig. 26, the inferred spin population is not significantly affected by changes in the waveform model for this event.

## 5. Impact of Sensitivity on Redshift Evolution Inference

As noted in Sec. II, one change in the sensitivity estimation procedure between this work and our previous study of GWTC-2 [11] is the use of injections that account for the effect of precession and as well as updates to our detection pipelines as detailed in [1]. Since precession was not included in the injections used in [11], the full spin distribution could not be reweighted to calculate the sensitivity via Equation A2, and thus, for the purposes of sensitivity estimation, an approximation was made that $S_{x,y} \in (-0.5, 0.5)$. Since we now use precessing injections, we do the reweighting procedure including the full spin distribution as a function of $\Lambda$. To test if this difference in our sensitivity estimation procedure is responsible for the change in the inferred redshift evolution, we repeat the population analysis reported in Sec. VI, using our updated sensitivity model, but only including events analyzed in the GWTC-2 populations study [11]. From this analysis, we infer $\kappa > 0$ at 97.6% credibility, as opposed to the 85% credibility reported in [11], indicating a much stronger preference for a merger rate increasing with redshift. We conclude that the differences between our current results for the evolution of the BBH merger rate and those reported in [11] are due to improvements to our sensitivity model rather than the presence of the additional events in GWTC-3.

In Fig. 27 we compare the redshift dependence of our current sensitivity model to that of the sensitivity model used in [11]. To make this comparison, we reweight the injections used in [11] to the same spin distribution assumed in that study, and assuming a fiducial PP and POWERLAW model for the mass and redshift distributions, respectively. We reweight the current injections to this same mass and redshift distribution, but reweight them to the median inferred spin distribution in [11], to mimic a astrophysically-realistic population. Both injection sets only cover the observing times of the O3a observing run. Taking the ratio of the corresponding sensitivities, we find our sensitivity has increased for low redshift events and decreased for high redshift events, relative to the sensitivity model used in [11]. We expect to see an increase in sensitivity between [11] and our current analysis due to updates to the detection pipelines. The relative decrease in sensitivity at higher redshifts indicates a bias in the previous sensitivity estimate, implying that the BBH merger rate at high redshift was underestimated in [11]. Accounting for the shift in sensitivity as a function of redshift causes a relative decrease in local BBH merger rate and a relative increase in high-redshift BBH merger rate, leading to a higher inferred value for $\kappa$.

One possible explanation for the shift in sensitivity is that the use of precession in the injections for sensitivity estimation caused a non-trivial change in the inferred sensitive hypervolume, given that we do observe precession in the BBH population. Our current detection pipelines use template banks that include only aligned-spin components; this can result in up to tens of percent reduced sensitivity to a population of BBHs with spin precession, depending on the degree of precession possible [330–332]. The farthest precessing sources, which, due to their distances, correspond to FARs closest to the detection threshold, are therefore the most susceptible to dropping below the detection threshold with our current pipelines, causing us to see a decrease in sensitivity to a population of BBHs with precession relative to a strictly non-precessing population.

Additionally, both the use of population-informed reweighting of the spin distribution to calculate sensitivity to a population and the incorporation of additional detection pipelines may have contributed to a more accurate estimate of our sensitivity across parameter space.

# Appendix D: Additional studies of the binary black hole distribution

## 1. Analyses from GWTC-2

We report updated Bayes factor comparisons for these various models in Table XV, showing that the BROKEN POWERLAW + PEAK model is slightly preferred over our fiducial POWER LAW + PEAK model. We highlight the key differences between the model priors for GWTC-2
<figure><img src="image_30.png" /><figcaption>FIG. 27. Comparison of our current BBH merger sensitivity estimate in the O3a observing run ($VT_{new}$) to that used in [11] ($VT_{old}$) as a function of redshift, for events with chirp masses between $20M_{\odot}$ and $50M_{\odot}$. Our current sensitivity model differs from what was used in [11] in two important ways: we use updated detection pipelines relative to those used in [11] and we use injections which include spin precession. Note the relative increase (decrease) in sensitivity at low (high) redshift. Computed by reweighting injections to a fiducial population for each of the two injection sets.</figcaption></figure>

<table><caption>TABLE XV. Bayes factors for each of the previously used phenomenological mass models relative to the model with highest marginal likelihood, POWER LAW + PEAK. The previous results from GWTC-2 are shown in the second column with the updated catalog results in the third column.</caption><thead><tr><th>Model</th><th>GWTC-2<br>log<sub>10</sub> B</th><th>GWTC-3<br>log<sub>10</sub> B</th></tr></thead><tbody><tr><td>POWER LAW + PEAK</td><td>0.0</td><td>0.0</td></tr><tr><td>BROKEN POWERLAW + PEAK</td><td>-0.11</td><td>-0.46</td></tr><tr><td>MULTI PEAK</td><td>-0.3</td><td>-0.22</td></tr><tr><td>BROKEN POWERLAW</td><td>-0.92</td><td>-2.0</td></tr></tbody></table>

compared to GWTC-3: the prior on $\beta_q$ was changed from $U(-4, 12)$ to $U(-2, 7)$, and the main population results now include an evolving redshift model where the prior on $\kappa$ is changed from 0 to $U(-10, 10)$.

In addition to these analyses, we used a variation of the MULTIPEAK model to study the feature in the mass distribution at $\sim 10M_{\odot}$. In GWTC-2 the prior on the mean of the peaks were $U(20, 50)$ and $U(50, 100)$ for the lower and upper mass peaks respectively. We modified these priors to be $U(5, 20)$ and $U(20, 100)$. This updated MULTIPEAK model is most preferred model with a $log_{10} B$ of 1.0 compared to the POWER LAW + PEAK model. This further supports our findings of the peak-like feature at $\sim 10M_{\odot}$ in the mass distribution.

## 2. Comprehensive BBH merger rates

In Table IV, we evaluate BBH merger rates over targeted mass subsets of the whole BBH space, using models specifically targeted to reproduce new features of the binary black hole mass distribution. For broader context, Table XVI also provides the corresponding merger rates in these intervals from all the models presented in this work.

## Appendix E: Population-weighted posteriors

With an increasing number of events, we can use the distribution of the population of compact binaries to inform our priors for parameter estimation. By reweighting the initial analysis of compact binaries with the population distribution we can obtain posterior distribution for the events in GWTC-3 with population-informed priors. Using our population analysis with models POWER LAW + PEAK and FLEXIBLE MIXTURES we provide population-weighted posteriors (Fig. 28) for $m_1$, $q$ and $\chi_{eff}$ for the BBHs population (69 events).

Some of our analyses will show apparent changes in the inferences about the mass ratio. These seemingly-substantive changes reflect the relatively weak constraints provided by the fiducial parameter inferences used as input and shown in black. Specifically, several low-amplitude or low-mass events have extremely weak constraints on mass ratios, with posterior support extending to $q < 0.4$. This extended feature reflects the prior distribution on component masses, conditioned on modest constraints on chirp mass. To be concrete, using the corresponding prior distributions for these events, conditioned on a suitable chirp mass interval, we often find a posterior distribution with comparable support for $q = 1$ (i.e., the Savage-Dickey estimate of the Bayes factor for unequal mass would be nearly unity).

Examining Fig. 28 in light of this caveat, we find that our population models and the fiducial model agree. For these events the population reweightings, as expected, strongly favor symmetric component masses (e.g., GW190503_185404, GW190720_000836, GW191127_050227). For a few binaries, however, the two population reweightings disagree. The most notable example is GW190513_205428, where the FLEXIBLE MIXTURES model pulls the posterior distribution to more symmetric component masses and a lower primary mass. Both population models also pull the majority of the posteriors closer to $\chi_{eff} \sim 0$. However, given the FLEXIBLE MIXTURES models spins as dependent on chirp mass the events with higher mass and higher spin do not drawn to $\chi_{eff} \sim 0$ as strongly as the POWER LAW + PEAK model (e.g., GW191109_010717) and in some cases the FLEXIBLE MIXTURES reweighted posterior move to higher $\chi_{eff}$ values (e.g., GW190706_222641).
<figure><img src="image_31.png" /></figure>
<figure><img src="image_32.png" /></figure>
<table><caption>TABLE XVI. Merger rates in Gpc⁻³ yr⁻¹ for black hole binaries, quoted at the 90% credible interval. Rates are given for three ranges of primary mass, $m_1$ as well as for the entire population. The PDB, MS, and BGP merger rates are derived assuming the merger rate does not increase with redshift, using a threshold FAR< 0.25 yr⁻¹ (Sec. IV). For FM, PS, and PP, merger rates are reported at $z = 0.2$, estimated using a threshold FAR< 1 yr⁻¹ (Sec. VI). The merged rates reported in the MERGED row are the union of the preceding three rows, which all account for distance-dependent merger rate and adopt a consistent threshold. The final row shows merger rates deduced from our analysis of GWTC-2 [11], which assumed a redshift-independent merger rate. Compare to Table IV.</caption><thead><tr><th></th><th>m<sub>1</sub> &isin; [5, 20]M<sub>&odot;</sub></th><th>m<sub>1</sub> &isin; [20, 50]M<sub>&odot;</sub></th><th>m<sub>1</sub> &isin; [50, 100]M<sub>&odot;</sub></th><th>All BBH</th></tr><tr><th></th><th>m<sub>2</sub> &isin; [5, 20]M<sub>&odot;</sub></th><th>m<sub>2</sub> &isin; [5, 50]M<sub>&odot;</sub></th><th>m<sub>2</sub> &isin; [5, 100]M<sub>&odot;</sub></th><th></th></tr></thead><tbody><tr><td>PDB (pair)</td><td>17<sup>+10</sup><sub>-6.0</sub></td><td>6.8<sup>+2.2</sup><sub>-1.7</sub></td><td>0.68<sup>+0.42</sup><sub>-0.29</sub></td><td>25<sup>+10</sup><sub>-7.0</sub></td></tr><tr><td>PDB (ind)</td><td>9.4<sup>+5.6</sup><sub>-3.7</sub></td><td>11<sup>+3.0</sup><sub>-2.0</sub></td><td>1.6<sup>+0.9</sup><sub>-0.7</sub></td><td>22<sup>+8.0</sup><sub>-6.0</sub></td></tr><tr><td>MS</td><td>30<sup>+23</sup><sub>-13</sub></td><td>6.6<sup>+2.9</sup><sub>-2.3</sub></td><td>0.73<sup>+0.87</sup><sub>-0.52</sub></td><td>37<sup>+24</sup><sub>-13</sub></td></tr><tr><td>BGP</td><td>20.0<sup>+11.0</sup><sub>-8.0</sub></td><td>6.3<sup>+3.0</sup><sub>-2.2</sub></td><td>0.75<sup>+1.1</sup><sub>-0.46</sub></td><td>33.0<sup>+16.0</sup><sub>-10.0</sub></td></tr><tr><td>PS</td><td>27<sup>+12</sup><sub>-8.8</sub></td><td>3.5<sup>+1.5</sup><sub>-1.1</sub></td><td>0.19<sup>+0.16</sup><sub>-0.09</sub></td><td>31<sup>+13</sup><sub>-9.2</sub></td></tr><tr><td>FM</td><td>21.1<sup>+11.6</sup><sub>-7.8</sub></td><td>4.3<sup>+2.0</sup><sub>-1.4</sub></td><td>0.2<sup>+0.2</sup><sub>-0.1</sub></td><td>26.5<sup>+11.7</sup><sub>-8.6</sub></td></tr><tr><td>PP</td><td>23.6<sup>+13.7</sup><sub>-9.0</sub></td><td>4.5<sup>+1.7</sup><sub>-1.3</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>28.3<sup>+13.9</sup><sub>-9.1</sub></td></tr><tr><td>MERGED</td><td>13.3 - 39</td><td>2.5 - 6.3</td><td>0.099 - 0.4</td><td>17.9 - 44</td></tr><tr><td>PP (O3a)</td><td>16.0<sup>+13.0</sup><sub>-7.7</sub></td><td>6.8<sup>+2.7</sup><sub>-1.9</sub></td><td>0.5<sup>+0.4</sup><sub>-0.3</sub></td><td>25.3<sup>+16.1</sup><sub>-9.9</sub></td></tr></tbody></table>

<figure><img src="image_33.png" /><figcaption>FIG. 28. Posterior distributions (black) for binary black hole events weighted by the population results from POWER LAW + PEAK (blue) and FLEXIBLE MIXTURES (orange)</figcaption></figure>

[1] R. Abbott, T. Abbott, S. Abraham, et al., GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Half of the Third Observing Run, Available as LIGO-P2000318 (2021).

[2] B. Abbott, R. Abbott, T. Abbott, et al., GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019).

[3] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, V. B. Adya, and et al., GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Available as arxiv:2108.01045, arXiv:2108.01045 (2021), arXiv:2108.01045 [gr-qc].

[4] R. Abbott, T. Abbott, S. Abraham, et al., GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run, Physical Review X 11, 021053 (2021), arXiv:2010.14527 [gr-qc].

[5] J. Aasi et al., Advanced LIGO, Class. Quant. Grav. 32,
074001 (2015).

[6] F. Acernese et al., Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. **32**, 024001 (2015).

[7] R. Essick and P. Landry, Discriminating between Neutron Stars and Black Holes with Imperfect Knowledge of the Maximum Neutron Star Mass, Astrophys. J. **904**, 80 (2020), arXiv:2007.01372 [astro-ph.HE].

[8] P. Landry, R. Essick, and K. Chatziiioannou, Non-parametric constraints on neutron star matter with existing and upcoming gravitational wave and pulsar observations, Phys. Rev. D **101**, 123007 (2020), arXiv:2003.04880 [astro-ph.HE].

[9] I. Legred, K. Chatziiioannou, R. Essick, S. Han, and P. Landry, Impact of the PSR J 0740 +6620 radius constraint on the properties of high-density matter, Phys. Rev. D **104**, 063003 (2021), arXiv:2106.05313 [astro-ph.HE].

[10] R. Abbott, T. Abbott, S. Abraham, et al., Observation of Gravitational Waves from Two Neutron Star-Black Hole Coalescences, Astrophys. J. Lett. **915**, L5 (2021), arXiv:2106.15163 [astro-ph.HE].

[11] R. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog, Astrophys. J. Lett. **913**, L7 (2021), arXiv:2010.14533 [astro-ph.HE].

[12] W. M. Farr, J. R. Gair, I. Mandel, and C. Cutler, Counting and confusion: Bayesian rate estimation with multiple populations, Phys. Rev. D **91**, 023005 (2015).

[13] R. Abbott, T. Abbott, S. Abraham, et al., Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. **116**, 061102 (2016), arXiv:1602.03837 [gr-qc].

[14] B. Abbott, R. Abbott, T. Abbott, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Binary black hole mergers in the first advanced ligo observing run, Phys. Rev. X **6**, 041015 (2016).

[15] B. P. Abbott et al. (LIGO Scientific, Virgo), GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence, Phys. Rev. Lett. **116**, 241103 (2016), arXiv:1606.04855 [gr-qc].

[16] B. Abbott, R. Abbott, T. Abbott, et al., GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2, Phys. Rev. Lett. **118**, 221101 (2017), arXiv:1706.01812 [gr-qc].

[17] B. Abbott, R. Abbott, T. Abbott, et al., GW170608: Observation of a 19 Solar-mass Binary Black Hole Coalescence, Astrophys. J. Lett. **851**, L35 (2017), arXiv:1711.05578 [astro-ph.HE].

[18] B. Abbott, R. Abbott, T. Abbott, et al., GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence, Phys. Rev. Lett. **119**, 141101 (2017), arXiv:1709.09660 [gr-qc].

[19] B. Abbott, R. Abbott, T. Abbott, et al. (LIGO Scientific, Virgo), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. **119**, 161101 (2017).

[20] R. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), GW190412: Observation of a binary-black-hole coalescence with asymmetric masses, Phys. Rev. D **102**, 043015 (2020).

[21] B. Abbott, R. Abbott, T. Abbott, et al., GW190425: Observation of a Compact Binary Coalescence with Total Mass ~3.4 M⊙, Astrophys. J. Lett. **892**, L3 (2020).

[22] R. Abbott et al. (LIGO Scientific, Virgo), GW190521: A Binary Black Hole Merger with a Total Mass of 150 M⊙, Phys. Rev. Lett. **125**, 101102 (2020), arXiv:2009.01075 [gr-qc].

[23] R. Abbott, T. Abbott, S. Abraham, et al., GW190814: Gravitational waves from the coalescence of a 23 solar mass black hole with a 2.6 solar mass compact object, Astrophys. J. Lett. **896**, L44 (2020).

[24] A. H. Nitz, C. D. Capano, S. Kumar, Y.-F. Wang, S. Kastha, M. Schäfer, R. Dhurkunde, and M. Cabero, 3-OGC: Catalog of gravitational waves from compact-binary mergers, Available as arxiv:2105.09151 (2021), arXiv:2105.09151 [astro-ph.HE].

[25] B. Abbott, R. Abbott, T. Abbott, et al., Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo, Astrophys. J. Lett. **882**, L24 (2019).

[26] P. Ajith et al., Inspiral-merger-ringdown waveforms for black-hole binaries with non-precessing spins, Phys. Rev. Lett. **106**, 241101 (2011), arXiv:0909.2867 [gr-qc].

[27] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, and et al., Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. **594**, A13 (2016), arXiv:1502.01589 [astro-ph.CO].

[28] C. D. Bailyn, R. K. Jain, P. Coppi, and J. A. Orosz, The Mass distribution of stellar black holes, Astrophys. J. **499**, 367 (1998), arXiv:astro-ph/9708032.

[29] W. M. Farr, N. Sravan, A. Cantrell, L. Kreidberg, C. D. Bailyn, I. Mandel, and V. Kalogera, The Mass Distribution of Stellar-Mass Black Holes, Astrophys. J. **741**, 103 (2011), arXiv:1011.1459 [astro-ph.GA].

[30] L. Kreidberg, C. D. Bailyn, W. M. Farr, and V. Kalogera, Mass Measurements of Black Holes in X-Ray Transients: Is There a Mass Gap?, Astrophys. J. **757**, 36 (2012), arXiv:1205.1805 [astro-ph.HE].

[31] F. Ozel, D. Psaltis, R. Narayan, and J. E. McClintock, The Black Hole Mass Distribution in the Galaxy, Astrophys. J. **725**, 1918 (2010), arXiv:1006.2834 [astro-ph.GA].

[32] C. L. Fryer, K. Belczynski, G. Wiktorowicz, M. Dominik, V. Kalogera, and D. E. Holz, Compact Remnant Mass Function: Dependence on the Explosion Mechanism and Metallicity, Astrophys. J. **749**, 91 (2012), arXiv:1110.1726 [astro-ph.SR].

[33] T. Liu, Y.-F. Wei, L. Xue, and M.-Y. Sun, Final Compact Remnants in Core-collapse Supernovae from 20 to 40 M⊙: The Lower Mass Gap, Astrophys. J. **908**, 106 (2021), arXiv:2011.14361 [astro-ph.HE].

[34] I. Mandel and B. Müller, Simple recipes for compact remnant masses and natal kicks, Mon. Not. Roy. Astron. Soc. **499**, 3214 (2020), arXiv:2006.08360 [astro-ph.HE].

[35] R. A. Patton, T. Sukhbold, and J. J. Eldridge, Comparing Compact Object Distributions from Mass- and Pre-supernova Core Structure-based Prescriptions, Available as arXiv:2106.05987 (2021), arXiv:2106.05978 [astro-ph.HE].

[36] M. Zevin, M. Spera, C. P. L. Berry, and V. Kalogera, Exploring the Lower Mass Gap and Unequal Mass Regime in Compact Binary Evolution, Astrophys. J. Lett. **899**, L1 (2020), arXiv:2006.14573 [astro-ph.HE].

[37] T. Jayasinghe, K. Z. Stanek, T. A. Thompson, C. S.
Kochanek, D. M. Rowan, P. J. Vallely, K. G. Strassmeier, M. Weber, J. T. Hinkle, F. J. Hambsch, D. V. Martin, J. L. Prieto, T. Pessi, D. Huber, K. Auchettl, L. A. Lopez, I. Ilyin, C. Badenes, A. W. Howard, H. Isaacson, and S. J. Murphy, A unicorn in monoceros: the 3 M☉ dark companion to the bright, nearby red giant V723 Mon is a non-interacting, mass-gap black hole candidate, Mon. Not. R. Ast. Soc. 504, 2577 (2021), arXiv:2101.02212 [astro-ph.SR].

[38] T. A. Thompson, C. S. Kochanek, K. Z. Stanek, C. Badenes, R. S. Post, T. Jayasinghe, D. W. Latham, A. Bieryla, G. A. Esquerdo, P. Berlind, M. L. Calkins, J. Tayar, L. Lindegren, J. A. Johnson, T. W. S. Holoien, K. Auchettl, and K. Covey, A noninteracting low-mass black hole-giant star binary system, Science 366, 637 (2019), arXiv:1806.02751 [astro-ph.HE].

[39] N. Farrow, X.-J. Zhu, and E. Thrane, The mass distribution of galactic double neutron stars, Astrophys. J. 876, 18 (2019).

[40] B. Kiziltan, A. Kottas, M. De Yoreo, and S. E. Thorsett, The Neutron Star Mass Distribution, arXiv e-prints, arXiv:1309.6635 (2013), arXiv:1309.6635 [astro-ph.SR].

[41] F. Özel and P. Freire, Masses, Radii, and the Equation of State of Neutron Stars, Ann. Rev. Astron. Astrophys. 54, 401 (2016), arXiv:1603.02698 [astro-ph.HE].

[42] T. Menne, awkde code, https://github.com/mennthor/awkde (2020), [Online; accessed Dec-2020].

[43] B. Wang and X. Wang, Bandwidth selection for weighted kernel density estimation (2011), arXiv:0709.1616 [stat.ME].

[44] J. Sadiq, T. Dent, and D. Wysocki, Flexible and Fast Estimation of Binary Merger Population Distributions with Adaptive KDE, Available as LIGO-T2100447 (2021), https://dcc.ligo.org/LIGO-T2100447/public.

[45] B. Edelman, Z. Doctor, J. Godfrey, and B. Farr, Ain't No Mountain High Enough: Semi-Parametric Modeling of LIGO-Virgos Binary Black Hole Mass Distribution, Available as arxiv:2109.06137, arXiv:2109.06137 (2021), arXiv:2109.06137 [astro-ph.HE].

[46] Y.-J. Li, Y.-Z. Wang, M.-Z. Han, S.-P. Tang, Q. Yuan, Y.-Z. Fan, and D.-M. Wei, A Flexible Gaussian Process Reconstruction Method and the Mass Function of the Coalescing Binary Black Hole Systems, Astrophys. J. 917, 33 (2021), arXiv:2104.02969 [astro-ph.HE].

[47] V. Tiwari and S. Fairhurst, The Emergence of Structure in the Binary Black Hole Mass Distribution, Astrophys. J. Lett. 913, L19 (2021), arXiv:2011.04502 [astro-ph.HE].

[48] D. Veske, I. Bartos, Z. Márka, and S. Márka, Characterizing the observation bias in gravitational-wave detections and finding structured population properties, Available as arxiv:2105.13983, arXiv:2105.13983 (2021), arXiv:2105.13983 [gr-qc].

[49] V. Tiwari, VAMANA: modeling binary black hole population with minimal assumptions, Class. Quant. Grav. 38, 155007 (2021), arXiv:2006.15047 [astro-ph.HE].

[50] C. Hanna et al., Fast evaluation of multi-detector consistency for real-time gravitational wave searches, Phys. Rev. D 101, 022003 (2020).

[51] C. Messick et al., Analysis Framework for the Prompt Discovery of Compact Binary Mergers in Gravitational-wave Data, Phys. Rev. D 95, 042001 (2017).

[52] S. Sachdev et al., The GstLAL Search Analysis Methods for Compact Binary Mergers in Advanced LIGO's Second and Advanced Virgo's First Observing Runs, arXiv eprints (2019), arXiv:1901.08580 [gr-qc].

[53] B. Allen, $\chi^2$ time-frequency discriminator for gravitational wave detection, Phys. Rev. D 71, 062001 (2005).

[54] B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E. Creighton, FIND-CHIRP: An Algorithm for detection of gravitational waves from inspiraling compact binaries, Phys. Rev. D 85, "122006" (2012).

[55] T. D. Canton et al., Implementing a search for aligned-spin neutron star-black hole systems with advanced ground based gravitational wave detectors, Phys. Rev. D 90, 082004 (2014).

[56] A. Nitz et al., (2019), gwastro/pycbc: Pycbc release v1.15.0.

[57] A. H. Nitz, T. Dent, T. D. Canton, S. Fairhurst, and D. A. Brown, Detecting binary compact-object mergers with gravitational waves: Understanding and Improving the sensitivity of the PyCBC search, Astrophys. J. 849, 118 (2017).

[58] S. A. Usman et al., The PyCBC search for gravitational waves from compact binary coalescence, Class. Quant. Grav. 33, 215004 (2016).

[59] T. Adams, D. Buskulic, V. Germain, G. M. Guidi, F. Marion, M. Montani, B. Mours, F. Piergiovanni, and G. Wang, Low-latency analysis pipeline for compact binary coalescences in the advanced gravitational wave detector era, Class. Quant. Grav. 33, 175012 (2016).

[60] S. Klimenko et al., Method for detection and reconstruction of gravitational wave transients with networks of advanced detectors, Phys. Rev. D 93, 042004 (2016).

[61] S. Klimenko and G. Mitselmakher, A wavelet method for detection of gravitational wave bursts, 8th Gravitational Wave Data Analysis Workshop (GWDAW 2003) Milwaukee, Wisconsin, December 17-20, 2003, Class. Quant. Grav. 21, S1819 (2004).

[62] S. M. Gaebel, J. Veitch, T. Dent, and W. M. Farr, Digging the population of compact binary mergers out of the noise, Mon. Not. R. Ast. Soc. 484, 4008 (2019), arXiv:1809.03815 [astro-ph.IM].

[63] S. Galaudage, C. Talbot, and E. Thrane, Gravitational-wave inference in the catalog era: Evolving priors and marginal events, Phys. Rev. D 102, 083026 (2020), arXiv:1912.09708 [astro-ph.HE].

[64] J. Roulet, T. Venumadhav, B. Zackay, L. Dai, and M. Zaldarriaga, Binary black hole mergers from LIGO/Virgo O1 and O2: Population inference combining confident and marginal events, Phys. Rev. D 102, 123022 (2020), arXiv:2008.07014 [astro-ph.HE].

[65] R. Magee et al., Sub-threshold Binary Neutron Star Search in Advanced LIGO's First Observing Run, Astrophys. J. 878, L17 (2019).

[66] A. H. Nitz, C. Capano, A. B. Nielsen, S. Reyes, R. White, D. A. Brown, and B. Krishnan, 1-OGC: The first open gravitational-wave catalog of binary mergers from analysis of public Advanced LIGO data, Astrophys. J. 872, 195 (2019).

[67] A. H. Nitz, T. Dent, G. S. Davies, S. Kumar, C. D. Capano, I. Harry, S. Mozzon, L. Nuttall, A. Lundgren, and M. Tápai, 2-OGC: Open Gravitational-wave Catalog of binary mergers from analysis of public Advanced LIGO and Virgo data, Astrophys. J. 891, 123 (2020).

[68] T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and
M. Zaldarriaga, New search pipeline for compact binary mergers: Results for binary black holes in the first observing run of Advanced LIGO, Phys. Rev. D **100**, 023011 (2019).

[69] T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, New Binary Black Hole Mergers in the Second Observing Run of Advanced LIGO and Advanced Virgo, Phys. Rev. D **101**, 083030 (2020).

[70] B. Zackay, L. Dai, T. Venumadhav, J. Roulet, and M. Zaldarriaga, Detecting Gravitational Waves With Disparate Detector Responses: Two New Binary Black Hole Mergers, Phys. Rev. D **101**, 083030 (2020).

[71] B. Zackay, T. Venumadhav, L. Dai, J. Roulet, and M. Zaldarriaga, Highly spinning and aligned binary black hole merger in the Advanced LIGO first observing run, Phys. Rev. D **100**, 023007 (2019).

[72] R. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo, SoftwareX **13**, 100658 (2021), arXiv:1912.11716 [gr-qc].

[73] A. Trovato et al. (Ligo Scientific, The Virgo), GWOSC: Gravitational Wave Open Science Center, Proceedings, The New Era of Multi-Messenger Astrophysics (ASTERICS 2019): Groningen, Netherlands, March 25-29, 2019, PoS Asterics2019, 082 (2020).

[74] R. Abbott, T. Abbott, S. Abraham, et al., Search for intermediate mass black hole binaries in the third observing run of Advanced LIGO and Advanced Virgo, Available as arxiv:2105.15120 (2021).

[75] G. Ashton and E. Thrane, The astrophysical odds of GW151216, Mon. Not. Roy. Astron. Soc. **498**, 1905 (2020), arXiv:2006.05039 [astro-ph.HE].

[76] G. Pratten and A. Vecchio, Assessing gravitational-wave binary black hole candidates with Bayesian odds, arXiv eprints (2020), arXiv:2008.00509 [gr-qc].

[77] J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, W. Del Pozzo, F. Feroz, J. Gair, C. J. Haster, V. Kalogera, T. Littenberg, I. Mandel, R. O'Shaughnessy, M. Pitkin, C. Rodriguez, C. Röver, T. Sidery, R. Smith, M. Van Der Sluys, A. Vecchio, W. Vousden, and L. Wade, Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library, Phys. Rev. D **91**, 042003 (2015).

[78] J. Lange, R. O'Shaughnessy, and M. Rizzo, Rapid and accurate parameter inference for coalescing precessing compact binaries, (2018), arxiv/1805.10457.

[79] D. Wysocki, R. O'Shaughnessy, J. Lange, and Y.-L. L. Fang, Accelerating parameter inference with graphics processing units, Phys. Rev. D **99** (2019).

[80] G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter, B. Goncharov, F. Hernandez Vivanco, J. Harms, M. E. Lower, G. D. Meadors, D. Melchor, E. Payne, M. D. Pitkin, J. Powell, N. Sarin, R. J. E. Smith, and E. Thrane, BILBY: A User-friendly Bayesian Inference Library for Gravitational-wave Astronomy, Astrophys. J. Supp. **241**, 27 (2019).

[81] I. M. Romero-Shaw, C. Talbot, S. Biscoveanu, V. D'Emilio, G. Ashton, C. P. L. Berry, S. Coughlin, S. Galaudage, C. Hoy, M. Hübner, K. S. Phukon, M. Pitkin, M. Rizzo, N. Sarin, R. Smith, S. Stevenson, A. Vajpeyi, M. Arène, K. Athar, S. Banagiri,

N. Bose, M. Carney, K. Chatziioannou, J. A. Clark, M. Colleoni, R. Cotesta, B. Edelman, H. Estellés, C. García-Quirós, A. Ghosh, R. Green, C. J. Haster, S. Husa, D. Keitel, A. X. Kim, F. Hernandez-Vivanco, I. Magaña Hernandez, C. Karathanasis, P. D. Lasky, N. De Lillo, M. E. Lower, D. Macleod, M. Mateu-Lucena, A. Miller, M. Millhouse, S. Morisaki, S. H. Oh, S. Ossokine, E. Payne, J. Powell, G. Pratten, M. Pürrer, A. Ramos-Buades, V. Raymond, E. Thrane, J. Veitch, D. Williams, M. J. Williams, and L. Xiao, Bayesian inference for compact binary coalesces with BILBY: validation and application to the first LIGO-Virgo gravitational-wave transient catalogue, Mon. Not. R. Ast. Soc. **499**, 3295 (2020), arXiv:2006.00714 [astro-ph.IM].

[82] B. Abbott, R. Abbott, T. Abbott, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Properties of the Binary Black Hole Merger GW150914, Phys. Rev. Lett. **116**, 241102 (2016).

[83] Y. Pan, A. Buonanno, A. Taracchini, L. E. Kidder, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Inspiral-merger-ringdown waveforms of spinning, precessing black-hole binaries in the effective-one-body formalism, Phys. Rev. D **89**, 084006 (2014), arXiv:1307.6232 [gr-qc].

[84] A. Taracchini, A. Buonanno, Y. Pan, T. Hinderer, M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, B. Szilágyi, N. W. Taylor, and A. Zenginoglu, Effective-one-body model for black-hole binaries with generic mass ratios and spins, Phys. Rev. D **89**, 061502 (2014), arXiv:1311.2544 [gr-qc].

[85] M. Hannam et al., Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms, Phys. Rev. Lett. **113**, 151101 (2014).

[86] A. Bohé et al., Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors, Phys. Rev. D **95**, 044028 (2017).

[87] S. Ossokine, A. Buonanno, S. Marsat, R. Cotesta, S. Babak, T. Dietrich, R. Haas, I. Hinder, H. P. Pfeiffer, M. Pürrer, C. J. Woodford, M. Boyle, L. E. Kidder, M. A. Scheel, and B. Szilágyi, Multipolar effective-one-body waveforms for precessing binary black holes: Construction and validation, Phys. Rev. D **102**, 044055 (2020), arXiv:2004.09442 [gr-qc].

[88] G. Pratten, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, M. Mateu-Lucena, R. Jaume, M. Haney, D. Keitel, J. E. Thompson, and S. Husa, Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes, Phys. Rev. D **103**, 104056 (2021), arXiv:2004.06503 [gr-qc].

[89] D. Estevez, P. Lagabbe, A. Masserot, L. Rolland, M. Seglar-Arroyo, and D. Verkindt, The Advanced Virgo photon calibrators, Class. Quant. Grav. **38**, 075007 (2021).

[90] L. Sun, E. Goetz, J. S. Kissel, J. Betzwieser, S. Karki, D. Bhattacharjee, P. B. Covas, L. E. H. Datrier, S. Kandhasamy, Y. K. Lecoeuche, G. Mendell, T. Mystery, E. Payne, R. L. Savage, A. Viets, M. Wade, A. J. Weinstein, S. Aston, C. Cahillane, J. C. Driggers, S. E. Dwyer, and A. Urban, Characterization of systematic error in Advanced LIGO calibration in the second half
of O3, Available as arxiv:2107.00129, arXiv:2107.00129 (2021), arXiv:2107.00129 [astro-ph.IM].

[91] L. Sun, E. Goetz, J. S. Kissel, J. Betzwieser, S. Karki, A. Viets, M. Wade, D. Bhattacharjee, V. Bossilkov, P. B. Covas, and et al., Characterization of systematic error in Advanced LIGO calibration, Class. Quant. Grav. 37, 225008 (2020), arXiv:2005.02531 [astro-ph.IM].

[92] I. Mandel, W. M. Farr, and J. R. Gair, Extracting distribution parameters from multiple uncertain observations with selection biases, Mon. Not. R. Ast. Soc. 486, 1086 (2019).

[93] E. Thrane and C. Talbot, An introduction to Bayesian inference in gravitational-wave astronomy: parameter estimation, model selection, and hierarchical models, PASA 36, E010 (2019).

[94] S. Vitale, One, No One, and One Hundred Thousand – Inferring the properties of a population in presence of selection effects, Available as arxiv:2007.05579 (2020), arXiv:2007.05579 [astro-ph.IM].

[95] T. J. Loredo, Accounting for Source Uncertainties in Analyses of Astronomical Survey Data, in American Institute of Physics Conference Series, Vol. 735 (2004) pp. 195–206.

[96] M. Fishbach, D. E. Holz, and W. M. Farr, Does the Black Hole Merger Rate Evolve with Redshift?, Astrophys. J. Lett. 863, L41 (2018), arXiv:1805.10270 [astro-ph.HE].

[97] C. Talbot, R. Smith, E. Thrane, and G. B. Poole, Parallelized inference for gravitational-wave astronomy, Phys. Rev. D 100, 043030 (2019).

[98] C. Talbot, Gwpopulation pipe, Zenodo, available as https://zenodo.org/record/5654673 10.5281/zenodo.5654673 (2021).

[99] D. Wysocki and R. O'Shaughnessy, Bayesian parametric population models (2017), bayesian-parametric-population-models.readthedocs.io.

[100] P. Landry and J. S. Read, The mass distribution of neutron stars in gravitational-wave binaries, Available as arxiv:2107.04559, arXiv:2107.04559 (2021), arXiv:2107.04559 [astro-ph.HE].

[101] B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell, Stan: A probabilistic programming language, Journal of Statistical Software 76, 10.18637/jss.v076.i01 (2017).

[102] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, Publ. Astron. Soc. Pac 125, 306 (2013).

[103] A. Riddell, A. Hartikainen, D. Lee, Riddell-Stan, et al., Stan-dev/pystan: V2.18.0.0 (2018).

[104] J. S. Speagle, dynasty: A Dynamic Nested Sampling Package for Estimating Bayesian Posteriors and Evidences, Mon. Not. R. Ast. Soc. 493, 3132 (2020).

[105] W. M. Farr, Accuracy Requirements for Empirically Measured Selection Functions, Research Notes of the American Astronomical Society 3, 66 (2019).

[106] J. Golomb and C. Talbot, Hierarchical Inference of Binary Neutron Star Mass Distribution and Equation of State with Gravitational Waves, arXiv e-prints, arXiv:2106.15745 (2021), arXiv:2106.15745 [astro-ph.HE].

[107] C. Talbot and E. Thrane, Measuring the binary black hole mass spectrum with an astrophysically motivated

parameterization, Astrophys. J. 856, 173 (2018).

[108] M. Fishbach, R. Essick, and D. E. Holz, Does Matter Matter? Using the Mass Distribution to Distinguish Neutron Stars and Black Holes, Astrophys. J. Lett. 899, L8 (2020), arXiv:2006.13178 [astro-ph.HE].

[109] A. M. Farah, M. Fishbach, R. Essick, D. E. Holz, and S. Galaudage, Bridging the gap: Categorizing gravitational-wave events at the transition between neutron stars and black holes, Available as arxiv:2111.03498 (2021), arXiv:2111.03498 [astro-ph.HE].

[110] Z. Doctor, D. Wysocki, R. O'Shaughnessy, D. E. Holz, and B. Farr, Black hole coagulation: Modeling hierarchical mergers in black hole populations, Astrophys. J. 893, 35 (2019).

[111] M. Fishbach and D. E. Holz, Picky Partners: The Pairing of Component Masses in Binary Black Hole Mergers, Astrophys. J. Lett. 891, L27 (2020), arXiv:1905.12669 [astro-ph.HE].

[112] B. Farr, E. Ochsner, W. M. Farr, and R. O'Shaughnessy, A more effective coordinate system for parameter estimation of precessing compact binaries from gravitational waves, Phys. Rev. D 90, 024018 (2014), arXiv:1404.7070 [gr-qc].

[113] D. Wysocki, J. Lange, and R. O'Shaughnessy, Reconstructing phenomenological distributions of compact binaries via gravitational wave observations, Phys. Rev. D 100, 043012 (2019), arXiv:1805.06442 [gr-qc].

[114] C. Talbot and E. Thrane, Determining the population properties of spinning black holes, Phys. Rev. D 96, 023012 (2017).

[115] S. Miller, T. A. Callister, and W. Farr, The Low Effective Spin of Binary Black Holes and Implications for Individual Gravitational-Wave Events, Astrophys. J. 895, 128 (2020).

[116] J. Roulet and M. Zaldarriaga, Constraints on binary black hole populations from LIGO-Virgo detections, Mon. Not. R. Ast. Soc. 484, 4216 (2019), arXiv:1806.10610 [astro-ph.HE].

[117] M. Fishbach and D. E. Holz, Where Are LIGO's Big Black Holes?, Astrophys. J. Lett. 851, L25 (2017).

[118] D. Foreman-Mackey, D. W. Hogg, and T. D. Morton, Exoplanet population inference and the abundance of earth analogs from noisy, incomplete catalogs, The Astrophysical Journal 795, 64 (2014).

[119] I. Mandel, W. M. Farr, A. Colonna, S. Stevenson, P. Tiño, and J. Veitch, Model-independent inference on compact-binary observations, Mon. Not. R. Ast. Soc. 465, 3254 (2017), arXiv:1608.08223 [astro-ph.HE].

[120] J. Salvatier, T. V. Wiecki, and C. Fonnesbeck, Probabilistic programming in python using PyMC3, PeerJ Computer Science 2, e55 (2016).

[121] S. J. Kapadia et al., A self-consistent method to estimate the rate of compact binary coalescences with a Poisson mixture model, Class. Quant. Grav. 37, 045007 (2020).

[122] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Multimesenger constraints on the neutron-star equation of state and the Hubble constant, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE].

[123] J.-L. Jiang, S.-P. Tang, Y.-Z. Wang, Y.-Z. Fan, and D.-M. Wei, PSR J0030+0451, GW170817, and the Nuclear Data: Joint Constraints on Equation of State and Bulk Properties of Neutron Stars, Astrophys. J. 892, 55
(2020), arXiv:1912.07467 [astro-ph.HE].

[124] Y. Lim and J. W. Holt, Bayesian modeling of the nuclear equation of state for neutron star tidal deformabilities and GW170817, *European Physical Journal A* **55**, 209 (2019), arXiv:1902.05502 [nucl-th].

[125] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, V. B. Adya, C. Affeldt, M. Agathos, and et al., Model comparison from LIGO-Virgo data on GW170817's binary components and consequences for the merger remnant, *Class. Quant. Grav.* **37**, 045006 (2020), arXiv:1908.01012 [gr-qc].

[126] B. Margalit and B. D. Metzger, Constraining the Maximum Mass of Neutron Stars from Multi-messenger Observations of GW170817, Astrophys. J. Lett. **850**, L19 (2017), arXiv:1710.05938 [astro-ph.HE].

[127] A. Nathanail, E. R. Most, and L. Rezzolla, GW170817 and GW190814: Tension on the Maximum Mass, Astrophys. J. Lett. **908**, L28 (2021), arXiv:2101.01735 [astro-ph.HE].

[128] L. Rezzolla, E. R. Most, and L. R. Weih, Using Gravitational-wave Observations and Quasi-universal Relations to Constrain the Maximum Mass of Neutron Stars, Astrophys. J. Lett. **852**, L25 (2018), arXiv:1711.00314 [astro-ph.HE].

[129] M. Ruiz, S. L. Shapiro, and A. Tsokaros, GW170817, general relativistic magnetohydrodynamic simulations, and the neutron star maximum mass, Phys. Rev. D **97**, 021501 (2018), arXiv:1711.00473 [astro-ph.HE].

[130] M. Shibata, E. Zhou, K. Kiuchi, and S. Fujibayashi, Constraint on the maximum mass of neutron stars using GW170817 event, Phys. Rev. D **100**, 023015 (2019), arXiv:1905.03656 [astro-ph.HE].

[131] S. Galaudage, C. Adamcewicz, X.-J. Zhu, S. Stevenson, and E. Thrane, Heavy Double Neutron Stars: Birth, Midlife, and Death, Astrophys. J. Lett. **909**, L19 (2021), arXiv:2011.01495 [astro-ph.HE].

[132] É. É. Flanagan and T. Hinderer, Constraining neutron-star tidal Love numbers with gravitational-wave detectors, Phys. Rev. D **77**, 021502 (2008), arXiv:0709.1915 [astro-ph].

[133] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, and et al., A Massive Pulsar in a Compact Relativistic Binary, Science **340**, 448 (2013), arXiv:1304.6875 [astro-ph.HE].

[134] H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. De morest, Z. Arzoumanian, H. Blumer, P. R. Brook, M. E. DeCesar, T. Dolch, J. A. Ellis, and et al., Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar, Nature Astronomy **4**, 72 (2020), arXiv:1904.06759 [astro-ph.HE].

[135] I. Bombaci, The maximum mass of a neutron star, Astron. Astrophys. **305**, 871 (1996).

[136] V. Kalogera and G. Baym, The Maximum Mass of a Neutron Star, Astrophys. J. Lett. **470**, L61 (1996), arXiv:astro-ph/9608059 [astro-ph].

[137] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, and et al., GW170817: Measurements of Neutron Star Radii and Equation of State, Phys. Rev. Lett. **121**, 161101 (2018), arXiv:1805.11581 [gr-qc].

[138] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Rapidly Rotating Neutron Stars in General Relativity: Realistic Equations of State, Astrophys. J. **424**, 823 (1994).

[139] C. Breu and L. Rezzolla, Maximum mass, moment of inertia and compactness of relativistic stars, Mon. Not. R. Ast. Soc. **459**, 646 (2016), arXiv:1601.06083 [gr-qc].

[140] F. Özel, D. Psaltis, R. Narayan, and A. Santos Vilolarreal, On the Mass Distribution and Birth Masses of Neutron Stars, Astrophys. J. **757**, 55 (2012), arXiv:1201.1006 [astro-ph.HE].

[141] J. Alsing, H. O. Silva, and E. Berti, Evidence for a maximum mass cut-off in the neutron star mass distribution and constraints on the equation of state, Mon. Not. R. Ast. Soc. **478**, 1377 (2018), arXiv:1709.07889 [astro-ph.HE].

[142] W. M. Farr and K. Chatziioannou, A Population-Informed Mass Estimate for Pulsar J0740+6620, Research Notes of the American Astronomical Society **4**, 65 (2020), arXiv:2005.00032 [astro-ph.GA].

[143] D.-S. Shao, S.-P. Tang, J.-L. Jiang, and Y.-Z. Fan, Maximum mass cutoff in the neutron star mass distribution and the prospect of forming supramassive objects in the double neutron star mergers, Phys. Rev. D **102**, 063006 (2020), arXiv:2009.04275 [astro-ph.HE].

[144] K. Chatziioannou and W. M. Farr, Inferring the maximum and minimum mass of merging neutron stars with gravitational waves, Phys. Rev. D **102**, 064063 (2020), arXiv:2005.00482 [astro-ph.HE].

[145] R. Essick, A. Farah, S. Galaudage, C. Talbot, M. Fishbach, E. Thrane, and D. E. Holz, Probing Extreme Gravitational-Wave Events with Coarse-Grained Likelihoods, arXiv e-prints, arXiv:2109.00418 (2021), arXiv:2109.00418 [astro-ph.HE].

[146] E. Wilson, Probable inference, the law of succession, and statistical inference, Journal of the American Statistical Association, **209** (1927).

[147] K. Belczynski, A. Heger, W. Gladysz, A. J. Ruiter, S. Woosley, G. Wiktorowicz, H. Y. Chen, T. Buk, R. O'Shaughnessy, D. E. Holz, C. L. Fryer, and E. Berti, The effect of pair-instability mass loss on black-hole mergers, Astron. Astrophys. **594**, A97 (2016), arXiv:1607.03116 [astro-ph.HE].

[148] W. A. Fowler and F. Hoyle, Neutrino Processes and Pair Formation in Massive Stars and Supernovae, Astrophys. J. Supp. **9**, 201 (1964).

[149] P. Marchant, M. Renzo, R. Farmer, K. M. W. Papas, R. E. Taam, S. E. de Mink, and V. Kalogera, Pulsational Pair-instability Supernovae in Very Close Binaries, Astrophys. J. **882**, 36 (2019), arXiv:1810.13412 [astro-ph.HE].

[150] M. Spera and M. Mapelli, Very massive stars, pair-instability supernovae and intermediate-mass black holes with the sevn code, Mon. Not. R. Ast. Soc. **470**, 4739 (2017).

[151] S. Stevenson, M. Sampson, J. Powell, A. Vigna-Gómez, C. J. Neijssel, D. Szécsi, and I. Mandel, The Impact of Pair-instability Mass Loss on the Binary Black Hole Mass Distribution, Astrophys. J. **882**, 121 (2019), arXiv:1904.02821 [astro-ph.HE].

[152] S. E. Woosley, Pulsational pair-instability supernovae, Astrophys. J. **836**, 244 (2017).

[153] S. E. Woosley and A. Heger, The Pair-instability Mass Gap for Black Holes, Astrophys. J. Lett. **912**, L31 (2021), arXiv:2103.07933 [astro-ph.SR].
[154] R. Abbott et al. (LIGO Scientific, Virgo), Properties and Astrophysical Implications of the 150 M⊙ Binary Black Hole Merger GW190521, Astrophys. J. Lett. 900, L13 (2020), arXiv:2009.01190 [astro-ph.HE].

[155] B. Edelman, Z. Doctor, and B. Farr, Poking Holes: Looking for Gaps in LIGO/Virgo's Black Hole Population, Astrophys. J. Lett. 913, L23 (2021), arXiv:2104.07783 [astro-ph.HE].

[156] M. Fishbach and D. E. Holz, Minding the Gap: GW190521 as a Straddling Binary, Astrophys. J. Lett. 904, L26 (2020), arXiv:2009.05472 [astro-ph.HE].

[157] A. H. Nitz and C. D. Capano, GW190521 May Be an Intermediate-mass Ratio Inspiral, Astrophys. J. Lett. 907, L9 (2021), arXiv:2010.12558 [astro-ph.HE].

[158] P. Madau and M. Dickinson, Cosmic Star-Formation History, Annual Review of Astronomy and Astrophysics 52, 415 (2014).

[159] F. Antonini and M. Gieles, Merger rate of black hole binaries from globular clusters: Theoretical error bars and comparison to gravitational wave data from GWTC-2, Phys. Rev. D 102, 123016 (2020), arXiv:2009.01861 [astro-ph.HE].

[160] M. Fishbach and V. Kalogera, The Time Delay Distribution and Formation Metallicity of LIGO-Virgos Binary Black Holes, Astrophys. J. Lett. 914, L30 (2021), arXiv:2105.06491 [astro-ph.HE].

[161] M. Mapelli and N. Giacobbo, The cosmic merger rate of neutron stars and black holes, Mon. Not. R. Ast. Soc. 479, 4391 (2018), arXiv:1806.04866 [astro-ph.HE].

[162] M. A. S. Martinez, G. Fragione, K. Kremer, S. Chatterjee, C. L. Rodriguez, J. Samsing, C. S. Ye, N. C. Weatherford, M. Zevin, S. Naoz, and F. A. Rasio, Black Hole Mergers from Hierarchical Triples in Dense Star Clusters, Astrophys. J. 903, 67 (2020), arXiv:2009.08468 [astro-ph.GA].

[163] R. O'Shaughnessy, V. Kalogera, and K. Belczynski, Binary Compact Object Coalescence Rates: The Role of Elliptical Galaxies, Astrophys. J. 716, 615 (2010), arXiv:0908.3635 [astro-ph.CO].

[164] C. L. Rodriguez and A. Loeb, Redshift Evolution of the Black Hole Merger Rate from Globular Clusters, Astrophys. J. Lett. 866, L5 (2018), arXiv:1809.01152 [astro-ph.HE].

[165] F. Santoliquido, M. Mapelli, Y. Bouffanais, N. Giacobbo, U. N. Di Carlo, S. Rastello, M. C. Artale, and A. Ballone, The Cosmic Merger Rate Density Evolution of Compact Binaries Formed in Young Star Clusters and in Isolated Binaries, Astrophys. J. 898, 152 (2020), arXiv:2004.09533 [astro-ph.HE].

[166] L. A. C. van Son, S. E. de Mink, T. Callister, S. Justham, M. Renzo, T. Wagg, F. S. Broekgaarden, F. Kummer, R. Pakmor, and I. Mandel, The redshift evolution of the binary black hole merger rate: a weighty matter, Available as arxiv:2110.01634, arXiv:2110.01634 (2021), arXiv:2110.01634 [astro-ph.HE].

[167] Y. Yang, I. Bartos, Z. Haiman, B. Kocsis, S. Márka, and H. Tagawa, Cosmic Evolution of Stellar-mass Black Hole Merger Rate in Active Galactic Nuclei, Astrophys. J. 896, 138 (2020), arXiv:2003.08564 [astro-ph.HE].

[168] N. Christensen, Stochastic gravitational wave backgrounds, Reports on Progress in Physics 82, 016903 (2019), arXiv:1811.08797 [gr-qc].

[169] T. Regimbau, The astrophysical gravitational wave

stochastic background, Research in Astronomy and Astrophysics 11, 369 (2011), arXiv:1101.2762 [astro-ph.CO].

[170] J. D. Romano and N. J. Cornish, Detection methods for stochastic gravitational-wave backgrounds: a unified treatment, Living Reviews in Relativity 20, 2 (2017), arXiv:1608.06889 [gr-qc].

[171] T. Callister, M. Fishbach, D. E. Holz, and W. M. Farr, Shouts and Murmurs: Combining Individual Gravitational-wave Sources with the Stochastic Background to Measure the History of Binary Black Hole Mergers, Astrophys. J. Lett. 896, L32 (2020), arXiv:2003.12152 [astro-ph.HE].

[172] R. Abbott, T. Abbott, S. Abraham, et al., Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo's third observing run, Phys. Rev. D 104, 022004 (2021), arXiv:2101.12130 [gr-qc].

[173] K. Belczynski, J. Klencki, C. E. Fields, A. Olejak, E. Berti, G. Meynet, C. L. Fryer, D. E. Holz, R. O'Shaughnessy, D. A. Brown, T. Bulik, S. C. Leung, K. Nomoto, P. Madau, R. Hirschi, E. Kaiser, S. Jones, S. Mondal, M. Chruslinska, P. Drozda, D. Gerosa, Z. Doctor, M. Giersz, S. Ekstrom, C. Georgy, A. Askar, V. Baibhav, D. Wysocki, T. Natan, W. M. Farr, G. Wiktorowicz, M. Coleman Miller, B. Farr, and J. P. Lasota, Evolutionary roads leading to low effective spins, high black hole masses, and O1/O2 rates for LIGO/Virgo binary black holes, Astron. Astrophys. 636, A104 (2020), arXiv:1706.07053 [astro-ph.HE].

[174] M. Mapelli, N. Giacobbo, F. Santoliquido, and M. C. Artale, The properties of merging black holes and neutron stars across cosmic time, Mon. Not. R. Ast. Soc. 487, 2 (2019), arXiv:1902.01419 [astro-ph.HE].

[175] C. J. Neijssel, A. Vigna-Gómez, S. Stevenson, J. W. Barrett, S. M. Gaebel, F. S. Broekgaarden, S. E. de Mink, D. Szécsi, S. Vinciguerra, and I. Mandel, The effect of the metallicity-specific star formation history on double compact object mergers, Mon. Not. R. Ast. Soc. 490, 3740 (2019).

[176] R. Farmer, M. Renzo, S. E. de Mink, P. Marchant, and S. Justham, Mind the Gap: The Location of the Lower Edge of the Pair-instability Supernova Black Hole Mass Gap, Astrophys. J. 887, 53 (2019), arXiv:1910.12874 [astro-ph.SR].

[177] T. Kinugawa, T. Nakamura, and H. Nakano, Formation of binary black holes similar to GW190521 with a total mass of ~150 M⊙ from Population III binary star evolution, Mon. Not. R. Ast. Soc. 501, L49 (2021), arXiv:2009.06922 [astro-ph.HE].

[178] J. S. Vink, E. R. Higgins, A. A. C. Sander, and G. N. Sabhahit, Maximum black hole mass across cosmic time, Mon. Not. R. Ast. Soc. 504, 146 (2021), arXiv:2010.11730 [astro-ph.HE].

[179] M. Fishbach, Z. Doctor, T. Callister, B. Edelman, J. Ye, R. Essick, W. M. Farr, B. Farr, and D. E. Holz, When are LIGO/Virgo's Big Black-Hole Mergers?, arXiv e-prints, arXiv:2101.07699 (2021), arXiv:2101.07699 [astro-ph.HE].

[180] W. M. Farr, S. Stevenson, M. Coleman Miller, I. Mandel, B. Farr, and A. Vecchio, Distinguishing Spin-Aligned and Isotropic Black Hole Populations With Gravitational Waves, Nature 548, 426 (2017).

[181] D. Gerosa, E. Berti, R. O'Shaughnessy, K. Belczynski,
M. Kesden, D. Wysocki, and W. Gladysz, Spin orientations of merging black holes formed from the evolution of stellar binaries, Phys. Rev. D **98**, 084036 (2018).

[182] D. Kushnir, M. Zaldarriaga, J. A. Kollmeier, and R. Waldman, GW150914: spin-based constraints on the merger time of the progenitor system, Mon. Not. R. Ast. Soc. **462**, 844 (2016), arXiv:1605.03839 [astro-ph.HE].

[183] I. Mandel and R. O'Shaughnessy, Compact binary coalescences in the band of ground-based gravitational-wave detectors, Class. Quant. Grav. **27**, 114007 (2010).

[184] R. O'Shaughnessy, D. Gerosa, and D. Wysocki, Inferences about Supernova Physics from Gravitational-Wave Measurements: GW151226 Spin Misalignment as an Indicator of Strong Black-Hole Natal Kicks, Phys. Rev. Lett. **119**, 011101 (2017).

[185] C. L. Rodriguez, M. Zevin, C. Pankow, V. Kalogera, and F. A. Rasio, Illuminating Black Hole Binary Formation Channels with Spins in Advanced LIGO, Astrophys. J. Lett. **832**, L2 (2016), arXiv:1609.05916 [astro-ph.HE].

[186] S. Stevenson, C. P. L. Berry, and I. Mandel, Hierarchical analysis of gravitational-wave measurements of binary black hole spin-orbit misalignments, Mon. Not. R. Ast. Soc. **471**, 2801 (2017), arXiv:1703.06873 [astro-ph.HE].

[187] S. Vitale, D. Gerosa, C.-J. Haster, K. Chatziioannou, and A. Zimmerman, Impact of bayesian priors on the characterization of binary black hole coalescences, Phys. Rev. Lett. **119**, 251103 (2017).

[188] J. Fuller and L. Ma, Most Black Holes Are Born Very Slowly Rotating, Astrophys. J. Lett. **881**, L1 (2019).

[189] J. Fuller, A. L. Piro, and A. S. Jermyn, Slowing the spins of stellar cores, Mon. Not. R. Ast. Soc. **485**, 3661 (2019).

[190] S. S. Bavera, T. Fragos, Y. Qin, E. Zapartas, C. J. Neijssel, I. Mandel, A. Batta, S. M. Gaebel, C. Kimball, and S. Stevenson, The origin of spin in binary black holes. Predicting the distributions of the main observables of Advanced LIGO, Astron. Astrophys. **635**, A97 (2020), arXiv:1906.12257 [astro-ph.HE].

[191] Y. Qin, P. Marchant, T. Fragos, G. Meynet, and V. Kalogera, On the Origin of Black Hole Spin in High-mass X-Ray Binaries, Astrophys. J. Lett. **870**, L18 (2019), arXiv:1810.13016 [astro-ph.SR].

[192] S. Biscoveanu, M. Isi, S. Vitale, and V. Varma, A new spin on LIGO-Virgo binary black holes, Available as arxiv:2007.09156 (2020), arXiv:2007.09156 [astro-ph.HE].

[193] V. Kalogera, Spin-Orbit Misalignment in Close Binaries with Two Compact Objects, Astrophys. J. **541**, 319 (2000), arXiv:astro-ph/9911417 [astro-ph].

[194] J. Roulet, H. S. Chia, S. Olsen, L. Dai, T. Venumadhav, B. Zackay, and M. Zaldarriaga, Distribution of effective spins and masses of binary black holes from the LIGO and Virgo O1-O3a observing runs, Phys. Rev. D **104**, 083010 (2021), arXiv:2105.10580 [astro-ph.HE].

[195] S. Galaudage, C. Talbot, T. Nagar, D. Jain, E. Thrane, and I. Mandel, Building better spin models for merging binary black holes: Evidence for non-spinning and rapidly spinning nearly aligned sub-populations, arXiv e-prints, arXiv:2109.02424 (2021), arXiv:2109.02424 [gr-qc].

[196] T. A. Callister, C.-J. Haster, K. K. Y. Ng, S. Vitale, and W. M. Farr, Who Ordered That? Unequal-Mass Binary Black Hole Mergers Have Larger Effective Spins, arXiv e-prints, arXiv:2106.00521 (2021), arXiv:2106.00521 [astro-ph.HE].

[197] D. Gerosa and M. Fishbach, Hierarchical mergers of stellar-mass black holes and their gravitational-wave signatures, Nature Astronomy **5**, 749 (2021), arXiv:2105.03439 [astro-ph.HE].

[198] V. Baibhav, D. Gerosa, E. Berti, K. W. K. Wong, T. Helfer, and M. Mould, The mass gap, the spin gap, and the origin of merging binary black holes, Phys. Rev. D **102**, 043002 (2020), arXiv:2004.00650 [astro-ph.HE].

[199] M. Fishbach, D. E. Holz, and B. Farr, Are LIGO's Black Holes Made From Smaller Black Holes?, Astrophys. J. Lett. **840**, L24 (2017), arXiv:1703.06869 [astro-ph.HE].

[200] D. Gerosa and E. Berti, Are merging black holes born from stellar collapse or previous mergers?, Phys. Rev. D **95**, 124046 (2017), arXiv:1703.06223 [gr-qc].

[201] B. McKernan, K. E. S. Ford, R. O'Shaughnessy, and D. Wysocki, Monte Carlo simulations of black hole mergers in AGN discs: Low $\chi_{eff}$ mergers and predictions for LIGO, Mon. Not. R. Ast. Soc. **494**, 1203 (2020), arXiv:1907.04356 [astro-ph.HE].

[202] Y. Huang, C.-J. Haster, J. Roulet, S. Vitale, A. Zimmerman, T. Venumadhav, B. Zackay, L. Dai, and M. Zaldarriaga, Source properties of the lowest signal-to-noise-ratio binary black hole detections, Phys. Rev. D **102**, 103024 (2020), arXiv:2003.04513 [gr-qc].

[203] J. W. Barrett, S. M. Gaebel, C. J. Neijssel, A. Vignagómez, S. Stevenson, C. P. L. Berry, W. M. Farr, and I. Mandel, Accuracy of inference on the physics of binary evolution from gravitational-wave observations, Mon. Not. R. Ast. Soc. **477**, 4685 (2018), arXiv:1711.06287 [astro-ph.HE].

[204] B. Farr, D. E. Holz, and W. M. Farr, Using Spin to Understand the Formation of LIGO and Virgo's Black Holes, Astrophys. J. Lett. **854**, L9 (2018).

[205] M. Fishbach, W. M. Farr, and D. E. Holz, The Most Massive Binary Black Hole Detections and the Identification of Population Outliers, Astrophys. J. Lett. **891**, L31 (2020), arXiv:1911.05882 [astro-ph.HE].

[206] D. Gerosa, E. Berti, R. O'Shaughnessy, K. Belczynski, M. Kesden, D. Wysocki, and W. Gladysz, Spin orientations of merging black holes formed from the evolution of stellar binaries, Phys. Rev. D **98**, 084036 (2018), arXiv:1808.02491 [astro-ph.HE].

[207] I. Mandel, B. Müller, J. Riley, S. E. de Mink, A. Vignagómez, and D. Chattopadhyay, Binary population synthesis with probabilistic remnant mass and kick prescriptions, Mon. Not. R. Ast. Soc. **500**, 1380 (2021), arXiv:2007.03890 [astro-ph.HE].

[208] I. M. Romero-Shaw, K. Kremer, P. D. Lasky, E. Thrane, and J. Samsing, Gravitational waves as a probe of globular cluster formation and evolution, Mon. Not. R. Ast. Soc. **506**, 2362 (2021), arXiv:2011.14541 [astro-ph.HE].

[209] S. Stevenson, F. Ohme, and S. Fairhurst, Distinguishing Compact Binary Population Synthesis Models Using Gravitational Wave Observations of Coalescing Binary Black Holes, Astrophys. J. **810**, 58 (2015), arXiv:1504.07802 [astro-ph.HE].

[210] V. Tiwari, S. Fairhurst, and M. Hannam, Constraining black-hole spins with gravitational wave observations, Astrophys. J. **868**, 140 (2018), arXiv:1809.01401 [gr-qc].

[211] K. W. K. Wong, K. Breivik, K. Kremer, and T. Callister, Joint constraints on the field-cluster mixing fraction, common envelope efficiency, and globular cluster radii from a population of binary hole mergers
via deep learning, Phys. Rev. D 103, 083021 (2021), arXiv:2011.03564 [astro-ph.HE].

[212] D. Wysocki, D. Gerosa, R. O'Shaughnessy, K. Belczynski, W. Gladysz, E. Berti, M. Kesden, and D. E. Holz, Explaining LIGO's observations via isolated binary evolution with natal kicks, Phys. Rev. D 97, 043014 (2018), arXiv:1709.01943 [astro-ph.HE].

[213] Y. Yang, I. Bartos, V. Gayathri, K. E. S. Ford, Z. Haiman, S. Klimenko, B. Kocsis, S. Márka, Z. Márka, B. McKernan, and R. O'Shaughnessy, Hierarchical Black Hole Mergers in Active Galactic Nuclei, Phys. Rev. Lett. 123, 181101 (2019), arXiv:1906.09281 [astro-ph.HE].

[214] M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez, F. Antonini, D. E. Holz, and C. Pankow, One Channel to Rule Them All? Constraining the Origins of Binary Black Holes Using Multiple Formation Pathways, Astrophys. J. 910, 152 (2021), arXiv:2011.10057 [astro-ph.HE].

[215] A. Askar, M. Arca Sedda, and M. Giersz, MOCCASURVEY Database I: Galactic globular clusters harbouring a black hole subsystem, Mon. Not. R. Ast. Soc. 478, 1844 (2018), arXiv:1802.05284 [astro-ph.GA].

[216] S. Banerjee, H. Baumgardt, and P. Kroupa, Stellar-mass black holes in star clusters: implications for gravitational wave radiation, Mon. Not. R. Ast. Soc. 402, 371 (2010), arXiv:0910.3954 [astro-ph.SR].

[217] J. Hong, E. Vesperini, A. Askar, M. Giersz, M. Szkudlarek, and T. Bulik, Binary black hole mergers from globular clusters: the impact of globular cluster properties, Mon. Not. R. Ast. Soc. 480, 5645 (2018), arXiv:1808.04514 [astro-ph.HE].

[218] S. R. Kulkarni, P. Hut, and S. McMillan, Stellar black holes in globular clusters, Nature 364, 421 (1993).

[219] R. M. O'Leary, R. O'Shaughnessy, and F. A. Rasio, Dynamical interactions and the black-hole merger rate of the Universe, Phys. Rev. D 76, 061504 (2007), arXiv:astro-ph/0701887 [astro-ph].

[220] S. F. Portegies Zwart and S. L. W. McMillan, Black Hole Mergers in the Universe, Astrophys. J. Lett. 528, L17 (2000), arXiv:astro-ph/9910061 [astro-ph].

[221] C. L. Rodriguez, S. Chatterjee, and F. A. Rasio, Binary black hole mergers from globular clusters: Masses, merger rates, and the impact of stellar evolution, Phys. Rev. D 93, 084029 (2016), arXiv:1602.02444 [astro-ph.HE].

[222] C. L. Rodriguez, M. Morscher, B. Pattabiraman, S. Chatterjee, C.-J. Haster, and F. A. Rasio, Binary Black Hole Mergers from Globular Clusters: Implications for Advanced LIGO, Phys. Rev. Lett. 115, 051101 (2015), arXiv:1505.00792 [astro-ph.HE].

[223] S. Sigurdsson and L. Hernquist, Primordial black holes in globular clusters, Nature 364, 423 (1993).

[224] C. L. Rodriguez, M. Zevin, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio, and C. S. Ye, Black holes: The next generation—repeated mergers in dense star clusters and their gravitational-wave properties, Phys. Rev. D 100, 043027 (2019), arXiv:1906.10260 [astro-ph.HE].

[225] S. Banerjee, Merger rate density of stellar-mass binary black holes from young massive clusters, open clusters, and isolated binaries: comparisons with LIGO-Virgo-KAGRA results, arXiv e-prints, arXiv:2108.04250 (2021), arXiv:2108.04250 [astro-ph.HE].

[226] S. Banerjee, Stellar-mass black holes in young massive and open stellar clusters - V. comparisons with LIGO-Virgo merger rate densities, Mon. Not. R. Ast. Soc. 503, 3371 (2021), arXiv:2011.07000 [astro-ph.HE].

[227] F. Santoliquido, M. Mapelli, N. Giacobbo, Y. Bouffanais, and M. C. Artale, The cosmic merger rate density of compact objects: impact of star formation, metallicity, initial mass function, and binary evolution, Mon. Not. R. Ast. Soc. 502, 4877 (2021), arXiv:2009.03911 [astro-ph.HE].

[228] F. Antonini and H. B. Perets, Secular Evolution of Compact Binaries near Massive Black Holes: Gravitational Wave Sources and Other Exotica, ApJ 757, 27 (2012), arXiv:1203.2938 [astro-ph.GA].

[229] F. Antonini and F. A. Rasio, Merging black hole binaries in galactic nuclei: Implications for advanced-ligo detections, Astrophys. J. 831, 187 (2016).

[230] M. Arca Sedda, M. Mapelli, M. Spera, M. Benacquista, and N. Giacobbo, Fingerprints of Binary Black Hole Formation Channels Encoded in the Mass and Spin of Merger Remnants, Astrophys. J. 894, 133 (2020), arXiv:2003.07409 [astro-ph.GA].

[231] G. Fragione, E. Grishin, N. W. C. Leigh, H. B. Perets, and R. Perna, Black hole and neutron star mergers in galactic nuclei, Mon. Not. R. Ast. Soc. 488, 47 (2019), arXiv:1811.10627 [astro-ph.GA].

[232] L. Gondán and B. Kocsis, High eccentricities and high masses characterize gravitational-wave captures in galactic nuclei as seen by Earth-based detectors, MNRAS 506, 1665 (2021), arXiv:2011.02507 [astro-ph.HE].

[233] B.-M. Hoang, S. Naoz, B. Kocsis, F. A. Rasio, and F. Dosopoulou, Black Hole Mergers in Galactic Nuclei Induced by the Eccentric Kozai-Lidov Effect, ApJ 856, 140 (2018), arXiv:1706.09896 [astro-ph.HE].

[234] M. Mapelli, M. Dall'Amico, Y. Bouffanais, N. Giacobbo, M. Arca Sedda, M. C. Artale, A. Ballone, U. N. Di Carlo, G. Iorio, F. Santoliquido, and S. Torniamenti, Hierarchical black hole mergers in young, globular and nuclear star clusters: the effect of metallicity, spin and cluster properties, Mon. Not. R. Ast. Soc. 505, 339 (2021), arXiv:2103.05016 [astro-ph.HE].

[235] K. E. S. Ford and B. McKernan, Binary Black Hole Merger Rates in AGN Disks versus Nuclear Star Clusters: Loud beats Quiet, arXiv e-prints, arXiv:2109.03212 (2021), arXiv:2109.03212 [astro-ph.HE].

[236] H. Tagawa, B. Kocsis, Z. Haiman, I. Bartos, K. Omukai, and J. Samsing, Mass-gap Mergers in Active Galactic Nuclei, Astrophys. J. 908, 194 (2021), arXiv:2012.00011 [astro-ph.HE].

[237] Y. Yang, I. Bartos, Z. Haiman, B. Kocsis, Z. Márka, N. C. Stone, and S. Márka, AGN Disks Harden the Mass Distribution of Stellar-mass Binary Black Hole Mergers, Astrophys. J. 876, 122 (2019), arXiv:1903.01405 [astro-ph.HE].

[238] K. Belczynski, D. E. Holz, T. Bulik, and R. O'Shaughnessy, The first gravitational-wave source from the isolated evolution of two stars in the 40-100 solar mass range, Nature 534, 512 (2016), arXiv:1602.04531 [astro-ph.HE].

[239] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O'Shaughnessy, Double Compact Objects. II. Cosmological Merger Rates, Astrophys. J. 779, 72 (2013).
[240] N. Giacobbo and M. Mapelli, The progenitors of compact-object binaries: impact of metallicity, common envelope and natal kicks, Mon. Not. R. Ast. Soc. 480, 2011 (2018), arXiv:1806.00001 [astro-ph.HE].

[241] G. Wiktorowicz, L. Wyrzykowski, M. Chruslinska, J. Klencki, K. A. Rybicki, and K. Belczynski, Populations of Stellar-mass Black Holes from Binary Systems, Astrophys. J. 885, 1 (2019), arXiv:1907.11431 [astro-ph.HE].

[242] S. E. de Mink and K. Belczynski, Merger Rates of Double Neutron Stars and Stellar Origin Black Holes: The Impact of Initial Conditions on Binary Evolution Predictions, Astrophys. J. 814, 58 (2015), arXiv:1506.03573 [astro-ph.HE].

[243] A. Olejak, K. Belczynski, and N. Ivanova, Impact of common envelope development criteria on the formation of LIGO/Virgo sources, Astron. Astrophys. 651, A100 (2021), arXiv:2102.05649 [astro-ph.HE].

[244] G. Costa, A. Bressan, M. Mapelli, P. Marigo, G. Iorio, and M. Spera, Formation of GW190521 from stellar evolution: the impact of the hydrogen-rich envelope, dredge-up, and $^{12}$C(α, γ)$^{16}$O rate on the pair-instability black hole mass gap, Mon. Not. R. Ast. Soc. 501, 4514 (2021), arXiv:2010.02242 [astro-ph.SR].

[245] D. Croon, S. D. McDermott, and J. Sakstein, New physics and the black hole mass gap, Phys. Rev. D 102, 115024 (2020).

[246] R. Farmer, M. Renzo, S. E. de Mink, M. Fishbach, and S. Justham, Constraints from Gravitational-wave Detectors of Binary Black Hole Mergers on the $^{12}$C(α, γ)$^{16}$O Rate, Astrophys. J. Lett. 902, L36 (2020), arXiv:2006.06678 [astro-ph.HE].

[247] E. Farrell, J. H. Groh, R. Hirschi, L. Murphy, E. Kaiser, S. Ekström, C. Georgy, and G. Meynet, Is GW190521 the merger of black holes from the first stellar generations?, Mon. Not. R. Ast. Soc. 502, L40 (2021), arXiv:2009.06585 [astro-ph.SR].

[248] K. Takahashi, The Low Detection Rate of Pair-instability Supernovae and the Effect of the Core Carbon Fraction, Astrophys. J. 863, 153 (2018), arXiv:1807.05373 [astro-ph.HE].

[249] F. Antonini, M. Gieles, and A. Gualandris, Black hole growth through hierarchical black hole mergers in dense star clusters: implications for gravitational wave detections, Mon. Not. R. Ast. Soc. 486, 5008 (2019), arXiv:1811.03640 [astro-ph.HE].

[250] G. Fragione, A. Loeb, and F. A. Rasio, On the Origin of GW190521-like Events from Repeated Black Hole Mergers in Star Clusters, Astrophys. J. Lett. 902, L26 (2020), arXiv:2009.05065 [astro-ph.GA].

[251] C. Kimball, C. Talbot, C. P. L. Berry, M. Carney, M. Zevin, E. Thrane, and V. Kalogera, Black Hole Genealogy: Identifying Hierarchical Mergers with Gravitational Waves, Astrophys. J. 900, 177 (2020), arXiv:2005.00023 [astro-ph.HE].

[252] A. Palmese and C. J. Conselice, GW190521 from the Merger of Ultradwarf Galaxies, Phys. Rev. Lett. 126, 181103 (2021), arXiv:2009.10688 [astro-ph.GA].

[253] H. Tagawa, Z. Haiman, I. Bartos, B. Kocsis, and K. Omukai, Signatures of hierarchical mergers in black hole spin and mass distribution, Mon. Not. R. Ast. Soc. 10.1093/mnras/stab2315 (2021), arXiv:2104.09510 [astro-ph.HE].

[254] U. D. Carlo et al., Merging black holes in young star clusters, Mon. Not. R. Ast. Soc. 487, 2947 (2019).

[255] U. N. Di Carlo, M. Mapelli, Y. Bouffanais, N. Giacobbo, F. Santoliquido, A. Bressan, M. Spera, and F. Haardt, Binary black holes in the pair instability mass gap, Mon. Not. R. Ast. Soc. 497, 1043 (2020), arXiv:1911.01434 [astro-ph.HE].

[256] K. Kremer, M. Spera, D. Becker, S. Chatterjee, U. N. Di Carlo, G. Fragione, C. L. Rodriguez, C. S. Ye, and F. A. Rasio, Populating the Upper Black Hole Mass Gap through Stellar Collisions in Young Star Clusters, Astrophys. J. 903, 45 (2020), arXiv:2006.10771 [astro-ph.HE].

[257] M. Renzo, M. Cantiello, B. D. Metzger, and Y. F. Jiang, The Stellar Merger Scenario for Black Holes in the Pair-instability Gap, Astrophys. J. Lett. 904, L13 (2020), arXiv:2010.00705 [astro-ph.SR].

[258] S. S. Bavera, T. Fragos, M. Zevin, C. P. L. Berry, P. Marchant, J. J. Andrews, S. Coughlin, A. Dotter, K. Kovlakas, D. Misra, J. G. Serra-Perez, Y. Qin, K. A. Rocha, J. Román-Garza, N. H. Tran, and E. Zapartas, The impact of mass-transfer physics on the observable properties of field binary black hole populations, Astron. Astrophys. 647, A153 (2021), arXiv:2010.16333 [astro-ph.HE].

[259] K. Inayoshi, Z. Haiman, and J. P. Ostriker, Hyper-Eddington accretion flows on to massive black holes, Mon. Not. R. Ast. Soc. 459, 3738 (2016), arXiv:1511.02116 [astro-ph.HE].

[260] Z. Roupas and D. Kazanas, Generation of massive stellar black holes by rapid gas accretion in primordial dense clusters, Astron. Astrophys. 632, L8 (2019), arXiv:1911.03915 [astro-ph.GA].

[261] M. Safarzadeh and Z. Haiman, Formation of GW190521 via Gas Accretion onto Population III Stellar Black Hole Remnants Born in High-redshift Minihalos, Astrophys. J. Lett. 903, L21 (2020), arXiv:2009.09320 [astro-ph.HE].

[262] L. A. C. van Son, S. E. De Mink, F. S. Broekgaarden, M. Renzo, S. Justham, E. Laplace, J. Morán-Fraile, D. D. Hendriks, and R. Farmer, Polluting the Pair-instability Mass Gap for Binary Black Holes through Super-Eddington Accretion in Isolated Binaries, Astrophys. J. 897, 100 (2020), arXiv:2004.05187 [astro-ph.HE].

[263] V. De Luca, V. Desjacques, G. Franciolini, P. Pani, and A. Riotto, GW190521 Mass Gap Event and the Primordial Black Hole Scenario, Phys. Rev. Lett. 126, 051101 (2021), arXiv:2009.01728 [astro-ph.CO].

[264] V. Baibhav, E. Berti, D. Gerosa, M. Mapelli, N. Giacobbo, Y. Bouffanais, and U. N. Di Carlo, Gravitational-wave detection rates for compact binaries formed in isolation: LIGO/Virgo O3 and beyond, Phys. Rev. D 100, 064060 (2019).

[265] M. Gallegos-Garcia, C. P. L. Berry, P. Marchant, and V. Kalogera, Binary Black Hole Formation with Detailed Modeling: Stable Mass Transfer Leads to Lower Merger Rates, arXiv e-prints, arXiv:2107.05702 (2021), arXiv:2107.05702 [astro-ph.HE].

[266] F. Santoliquido, M. Mapelli, Y. Bouffanais, N. Giacobbo, U. N. Di Carlo, S. Rastello, M. C. Artale, and A. Ballone, The cosmic merger rate density evolution of compact binaries formed in young star clusters and in isolated binaries, arXiv e-prints, arXiv:2004.09533 (2020).
[267] G. Fragione and B. Kocsis, Black Hole Mergers from an Evolving Population of Globular Clusters, Phys. Rev. Lett. 121, 161103 (2018), arXiv:1806.02351 [astro-ph.GA].

[268] H. Tagawa, Z. Haiman, I. Bartos, and B. Kocsis, Spin Evolution of Stellar-mass Black Hole Binaries in Active Galactic Nuclei, Astrophys. J. 899, 26 (2020), arXiv:2004.11914 [astro-ph.HE].

[269] S. Albrecht, S. Reffert, I. A. G. Snellen, and J. N. Winn, Misaligned spin and orbital axes cause the anomalous precession of DI Herculis, Nature 461, 373 (2009), arXiv:0909.2861 [astro-ph.SR].

[270] S. Albrecht, J. Setiawan, G. Torres, D. C. Fabrycky, and J. N. Winn, The BANANA Project. IV. Two Aligned Stellar Rotation Axes in the Young Eccentric Binary System EP Crucis: Primordial Orientation and Tidal Alignment, Astrophys. J. 767, 32 (2013), arXiv:1211.7065 [astro-ph.SR].

[271] S. Albrecht, J. N. Winn, J. A. Carter, I. A. G. Snellen, and E. J. W. de Mooij, The Banana Project. III. Spin-Orbit Alignment in the Long-period Eclipsing Binary NY Cephei, Astrophys. J. 726, 68 (2011), arXiv:1011.0425 [astro-ph.SR].

[272] S. Albrecht, J. N. Winn, G. Torres, D. C. Fabrycky, J. Setiawan, M. Gillon, E. Jehin, A. Triaud, D. Queloz, I. Snellen, and P. Eggleton, The BANANA Project. V. Misaligned and Precessing Stellar Rotation Axes in CV Velorum, Astrophys. J. 785, 83 (2014), arXiv:1403.0583 [astro-ph.SR].

[273] F. Antonini, C. L. Rodriguez, C. Petrovich, and C. L. Fischer, Precessional dynamics of black hole triples: binary mergers with near-zero effective spin, Mon. Not. R. Ast. Soc. Lett. 480, L58 (2018).

[274] A. S. Hamers, G. Fragione, P. Neunteufel, and B. Kocsis, First- and second-generation black hole and neutron star mergers in 2+2 quadruples: population statistics, MNRAS 506, 5345 (2021), arXiv:2103.03782 [astro-ph.HE].

[275] B. Liu, D. Lai, and Y.-H. Wang, Black Hole and Neutron Star Binary Mergers in Triple Systems. II. Merger Eccentricity and Spin-Orbit Misalignment, Astrophys. J. 881, 41 (2019).

[276] C. L. Rodriguez and F. Antonini, A Triple Origin for the Heavy and Low-spin Binary Black Holes Detected by LIGO/VIRGO, Astrophys. J. 863, 7 (2018).

[277] C. Petrovich and F. Antonini, Greatly Enhanced Merger Rates of Compact-object Binaries in Non-spherical Nuclear Star Clusters, Astrophys. J. 846, 146 (2017), arXiv:1705.05848 [astro-ph.HE].

[278] Y. Su, D. Lai, and B. Liu, Spin-orbit misalignments in tertiary-induced binary black-hole mergers: Theoretical analysis, Phys. Rev. D 103, 063040 (2021), arXiv:2010.11951 [gr-qc].

[279] H. Yu, S. Ma, M. Giesler, and Y. Chen, Spin and eccentricity evolution in triple systems: From the Lidov-Kozai interaction to the final merger of the inner binary, Phys. Rev. D 102, 123009 (2020), arXiv:2007.12978 [gr-qc].

[280] J. Stegmann and F. Antonini, Flipping spins in mass transferring binaries and origin of spin-orbit misalignment in binary black holes, Phys. Rev. D 103, 063007 (2021), arXiv:2012.06329 [astro-ph.HE].

[281] C. Chan, B. Müller, and A. Heger, The impact of falloff on the compact remnants and chemical yields of

core-collapse supernovae, Mon. Not. R. Ast. Soc. 495, 3751 (2020), arXiv:2003.04320 [astro-ph.SR].

[282] W. M. Farr, K. Kremer, M. Lyutikov, and V. Kalogera, Spin Tilts in the Double Pulsar Reveal Supernova Spin Angular-momentum Production, Astrophys. J. 742, 81 (2011), arXiv:1104.5001 [astro-ph.HE].

[283] I. Mandel and T. Fragos, An alternative interpretation of GW190412 as a binary black hole merger with a rapidly spinning secondary, Astrophys. J. Lett. 895, L28 (2020), arXiv:2004.09288 [astro-ph.HE].

[284] Y. Qin, T. Fragos, G. Meynet, J. Andrews, M. Sørensen, and H. F. Song, The spin of the second-born black hole in coalescing binary black holes, Astron. Astrophys. 616, A28 (2018), arXiv:1802.05738 [astro-ph.SR].

[285] I. Mandel and S. E. de Mink, Merging binary black holes formed through chemically homogeneous evolution in short-period stellar binaries, Mon. Not. R. Ast. Soc. 458, 2634 (2016).

[286] P. Marchant, N. Langer, P. Podsiadlowski, T. M. Tauris, and T. J. Moriya, A new route towards merging massive black holes, Astron. Astrophys. 588, A50 (2016).

[287] T. Chiba and S. Yokoyama, Spin distribution of primordial black holes, Progress of Theoretical and Experimental Physics 2017, 083E01 (2017), arXiv:1704.06573 [gr-qc].

[288] V. De Luca, V. Desjacques, G. Franciolini, A. Malhotra, and A. Riotto, The initial spin probability distribution of primordial black holes, J. Cosmology Astropart. Phys. 2019, 018 (2019), arXiv:1903.01179 [astro-ph.CO].

[289] M. Safarzadeh, W. M. Farr, and E. Ramirez-Ruiz, A Trend in the Effective Spin Distribution of LIGO Binary Black Holes with Mass, Astrophys. J. 894, 129 (2020), arXiv:2001.06490 [gr-qc].

[290] C. Hoy, S. Fairhurst, M. Hannam, and V. Tiwari, Understanding how fast black holes spin by analysing data from the second gravitational-wave catalogue, arXiv e-prints (2021), arXiv:2110.13542 [gr-qc].

[291] I. M. Romero-Shaw, N. Farrow, S. Stevenson, E. Thrane, and X.-J. Zhu, On the origin of GW190425, Mon. Not. R. Ast. Soc. 496, L64 (2020), arXiv:2001.06492 [astro-ph.HE].

[292] M. Safarzadeh, E. Ramirez-Ruiz, and E. Berger, Does GW190425 Require an Alternative Formation Pathway than a Fast-merging Channel?, Astrophys. J. 900, 13 (2020), arXiv:2001.04502 [astro-ph.HE].

[293] A. Heger, C. L. Fryer, S. E. Woosley, N. Langer, and D. H. Hartmann, How Massive Single Stars End their Life, Astrophys. J. 591, 288 (2003).

[294] S. E. Woosley, A. Heger, and T. A. Weaver, The evolution and explosion of massive stars, Rev. of Mod. Phys. 74, 1015 (2002).

[295] M. W. Coughlin, T. Dietrich, B. Margalit, and B. D. Metzger, Multimessenger Bayesian parameter inference of a binary neutron star merger, Mon. Not. R. Ast. Soc. 489, L91 (2019), arXiv:1812.04803 [astro-ph.HE].

[296] D. Radice, A. Perego, K. Hotokezaka, S. A. Fromm, S. Bernuzzi, and L. F. Roberts, Binary Neutron Star Mergers: Mass Ejection, Electromagnetic Counterparts, and Nucleosynthesis, Astrophys. J. 869, 130 (2018), arXiv:1809.11161 [astro-ph.HE].

[297] J. M. Lattimer and D. N. Schramm, Black-Hole-Neutron-Star Collisions, Astrophys. J. Lett. 192, L145 (1974).
[298] K. Kiuchi, K. Kyutoku, M. Shibata, and K. Taniguchi, Revisiting the Lower Bound on Tidal Deformability Derived by AT 2017gfo, Astrophys. J. Lett. 876, L31 (2019), arXiv:1903.01466 [astro-ph.HE].

[299] F. Foucart, T. Hinderer, and S. Nissanke, Remnant baryon mass in neutron star-black hole mergers: Predictions for binary neutron star mimickers and rapidly spinning black holes, Phys. Rev. D 98, 081501 (2018), arXiv:1807.00011 [astro-ph.HE].

[300] C. J. Krüger and F. Foucart, Estimates for disk and ejecta masses produced in compact binary mergers, Phys. Rev. D 101, 103002 (2020), arXiv:2002.07728 [astro-ph.HE].

[301] H. T. Janka, T. Eberl, M. Ruffert, and C. L. Fryer, Black Hole-Neutron Star Mergers as Central Engines of Gamma-Ray Bursts, Astrophys. J. Lett. 527, L39 (1999), arXiv:astro-ph/9908290 [astro-ph].

[302] R. Mochkovitch, M. Hernanz, J. Isern, and X. Martin, Gamma-ray bursts as collimated jets from neutron star/black hole mergers, Nature 361, 236 (1993).

[303] S. E. Woosley, Gamma-Ray Bursts from Stellar Mass Accretion Disks around Black Holes, Astrophys. J. 405, 273 (1993).

[304] B. Margalit and B. D. Metzger, The Multi-messenger Matrix: The Future of Neutron Star Merger Constraints on the Nuclear Equation of State, Astrophys. J. Lett. 880, L15 (2019), arXiv:1904.11995 [astro-ph.HE].

[305] B. Liu and D. Lai, Hierarchical black hole mergers in multiple systems: constrain the formation of GW190412-, GW190814-, and GW190521-like events, Mon. Not. R. Ast. Soc. 502, 2049 (2021), arXiv:2009.10068 [astro-ph.HE].

[306] W. Lu, P. Beniamini, and C. Bonnerot, On the formation of GW190814, Mon. Not. R. Ast. Soc. 500, 1817 (2021), arXiv:2009.10082 [astro-ph.HE].

[307] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, et al., Search for the isotropic stochastic background using data from Advanced LIGO's second observing run, prd 100, 061101 (2019), arXiv:1903.02886 [gr-qc].

[308] E. Vangioni, S. Goriely, F. Daigne, P. François, and K. Belczynski, Cosmic neutron-star merger rate and gravitational waves constrained by the r-process nucleosynthesis, Mon. Not. R. Ast. Soc. 455, 17 (2016), arXiv:1501.01115 [astro-ph.HE].

[309] E. Vangioni and K. A. Olive, The cosmic evolution of magnesium isotopes, Mon. Not. R. Ast. Soc. 484, 3561 (2019), arXiv:1809.10514 [astro-ph.GA].

[310] E. Vangioni, K. A. Olive, T. Prestegard, J. Silk, P. Peitjean, and V. Mandic, The impact of star formation and gamma-ray burst rates at high redshift on cosmic chemical evolution and reionization, Mon. Not. R. Ast. Soc. 447, 2575 (2015), arXiv:1409.2462 [astro-ph.GA].

[311] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O'Shaughnessy, Double Compact Objects. I. The Significance of the Common Envelope on Merger Rates, apj 759, 52 (2012), arXiv:1202.4901 [astro-ph.HE].

[312] M. Chruslinska, G. Nelemans, and K. Belczynski, The influence of the distribution of cosmic star formation at different metallicities on the properties of merging double compact objects, Mon. Not. R. Ast. Soc. 482, 5012 (2019), arXiv:1811.03565 [astro-ph.HE].

[313] E. Thrane and J. D. Romano, Sensitivity curves for

searches for gravitational-wave backgrounds, prd 88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM].

[314] W. M. Farr, M. Fishbach, J. Ye, and D. E. Holz, A Future Percent-level Measurement of the Hubble Expansion at Redshift 0.8 with Advanced LIGO, Astrophys. J. Lett. 883, L42 (2019), arXiv:1908.09084 [astro-ph.CO].

[315] S. R. Taylor, J. R. Gair, and I. Mandel, Cosmology using advanced gravitational-wave detectors alone, Phys. Rev. D 85, 023535 (2012), arXiv:1108.5161 [gr-qc].

[316] R. Abbott, T. Abbott, S. Abraham, et al., Constraints on the cosmic expansion history from the third LIGO-Virgo-KAGRA Gravitational-Wave Transient Catalog, Available as LIGO-P2100185 (2021).

[317] R. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Reviews in Relativity 23, 3 (2020).

[318] R. Abbott, T. Abbott, S. Abraham, et al., GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run - O3 search sensitivity estimates, Zenodo, available as https://zenodo.org/deposit/5546676 10.5281/zenodo.5546675 (2021).

[319] V. Tiwari, Estimation of the sensitive volume for gravitational-wave source populations using weighted Monte Carlo integration, Class. Quant. Grav. 35, 145009 (2018), arXiv:1712.00482 [astro-ph.HE].

[320] R. Abbott, T. Abbott, S. Abraham, et al., GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run - O1+O2+O3 search sensitivity estimates, Zenodo, available as https://zenodo.org/deposit/5636816 10.5281/zenodo.5636816 (2021).

[321] D. Wysocki and R. O'Shaughnessy, Popmodels o3a aps april 2021 presentation, https://gitlab.com/dwysocki/pop-models-o3a-aps-april-2021 (2021).

[322] C. Cutler and E. E. Flanagan, Gravitational waves from merging compact binaries: How accurately can one extract the binary's parameters from the inspiral waveform?, Physical Review D 49, 2658 (1994), publisher: American Physical Society.

[323] M. Hannam, D. A. Brown, S. Fairhurst, C. L. Fryer, and I. W. Harry, WHEN CAN GRAVITATIONAL-WAVE OBSERVATIONS DISTINGUISH BETWEEN BLACK HOLES AND NEUTRON STARS?, The Astrophysical Journal 766, L14 (2013), publisher: American Astronomical Society.

[324] C. P. L. Berry, I. Mandel, H. Middleton, L. P. Singer, A. L. Urban, A. Vecchio, S. Vitale, K. Cannon, B. Farr, W. M. Farr, P. B. Graff, C. Hanna, C.-J. Haster, S. Mohapatra, C. Pankow, L. R. Price, T. Sidery, and J. Veitch, PARAMETER ESTIMATION FOR BINARY NEUTRON-STAR COALESCENCES WITH REALISTIC NOISE DURING THE ADVANCED LIGO ERA, The Astrophysical Journal 804, 114 (2015), publisher: American Astronomical Society.

[325] B. Farr, C. P. L. Berry, W. M. Farr, C.-J. Haster, H. Middleton, K. Cannon, P. B. Graff, C. Hanna, I. Mandel, C. Pankow, L. R. Price, T. Sidery, L. P. Singer, A. L. Urban, A. Vecchio, J. Veitch, and S. Vitale, PARAMETER ESTIMATION ON GRAVITATIONAL WAVES FROM NEUTRON-STAR BINARIES WITH
SPINNING COMPONENTS, The Astrophysical Journal 825, 116 (2016), publisher: American Astronomical Society.

[326] K. K. Ng, S. Vitale, A. Zimmerman, K. Chatzioioannou, D. Gerosa, and C.-J. Haster, Gravitational-wave astrophysics with effective-spin measurements: Asymmetries and selection biases, Phys. Rev. D 98, 083007 (2018).

[327] D. Wysocki, R. O'Shaughnessy, L. Wade, and J. Lange, Inferring the neutron star equation of state simultaneously with the population of merging neutron stars, arXiv e-prints, arXiv:2001.01747 (2020), arXiv:2001.01747 [gr-qc].

[328] S. Biscoveanu, C. Talbot, and S. Vitale, Sources of systematic error in gravitational-wave measurements of the binary neutron star mass distribution, arXiv:2111.13619 [astro-ph, physics:gr-qc] (2021), arXiv: 2111.13619.

[329] R. Abbott, T. D. Abbott, S. Abraham, et al, Ligo Scientific Collaboration, VIRGO Collaboration, and KAGRA Collaboration, Observation of Gravitational Waves from Two Neutron Star-Black Hole Coalescences, Astrophys. J. Lett. 915, L5 (2021), arXiv:2106.15163 [astro-ph.HE].

[330] P. Ajith, Addressing the spin question in gravitational-wave searches: Waveform templates for inspiralling compact binaries with nonprecessing spins, Phys. Rev. D 84, 084037 (2011), arXiv:1107.1267 [gr-qc].

[331] P. Ajith, N. Fotopoulos, S. Privitera, A. Neunzert, N. Mazumder, and A. J. Weinstein, Effectual template bank for the detection of gravitational waves from inspiralling compact binaries with generic spins, Phys. Rev. D 89, 084041 (2014).

[332] I. Harry, S. Privitera, A. Bohé, and A. Buonanno, Searching for gravitational waves from compact binaries with precessing spins, Physical Review D 94, 10.1103/physrevd.94.024012 (2016).