# Searching for binary black hole sub-populations in gravitational wave data
using binned Gaussian processes

Anarya Ray<EMAIL_ADDRESS>University of Wisconsin-Milwaukee, Milwaukee, WI
53201, USA  Ignacio Magaña Hernandez<EMAIL_ADDRESS>McWilliams Center
for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh,
PA 15213, USA  Katelyn Breivik McWilliams Center for Cosmology, Department of
Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA  Jolien
Creighton University of Wisconsin-Milwaukee, Milwaukee, WI 53201, USA

###### Abstract

Astrophysically motivated population models for binary black hole observables
are often insufficient to capture the imprints of multiple formation channels.
This is mainly due to the strongly parametrized nature of such investigations.
Using a non-parametric model for the joint population-level distributions of
binary black hole component masses and effective inspiral spins, we find hints
of multiple subpopulations in the third gravitational wave transient catalog.
The higher (more positive) spin subpopulation is found to have a mass-spectrum
without any feature at the $30-40M_{\odot}$, which is consistent with the
predictions of isolated stellar binary evolution, simulations for which place
the pile up due to pulsational pair-instability supernovae near $50M_{\odot}$
or higher. The other sub-population with effective spins closer to zero shows
a feature at $30-40M_{\odot}$ and is consistent with binary black holes formed
dynamically in globular clusters, which are expected to peak around
$30M_{\odot}$. We also compute merger rates for these two subpopulations and
find that they are consistent with the theoretical predictions of the
corresponding formation channels. We validate our results by checking their
robustness against variations of several model configurations and by analyzing
large simulated catalogs with the same model.

##  1 Introduction

The formation of compact object binaries that merge within a Hubble time via
gravitational wave (GW) emission, particularly stellar-mass binary black holes
(BBHs), is poorly understood yet critically important to several fields of
modern astrophysics (Mandel & Farmer, 2022; Mapelli, 2020a, b). The proposed
formation scenarios for such BBH systems can be broadly classified into three
categories: isolated evolution of massive stellar binaries that undergo
orbital hardening, via either common-envelope, stable mass transfer, or
chemical mixing (e.g., Postnov & Yungelson, 2006; Paczynski, 1976; van den
Heuvel et al., 2017; Marchant et al., 2016; Mandel & de Mink, 2016), dynamical
assembly assisted by either a tertiary companion, multiple exchanges in dense
clusters, or gas-assisted migration (e.g., Wen, 2003; Antonini & Perets, 2012;
Benacquista & Downing, 2013; Bartos et al., 2017), and hierarchical mergers in
star clusters and AGN disks (Gerosa & Fishbach, 2021). These formation
channels are usually characterized by several unknown parameters and initial
conditions which can uniquely shape and correlate the observable features in
an astrophysical BBH population, and hence be constrained from a catalog of GW
observations.

At the end of its third observing run, the LIGO-Virgo-KAGRA detector network
(LVK, Aasi et al., 2015; Acernese et al., 2015; Akutsu et al., 2021) unveiled
a catalog of GW observations (GWTC-3, Abbott et al., 2023a) around 70 of which
were confidently111With a false alarm rate of less than 1 per year identified
to be BBH signals. Studying the ensemble properties of these BBHs has offered
a plethora of new insights into the underlying astrophysics of their
progenitor systems and environments. With the LVK’s ongoing fourth observing
run promising a large number of new detections, the search for the imprints of
BBH formation mechanisms on the observable population properties of BBHs will
continue to enable exciting new discoveries. Traditional approaches towards
conducting such investigations have relied on the theoretical predictions of
one or more formation models to construct a fiducial form of the distribution
function of BBH parameters. The hyperparameters characterizing this
distribution function encode information on the astrophysical processes
underlying the corresponding formation scenarios and are inferred from a
collection of GW observations (e.g. Fishbach & Holz, 2017; Edelman et al.,
2021; Kimball et al., 2021).

These parametric population models, while useful in extracting information on
a restrictive set of astrophysical assumptions regarding BBH formation, are
not flexible enough to capture all of the physics that underlies BBH formation
and are also susceptible to model-induced biases. For example, BBHs formed in
isolation are expected to be distributed as a Power-law in component masses
extending up to a maximum range followed by a steep fall-off known as the
pair-instability gap and a small bump right before the truncation arising from
the pile-up due to pulsational pair-instability supernovae (PPISN, Fishbach &
Holz, 2017; Woosley, 2017; Talbot & Thrane, 2018). However, when this
Powerlaw+Peak model for the BBH mass distribution was constrained using GWTC-3
data (Abbott et al., 2023b), several concerns regarding our understanding of
massive stellar interiors came to light. Firstly, the location of the PPISN
feature was found to be within the $30-40M_{\odot}$ mass range, leading to
tension with the predictions of stellar evolution simulations which place it
above $50M_{\odot}$ (see eg, Hendriks et al., 2023, and references therein).
Secondly, flexible alternatives to the Powerlaw+Peak model have found hints of
additional features in the BBH mass spectrum (Edelman et al., 2023; Toubiana
et al., 2023; Abbott et al., 2023b; Callister & Farr, 2023; Tiwari, 2021),
which might be indicative of additional physics at play, beyond the
restrictions built into the underlying assumptions of Powerlaw+Peak.

A plausible explanation for both of these issues could be the existence of
multiple sub-populations of BBHs in GWTC-3 data, each corresponding to
different formation channels. The peak at the $30-40M_{\odot}$ range, for
example, can be explained by a sub-population of dynamically formed BBHs in
globular clusters, which are expected to peak around the same mass range
(Antonini et al., 2023; Wong et al., 2021). Given that the distribution of the
effective inspiral spin parameters for dynamically formed BBHs in globular
clusters is expected to be more symmetric about zero than that of BBHs formed
in isolation (Mapelli et al., 2022; Chattopadhyay et al., 2023; Rodriguez et
al., 2018, 2022), observational evidence for these two sub-populations
coexisting together can manifest in the form of population-level correlations
between BBH masses and spins, and hence be searched for in current and future
GW catalogs (Baibhav et al., 2023). However, several other astrophysical
processes such as mass-ratio reversal in isolated binaries (Broekgaarden et
al., 2022) can introduce mass-spin correlations in the BBH population. Hence,
the use of flexible population models that can account for a wide variety of
known, as well as previously unmodelled astrophysical scenarios, is crucial to
these kinds of investigations.

In this letter, using a non-parametric model for the joint population-level
distribution of BBH component masses and effective inspiral spins based on
binned Gaussian processes (BGP, Ray et al., 2023; Mohite, 2022; Mandel et al.,
2017), we find hints of multiple subpopulations of BBHs in GWTC-3 data, one
consistent with the predictions of isolated BBH formation, and the other
dynamical. In particular, we find two distinct shapes in the BBH mass
spectrum, one of which has a powerlaw-like fall-off with no feature in the
$30-40M_{\odot}$ region and is associated with higher (more positive) values
of effective inspiral spins. The other one is found to have a bump at the
$30-40M_{\odot}$ and is associated with smaller values of effective inspiral
spins. Furthermore, BBHs with at least one component in the $30-40M_{\odot}$
range are found to have a more symmetric distribution of effective inspiral
spins about zero, as compared to that of the complementary case which is
characterized by a spin distribution skewed towards more positive values. We
also constrain the total merger rate of BBHs within different mass and
effective spin ranges and find that they are consistent with the theoretical
predictions of the corresponding formation models.

Given the size of contemporary datasets, we restrict the dimensionality of our
population model to account for the distributions of only the best-measured
among BBH observables that are relevant to the astrophysics we aim to extract.
In particular, we only model the joint distribution of source frame component
masses and effective inspiral spins using BGPs, and fix that of other BBH
observables to be of much simpler and less flexible functional forms. While
this restriction amounts to ignoring higher dimensional correlations in the
BBH population that can in principle bias our conclusions, we validate our
analysis using large simulated catalogs and demonstrate that the ignored
correlations cannot manifest into the trends we find in the joint mass and
effective spin distributions. Even though a higher dimensional BGP on masses,
redshifts, and all of the spin parameters (component spin magnitudes and
tilts) would be the most generalized and robust approach to this
investigation, we expect such a model to yield uninformative constraints given
GWTC-3-sized datasets. However, as catalogs continue to grow, we aim to
implement higher dimensional BGP models for obtaining a more generic and
unbiased characterization of the various sub-populations of BBHs.

In a previous study, Godfrey et al. (2023) had also looked at mass-spin
correlations in the BBH population using a mixture of semi-parametric models
and had come to conclusions that are broadly consistent with ours except for a
few subtle distinctions. They identify a sub-population of binaries with a
mass distribution that demonstrates a strong peak near $10M_{\odot}$ followed
by a sharp fall-off, and a spin distribution that slightly prefers aligned
components over isotropic ones. They identify a second sub-population that
demonstrates a peak in the 30-40$M_{\odot}$ range and has a spin distribution
consistent with that of the other subpopulation. In addition to component
spins, they infer an effective spin distribution that is completely consistent
between the two subpopulations they identify. On the other hand, our study
identifies two subpopulations that correspond to significantly different
effective spin distributions. The astrophysical implications of Godfrey et al.
(2023) are therefore slightly different from ours even though both are broadly
in agreement regarding the existence of a subpopulation in the mass-spectrum
that falls off steadily beyond 15$M_{\odot}$. These differences might be
attributed to the differences in our modeling choices which are discussed
further in Sec. 2.

This letter is organized as follows. In Sec. 2 we describe in detail our
flexible population model and the associated inference framework. In Sec. 3,
we summarize our results for GWTC-3 and describe the data used to obtain it.
In Sec. 4, we validate our model by analyzing large simulated catalogs of five
different BBH populations and demonstrate that our results for real data are
robust against known sources of systematic biases. In Sec. 5, we discuss the
astrophysical implications of our findings. Lastly, In Sec. 6, we conclude
with a summary of our model, findings, and future investigations.

##  2 Population model and hierarchical inference

To search for multiple sub-populations in the BBH mass spectrum that
potentially correspond to different distributions of effective inspiral spin,
we construct a flexible model for the population properties of BBH
observables, agnostic of astrophysical predictions. Specifically, we model the
joint distribution of BBH component masses $(m_{1},m_{2})$, and effective
inspiral spins $(\chi_{\text{eff}})$ as a piece-wise binned function, such as:

$\frac{dN}{dm_{1}dm_{2}dzd\chi_{\rm{eff}}}(m_{1},m_{2},z,\chi_{\rm{eff}}|\vec{n},\kappa)=\sum_{\gamma}\frac{n^{\gamma}}{m_{1}m_{2}}\frac{dV}{dz}T_{\rm{obs}}(1+z)^{\kappa-1}\times\begin{cases}1&\textrm{if
}(m_{1},m_{2},\chi_{\rm{eff}})\in\gamma^{th}\textrm{ bin}\\\
0&\textrm{otherwise}\end{cases}$ (1)

where, the LHS is the number of CBCs per component masses, effective spin and
redshift, $n^{\gamma}$ is the merger rate density per comoving volume, source-
frame time, log component mass, and effective spin in the $\gamma^{th}$ bin,
and $\kappa$ is the parameter that controls the redshift evolution of the
overall merger rate. Simultaneously inferring the merger rate density in each
bin from multiple uncertain measurements amounts to learning the shape of the
population distribution from GW data up to the resolution limit imposed by our
choice of binning.

While Eq. (1) is flexible enough to capture any correlation between the mass
and effective spin populations, it is unable to inform on the existence of
mass-redshift and spin-redshift correlations. Even though contemporary
datasets are mostly uninformative on the existence of mass-redshift
correlations (Abbott et al., 2023b), particularly in the context of binned
models (Ray et al., 2023), the studies of Heinzel et al. (2023) and Biscoveanu
et al. (2022) have revealed strong evidence for spin-redshift correlations in
the same. Biscoveanu et al. (2022) has further shown that certain population
models designed to search only for a mass-spin correlation can falsely infer
their existence from a simulated catalog with spin-redshift correlations. We
show is Sec. 4 that our results are robust against such model-induced
systematics by testing our model on large simulated catalogs characterized by
various correlated populations. See Rinaldi et al. (2023); Karathanasis et al.
(2023) for more investigations on mass-redshift correlations using alternative
models.

To infer the merger rate density in each bin from multiple uncertain
measurements of BBH masses and effective spins we implement the framework of
Bayesian hierarchical inference (Thrane & Talbot, 2019; Mandel et al., 2019;
Loredo, 2004; Wysocki et al., 2019). By modeling the occurrence of BBHs as an
inhomogeneous Poisson process, we construct the likelihood function of the
rate densities from a collection of single event measurements and use a
simulated population of detectable BBHs to account for Malmquist biases in the
inferred distributions. In the context of our population model, the likelihood
function takes the following form:

$\log
p(\vec{d}|\vec{n})=-\sum_{\gamma}n^{\gamma}\left<VT\right>^{\gamma}+\sum_{i}\log\left(\sum_{\gamma}n^{\gamma}w^{\gamma}(d_{i})\right)$
(2)

where $w^{\gamma}(d_{i})$ is the posterior support of the $\gamma^{th}$ bin
from the $i^{th}$ event and $\left<VT\right>^{\gamma}$ is the sensitive time-
volume within which BBHs in the $\gamma^{th}$ bin are expected to be
detectable (Ray et al., 2023).

To regularize the inferred shape of the population distribution over sparse
regions of parameter space we choose the prior on logarithmic rate densities
to be a Gaussian process (GP) with an exponential quadratic kernel. The means,
correlation lengths, and covariance amplitudes characterizing the GP are
themselves modeled using normal, log-normal, and half-normal priors
respectively (Ray et al., 2023). Using Hamiltonian Monte Carlo (HMC) sampling
implemented through the No U-Turn Sampler (Neal, 2011; Homan & Gelman, 2014),
we infer the GP hyper-parameters along with the rate densities characterizing
our population model from their joint posterior distribution. The code
developed to implement this analysis is publicly available as the python
package gppop222https://github.com/AnaryaRay1/gppop/tree/spin-dev, which in
turn relies on the Pymc package (Oriol et al., 2023) for conducting the HMC
sampling.

Previous studies by Farr (2019); Essick & Farr (2022) have demonstrated that
the convergence of the Monte Carlo sums used to compute the posterior weights
and detectable time-volumes of Eq. (2) needs to be ensured in order to avoid
biases in the corresponding hierarchical inference. Following the
implementation of their convergence criteria in the context of the binned
model, as derived by Ray et al. (2023), we verify that all of our inferred
hyper-parameter samples support regions of parameter space that have enough
effectively independent Monte Carlo samples, both for event-specific
posteriors and detectable injections to avoid biases resulting from improperly
converged estimates of posterior weights and detectable time-volumes.

Even though we do not fit for the redshift evolution parameter $\kappa$ and
instead fix it to be $\kappa=2.9$, we show in appendix A that varying $\kappa$
across its measured confidence interval leaves our results unchanged. The
reason for this choice is as follows. For fixed $\kappa$, the posterior
weights and sensitive time-volumes of Eq. (2) are pre-computable leading to
computationally cheap inference of the joint population posterior
distribution. Inferring $\kappa$ amounts to re-computing the weights and time-
volumes for every step of the sampling which is intractable for our current
Central Processing Unit (CPU) based implementation. While parallelizing these
computations on Graphics Processing Units (GPUs) offers a potential solution
to this problem, the associated code developments are part of an ongoing study
and beyond the scope of this work given the insensitivity of our conclusions
to the value of $\kappa$. We however note that our current implementation is
much more scalable on CPUs than the earlier ones of Ray et al. (2023). See
appendix C for more details.

We further note that while previous attempts (Godfrey et al., 2023; Li et al.,
2023) at obtaining data-driven constraints on mass-spin correlations in the
BBH population have accounted for more spin parameters such as component spin
magnitudes and orientations, most of these studies have imposed several
restrictions on the nature of mass-spin correlations, in contrast to our
flexible model. In particular, they have used mixture models wherein, each
mixture component is assumed to have an uncorrelated mass-spin distribution.
Therefore, within their framework, mass-spin correlations can only manifest in
the form of multiple sub-populations coexisting together. Hence, mass-spin
correlations arising from the astrophysical processes within a single
formation channel (such as mass-ratio reversal during stable mass transfer
introducing a spin-mass-ratio anti-correlation during field formation
Broekgaarden et al., 2022), can potentially bias the astrophysical conclusions
of such models. On the other hand, we impose no similar restrictions on our
joint mass-effective-spin distribution and reconstruct its functional form
directly from the data, up to the resolution limit imposed by our choice of
binning.

An additional source of systematics in our population model involves the
potential sensitivity of the inferred shapes of the population distribution to
binning choices. Similar to the example presented in Ray et al. (2023) we show
in appendix B that our results are stable against multiple binning choices by
varying the mass and $\chi_{\rm{eff}}$ bin locations and demonstrating that
the trends we infer are robust against this variation.

##  3 Results

We re-analyze public LVK data (Abbott et al., 2023c) comprising all the BBH
events that were observed through GWTC-3 with a FAR of less than 1 per year.
Following previous works, we exclude the outlier event GW190814 given the
uncertainty regarding its system of origin, which leaves us with a set of 69
high confidence BBH observations (Abbott et al., 2021, 2020, 2023b; Essick et
al., 2022). For each of these events, we use $(m_{1},m_{2},\chi_{\rm{eff}})$
samples publicly released by the LVK (Abbott et al., 2023c) to compute the
posterior weights in Eq. (2). Specifically, following (Abbott et al., 2023b),
we convert detector frame mass and luminosity distance samples to source frame
by assuming a particular cosmological model, which we choose to be Planck 2015
(Ade et al., 2016). For further details of the single-event PE analyses see
(Abbott et al., 2023b). We also use LVK’s publicly released set of detectable
injections (Abbott et al., 2023c) to compute the detectable time-volumes of
Eq. (2).

Figure 1:  Conditional distributions of BBH component masses and spins
obtained from GWTC-3. The parametric inferences of Abbott et al. (2023b) that
used the Powerlaw+Peak and Gaussian effective spin models are overplotted for
comparison. The top panels show the primary (left) and secondary (right) mass
distributions conditioned on two different ranges of effective inspiral spins,
namely $0.05<\chi_{\rm{eff}}<0.1$, and $0.1<\chi_{\rm{eff}}<0.2$. The bottom
panels show the density (left) and cumulative densitiy (right) functions
corresponding to the distribution of effective inspiral spins, conditioned on
two different mass ranges, namely $\\{m_{1}\in(30M_{\odot},40M_{\odot})\text{
or }m_{2}\in(30M_{\odot},40M_{\odot})\\}$, and its complement. The blue dot on
the lower right figure is the point (0,0.5).

As seen in Figure 1, we find trends indicative of at least two subpopulations
in the BBH mass-spectrum that differ by 90% in shape for the $30-40M_{\odot}$
bin and are each conditional on different ranges of effective inspiral spin.
In particular, we find that for $0.1<\chi_{\rm{eff}}<0.2$, the primary and
secondary mass distributions resemble a power-law like fall off with no
feature in the $30-40M_{\odot}$ range. On the other hand, for
$0.05<\chi_{\rm{eff}}<0.1$, we find that the primary and secondary mass-
distributions both show an excess of BBHs in the $30-40M_{\odot}$ range.
Similarly, constraints on the conditional spin distribution depict that
binaries with one or both components in the $30-40M_{\odot}$ mass range have a
more symmetric effective spin distribution about zero as compared to other
binaries which are found to have a skewed distribution of $\chi_{\rm{eff}}$
that prefers more positive values. We also compute the fraction of events with
$\chi_{\text{eff}}<0$ to be $21^{+43}_{-18}\%$ for binaries with at least one
component in the the $30-40M_{\odot}$ mass range and $10^{+22}_{-7}\%$ for all
other binaries.

##  4 Validation Study

Figure 2: Conditional mass-spin distributions for the three simulated
catalogs. The uncorrelated population is shown on the left, the one with
$\chi_{eff}\mbox{-}q$ correlations in the middle and the one with spin-
redshift correlations on the right. Note that for the right panels, the
reconstructed $\chi_{eff}$ distributions are not expected to match the
injected population since our model assumes the joint $\chi_{eff}-z$
distribution to be separable, i.e.,
$p(\chi_{\text{eff}},z)=p(\chi_{\text{eff}})p(z)$.

To validate the robustness of our results obtained from real data we test our
model on large simulated catalogs. These datasets are generated by drawing BBH
parameters from a known population distribution and injecting the
corresponding signals into detector noise realizations. We use Bayesian
parameter estimation on the injected signals using the BILBY package (Ashton
et al., 2019; Romero-Shaw et al., 2020), its DYNESTY-based nested sampler
(Speagle, 2020), and the aligned spin waveform approximant known as
IMPRphenomD (Husa et al., 2016; Khan et al., 2016). We use LVK’s noise power
spectral densities that were measured during the first three months of its
third observing run (Collaboration & orthers, 2022) and draw 276 events for
each simulated population which is roughly 4 times the size of the GWTC-3 BBH
catalog. We use sensitivity injections found during the third observing run as
publicly released by the LVK (Abbott et al., 2023c) to compute the detectable
time-volumes required for analyzing these simulated populations.

###  4.1 Testing for false positives

We first simulate three BBH populations each of which is characterized by a
Powerlaw + Peak model for the marginal primary mass-distribution (Talbot &
Thrane, 2018), a mass-ratio dependent pairing function (Fishbach & Holz,
2020), a power-law in one plus redshift for the evolution of the merger rate
(Fishbach et al., 2018), and a Gaussian distribution for effective inspiral
spins. In the first population, we impose no intrinsic correlations between
masses, spins, and redshifts. In the second one, we introduce an anti-
correlation between mass-ratio and effective inspiral spins by making the mean
of the $\chi_{\rm{eff}}$ Gaussian a linear function of mass-ratio with a
negative slope (Callister et al., 2021). In the third population we allow the
$\chi_{\rm{eff}}$ distribution to broaden with redshift, similar to the models
of Biscoveanu et al. (2022).

All three population distributions can be expressed in the following form:

$\frac{dN}{dm_{1}dqdzd\chi_{\rm{eff}}}=N_{0}\left(\lambda\frac{(1-\alpha)m_{1}^{-\alpha}}{m_{\rm{max}}^{1-\alpha}-m_{\rm{min}}^{1-\alpha}}+(1-\lambda)\frac{e^{-\frac{1}{2}\left(\frac{m_{1}-\mu_{m}}{\sigma_{m}}\right)^{2}}}{\sqrt{2\pi}\sigma_{m}}\right)\frac{(1+\beta)q^{\beta}}{1-\left(\frac{m_{1}}{m_{\rm{min}}}\right)^{\beta+1}}\frac{dV}{dz}(1+z)^{\kappa-1}\frac{e^{-\frac{1}{2}\left(\frac{\chi_{\rm{eff}}-\mu_{\chi}+\delta\mu_{\chi}(q-0.5)}{\sigma_{\chi}+\delta\sigma_{\chi}(z-0.001)}\right)^{2}}}{\sqrt{2\pi}(\sigma_{\chi}+\delta\sigma_{\chi}(z-0.001))}$
(3)

where $q=\frac{m_{2}}{m_{1}}$ is the mass-ratio and $N_{0}$ the total number
of BBHs that have occured through out the observation time. The hyper-
parameters
$(m_{\rm{min}},m_{\rm{max}},\lambda,\alpha,\beta,\mu_{m},\sigma_{m},\kappa,\mu_{\chi},\delta\mu_{\chi},\sigma_{\chi},\delta\sigma_{\chi})$
corresponding to the three simulated populations are listed in Table 1.

The reasons for choosing these three populations to conduct our validation
study are as follows. Using the un-correlated mock dataset, we demonstrate
that the kind of mass-spin correlations we are seeing in real data are not
spuriously manifesting from any artifacts built into our model, nor from the
correlations between the measurement uncertainties of the BBH observables of
interest (left panels of Fig. 2). Using the simulated population with
intrinsic $q\mbox{-}\chi_{\rm{eff}}$ correlation, we demonstrate that such a
correlation cannot manifest into multiple subpopulations in the marginal
$m_{1}\mbox{-}\chi_{\rm{eff}}$ plane (middle panels of Fig. 2). Finally, using
the population with intrinsic spin-redshift correlations we validate that our
results for GWTC-3 are not indicative of spin-redshift correlations but
instead of new astrophysics underlying BBH formation (right panels of Fig. 2).

We note that for the simulated dataset with intrinsic spin-redshift
correlation, the reconstructed $\chi_{\text{eff}}$ distributions are not
expected to match the injected population since our model assumes the joint
$\chi_{\text{eff}}\mbox{-}z$ distribution to be separable, i.e.,
$p(\chi_{\text{eff}},z)=p(\chi_{\text{eff}})p(z)$. The point of this
validation study is to instead show that the ignored spin-redshift
correlations cannot manifest into the kind of mass-spin correlations we are
seeing, in the context of our binned model. On the other hand, for the
simulations with intrinsic spin-mass ratio correlation, we choose a higher
signal to noise ratio (SNR) threshold than the other two to make sure that the
poor measurability of both mass ratio and effective spin at low SNR does not
bias the inferred shapes.

Parameter | Uncorrelated | Spin-mass-ratio | Spin-redshift
---|---|---|---
$m_{\rm{min}}/M_{\odot}$ | 6.0 | 6.0 | 6.0
$m_{\rm{max}}/M_{\odot}$ | 70.0 | 70.0 | 70.0
$\lambda$ | 0.04 | 0.04 | 0.04
$\alpha$ | 4.0 | 4.0 | 4.0
$\beta$ | 1.5 | 1.5 | 1.5
$\mu_{m}/M_{\odot}$ | 35.0 | 35.0 | 35.0
$\sigma_{m}/M_{\odot}$ | 4.0 | 4.0 | 4.0
$\kappa$ | 3.0 | 3.0 | 3.0
$\mu_{\chi}$ | 0.2 | 0.2 | 0.2
$\delta\mu_{\chi}$ | 0.0 | $-0.5$ | 0.0
$\sigma_{\chi}$ | $10^{-1.8}$ | $10^{-1.8}$ | $10^{-0.85}$
$\delta\sigma_{\chi}$ | 0.0 | 0.0 | $10^{-0.88}$
Table 1: True values for the hyper-parameters characterizing the underlying
population in Eq. (3) for the three simulated catalogs. The parameters
controlling the correlations are chosen to have best-fit values from Abbott et
al. (2023b)

###  4.2 Recovering different kinds of mass-spin correlations

In the previous section, we had shown that various forms of
$\chi_{\text{eff}}$ correlations which have already been found in GWTC-3,
cannot manifest into the trends we have found in the same dataset. In this
section we simulate two additional populations with intrinsic mass-spin
correlations, similar to the trends we have found in GWTC-3, and show that our
model can accurately recover injected conditional distributions of mass given
effective spin and vice-versa. For the primary mass distribution, we choose a
mixture of two subpopulations, one of which is a power-law without any peak
and the other a Powerlaw+Peak with a very shallow spectral index and a large
peak fraction. Each of these subpopulations is associated with a Gaussian
$\chi_{\text{eff}}$ distribution whose mean and standard deviations are
different from that of the other subpopulations. The mass-ratio distribution
and redshift evolution of the merger rate are both chosen to be the same as
that of the mock populations in the previous section. These two population
distributions can be expressed in the following form:

$\frac{dN}{dm_{1}dqdzd\chi_{\rm{eff}}}=\sum_{i=1}^{2}f_{i}N_{0,i}\left(\lambda_{i}\frac{(1-\alpha)m_{1}^{-\alpha_{i}}}{m_{\rm{max}}^{1-\alpha_{i}}-m_{\rm{min}}^{1-\alpha_{1}}}+(1-\lambda_{i})\frac{e^{-\frac{1}{2}\left(\frac{m_{1}-\mu_{m,i}}{\sigma_{m,i}}\right)^{2}}}{\sqrt{2\pi}\sigma_{m,i}}\right)\frac{(1+\beta)q^{\beta}}{1-\left(\frac{m_{1}}{m_{\rm{min}}}\right)^{\beta+1}}\frac{dV}{dz}(1+z)^{\kappa-1}\frac{e^{-\frac{1}{2}\left(\frac{\chi_{\rm{eff}}-\mu_{\chi,i}}{\sigma_{\chi,i}}\right)^{2}}}{\sqrt{2\pi}(\sigma_{\chi,i})}$
(4)

where the hyper-parameters corresponding to the two simulated populations are
listed in Table 2 We set the mean of the effective-spin Gaussian associated
with the non-zero peak mass distribution to be 0.5 in one of the simulated
datasets (peak+aligned) and 0 in the other (peak+isotropic, table 2). The
subpopulation without the peak is chosen to have an effective spin
distribution with mean 0.2 for both cases. We analyze both mock catalogs with
our model and show that it can accurately reconstruct the injected
distribution (Figure 3). In other words, for GWTC-3, had the subpopulation
with the $30-40M_{\odot}$ feature been associated with more positively aligned
spin distributions, our model would have recovered it correctly.

Figure 3: Inferred conditional distributions for the two simulated catalogs with intrinsic mass-spin correlations resulting from the mixture of two subpopulations. The top and bottom pannel represrents the peak+aligned and peak+isotropic populations respectively from table 2 Parameter | Peak+aligned | Peak+isotropic
---|---|---
$m_{\rm{min}}/M_{\odot}$ | 8.0 | 8.0
$m_{\rm{max}}/M_{\odot}$ | 70.0 | 70.0
$\lambda_{1}$ | 0 | 0
$\lambda_{2}$ | 0.4 | 0.4
$\alpha_{1}$ | 3 | 3
$\alpha_{2}$ | 0.5 | 0.5
$\beta$ | 1.5 | 1.5
$\mu_{m}/M_{\odot}$ | 35.0 | 35.0
$\sigma_{m}/M_{\odot}$ | 3.0 | 3.0
$\kappa$ | 3.0 | 3.0
$\mu_{\chi,1}$ | 0.2 | 0.2
$\mu_{\chi,2}$ | 0.5 | 0
$\sigma_{\chi,1}$ | $10^{-1.1}$ | $10^{-1.1}$
$\sigma_{\chi,2}$ | $10^{-0.7}$ | $10^{-0.8}$
$f{1}$ | 0.6 | 0.6
$f{2}$ | 0.4 | 0.4
Table 2: True values for the hyper-parameters characterizing the underlying
population in Eq. (4) for the three simulated catalogs. The parameters
controlling the correlations are chosen to have best-fit values from Abbott et
al. (2023b)

##  5 Astrophysical Implications

We have found evidence for a subpopulation of BBH mergers in the mass range
$30-40\ M_{\odot}$ with preference for a symmetric $\chi_{\rm{eff}}$
population distribution. As explored by different studies, the $30-40\
M_{\odot}$ feature seen in the GWTC-3 BBH inferred mass spectrum with models
such as Powerlaw + Peak, cannot be explained by the PISN/PPISN process as this
feature is expected to be above the $50M_{\odot}$ mass range (Hendriks et al.,
2023; Woosley, 2017; Thrane & Talbot, 2019; Golomb et al., 2023). In the rest
of this section we consider different proposed formation channels which could
be responsible for this subpopulation.

Isolated binary evolution formation channels generally rely on Roche-overflow
mass transfer or close initial binary separations to produce BBHs that are
close enough to merge within a Hubble time. In the case of Roche overflow, the
mass transfer can proceed stably in all cases potentially leading to mass
ratio reversals (e.g. Gallegos-Garcia et al., 2021; Broekgaarden et al., 2022;
van Son et al., 2022; Bavera et al., 2023) or unstably leading to common
envelope evolution (e.g. Mapelli, 2020a; Zevin et al., 2021; Wong et al.,
2021; Belczynski et al., 2022; van Son et al., 2022). In the case of close
initial binary separations, the stellar components are expected to rotate with
periods synchronized to the orbital period such that the stars evolve in a
chemically homogeneous fashion and thus do not experience significant radius
increases over their lifetimes (e.g. Mandel & de Mink, 2016; Marchant et al.,
2016; de Mink & Mandel, 2016). While common envelope evolution is primarily
expected to produce BBHs with masses $<30M_{\odot}$, both stable mass transfer
and chemically homogeneous evolution could produce BBHS in this mass range.
However, due to the coevolution of the stars in the binary, our results
suggest that most of the merger rate contribution from either of these
channels are likely not significantly contributing to the $30-40\ M_{\odot}$
merger rate due to the preference for a symmetric $\chi_{\rm{eff}}$ population
distribution (e.g, Bavera et al., 2020; Briel et al., 2023).

Dynamically formed BBHs can originate from a wide variety of astrophysical
environments from triple star systems (e.g., Antonini et al., 2017; Martinez
et al., 2022) to young stellar clusters (e.g., Di Carlo et al., 2019, 2020),
globular clusters (e.g., Rodriguez et al., 2018; Mapelli et al., 2022), and
nuclear star clusters (e.g. Antonini & Perets, 2012; Petrovich & Antonini,
2017), as well as the disks of active galactic nuclei (AGN) (e.g. Bartos et
al., 2017; Stone et al., 2017; McKernan et al., 2020). While each of these
environments are capable of producing masses near $30-40\ M_{\odot}$, only
globular clusters are predicted to produce BBHs with a mass distributions that
peaks near $30-40\ M_{\odot}$ _and_ a $\chi_{\rm{eff}}$ distribution that is
symmetric about 0. In the case of nuclear star clusters the mass peaks toward
lower values, while in the case of young stellar clusters the mass peaks
toward higher values (Cheng et al., 2023; Mapelli et al., 2022). In the case
of BBH formation in AGN disks, the mass and spin distributions are heavily
influenced by the existence (or lack) of migration traps which predominantly
produce low-mass-ratio mergers (e.g., Bellovary et al., 2016; McKernan et al.,
2020).

We can isolate the $30-40\ M_{\odot}$ subpopulation and compute the measured
local merger rate $R_{0}^{\rm{peak}}$ (at $z=0$) under various assumptions.
First, we compute the contribution by marginalizing over the remaining mass
bins and across the entire $\chi_{\rm{eff}}$ range. Second, we can compute
this contribution $R_{0}^{\rm{low}\mathrm{\chi}}$ over the range
$0.05<\chi_{\rm{eff}}<0.1$ only and also provide an estimate for the local
merger rate $R_{0}^{\rm{other}}$ in the complementary mass and
$\chi_{\rm{eff}}$ bins, a proxy, for the merger rate due to field binaries and
other formation channels which we neglect. In the first case we find,
$R_{0}^{peak}=2.04_{-0.73}^{+0.79}\ \rm{Gpc^{-3}yr^{-1}}$ with the merger rate
of other binaries to be $R_{0}^{\rm{other}}=22.4_{-9.4}^{+10.2}\
\rm{Gpc^{-3}yr^{-1}}$. Similarly, for the second more restrictive case we
consider, we find, $R_{0}^{\rm{low}\mathrm{\chi}}=1.09^{+0.81}_{-0.95}\
\rm{Gpc^{-3}yr^{-1}}$, while for the remaining binaries we find,
$R_{0}^{\rm{other}}=23.27^{+10.51}_{-9.50}$. We measure a total merger rate of
$R_{0}=24.4_{-9.2}^{+10.7}\ \rm{Gpc^{-3}yr^{-1}}$ for reference. From the
estimated local merger rates we estimate the percentage of binaries
contributing to the $30-40\ M_{\odot}$ feature to be in the range of $1-8\%$.

Under the assumption that the $30-40\ M_{\odot}$ and
$0.05<\chi_{\rm{eff}}<0.1$ regions contain dominant contributions from the
dynamical channel, our results indicate that we find merger rates in the range
$R_{0}^{\rm{dyn}}=0.1-2.8\ \rm{Gpc^{-3}yr^{-1}}$ for these subpopulations. We
also find merger rates in the range $R_{0}^{\rm{other}}=13-34\
\rm{Gpc^{-3}yr^{-1}}$ for the remaining contribution due to other formation
channels. If we assume that this subpopulation is due entirely to mergers
occuring in GCs, we expect to have merger rates in the range $1-4\
\rm{Gpc^{-3}yr^{-1}}$ and hence find consistency with theoretical models
(Mapelli et al., 2022; Rodriguez et al., 2018). However, we find a preference
towards the lower and more pessimistic end of the rate predictions. We note
that our measured rates are also broadly consistent with other dynamical
formation environments such as young star clusters and nuclear star clusters.
However, the mass feature being at $30-40\ M_{\odot}$, makes GCs most
consistent with our findings (Mapelli et al., 2022).

##  6 Discussion and Conclusion

In this work, for the first time, we have characterized the joint distribution
of BBH component masses and effective inspiral spins using a highly data-
driven BGP-based population model. Using GWTC-3 data, we have found hints of
multiple sub-populations in the BBH mass spectrum that are associated with
different ranges of effective inspiral spins. The high (more positive) spin
sub-population is found to have a powerlaw-like shape with no feature at the
$30-40M_{\odot}$ range, which is consistent with the isolated binary formation
channel and with the predictions of stellar evolution models that place the
PPISN feature at $50M_{\odot}$ or higher (Hendriks et al., 2023). The low
(closer to zero) spin sub-population demonstrates a feature at the
$30-40M_{\odot}$ range and is consistent with the expectations for BBHs formed
dynamically in globular clusters, which are thought to peak around
$30M_{\odot}$. These implications are corroborated by the conditional
effective spin distributions that we have inferred. We have found that BBHs in
the $30-40M_{\odot}$ range have an effecetive spin distribution more symmetric
about zero as compared to that of other BBHs which correspond to a more
positively skewed effective spin distribution. We have also computed the
combined merger rates in the $30-40M_{\odot}$ and $0.05<\chi_{\rm{eff}}<0.1$
ranges and have found that they are consistent with the theoretical
predictions for dynamically formed BBHs in globular clusters.

Our conclusions align broadly with the findings of Godfrey et al. (2023) who
employ flexible mixture models to search for BBH subpopulations in GWTC-3
data, except a few subtle distinctions. We both identify two subpopulations in
the mass spectrum one of which peaks near $10M_{\odot}$ and falls off steeply
afterward, and the other one with a feature near the $30-40M_{\odot}$ range.
However, Godfrey et al. (2023) find that the spin distributions of these two
subpopulations are mostly consistent with each other, particularly the
effective-spin distributions which they find to be almost completely
overlapping. On the other hand, we find that the subpopulation that is
comprised mostly of BBHs with at least one component in the $30-40M_{\odot}$
range has a significantly different effective spin distribution from the
complementary subpopulation of binaries and hence arrive at slightly different
astrophysical interpretations.

We have verified that our results are not susceptible to the inflexibilities
that remain in our highly data-driven model. We have shown that our findings
from GWTC-3 are robust against variations of binning choices as well as that
of the redshift evolution parameter $\kappa$. Using large simulated catalogs
corresponding to three different populations of BBHs, we have shown that the
mass-spin correlations we are seeing cannot manifest spuriously due to
artifacts built into our model assumptions, nor from any other kinds of
correlations that are known to exist in the astrophysical BBH population and
yet have been ignored in our simplified model. Additionally, using two
simulated populations of BBHs with intrinsic mass-spin correlations, we have
shown that if the subpopulation contributing dominantly to the
$30-40M_{\odot}$ range had instead been associated with a more positively
skewed effective spin distribution, our model would have accurately identified
that trend.

At the end of LVK’s ongoing fourth observing run, we expect our model to yield
tighter constraints on the nature and existence of subpopulations in the BBH
mass-spin distributions. As the size of GW catalogs continue to grow, we aim
to generalize our model by constructing higher-dimensional BGPs that will
self-consistently inform on mass-spin, spin-redshift and mass-redshift
correlations in the BBH population and thereby constrain the underlying
astrophysics of several BBH formation channels, beyond the limitations of
parametric population modeling.

##  7 Acknowledgements

The authors would like to thank Thomas Callister, Vicky Kalogera, and
Soumendra Roy for useful discussions and suggestions. This work was supported
by the National Science Foundation award PHY-2207728. IMH is supported by a
McWilliams postdoctoral fellowship at Carnegie Mellon University. The authors
are grateful for computational resources provided by the LIGO Laboratory and
supported by National Science Foundation Grants PHY-0757058 and PHY-0823459,
PHY-0823459, PHY-1626190, and PHY-1700765. This material is based upon work
supported by NSF’s LIGO Laboratory which is a major facility fully funded by
the National Science Foundation

##  Appendix A Variation of the redshift evolution parameter

In this appendix we demonstrate that our results for GWTC-3 are robust against
variations in the redshift evolution parameter $\kappa$, which we had
originally fixed to be the median value found by Abbott et al. (2023b) using
the Powerlaw+Peak model from the same dataset $(\kappa=2.9)$. Here we compare
our results for the the joint mass-spin distributions for two additional
choices of $\kappa$ which are taken to be the 90% lower $(\kappa=1.1)$ and
upper $(\kappa=4.6)$ bounds on it as reported by Abbott et al. (2023b).

Figure 4: Variation of posterior weights ($w^{\gamma}_{i}$, left) and
detectable time-volumes ($\rm{VT}^{\gamma}_{i}$,right) of Eq. (2) in each bin
with $\kappa$.

In figure 4, we first show how the posterior weights ($w^{\gamma}_{i}$) and
detectable time-volumes ($\rm{VT}^{\gamma}_{i}$) of Eq. (2) required to
compute our population likelihood vary with $\kappa$. In Figures. 5, and 6, we
show that the trends we find in the astrophysical mass and spin distributions
of GWTC-3 BBHs are robust against variations in $\kappa$ with the upper and
lower bounds of its measured values that were reported by Abbott et al.
(2023b). We therefore conclude that our conclusions regarding the mass-spin
sub-populations are not biased by our restrictions on $\kappa$.

Figure 5:  GWTC-3 results for $\kappa=1.1$

Figure 6:  GWTC-3 results for $\kappa=4.4$

##  Appendix B Variation of binning choices

In this appendix, we show that the trends we find in the mass-spin population
of GWTC-3 BBHs are robust against changes in choices of binning, for both
masses and effective inspiral spins. We change both the location and number of
bins which is done in two ways. In the first case we change the mass bins to
while keeping the same $\chi_{\rm{eff}}$ bins as in Sec. 3. In the second case
we change the $\chi_{\rm{eff}}$ bins while keeping the mass bins the same as
they were in Sec. 3.

Figure 7: Posterior distributions of the GP lengthscales along the mass (per
log mass bin center,left) and spin (per spin bin center right) dimensions.

.

In both scenarios, we recover the same trends as in Sec. 3. This is expected
given that the GP hyper-parameters such as the correlation lengths are
themselves inferred from the data, with their priors determined by the
distribution of bin centers (Ray et al., 2023). In Figure 7 we show that the
posterior distributions of the GP length scales along both the mass and spin
dimensions remain consistent for all three choices of binning. This is
indicative of the GP correlating more or less bins over the same range of
parameter space depending on how many bins reside within that range,
ultimately leading to consistency in the inferred shapes of the population
disribution (Ray et al., 2023), as shown in Figures 8 and 9.

Figure 8:  GWTC-3 results for changed mass bins

Figure 9:  GWTC-3 results for changed $\chi_{\rm{eff}}$ bins

We note that increasing the number of bins makes the constraints broader,
which is to be expected given that we are trying to fit a larger number of
hyper-paramaters from the same dataset. Furthermore, the increased resolution
does not reveal any additional features in the dataset which justifies our
choice of binning used in Sec. 3.

##  Appendix C Scalability

Previous implementations of the BGP method had scalability issues which
limited the resolution that could be achieved within a tractable inference run
(Ray et al., 2023). This was mostly due to the cubic time-complexity of GP
engines which in our case previously translated to
$O(\prod_{\theta}N_{bins,\theta}^{3})$, where $N_{bins,\theta}$ are the number
of bins along dimension $\theta$, and was bound to become prohibitive with the
addition of more parameters or an increase in resolution. Specifically this
particular form of complexity arises from the Cholesky decomposition of the
GP’s covariance kernel, which takes the following form:

$\displaystyle K^{\gamma\gamma^{\prime}}(\lambda,\sigma)$ $\displaystyle=$
$\displaystyle\sigma^{2}e^{-\sum_{\theta}\frac{(c^{\gamma}_{\theta}-c_{\theta}^{\gamma^{\prime}})^{2}}{2\lambda_{\theta}^{2}}}$
(C1) $\displaystyle\mathbf{K}$ $\displaystyle=$
$\displaystyle\mathbf{L}\mathbf{L}^{T},$ (C2)

where $L$ is a lower triangular matrix. Given that
$\gamma\in\\{1,2,...,\prod_{\theta}N_{bins,\theta}\\}$, and the cubic
complexity of Cholesky decomposition algorithms, the total computational cost
is expected to scale as $O(\prod_{\theta}N_{bins,\theta}^{3})$.

Fortunately, the structure of the kernel in Eq. (C1) provides a nice work
around to this problem. For example, the matrix in Eq. (C1) can be expressed
as a Kronecker product of smaller matrices (one for each BBH parameter):

$\displaystyle\mathbf{K}$ $\displaystyle=$
$\displaystyle\bigotimes_{\theta}\mathbf{K}_{\theta}$ (C3) $\displaystyle
K_{\theta}^{\gamma\gamma^{\prime}}$ $\displaystyle=$
$\displaystyle\sigma^{\frac{2}{n_{\theta}}}e^{-\frac{(c^{\gamma}_{\theta}-c_{\theta}^{\gamma^{\prime}})^{2}}{2\lambda_{\theta}^{2}}},$
(C4)

where $n_{\theta}$ is the total number of BBH parameters considered in the
population model. Exploiting the fact that the Cholesky decomposition of
Kronecker product is the Kornecker product of Cholesky decompositions, as in:

$\displaystyle\mathbf{L}$ $\displaystyle=$
$\displaystyle\bigotimes_{\theta}\mathbf{L}_{\theta}$ (C5)
$\displaystyle\mathbf{K}_{\theta}$ $\displaystyle=$
$\displaystyle\mathbf{L}_{\theta}\mathbf{L}^{T}_{\theta},$ (C6)

we are able to implement an inference whose cost scales as
$O(\sum_{\theta}N_{bins,\theta}^{3})$. This leads to higher dimensional
generalizations of our models straight forward to implement in a scalable
manner. Similarly, an increase in bin resolution along a particular dimension
can also be implemented tractably, upto the scenario wherein we have around
50+ bins along each dimension, a level of resolution that is likely
unnecessary for contemporary and near-future datasets.

## References

  * Aasi et al. (2015) Aasi, J., et al. 2015, Class. Quant. Grav., 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2020) Abbott, R., et al. 2020, Astrophys. J. Lett., 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abbott et al. (2021) —. 2021, Astrophys. J. Lett., 913, L7, doi: 10.3847/2041-8213/abe949
  * Abbott et al. (2023a) —. 2023a, Phys. Rev. X, 13, 041039, doi: 10.1103/PhysRevX.13.041039
  * Abbott et al. (2023b) —. 2023b, Phys. Rev. X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * Abbott et al. (2023c) —. 2023c, Open data from the third observing run of LIGO, Virgo, KAGRA and GEO.  https://arxiv.org/abs/2302.03676
  * Acernese et al. (2015) Acernese, F., et al. 2015, Class. Quant. Grav., 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Ade et al. (2016) Ade, P. A. R., et al. 2016, Astron. Astrophys., 594, A13, doi: 10.1051/0004-6361/201525830
  * Akutsu et al. (2021) Akutsu, T., et al. 2021, PTEP, 2021, 05A102, doi: 10.1093/ptep/ptab018
  * Antonini et al. (2023) Antonini, F., Gieles, M., Dosopoulou, F., & Chattopadhyay, D. 2023, Mon. Not. Roy. Astron. Soc., 522, 466, doi: 10.1093/mnras/stad972
  * Antonini & Perets (2012) Antonini, F., & Perets, H. B. 2012, Astrophys. J., 757, 27, doi: 10.1088/0004-637X/757/1/27
  * Antonini et al. (2017) Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77, doi: 10.3847/1538-4357/aa6f5e
  * Ashton et al. (2019) Ashton, G., et al. 2019, Astrophys. J. Suppl., 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Baibhav et al. (2023) Baibhav, V., Doctor, Z., & Kalogera, V. 2023, Astrophys. J., 946, 50, doi: 10.3847/1538-4357/acbf4c
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, Astrophys. J., 835, 165, doi: 10.3847/1538-4357/835/2/165
  * Bavera et al. (2020) Bavera, S. S., Fragos, T., Qin, Y., et al. 2020, A&A, 635, A97, doi: 10.1051/0004-6361/201936204
  * Bavera et al. (2023) Bavera, S. S., Fragos, T., Zapartas, E., et al. 2023, Nature Astronomy, 7, 1090, doi: 10.1038/s41550-023-02018-5
  * Belczynski et al. (2022) Belczynski, K., Doctor, Z., Zevin, M., et al. 2022, ApJ, 935, 126, doi: 10.3847/1538-4357/ac8167
  * Bellovary et al. (2016) Bellovary, J. M., Mac Low, M.-M., McKernan, B., & Ford, K. E. S. 2016, ApJ, 819, L17, doi: 10.3847/2041-8205/819/2/L17
  * Benacquista & Downing (2013) Benacquista, M. J., & Downing, J. M. B. 2013, Living Rev. Rel., 16, 4, doi: 10.12942/lrr-2013-4
  * Biscoveanu et al. (2022) Biscoveanu, S., Callister, T. A., Haster, C.-J., et al. 2022, Astrophys. J. Lett., 932, L19, doi: 10.3847/2041-8213/ac71a8
  * Briel et al. (2023) Briel, M. M., Stevance, H. F., & Eldridge, J. J. 2023, MNRAS, 520, 5724, doi: 10.1093/mnras/stad399
  * Broekgaarden et al. (2022) Broekgaarden, F. S., Stevenson, S., & Thrane, E. 2022, Astrophys. J., 938, 45, doi: 10.3847/1538-4357/ac8879
  * Callister & Farr (2023) Callister, T. A., & Farr, W. M. 2023.  https://arxiv.org/abs/2302.07289
  * Callister et al. (2021) Callister, T. A., Haster, C.-J., Ng, K. K. Y., Vitale, S., & Farr, W. M. 2021, Astrophys. J. Lett., 922, L5, doi: 10.3847/2041-8213/ac2ccc
  * Chattopadhyay et al. (2023) Chattopadhyay, D., Stegmann, J., Antonini, F., Barber, J., & Romero-Shaw, I. M. 2023.  https://arxiv.org/abs/2308.10884
  * Cheng et al. (2023) Cheng, A. Q., Zevin, M., & Vitale, S. 2023, Astrophys. J., 955, 127, doi: 10.3847/1538-4357/aced98
  * Collaboration & orthers (2022) Collaboration, T. L. S., & orthers. 2022, LIGO-T2000012-v2: Noise curves used for Simulations in the update of the Observing Scenarios Paper — dcc.ligo.org, https://dcc.ligo.org/LIGO-T2000012/public
  * de Mink & Mandel (2016) de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545, doi: 10.1093/mnras/stw1219
  * Di Carlo et al. (2019) Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947, doi: 10.1093/mnras/stz1453
  * Di Carlo et al. (2020) Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020, MNRAS, 498, 495, doi: 10.1093/mnras/staa2286
  * Edelman et al. (2021) Edelman, B., Doctor, Z., & Farr, B. 2021, Astrophys. J. Lett., 913, L23, doi: 10.3847/2041-8213/abfdb3
  * Edelman et al. (2023) Edelman, B., Farr, B., & Doctor, Z. 2023, Astrophys. J., 946, 16, doi: 10.3847/1538-4357/acb5ed
  * Essick et al. (2022) Essick, R., Farah, A., Galaudage, S., et al. 2022, The Astrophysical Journal, 926, 34, doi: 10.3847/1538-4357/ac3978
  * Essick & Farr (2022) Essick, R., & Farr, W. 2022, Precision Requirements for Monte Carlo Sums within Hierarchical Bayesian Inference.  https://arxiv.org/abs/2204.00461
  * Farr (2019) Farr, W. M. 2019, Research Notes of the AAS, 3, 66, doi: 10.3847/2515-5172/ab1d5f
  * Fishbach & Holz (2017) Fishbach, M., & Holz, D. E. 2017, Astrophys. J. Lett., 851, L25, doi: 10.3847/2041-8213/aa9bf6
  * Fishbach & Holz (2020) —. 2020, Astrophys. J. Lett., 891, L27, doi: 10.3847/2041-8213/ab7247
  * Fishbach et al. (2018) Fishbach, M., Holz, D. E., & Farr, W. M. 2018, The Astrophysical Journal Letters, 863, L41, doi: 10.3847/2041-8213/aad800
  * Gallegos-Garcia et al. (2021) Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, ApJ, 922, 110, doi: 10.3847/1538-4357/ac2610
  * Gerosa & Fishbach (2021) Gerosa, D., & Fishbach, M. 2021, Nature Astron., 5, 749, doi: 10.1038/s41550-021-01398-w
  * Godfrey et al. (2023) Godfrey, J., Edelman, B., & Farr, B. 2023.  https://arxiv.org/abs/2304.01288
  * Golomb et al. (2023) Golomb, J., Isi, M., & Farr, W. 2023.  https://arxiv.org/abs/2312.03973
  * Heinzel et al. (2023) Heinzel, J., Biscoveanu, S., & Vitale, S. 2023.  https://arxiv.org/abs/2312.00993
  * Hendriks et al. (2023) Hendriks, D. D., van Son, L. A. C., Renzo, M., Izzard, R. G., & Farmer, R. 2023.  https://arxiv.org/abs/2309.09339
  * Homan & Gelman (2014) Homan, M. D., & Gelman, A. 2014, J. Mach. Learn. Res., 15, 1593–1623
  * Husa et al. (2016) Husa, S., Khan, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044006, doi: 10.1103/PhysRevD.93.044006
  * Karathanasis et al. (2023) Karathanasis, C., Mukherjee, S., & Mastrogiovanni, S. 2023, Mon. Not. Roy. Astron. Soc., 523, 4539, doi: 10.1093/mnras/stad1373
  * Khan et al. (2016) Khan, S., Husa, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044007, doi: 10.1103/PhysRevD.93.044007
  * Kimball et al. (2021) Kimball, C., et al. 2021, Astrophys. J. Lett., 915, L35, doi: 10.3847/2041-8213/ac0aef
  * Li et al. (2023) Li, Y.-J., Wang, Y.-Z., Tang, S.-P., & Fan, Y.-Z. 2023.  https://arxiv.org/abs/2303.02973
  * Loredo (2004) Loredo, T. J. 2004, in AIP Conference Proceedings (AIP), doi: 10.1063/1.1835214
  * Mandel & de Mink (2016) Mandel, I., & de Mink, S. E. 2016, Mon. Not. Roy. Astron. Soc., 458, 2634, doi: 10.1093/mnras/stw379
  * Mandel & Farmer (2022) Mandel, I., & Farmer, A. 2022, Phys. Rept., 955, 1, doi: 10.1016/j.physrep.2022.01.003
  * Mandel et al. (2017) Mandel, I., Farr, W. M., Colonna, A., et al. 2017, Mon. Not. Roy. Astron. Soc., 465, 3254, doi: 10.1093/mnras/stw2883
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, Mon. Not. Roy. Astron. Soc., 486, 1086, doi: 10.1093/mnras/stz896
  * Mapelli (2020a) Mapelli, M. 2020a, Proc. Int. Sch. Phys. Fermi, 200, 87, doi: 10.3254/ENFI200005
  * Mapelli (2020b) —. 2020b, Front. Astron. Space Sci., 7, 38, doi: 10.3389/fspas.2020.00038
  * Mapelli et al. (2022) Mapelli, M., Bouffanais, Y., Santoliquido, F., Arca Sedda, M., & Artale, M. C. 2022, MNRAS, 511, 5797, doi: 10.1093/mnras/stac422
  * Mapelli et al. (2022) Mapelli, M., Bouffanais, Y., Santoliquido, F., Sedda, M. A., & Artale, M. C. 2022, Mon. Not. Roy. Astron. Soc., 511, 5797, doi: 10.1093/mnras/stac422
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, Astron. Astrophys., 588, A50, doi: 10.1051/0004-6361/201628133
  * Martinez et al. (2022) Martinez, M. A. S., Rodriguez, C. L., & Fragione, G. 2022, ApJ, 937, 78, doi: 10.3847/1538-4357/ac8d55
  * McKernan et al. (2020) McKernan, B., Ford, K. E. S., O’Shaugnessy, R., & Wysocki, D. 2020, MNRAS, 494, 1203, doi: 10.1093/mnras/staa740
  * Mohite (2022) Mohite, S. 2022, PhD thesis, Wisconsin U., Milwaukee
  * Neal (2011) Neal, R. M. 2011, MCMC Using Hamiltonian Dynamics, ed. S. Brooks, A. Gelman, G. Jones, & X.-L. Meng (Chapman and Hall/CRC), doi: 10.1201/b10905
  * Oriol et al. (2023) Oriol, A.-P., Virgile, A., Colin, C., et al. 2023, PeerJ Computer Science, 9, e1516, doi: 10.7717/peerj-cs.1516
  * Paczynski (1976) Paczynski, B. 1976, in Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, Vol. 73, 75
  * Petrovich & Antonini (2017) Petrovich, C., & Antonini, F. 2017, ApJ, 846, 146, doi: 10.3847/1538-4357/aa8628
  * Postnov & Yungelson (2006) Postnov, K., & Yungelson, L. 2006, Living Rev. Rel., 9, 6, doi: 10.12942/lrr-2006-6
  * Ray et al. (2023) Ray, A., Magaña Hernandez, I., Mohite, S., Creighton, J., & Kapadia, S. 2023, Astrophys. J., 957, 37, doi: 10.3847/1538-4357/acf452
  * Rinaldi et al. (2023) Rinaldi, S., Del Pozzo, W., Mapelli, M., Medina, A. L., & Dent, T. 2023.  https://arxiv.org/abs/2310.03074
  * Rodriguez et al. (2018) Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., et al. 2018, Phys. Rev. D, 98, 123005, doi: 10.1103/PhysRevD.98.123005
  * Rodriguez et al. (2022) Rodriguez, C. L., et al. 2022, Astrophys. J. Supp., 258, 22, doi: 10.3847/1538-4365/ac2edf
  * Romero-Shaw et al. (2020) Romero-Shaw, I. M., et al. 2020, Mon. Not. Roy. Astron. Soc., 499, 3295, doi: 10.1093/mnras/staa2850
  * Speagle (2020) Speagle, J. S. 2020, Mon. Not. Roy. Astron. Soc., 493, 3132, doi: 10.1093/mnras/staa278
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946, doi: 10.1093/mnras/stw2260
  * Talbot & Thrane (2018) Talbot, C., & Thrane, E. 2018, The Astrophysical Journal, 856, 173, doi: 10.3847/1538-4357/aab34c
  * Thrane & Talbot (2019) Thrane, E., & Talbot, C. 2019, Publ. Astron. Soc. Austral., 36, e010, doi: 10.1017/pasa.2019.2
  * Tiwari (2021) Tiwari, V. 2021, Class. Quant. Grav., 38, 155007, doi: 10.1088/1361-6382/ac0b54
  * Toubiana et al. (2023) Toubiana, A., Katz, M. L., & Gair, J. R. 2023, Mon. Not. Roy. Astron. Soc., 524, 5844, doi: 10.1093/mnras/stad2215
  * van den Heuvel et al. (2017) van den Heuvel, E. P. J., Portegies Zwart, S. F., & de Mink, S. E. 2017, Mon. Not. Roy. Astron. Soc., 471, 4256, doi: 10.1093/mnras/stx1430
  * van Son et al. (2022) van Son, L. A. C., de Mink, S. E., Chruslinska, M., et al. 2022, doi: 10.3847/1538-4357/acbf51
  * van Son et al. (2022) van Son, L. A. C., de Mink, S. E., Renzo, M., et al. 2022, ApJ, 940, 184, doi: 10.3847/1538-4357/ac9b0a
  * Wen (2003) Wen, L. 2003, Astrophys. J., 598, 419, doi: 10.1086/378794
  * Wong et al. (2021) Wong, K. W. K., Breivik, K., Kremer, K., & Callister, T. 2021, Phys. Rev. D, 103, 083021, doi: 10.1103/PhysRevD.103.083021
  * Woosley (2017) Woosley, S. E. 2017, The Astrophysical Journal, 836, 244, doi: 10.3847/1538-4357/836/2/244
  * Wysocki et al. (2019) Wysocki, D., Lange, J., & O’Shaughnessy, R. 2019, Phys. Rev. D, 100, 043012, doi: 10.1103/PhysRevD.100.043012
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, Astrophys. J., 910, 152, doi: 10.3847/1538-4357/abe40e