# Physical Models for the Astrophysical Population of Black Holes: Application
to the Bump in the Mass Distribution of Gravitational Wave Sources

Jacob Golomb<EMAIL_ADDRESS>Department of Physics, California Institute
of Technology, Pasadena, California 91125, USA  LIGO Laboratory, California
Institute of Technology, Pasadena, California 91125, USA  Maximiliano Isi
<EMAIL_ADDRESS>Center for Computational Astrophysics, Flatiron
Institute, 162 5th Ave, New York, NY 10010, USA  Will M. Farr
<EMAIL_ADDRESS>Center for Computational Astrophysics, Flatiron
Institute, 162 5th Ave, New York, NY 10010, USA  Department of Physics and
Astronomy, Stony Brook University, Stony Brook NY 11794, USA

###### Abstract

Gravitational wave observations of binary black holes have revealed unexpected
structure in the black hole mass distribution. Previous studies employ
physically-motivated phenomenological models and infer the parameters that
control the features of the mass distribution that are allowed in their model,
associating the constraints on those parameters with their physical
motivations a posteriori. In this work, we take an alternative approach in
which we introduce a model parameterizing the underlying stellar and core-
collapse physics and obtaining the remnant black hole distribution as a
derived byproduct. In doing so, we constrain the stellar physics necessary to
explain the astrophysical distribution of black hole properties under a given
model. We apply this to the mapping between stellar core mass and remnant
black hole mass, accounting for mass loss due to the pulsational pair
instability supernova (PPISN) process, previously proposed as an explanation
for the observed excess of black holes at ${\sim}35\,M_{\odot}$. Placing
constraints on the reaction rates necessary to explain the PPISN parameters,
we conclude that the peak observed at ${\sim}35\,M_{\odot}$ is highly unlikely
to be a signature from the PPISN process. Allowing the parameters of the core-
remnant mass relationship to evolve with redshift permits correlated and
physically reasonable changes to features in the mass function. We find that
the current data are consistent with no redshift evolution in the core-remnant
mass relationship, but place only weak constraints on the change of these
parameters. This procedure can be applied to modeling any physical process
underlying the astrophysical distribution.

##  1 Introduction

Observations of gravitational waves from binary-black-hole (BBH) and binary-
neutron-star (BNS) mergers with the LIGO-Virgo-KAGRA detector network have
provided otherwise inaccessible information on the properties of those compact
objects (Abbott et al., 2018; Aasi et al., 2015; Acernese et al., 2015). While
individual events offer a glimpse into the details of a particular black hole
(BH) or neutron star, studying observations collectively on a population-level
allows us to draw inferences about stellar astrophysics, the formation
channels of BHs and neutron stars, the overall rates at which BBH/BNS mergers
occur in the universe, and cosmology (LIGO Scientific Collaboration et al.,
2023a, b, c; LIGO Scientific Collaboration & Virgo Collaboration, 2021), as
well as tests of general relativity (Abbott et al., 2021; Payne et al., 2023).

Models adopted for population inference of BBHs tend to take one of two major
approaches. The first are so-called parametric methods, in which a
phenomenological model is constructed using relatively few parameters, with
these parameters directly controlling well-defined features encoded in the
model. This commonly involves assuming a functional form for the global
structure of the distribution (e.g., a truncated power law for the mass
distribution), enhanced by features such as a bump, dip, or break, and jointly
fitting for the properties of the global structure and additional features
(Talbot & Thrane, 2018; LIGO Scientific Collaboration & Virgo Collaboration,
2021; LIGO Scientific Collaboration et al., 2023a; Kovetz et al., 2017;
Fishbach & Holz, 2017). The other popular approach consists of data-driven
methods (sometimes called “non-parametric” methods, in spite of their
overabundance of parameters) in which a flexible model is allowed to fit
nearly-arbitrary shapes to the distribution. Such fits have been achieved with
tools such as splines, histograms, and Gaussian process regression (Callister
& Farr, 2023; Golomb & Talbot, 2023; Mandel et al., 2017; Ray et al., 2023;
Edelman et al., 2022; Tiwari & Fairhurst, 2021).

While the latter method is more general and can capture features not
explicitly defined in a model, the former offers the ability to encode
signatures from expected physical processes in parameters controlling the
shape of the distribution, making it possible to interpret the constraints in
terms of the underlying physical motivations. Nevertheless, interpreting the
constraints on these parameters as constraints on the underlying physics can
be difficult, as there are often unmodeled assumptions as to how the
underlying physical processes translate into resulting distribution that is
being modeled.

In this work, we provide an alternative approach to prescribing a parametric
population model: instead of modeling the BBH distribution directly, we
introduce parameters to describe the underlying astrophysics and the
associated mapping to remnant BH properties and we derive the resulting
population distribution as a result of these underlying parameters and its
associated mapping. Subject to model assumptions, we more directly infer the
physics as informed by BBH properties in a way that avoids strong
phenomenological approximations in the BBH population model itself. This is
related to the “backpropagation” approach in Wong et al. (2023) and the
related method in Andrews et al. (2021).

Previous parametric approaches to modeling the BBH mass distribution have been
particularly useful to place constraints with relatively little data when
strong assumptions about the structure of the mass distribution are warranted.
For example, parametric population analyses of the first catalogs of BBH
events have revealed that the mass distribution is well-described by a
truncated power law that peaks at ${\sim}8\,M_{\odot}$, decaying to high
masses, and featuring an overdensity (modeled as a Gaussian bump) at
${\sim}\,35\,M_{\odot}$ (LIGO Scientific Collaboration et al., 2023a; LIGO
Scientific Collaboration & Virgo Collaboration, 2021).

We apply our novel approach to explain this overdensity in the mass
distribution at ${\sim}35\,M_{\odot}$, which may or may not be accompanied by
a subsequent dip (Talbot & Thrane, 2018; Edelman et al., 2022). The original
motivation for allowing for this feature in the mass distribution was to
capture the expected “pile-up” of BHs that resulted from progenitors that had
undergone pair instability pulsations (Talbot & Thrane, 2018). This pile-up
results from the remnants of stars with zero-age main-sequence (ZAMS) masses
between ${\sim}60M_{\odot}$ through ${\sim}\,140M_{\odot}$ (Rahman et al.,
2022; Woosley & Heger, 2021; Woosley, 2017). As stars in this mass range
evolve past the helium burning phase, their cores reach temperatures hot
enough to produce electron-positron pairs, reducing radiation pressure and
softening the star’s equation of state, which causes the star to contract
(Fowler & Hoyle, 1964; Rakavy & Shaviv, 1967). This contraction can
explosively ignite oxygen in the core, driving an outward pulse that removes
varying amounts of mass from the star. The star then returns to a contracting
phase and the cycle can repeat itself until enough mass has been shed that it
evolves toward normal core collapse. For stars with ZAMS masses above this
range, the contraction-driven core ignition drives a shock that completely
disrupts the star, leaving behind no remnant. Since the pulsational pair
instability process produces a small range of remnant BH masses from stars
from a wide range of ZAMS masses, it is expected that the mass distribution
will exhibit the bump due to this pileup (sometimes referred to as the
“pulsational-pair-instability supernova (PPISN) graveyard”) followed by a
suppression of sources, known as the upper mass gap (Woosley, 2017, 2019;
Woosley & Heger, 2021; Farmer et al., 2019).

Simulations of stellar evolution (e.g., Farmer et al., 2019; Mehta et al.,
2022; Farag et al., 2022) have explored the relationship between initial
stellar mass (in particular the mass of the Carbon-Oxygen (CO) core, $M_{\rm
CO}$) and the final BH mass ($M_{\rm BH}$) after core collapse. They have also
quantified the dependence of the location of the lower edge of this mass gap
and its associated mass range on other physical parameters such as reaction
rates, metallicity, and details of assumed neutrino physics. Previous studies
have used this relationship to place constraints on the astrophysical
properties of the pulsational pair instability process, assuming the most
massive sources observed through LIGO are below the upper mass gap (Farmer et
al., 2020; Mehta et al., 2022; Farag et al., 2022; Stevenson et al., 2019).
Baxter et al. (2021) instead infers the population of BHs coming from the
first-generation (1G) subpopulation below the upper mass gap along with the
subpopulation of higher-generation (2G+) BHs (i.e., BHs that are themselves
the product of past mergers) whose masses can lie within the upper mass gap.

In this work, we implement a simple model for the initial mass function of
stellar CO cores and the associated map from CO core mass to remnant mass,
motivated by the relationship found from the simulations in Farmer et al.
(2019, 2020). Rather than informing our model with individual massive sources,
we construct a full population model for the mass distribution, including a
subpopulation in the upper mass gap due to higher-generation mergers. Using
data from gravitational wave events from the third Gravitational Wave
Transient Catalog (GWTC-3, LIGO Scientific Collaboration et al., 2023c), we
infer the shape of the initial mass function of CO cores, the associated
mapping to the remnant BH distribution, and the relative contribution of
sources formed through 1G mergers.

Having a physically-motivated model for the BH mass distribution facilitates
extensions that incorporate richer physics. As an example, here we allow the
underlying physics to evolve with redshift, as may be expected from cosmic
history considerations. Such evolution in the underlying physical parameters
captures the correlated changes in shape and height of the bump that must
occur in the presence of changing progenitor to remnant mass relationships. It
may be possible to use shape measurements to calibrate changes in the mass
scale of the bump with redshift to reduce or eliminate systematic
uncertianities in cosmological parameter inference from the BBH mass function,
sometimes called the “spectral siren” method (Farr et al., 2019; Ezquiaga &
Holz, 2022).

While here we apply this approach to modeling the PPISN process underlying the
astrophysical BBH mass distribution, it can more generally be used as a model
to place constraints on any relationship between progenitor mass and remnant
BH mass as informed by gravitational wave observations. The model introduced
here can readily be applied to any process with accelerating mass loss as a
function of progenitor mass, but this method can be useful for inferring the
physics of any arbitrary relationship underlying an observable distribution
associated with BBHs (cf. Fishbach & van Son, 2023, for a related approach
applied to inferring the delay times between binary formation and merger).

We begin with an overview of hierarchical Bayesian inference in Sec. 2. We
then outline the models with and without evolution with redshift in Sec. 3. In
Sec. 4, we present results for both model configurations, using data from the
third oberseving run (O3) of LIGO-Virgo. We offer interpretations of our
results in Sec. 5 and provide concluding remarks in Sec. 6.

##  2 Hierarchical Bayesian Inference

We conduct our inference on the population parameters $\Lambda$ with a
hierarchical Bayesian framework, in which we inform our population model with
a catalog $N_{\rm det}$ events, to compute the likelihood (see, e.g., Mandel
et al. (2019); Thrane & Talbot (2019)):

$\mathcal{L}(\\{d\\}|\Lambda)\propto\frac{K(\Lambda)^{N_{d}}e^{-K(\Lambda)}}{p_{\rm
det}(\Lambda)^{N_{d}}}\prod_{i=1}^{N_{\rm
det}}\int\mathcal{L}(d_{i}|\theta)\pi(\theta|\Lambda)d\theta$ (1)

where $\mathcal{L}(d_{i}|\theta)$ is the likelihood of the data for the $i$th
event, given physical parameters $\theta$ (i.e., masses, distances), and
$\pi(\theta|\Lambda)$ is our population model with a predicted number of
detections $K$. The $p_{\rm det}(\Lambda)$ prefactor accounts for the
selection effects associated with observing a catalog biased toward sources
with parameters that favor detectability (i.e., the Malmquist bias). See
Appendix A for details on this likelihood.

Following the approach in Farr (2019) and Tiwari (2018), we compute $p_{\rm
det}(\Lambda)$ with injections of sources from a fiducial population in
detector noise, and assigning weights to each of the sources that pass our
detection threshold. These sensitivity injections are from the O3 injection
set released in LIGO Scientific Collaboration et al. (2023d). We compute the
per-event population evidence (the integral in Eq. (1)) by reweighting samples
from individual event posterior distributions and dividing by the event-
specific sampling priors. Since our population model is written only in terms
of masses and distances, we effectively adopt the prior from parameter
estimation for the spin parameters (isotropic in direction and uniform in spin
magnitude).

For our analyses involving third observing run (O3) data, we obtain posterior
samples for each event from LIGO Scientific Collaboration et al. (2021), using
the same BBH events from O3 as in LIGO Scientific Collaboration et al.
(2023a). This results in 59 events meeting the False Alarm Rate (FAR)
threshold of 1 per year. Throughout this work we assume the best-fit
cosmological parameters from the Planck 2018 release (Planck Collaboration,
2020).

We sample the population posterior using the No-U-Turn-Sampler (NUTS) in
Numpyro (Phan et al., 2019; Bingham et al., 2019), and we write the functions
for computing Eq. (1) in jax (Bradbury et al., 2018) to take advantage of
automatic differentiation when sampling with Hamiltonian Monte Carlo (Duane et
al., 1987).

We do not enforce the convergence conditions from LIGO Scientific
Collaboration et al. (2023a) for the Monte Carlo integrals in our likelihood,
but we confirm that all points in our posterior have a reasonable enough
number of effective per event posterior samples and injections for good
convergence.

##  3 Mass Distribution of Black Holes From CO Core IMF

###  3.1 Mass Distribution Model

We begin to construct our mass distribution model by assuming a functional
form for the initial mass function (IMF) of stellar CO cores. Surveys have
shown that the stellar IMF on the main sequence can be well-modeled as a
featureless power law at high masses, with a power law index of approximately
$-2.3$ (see, e.g., Kroupa, 2001; Salpeter, 1955; Kroupa & Jerabkova, 2019).
Recent studies have shown through simulations and analytic approximations that
the there may be an approximately linear relationship between a high-mass
star’s zero-age main-sequence (ZAMS) mass and the mass of its core before
undergoing supernova, although this is uncertain (Woosley, 2019; Belczynski et
al., 2016). We therefore assume that the IMF of CO cores can also be modeled
with a power law to good approximation, but allow for a break at
$20\,M_{\odot}$ for additional flexibility (see, e.g., Schneider et al.,
2018). Even if this relationship has nontrivial nonlinearities, modelling the
shape of the broken power law should capture the dominant resolvable structure
of the distribution.

We express the distribution CO core masses, $M_{\rm CO}$, as

$\frac{dN}{dM_{\rm CO}}(M_{\rm CO})\propto\begin{cases}\left(\frac{M_{\rm
CO}}{20M_{\odot}}\right)^{-a}&\text{if }M_{\rm CO}<20M_{\odot}\,,\\\
\left(\frac{M_{\rm CO}}{20M_{\odot}}\right)^{-b}&\text{if }M_{\rm
CO}>20M_{\odot}\,.\end{cases}$ (2)

In order to obtain the resulting BH mass distribution from the CO core mass
distribution, we require a mapping between $M_{\rm CO}$ for a star and the
mass of its remnant after undergoing core collapse. From simulations with the
stellar evolution code MESA (Paxton et al., 2019), Farmer et al. (2019) find
that $M_{\rm CO}$ prior to core collapse is in fact the dominant variable
determining the remnant mass post core collapse. Figure 4 in Farmer et al.
(2019) shows the resulting $M_{\rm BH}$ vs $M_{\rm CO}$ relationships obtained
for a range of choices of input physics and metallicity. The authors note that
for a given choice of metallicity, this relationship is well-modeled by a
piecewise map: a linear relationship, turning over to a quadratic at the CO
core mass at which pulsations begin to remove notable mass, followed by a
decay to $M_{\rm BH}=0\,M_{\odot}$, corresponding to the mass at which
pulsational pair instability fully disrupts the star, leaving no remnant. This
general trend has been confirmed by other simulation-based studies (e.g.,
Mehta et al., 2022; Woosley, 2017).

We express this piecewise mapping through a functional form $\bar{M}_{\rm
BH}(M_{\rm CO}|M_{\rm PISN},M_{\rm BHmax})$ given by

$\bar{M}_{\rm BH}(M_{\rm CO}|M_{\rm PISN},M_{\rm BHmax})=\begin{cases}M_{\rm
CO}&\text{if }M_{\rm CO}<M_{\rm PPISN}\,,\\\ M_{\rm BHmax}+\frac{(M_{\rm
CO}-2M_{\rm BHmax}+M_{\rm PPISN})^{2}}{4(M_{\rm PPISN}-M_{\rm
BHmax})}&\text{if }M_{\rm PPISN}<M_{\rm CO}<2M_{\rm BHmax}-M_{\rm PPISN}\,,\\\
0&\text{otherwise}.\end{cases}$ (3)

This captures the linear relationship at lower masses, for which the relative
amount of mass loss in core collapse is independent of $M_{\rm CO}$; it then
transitions to a quadratic function at the onset of pulsations, $M_{\rm
PPISN}$, with the range $M_{\rm PPISN}$ to $M_{\rm BH,max}$ approximately
representing the range of final BH masses that “pile up” to form a feature in
the BH mass distribution. This pile-up is usually expressed as a relative
overdensity through a Gaussian bump (Talbot & Thrane, 2018; LIGO Scientific
Collaboration & Virgo Collaboration, 2021; LIGO Scientific Collaboration et
al., 2023a).

When $M_{\rm CO}$ reaches $M_{\rm PPISN}$, the $M_{\rm CO}$ to $M_{\rm BH}$
mapping transitions from its linear relationship to a nonlinear one as the
pulsation process causes mass loss whose efficiency increases with the star’s
mass (Marchant et al., 2019; Farmer et al., 2019; Woosley, 2017; Woosley &
Heger, 2021). The form of the quadratic function in Eq. (3) puts the peak
value at $M_{\rm BH,max}$ and enforces that the transition and its derivative
be continuous at $M_{\rm PPISN}$. This results in a BH mass distribution in
which remnants between roughly $M_{\rm PPISN}$ and $M_{\rm BH,max}$ can map
back to a wider range of progenitor masses, and each BH mass bin $dM_{\rm BH}$
in this range contains more systems than it otherwise would had the $M_{\rm
CO}-M_{\rm BH}$ relationship continued to be linear.

As shown in Farmer et al. (2019), the map from $M_{\rm CO}$ to $M_{\rm BH}$ is
sensitive to unknown physics affecting the core collapse and stellar evolution
process, even given a fixed $M_{\rm PPISN}$. We therefore do not know with
certainty of a one-to-one map of $M_{\rm CO}$ to $M_{\rm BH}$ and instead
treat this as a probabilistic model, effectively capturing the uncertainty in
physical parameters ignored in our mapping. Even if we knew global physical
parameters for the core collapse process (e.g., reaction rates) with
certainty, a given $M_{\rm CO}$ will always have a range of possible
associated remnant masses due to factors such as the unknown metallicity at
formation (see Sec. 3.4). Treating the $M_{\rm CO}$ to $M_{\rm BH}$ map
probabilistically allows us to capture some of this uncertainty, for example,
if sources from multiple metallicities are present in our dataset, spreading
out this mapping. We simulate such uncertainty in the $M_{\rm CO}-M_{\rm BH}$
mapping by treating the natural logarithm of the remnant mass as a realization
from a Gaussian distribution, with standard deviation $\sigma$, namely

$p\left(\ln(M_{\rm BH})\mid\bar{M}_{\rm
BH},\sigma\right)=\mathcal{N}\left(\ln[\bar{M}_{\rm BH}(M_{\rm
CO})],\sigma\right),$ (4)

where $\bar{M}_{\rm BH}(M_{CO})$ is given in Eq. (3). Since Eq. (4) specifies
that the logarithm of $M_{\rm BH}$ values are normally distributed around
$M_{\rm CO}$ with standard deviation $\sigma$, the uncertainty on the physical
value of the mass $M_{\rm BH}$ will grow with $M_{\rm CO}$ for fixed $\sigma$.
This is intended to be consistent with the trend in the spread around the
$M_{\rm CO}-M_{\rm BH}$ map seen in Fig. 4 of Farmer et al. (2019), in which
higher values of $M_{\rm CO}$ can create a larger range of remnant $M_{\rm
BH}$ values.

Any confident measurement of a nonzero value of $\sigma$ would mean there is
variation in the $M_{\rm CO}-M_{\rm BH}$ mapping. This could originate in a
number of factors, including physical properties affecting stellar evolution
manifesting differently between events in the dataset. For example, since
metallicity is expected to have a slight effect on the remnant mass given an
initial CO core mass (Farmer et al., 2019), resolvable contributions from
sources with differing birth metallicities in our dataset would result in a
preferentially nonzero value for $\sigma$.

Finally, to get the resulting mass distribution for first generation BHs
${dN}/{dM_{1G}}$, we can simply integrate over progenitor CO core masses,

$\begin{split}\frac{dN}{dM_{\rm 1G}}&\left(M_{\rm BH}\mid a,b,M_{\rm
PPISN},M_{\rm BH,max},\sigma\right)=\\\ &\int dM_{\rm CO}\frac{dN}{dM_{\rm
CO}}\,p(M_{\rm BH}\mid M_{\rm CO}),\end{split}$ (5)

where we implicitly include the Jacobian from computing $p(M_{BH}|M_{CO})$ in
logarithmic mass based on Eq. (4).

Figure 1: Obtaining a 1G BH mass distribution from the CO core IMF. The CO IMF
(bottom panel) gets weighted by the $M_{\rm CO}$ to $M_{\rm BH}$ mapping (top
right panel), resulting in a distribution for ${dN}/{dM_{1G}}$ of BHs in 1G
systems. We label parameters underlying the 1G BBH mass distribution: $a$ and
$b$ are the low-mass and high-mass spectral indices of the CO IMF; $M_{\rm
PPISN}$ and $M_{\rm BH,max}$ (dashed lines) control the start of the nonlinear
mapping and maximum of the quadratic part of the mapping, respectively. The
blue shaded region is the 90% credible region of the lognormal mapping for a
single realization (see Eq 5 and preceding equations for the functional form).
We vary these parameters in our fit to the LIGO-Virgo data, together with the
parameters for the 2G population.

In Fig. 1, we show how the distribution of $M_{\rm BH}$ is derived from an
initial distribution of CO core masses $M_{\rm CO}$ according to our model
with some fiducial values. Each BH in the 1G population is assumed to come
from a progenitor from the $dN/dM_{\rm CO}$ distribution (bottom panel), which
is mapped to a remnant BH mass through the $M_{\rm CO}-M_{\rm BH}$
relationship (upper right panel). Finally, the resulting BH mass distribution
$dN/dM_{\rm BH}$ is obtained by integrating this distribution in the upper
right panel across $M_{\rm CO}$, weighted by $dN/dM_{\rm CO}$ (Eq. 5). This
differs from the procedure in Baxter et al. (2021), as we directly infer the
$dN/dM_{\rm CO}$ and $p(M_{\rm BH}\mid M_{\rm CO})$ distributions, which
uniquely specify $dN/dM_{1G}$ rather than inferring a phenomenological
representation of a resulting $dN/dM_{1G}$ distribution.

###  3.2 Full Mass Distribution

The model outlined in the previous section is only directly applicable to 1G
mergers, in which $M_{\rm BH}$ is the remnant mass from the core collapse of
one of the sources from $dN/dM_{\rm CO}$. Realistically, a catalog of observed
gravitational wave sources should also contain higher-generation merger
events—namely, events that involve BHs who are themselves remnants of previous
BH mergers—although it is commonly assumed that these systems will only
subdominantly contribute to the inferred mass distribution (Miller & Hamilton,
2002; Rodriguez et al., 2019; Kimball et al., 2021; Gerosa & Fishbach, 2021).
As the component masses of these events will be approximately (slightly less
than) the sum of the masses of the BHs from its previous mergers, the masses
in this population can exceed $M_{\rm BH,max}$ and will not follow the same
distribution as the 1G BHs.

The details of the 2G distribution depend on unknown factors that make it
difficult to prescribe a specific functional form (Rodriguez et al., 2019;
Kimball et al., 2021; Doctor et al., 2020). In order to capture these events
in a relatively agnostic manner, we enhance our model with a power law tail
with a spectral index $c$ that smoothly turns on just below $M_{\rm BH,max}$,
and has a height $f_{pl}$ relative to $dN/dM_{1G}$ at $M_{\rm BH,max}$ (see
the bottom right panel in Fig. 2). We express the full mass distribution as:

$\frac{dN}{dm}=\frac{dN}{dM_{1G}}+\delta(m\mid M_{\rm
BH,max})\,f_{pl}\left.\frac{dN}{dM_{1G}}\right|_{M_{\rm
BH,max}}\left(\frac{m}{M_{\rm BH,max}}\right)^{-c}\,,$ (6)

where $dN/dM_{1G}$ is given in Eq. (5) and $\delta(m)$ is an exponential
tapering function that smoothly turns on to $M_{\rm BH,max}$; the parameter
$f_{pl}$ controls the relative height between the peak of the 2G power law and
the 1G mass distribution at $M_{\rm BH,max}$. By adopting this two-component
model, we can prevent 2G sources from biasing the inference of the parameters
of the ${dN}/{dM_{1G}}$ distribution. This assumes that the 2G sources have a
minimal contribution to the mass distribution below ${\sim}M_{\rm BH,max}$,
consistent with the conclusions from, e.g., Fishbach et al. (2022).

We model both component masses as coming from the same mass distribution
${dN}/{dM_{\rm BH}}$ and include a pairing function with power law slope
$\beta$ to get the full mass distribution:

$\frac{dN}{dm_{1}dm_{2}}(m_{1},m_{2})=\left(m_{1}+m_{2}\right)^{\beta}\frac{dN}{dm_{1}}\frac{dN}{dm_{2}}\,,$
(7)

where each $dN/dm_{1/2}$ factor corresponds to a density as in Eq. (6). The
first factor in Eq. (7) constitutes the pairing function, by which the
component masses do not only inform the mass distribution independently but
also by how they pair together to form a total mass (Fishbach & Holz, 2020;
Farah et al., 2023b). The parameter $\beta$ is the exponent on the total mass,
such that positive (negative) values for $\beta$ mean that masses pair up to
preferentially form systems of higher (lower) total mass. We choose this form
of the pairing function, first suggested in Fishbach & Holz (2020), to permit
the possibility of _breaking_ factorization symmetry, so that when $\beta\neq
0$ the joint mass function is not the product of a function of $m_{1}$ and a
function of $m_{2}$; many of the models highlighted in LIGO Scientific
Collaboration et al. (2023a) are forced to be symmetric in this sense.

Fig. 2 shows how the mass distribution $dN/dm$ changes as a function of
various population hyperparameters. The top left panel shows how $M_{\rm
PPISN}$ predictably controls the location of the peak in the mass distribution
by changing where the $M_{\rm CO}-M_{\rm BH}$ mapping becomes nonlinear. The
upper right panel shows how this peak moves to higher $M_{\rm BH}$ and gets
wider as $M_{\rm BH,max}$ increases for fixed $M_{\rm PPISN}$. This is due to
the fact that the difference between $M_{\rm PPISN}$ and $M_{\rm BH,max}$
controls the range of sources that end up in the PPISN pileup, which results
in a less sharp cutoff of the pile-up as $M_{\rm BH,max}$ increases (see, for
example, the top right panel of Fig. 1 for how the quadratic turnover between
$M_{\rm PPISN}$ and $M_{\rm BH,max}$ causes the resulting pile-up in the BH
mass distribution). The bottom left panel of Fig. 2 shows the effect of
varying $\sigma$ on the resulting mass distribution. As $\sigma$ increases,
the remnant masses for a given core mass scatter more broadly, smoothing the
remnant mass distribution and softening the peak. Together, these physical
parameters govern the location, strength, and width of the peak in the BH mass
distribution, as well as the strength of its cutoff. The final panel of Fig. 2
demonstrates the increasing contribution of the high-mass power law tail when
raising $f_{pl}$.

Figure 2: BH mass spectrum following our model in Eq. (7) for different
choices of (clockwise from top left) $M_{\rm PPISN}$, $M_{\rm BH,max}$,
$f_{\rm pl}$, and $\sigma$. For the top left figure, we consider constant
difference between $M_{\rm PPISN}$ and $M_{\rm BH,max}$. Unless being varied,
we assume the following fiducial parameters: $a=2,b=1,c=2.5,M_{\rm
PPISN}=38M_{\odot},M_{\rm BH,max}=41M_{\odot},\sigma=0.05$, and $f_{pl}=0.04$.

###  3.3 Redshift Model

Studies of cosmic star formation history with astronomical surveys show that
the star formation rate increases to a redshift of $z\approx 2$, then smoothly
decays at high redshifts, which is well-modeled by a smoothly broken power law
(Madau & Dickinson, 2014; Vangioni et al., 2015; Ghirlanda et al., 2016). When
convolved with a reasonable delay-time distribution, this also gives rise to a
smoothly broken power law for the merger rate $R(z)$ (Fishbach et al., 2018),
i.e., the number of mergers per comoving volume ($V_{c}$) per time interval in
the source frame ($t_{s}$). Accordingly, we assume a redshift distribution
such that

$\frac{dN}{dV_{c}dt_{s}}(z)\equiv
R(z)\propto\frac{\left(1+z\right)^{\lambda}}{1+\left(\frac{1+z}{1+z_{\rm
peak}}\right)^{\kappa}}\,,$ (8)

where, $\lambda$ controls the low-redshift merger rate, estimated to be
$\lambda\approx 3$ in current LIGO/Virgo studies (LIGO Scientific
Collaboration et al., 2023a); meanwhile, the parameter $z_{\rm peak}$ controls
the redshift at which the merger rate peaks and the slope becomes negative.

Both the peak redshift and the high-redshift merger rate are expected to be
directly informed from detections beyond the horizon of current ground-based
detectors or at redshifts where detections are scarce (Callister & Farr, 2023;
Vitale et al., 2019; LIGO Scientific Collaboration et al., 2023a). However,
combining upper-limits from stochastic gravitational wave searches with
population inference studies can place limits on these parameters with current
detections (Callister et al., 2020). Future observations with 3G detectors
will allow us to significantly constrain the merger rate history across cosmic
time using direct detections of BBH mergers at nearly all relevant redshifts
(Ng et al., 2021; Mancarella et al., 2023).

The merger rate as a function of redshift in the detector frame is expressed
simply in terms of mergers per redshift $z$ per detector-frame time $t_{\rm
det}$:

$\frac{dN}{dzdt_{\rm
det}}(z)=\frac{dN}{dV_{c}\,dt_{s}}\frac{dV_{c}}{dz}\frac{dt_{s}}{dt_{\rm
det}}=R(z)\frac{dV_{c}}{dz}\frac{1}{1+z}\,,$ (9)

where $R(z)$ is as in Eq. (8), and $dV_{c}/dz$ is the differential comoving
volume per redshift bin as determined by cosmology.

###  3.4 Allowing the Mass Spectrum to Evolve with Redshift

Studies of stellar evolution predict that stars formed in lower-metallicity
environments can reach higher remnant BH masses before hitting the PISN
cutoff. This is generally attributed to the ability for metal-rich stellar
winds to carry off significant mass, resulting in lower remnant BH masses
after undergoing pulsations (Van Son et al., 2022; Marchant et al., 2019;
Farmer et al., 2019). While we do not get direct information about the
progenitor metallicities of gravitational-wave sources from the observed data,
we can use known correlations between metallicity and observables in
gravitational-wave data to look for this evolution. Redshift and metallicity
are anticorrelated: stars formed earlier in the universe (i.e., at higher
redshift) are relatively metal-poor when compared to those formed more
recently (at lower redshift), due to the need for the existence of pre-
existing stars to deposit metals into the interstellar medium in order to
birth further generations of stars with higher metallicities (Maiolino et al.,
2008; Belczynski et al., 2016).

Several previous studies have used this trend as motivation to search for
redshift-dependence in the observed BBH mass distribution.111Van Son et al.
(2022) proposes that differing delay time distributions between the high and
low mass portions of the mass distribution may also result in an evolving mass
distribution. Unlike evolution due to birth metalicity, this trend would not
be tracked by the evolution in our model. These studies have typically adopted
phenomenological approaches to modeling this effect, directly encoding
redshift dependence in the location of features in the BBH mass distribution.
Such features include the location of the Gaussian peak and the truncation
point of the mass distribution, allowing these features to vary, for example,
linearly with redshift or with some function of expected metallicity at a
particular redshift (Safarzadeh & Farr, 2019; Fishbach et al., 2018).

We can leverage the physical framework we introduced in Sec. 3.2 to model the
redshift dependence in the mapping from $M_{\rm CO}$ to $M_{\rm BH}$. This
allows us to treat the redshift evolution in the observed BBH mass
distribution as a derived byproduct from an astrophysical process expected to
evolve with redshift, rather than encoding the redshift dependence in the BBH
mass distribution directly.

We express this evolution in term of a linear expansion for the location of
the $M_{\rm PPISN}$ turnover in the mass distribution:

$M_{\rm PPISN}(z)=M_{\rm PPISN}(z=0)+\dot{M}_{\rm
PPISN}\left(1-\frac{1}{1+z}\right),$ (10)

where $M_{\rm PPISN}$ and $\dot{M}_{\rm PPISN}$ are free parameters which we
can interpret as the PPISN turnover location at $z=0$ and the change in this
location over a Hubble time, respectively. This is an indirect model of the
evolution of this feature from high metallicity (late universe, $z=0$) to low-
metallicity (early universe, $z=\infty$) environments. One could alternatively
construct a more explicit model for $M_{\rm PPISN}$ as a function of
metallicity and then metallicity as a function of redshift.

Fig. 3 shows how our model for $dN/dM_{1G}$ appears for different values of
redshift for two choices of $\dot{M}_{\rm PPISN}$. To be consistent with
predictions from stellar models, as described above, we expect a positive
value for $\dot{M}_{\rm PPISN}$ such that the turnover to PPISN occurs at
higher masses at higher redshifts; equivalently, this means we expect the bump
in the mass distribution moves toward higher mass at higher redshifts.

We treat the difference between $M_{\rm PPISN}$ and $M_{\rm BH,max}$ as
constant across redshift, such that $M_{\rm PPISN}$ and $M_{\rm BH,max}$ vary
together with redshift; this could be generalized if needed. As with the model
introduced in the previous section, we include a high-mass power law tail and
a pairing function for the total mass. We use the same $R(z)$ model from Sec.
3.3.

Figure 3: Model for a redshift-dependent mass distribution evaluated at
selected redshifts assuming $\dot{M}_{\rm PISN}=5\,M_{\odot}$ in Eq. (10). All
other parameters are the same as the fiducial set in Fig. 2.

##  4 Results

###  4.1 Non-Evolving Mass Distribution

Adopting the non-evolving mass and redshift models introduced in the previous
section, we infer the corresponding hyperparameters using the hierarchical
Bayesian inference approach described in Sec. 2 and the priors in Table 1.

Parameter | Prior | Description
---|---|---
$a$ | $\mathcal{N}(2.35,2)[-1.65,6.35]$ | Power law index of low-mass CO IMF
$b$ | $\mathcal{N}(1.9,2)[-2.1,5.9]$ | Power law index of high-mass CO IMF
$c$ | $\mathcal{N}(4,2)[0,8]$ | Power law index of 2G high-mass tail
$M_{\rm PPISN}$ | $\mathcal{N}(35,5)[20,50]$ | Mass at start of PPISN process $[M_{\odot}]$
$M_{\rm BH,max}-M_{\rm PPISN}$ | $\mathcal{N}(3,2)[0.5,7]$ | Maximum remnant mass produced by the 1G channel, relative to $M_{\rm PPISN}$ $[M_{\odot}]$
$\sigma$ | $\mathcal{N}(0.1,0.1)[0.05,]$ | Width of lognormal distribution for $M_{\rm CO}$ to $M_{\rm BH}$ mapping
$\beta$ | $\mathcal{N}(0,2)$ | Exponent on total mass pairing function
$\log(f_{pl})$ | $\mathcal{U}[0.01,0.5]$ | Log of relative height between the start of the 2G powerlaw and end of ${dN}/{dM_{1G}}$
$\lambda$ | $\mathcal{N}(2.7,2)[-1.3,6.7]$ | Exponent controlling $R(z)$ at low redshift
$\kappa-\lambda$ | $\mathcal{N}(2.9,2)[1,6.9]$ | Exponent controlling $R(z)$ at high redshift
$z_{\rm peak}$ | $\mathcal{N}(1.9,1)[0,3.9]$ | Redshift at peak $R(z)$
$\dot{M}_{\rm PPISN}$ | $\mathcal{U}[-2,8]$ | Difference in $M_{\rm PPISN}(z)$ over Hubble time $[M_{\odot}]$
Table 1: Priors used in this work. $\mathcal{U}$ is a uniform distribution and
$\mathcal{N}$ is a Gaussian distribution with mean and standard deviation
specified in the parentheses. Numbers in square brackets are upper and lower
bounds of the prior.

We plot draws from the mass distribution posterior in Fig. 4, which shows the
inferred decaying power law shape of the mass distribution and the feature
just below ${\sim}40M_{\odot}$. Turning attention to the parameters that most
directly control the location and strength of the peak in the mass
distribution, Fig. 5 shows the posterior distributions for $\sigma$, $M_{\rm
PPISN}$, and $M_{\rm BH,max}$. The recovered distribution for $\sigma$ shows
notable preference for low values, converging toward the lower bound of the
prior ($\sigma$ = 0.05). This means that the data are consistent with little
to no scatter around the $M_{\rm CO}-M_{\rm BH}$ mapping, while ruling out
high values of $\sigma$ that would over-smoothen the peak in the mass
distribution (cf. Fig. 2, bottom panel). By the same token, the strong support
for low values of $\sigma$ indicates that the data do allow for a relatively
sharp cutoff in the peak; this is such that a suppression of the high end of
the peak need not be compensated by a higher rate in the start of the 2G tail.

With the data supporting a peak at ${\sim}35M_{\odot}$, we measure a
correlation between $M_{\rm PPISN}$ and $M_{\rm BH,max}$, driven by the
constraint that $M_{\rm BH,max}>M_{\rm PPISN}$. We note that the posteriors
for $M_{\rm BH,max}$ and $M_{\rm PPISN}$ are different from the prior,
indicating that the data are informing both the location and the width of the
bump.

Figure 4: Draws from the non-evolving mass distribution posterior, evaluated
at $z=0$.

Taking the model at face value, we infer the location of $M_{\rm BH,max}$ to
be much lower than where stellar nucleosynthesis simulations generally predict
the upper mass gap to begin. For example, Farmer et al. (2019, 2020) finds the
lower edge of the PISN mass gap to range between ${\sim}45{-}50\,M_{\odot}$,
when varying the CO reaction rate within its 1$\sigma$ uncertainty (with
standard deviation $\sigma_{C12}$) with respect to the distribution of
reaction rates given in STARLIB (Sallaska et al., 2013). Using their fit to
the start of the mass gap as a function of $\sigma_{C12}$, and extrapolating
down to our inferred values of $M_{\rm BH,max}$, we infer
$\sigma_{C12}=4.8^{+3.1}_{-2.8}$ at 90% credible levels.222For reaction rates
this high, the fraction of carbon in the core is too low ($X_{\rm C}\ll
10^{-3}$) to be considered a CO core. To be consistent with the assumptions
stated in Farmer et al. (2019), $M_{\rm CO}$ can instead be interpreted as the
mass within the convective zone during helium burning.

Although the simulation coverage is sparse at these masses and therefore these
constraints are largely extrapolation-driven, the anomalous value inferred for
this parameter casts doubt on this PPISN model as a an explanation for the
${\sim}35\,M_{\odot}$ peak in question. See Sec. 5 for further discussion.

In Fig. 6 we plot the inferred inferred $M_{\rm CO}-M_{\rm BH}$ mapping, the
$M_{\rm CO}$ IMF, and 1G $M_{\rm BH}$ mass distribution. We find that the
$M_{\rm CO}$ IMF steeply decreases before the break at $M_{\rm
CO}=20M_{\odot}$ and then becomes shallower or flattens out. The $M_{\rm
CO}-M_{\rm BH}$ mapping is very uncertain, but turns over and creates a peak
at the expected location of $\sim 35M_{\odot}$.

Figure 5: Posterior for selected mass distribution parameters using the model
outlined in Sec. 3.2. Dark- and light-green shaded regions are the $1\sigma$
and $2\sigma$ contours, enclosing 39% and 86% of the probability respectively.
Prior distribution is shaded grey for reference. We find that widening the
prior does not meaningfully increase the posterior support for the PPISN
feature at higher masses. Figure 6: Representation of the $M_{\rm CO}-M_{\rm
BH}$ relationship of the non-evolving $\frac{dN}{dM_{1G}}$ model in Sec. 3.2.
(Bottom panel) Inferred distribution for the initial mass function of $M_{\rm
CO}$ in merging binaries. (Top right panel) 1-$\sigma$ credible region of the
95th percentile (higher blue shaded region) and 5th percentile (lower blue
shaded region) for the inferred $M_{\rm CO}-M_{\rm BH}$ mapping. Dotted line
is a single representative draw from the posterior.

###  4.2 Evolving Mass Distribution

By adopting the more general model from Sec. 3.4 we can relax some of the
assumptions made in the previous sections and now infer the mass distribution
in the presence of an $M_{\rm PPISN}$ that evolves with redshift. In Fig. 7,
we present the posterior probability density on several mass distribution
parameters from this model. Most of the events in the O3 catalog lie at
relatively low redshift, and therefore do not provide good coverage across
redshift scales to inform $\dot{M}_{\rm PPISN}$ meaningfully. Due to these
poor constraints on $\dot{M}_{\rm PPISN}$, we find that the inferred
distributions for $M_{\rm PPISN}$ and $M_{\rm BH,max}$ are consistent with
those obtained when using the non-evolving model. In other words, the feature
at $M_{\rm PPISN}$ is being informed by structure in the data that does not
appear to need to vary with redshift. The resulting mass distribution is
consistent with that from the non-evolving model (plotted below in Fig. 10).

When we extend the prior on $\dot{M}_{\rm PPISN}$ considerably, we find that
we only rule out redshift evolution of the intrinsic mass function for very
large values of $\dot{M}_{\rm PPISN}$. In Fig. 8, we show the posterior
distribution for relevant parameters: constraints on $\dot{M}_{\rm PPISN}$ are
broad, encapsulating a 90%-credible region from $\dot{M}_{\rm PPISN}$ from
$-20\,M_{\odot}$ to $36\,M_{\odot}$. This range is much broader than a
physical prior should allow for these parameters, as metallicity should not
affect $M_{\rm PPISN}$ to this degree (Farmer et al., 2019). While this
demonstrates that data rule out extreme values of $\dot{M}_{\rm PPISN}$, we
cannot currently place constraints within a narrower, more physically-relevant
prior range; future observations may change this. The strong anticorrelation
between $\dot{M}_{\rm PPISN}$ and $M_{\rm PPISN}$ likely indicates that we are
observing the peak from sources in a small range of redshifts, and we
therefore cannot constrain both free parameters in Eq. (10). In order to break
this degeneracy and get constraints on $\dot{M}_{\rm PPISN}$, we would need
additional observations across redshifts. For the remainder of this section,
we present results using the narrower prior for $\dot{M}_{\rm PPISN}$ (in
Table 1).

Figure 7: Posterior for selected mass distribution population parameters from
the evolving mass distribution model in Sec. 3.4 (magenta). One-dimensional
posteriors from the non-evolving mass model overplotted in green for
reference. Figure 8: Inferred mass distribution parameters using the evolving
mass model in Sec. 3.4, but adopting a wider prior on $\dot{M}_{\rm PPISN}$.
Prior distribution is shaded grey for reference. Figure 9: Redshift
distribution parameters inferred with the evolving mass distribution model in
Sec. 3.4. Prior distribution is shaded grey for reference.

Also of note is the similarity between the distribution of $\sigma$ obtained
with this model and that obtained with the non-evolving mass model (see the
comparison in Fig. 7). We discuss the implications of this in Sec. 5.

We present the distribution for redshift parameters, inferred jointly with the
evolving mass distribution, in Fig. 9. The parameter best constrained is
$\lambda$, which controls the evolution of the low-redshift merger rate. We
infer $\lambda=4.3^{+1.6}_{-1.6}$, preferring a merger rate that evolves
steeper than the low-redshift star formation rate ($\lambda\sim 2.7$).
However, the evolution of the merger rate is still consistent with that
implied by the star formation rate along with a short delay-time distribution.
Narrower constraints on this parameter may reveal information on different
formation channels contributing to the observed catalog of BBHs.

Additionally, the posterior distribution for $z_{\rm peak}$ is shifted
slightly to the right of the prior, meaning that we are able to begin to place
very conservative lower limits on $z_{\rm peak}$ due to the lack of a visible
start of a turnover in the inferred $R(z)$ distribution. The lack of support
at the tails towards higher $z_{\rm peak}$ is not due to information gained
from the data, but rather from the prior (see Table 1). These constraints
appear despite having only very little high-redshift information and are
therefore very weak limits.

Given the similarities between the inferred distributions with and without
redshift evolution in the $M_{\rm CO}$ to $M_{\rm BH,max}$ map, we do not
currently obtain improved constraints on physical parameters of interest when
modifying the model in this way. For example, we infer
$\sigma_{C12}=5.2^{+3.4}_{-2.9}$, which is similar to what we reported in Sec.
4.1. Future detections at higher redshift may further inform these aspects of
the model.

##  5 Interpretation of Results

By adopting the model introduced in this work, we can draw conclusions from
the inferred underlying physics represented in our models and explore how the
population model compares to those reported in other works.

Figure 10: Mass spectrum derived from our model (green) with no evolution of
the mass distribution with redshift compared to the Powerlaw + Peak mass
spectrum informed by the same events (blue). We do not include redshift
evolution of the mass distribution in this comparison as the models in LIGO
Scientific Collaboration et al. (2023a) do not include mass-redshift
correlations. Powerlaw + Peak results were obtained using GWPopulation (Talbot
et al., 2019). (Top) Comparison of the primary mass distributions. (Bottom)
Comparison of the mass ratio distributions.

###  5.1 Global Shape of the Mass Distribution

For comparison, we obtain results using the same set of O3 events, adopting
the Powerlaw + Peak mass distribution model, a flat spin magnitude and tilt
model, and a broken power law redshift distribution as implemented in
gwpopulation (Talbot et al., 2019; Talbot & Thrane, 2018; Fishbach et al.,
2018). Qualitatively, we infer a mass distribution (marignalized over $q$)
consistent with the Powerlaw + Peak model, with major features such as the
slope at higher BH mass as well as the bump location showing good agreement in
Fig. 10. This indicates that this overdensity is a confident feature in the
data whose location and prominence is not affected by systematic differences
between these two models. This is reinforced by several other works, which
find that models must include such a feature in order to faithfully capture
the observed mass spectrum (LIGO Scientific Collaboration et al., 2023a; Farah
et al., 2023a; Edelman et al., 2022; Callister & Farr, 2023).

A notable difference is in the merger rate relative to Powerlaw + Peak. While
the 90% credible regions overlap in Fig. 10, their relative heights show that
our model tends to prefer a higher merger rate than what is predicted by
Powerlaw + Peak, particularly at masses above ${\sim}15M_{\odot}$.

In our model, $m_{1}$ and $m_{2}$ both directly inform the physical mass
distribution model (along with the pairing function, see Eq. 7). This is in
contrast to the Powerlaw+Peak model which has separate distributions for
$p(m_{1})$ and $p(q|m_{1})$, such that $m_{2}$ does not directly inform
$dN/dm_{1}$. As demonstrated in Farah et al. (2023b), this makes the mass
distribution feature a peak in the joint $m_{1}-m_{2}$ space, rather than in
the marginal $m_{1}$ distribution. Furthermore, in the marginal mass ratio
distribution (bottom panel of Fig. 10), we see that our model prefers a much
flatter distribution in $q$ than what is preferred by Powerlaw + Peak, which
explicitly models $p(q|m_{1})$ as a power law. This flat mass ratio
distribution is consistent with what is found in Fishbach & Holz (2020) when
adopting a pairing function that is a power law in $q$. Given that we infer
$\beta$ with a preference for negative values, we find that BHs tend to pair
up in binaries that favor lower total masses; this may cause a relative lack
of support at higher masses in the mass distribution which models both
${dN}/{dm_{1}}$ and ${dN}/{dm_{2}}$. We also note that the inferred local rate
$R(z=0)$ is consistent between models. The Powerlaw + Peak model fits the
underlying distribution (i.e., not including the bump) with single power law,
limiting the possible morphologies. We have checked that allowing the
underlying power law in Powerlaw + Peak to include a break does not resolve
the discrepancy.

The distribution for $p(M_{\rm CO})$ we infer (see Fig. 6) disagrees with what
one may expect from a CO IMF resulting from the ZAMS mass IMF assuming a
linear relationship between ZAMS mass and $M_{\rm CO}$. While the distribution
is consistent with a decaying power law for low masses, the distribution
appears to flatten out above our break point of $20\,M_{\odot}$. This trend is
not strongly correlated with the $M_{\rm PPISN}$ and $M_{\rm BH,max}$ we
infer.

Comparing our results to those obtained in Baxter et al. (2021), we find
strong tension with the maximum BH mass in the 1G channel (the start of the
upper mass gap). Motivated by stellar evolution simulations to model the 1G BH
mass distribution with a phenomenological approximation to the shape and
location of a overdensity due to PPISN pileup, Baxter et al. (2021) finds the
PPISN feature and corresponding start of the upper mass gap to be at
${\sim}46M_{\odot}$, in very good agreement with predictions from typical
values of the $\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate. Notably,
Baxter et al. (2021) does not find the feature at ${\sim}35\,M_{\odot}$ we
find and is consistently found in the literature.

###  5.2 Evolution of Mass Distribution with Redshift

Our finding that $\dot{M}_{\rm PPISN}$ is consistent with zero agrees with
other studies that do not find strong preference for evolution of the BBH mass
distribution with redshift. For example, Fishbach et al. (2021) models the
mass distribution as a broken power law where the mass at which the power law
breaks is allowed to vary with redshift. While this is a very different model,
it should qualitatively reproduce some of the features of our model,
particularly at the ${\sim}35\,M_{\odot}$ feature (see Fig. 2). We therefore
expect that if Fishbach et al. (2021) had found strong preferences for an
evolving mass distribution, we would confidently find $\dot{M}_{\rm PPISN}>0$.
We also agree that the data are still consistent with a mass distribution that
has some evolution with redshift, but again we do not have positive evidence
that this is the preferred scenario.

Karathanasis et al. (2023) also looks for evolution of the mass distribution
with redshift. The authors allow the Gaussian bump in a Powerlaw + Peak-like
model to vary with redshift, where the placement of this peak at a given
redshift is determined by the delay time distribution and a jointly-inferred
model for the evolution of (birth) metallicity with redshift. The value they
find for the lower edge of the upper mass gap of ${\sim}44M\,_{\odot}$ is
nominally in better agreement with the prediction from stellar physics models.
However, this value is cited at low metallicity, and they also find there must
be a very strong evolution of this mass scale with metallicity. Extrapolating
their results to the local universe, they find that the upper mass gap at
$z=0$ starts at ${\sim}30\,M_{\odot}$, which is closer to the corresponding
value we obtain for the start of the $M_{\rm CO}{-}M_{\rm BH}$ turnover. This
result seems in tension with theoretical predictions given how small of an
effect metallicity is expected to have on $M_{\rm PPISN}$. There are
unexplained differences in our results, however, as such a strong evolution of
$M_{\rm PPISN}$ with metallicity should mean that we would infer a positive
$\dot{M}_{\rm PPISN}$, assuming delays do not mix events from many different
birth metallicities into similar merger times.

If metallicity evolution effects were causing some of the support for nonzero
values of $\sigma$ obtained in the non-evolving mass model (i.e., from scatter
in the $M_{\rm CO}-M_{\rm BH}$ relation), we would expect $\sigma$ to be
constrained closer to zero with the evolving model, as some of that scatter
would have been absorbed by the redshift evolution. Given that this is not the
case, we conclude that either (1) birth metallicity effects fundamentally have
a subdominant impact on the $M_{\rm CO}{-}M_{\rm BH}$ relationship compared to
other physical parameters that vary between BBH systems, or that (2) the birth
metallicities of the systems in our catalog are not strongly correlated with
the redshifts at which they merge. The latter scenario could result from the
delay time distribution between formation and merger having enough support in
the long-delay tails such that we cannot yet discern a strong correlation
between birth time and merger redshift for systems merging at redshifts of
$z\lesssim 1$.

###  5.3 Physical Interpretations

####  5.3.1 PPISN Process

We can take advantage of the underlying physical parameterization of our model
and offer corresponding interpretations of the implied stellar physics,
assuming $M_{\rm BH,max}$ corresponds to the maximum 1G BH mass as determined
by the PPISN process at a given redshift/metallicity. As found in several
simulation-based studies, the $\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate
is the most important tunable parameter for determining the maximum BH mass
(Farmer et al., 2020; Mehta et al., 2022; Farag et al., 2022). We therefore
consider the $\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate to be the
dominant unknown physical parameter leading to the measured value of $M_{\rm
BH,max}$. Under this assumption, Fig. 11 shows the fit for maximum BH mass as
a function of $\sigma_{C12}$ (each corresponding to a different value of the
$\rm{}^{12}C(\alpha,\gamma)^{16}O$ rate) from the MESA simulations in Farmer
et al. (2020) (reproduced from their data release). While their simulations
only cover the range $-3<\sigma_{C12}<3$, there is a clear trend that
$\sigma_{C12}$ must rise very steeply to reach values of maximum BH mass below
${\sim}45\,M_{\odot}$.333The physical reason for this trend is that during
contraction of the stellar core, hydrostatic equilibrium can be maintained by
convective carbon burning. Higher $\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction
rates lead to cores of lower carbon fractions, $X_{C}$. When the core gets hot
enough to produce electron-positron pairs, the equation of state softens,
leading to a contraction. With little carbon present to provide convective-
driven pressure to stabilize the star, contraction can continue until it
drives thermonuclear ignition of oxygen. This explosive process leads to an
outward-moving shock, removing mass from the star when the shock reaches the
surface with enough velocity. Once this shock breaks through the surface of
the star, contraction begins again. This sequence of pulsations continues
until oxygen in the core is depleted, core elements burn through the normal
pre-SN process, and the star undergoes normal core collapse. If the
${}^{12}C(\alpha,\gamma)^{16}O$ reaction rate is lower, then relatively more
carbon is present and able to burn convectively, counteracting the
contractions in a stable manner. As the carbon fraction gets higher, the star
is able to remain stable against pair-production contractions and stably burn
through the core oxygen (Farmer et al., 2020; Woosley, 2017, 2019)

With the core temperature strongly increasing with stellar mass, there exists
a core mass at which the softening of the equation of state is too extreme to
be resisted by available sources of outward pressure. Since stable outward
pressure support at this stage is largely provided by shell carbon-burning,
the lower carbon fraction, $X_{C}$, makes it now easier for a given
contraction to compress and fully ignite the oxygen core, driving a subsequent
pulsation so powerful that further pair production in the core cannot re-
soften the equation of state fast enough to return it to a contraction phase.
This is basically equivalent to a single pulsation during the pulsational
pair-instability process blowing away the total mass of the star (Woosley,
2017, 2019; Woosley & Heger, 2021). A lower carbon fraction results in this
full disruption of the star (PISN) occurring at lower masses, controlling
where the $M_{\rm CO}-M_{\rm BH}$ map decays to zero after $M_{\rm BH,max}$.

####  5.3.2 Implications for PPISN Physics

Each value of $\sigma_{C12}$ corresponds to a different
$\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate, where $\sigma_{C12}$ is the
number of standard deviations from the median of of the distribution of
reaction rates given in STARLIB (Sallaska et al., 2013), which is adopted from
Kunz et al. (2002) and is the distribution of reaction rates used in the MESA
simulations in Farmer et al. (2020). After fitting for the $\sigma_{C12}$
values necessary to give the our inferred $M_{\rm BH,max}$ (via the top of
Fig. 11), we use the method suggested in Farmer et al. (2020) to approximate
these reaction rates to S-factors for the $\rm{}^{12}C(\alpha,\gamma)^{16}O$
reaction.444While the total reaction rate is given by the cross-section, it is
common in the stellar nuclear physics literature to characterize it in terms
of the so-called “astrophysical S-factor”. This factor can be thought of as
the part of the cross section that is given by the matrix element for the
nuclear reaction itself, ignoring Coulomb repulsion effects between nuclei. At
the (relatively low) typical energies for this reaction in the relevant phase
of stellar evolution (E ${\sim}300\rm keV$), the effects from the Coulomb
repulsion between the positively-charged ${}^{12}C$ and $\alpha$ nuclei
exponentially suppresses the cross-section. It is therefore sometimes
convenient to express the cross-section in terms of the product of two
distinct factors: the S-factor for the reaction, which varies slowly with
energy, and the Gamow factor which characterizes the tunneling probability
through the Coulomb barrier and gives exponential suppression of the cross
section at low energies (see Kippenhahn & Weigert (1994)) We arrive at a value
of the $S$-factor at $300\,{\rm keV}$ of
$S_{300}=932^{+1929}_{-581}\,\mathrm{keV\cdot barn}$.

In Fig. 11, we compare this estimate of $S_{300}$ to that which is reported in
deBoer et al. (2017), in which the authors aggregate a large catalog of
experimental studies related to the $\rm{}^{12}C(\alpha,\gamma)^{16}O$ cross-
section and use numerical methods to extrapolate to an estimated low-energy
$S$-factor. That value is in strong tension with respect to our estimate: our
estimate rules out the value from deBoer et al. (2017) at $>99.9\%$
credibility.

This is in contrast with Farag et al. (2022) and Farmer et al. (2020), which
present inferred values or lower limits of $S_{300}$ which are consistent with
the predicted value. However, these studies draw these conclusions assuming
that $S_{300}$ can be inferred from the most massive BHs in the observed
catalog, assuming that individual high-mass events are 1G BBH mergers. In
contrast, we infer the entire astrophysical mass distribution, with a
subpopulation of 1G BHs that cuts off at a mass determined by the inferred
location and strength of the turnover in the $M_{\rm CO}-M_{\rm BH}$ mapping.
We furthermore allow for high-mass events to belong to a separate 2G BBH
subpopulation that contaminates the 1G upper mass gap, inferring $M_{\rm
BH,max}$ to be where the 1G BBH subpopulation ends, therefore leading to a
lower inferred lower edge of the upper mass gap. Given the confidence with
which the values of $S_{300}$ reported in the nuclear physics literature lie
below what we report, we conclude that one of the assumptions we make when
calculating $S_{300}$ from our inferred mass distribution is incorrect.

Figure 11:  (Top) Relationship between $\sigma_{C12}$ and lower edge of upper
mass gap, reproduced from data release of Farmer et al. (2020). (Center)
Posterior distribution of $\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate, in
terms of standard deviations away from the median reaction rate in STARLIB
(Sallaska et al., 2013); computed by evaluating the fit in the top panel for
the samples of $M_{\rm BH,max}$ in the posterior in Fig 7. (Bottom) Inferred
distribution of $S_{300}$, extrapolated from distribution of $\sigma_{C12}$ as
calculated from the fit in the top panel. Constraints on $S_{300}$ from deBoer
et al. (2017) plotted for comparison, showing tension with the values implied
from our results.

One assumption in this calculation is that the
$\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate is the dominant physical
parameter controlling the location of $M_{\rm BH,max}$. While neglecting the
uncertainty from other nuclear reaction rates certainly causes some sort of
bias when we infer $S_{300}$, such a bias is likely too weak to result in the
tension we see, as multiple studies find that varying the
$\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate has a much stronger effect on
the location of the start of the upper mass gap than other relevant reaction
rates (i.e. 3$\alpha$, ${}^{12}C+^{12}C$, ${}^{12}C+^{16}O$, and
${}^{16}O+^{16}O$) (Farmer et al., 2020; Farag et al., 2022). The caveat here
again is that these simulations only go up to a maximum reaction rate of
$\sigma_{C12}=3$ with respect to the reaction rate distribution in STARLIB
Sallaska et al. (2013), so we cannot confirm the effect of varying these other
reactions rates when $\sigma_{C12}$ is high (let alone the predicted lower
edge of the upper mass gap for these reaction rates, for which we must
extrapolate beyond the simulation space of reaction rates). While it is
possible that one of these other reaction rates can be varied within their
uncertainties to allow us to infer a lower $\sigma_{C12}$, it would have to
change the location of the upper mass gap substantially to agree with the
data; again, studies currently in the literature indicate that these reaction
rates do not have such a strong effect (Farmer et al., 2019; Farag et al.,
2022; Woosley & Heger, 2021).

Another assumption is that the simulations which inspire our model accurately
represent the $M_{\rm CO}-M_{\rm BH}$ mapping for a fixed
$\rm{}^{12}C(\alpha,\gamma)^{16}O$ reaction rate. Mass-loss prescriptions and
temporal resolutions of the simulations (particularly given the high
$\sigma_{C12}$ we wish to probe) may introduce an unknown bias into the
assumptions we make when associating the inferred $M_{\rm CO}-M_{\rm BH}$
mapping with an underlying physical parameter (Mehta et al., 2022; Farag et
al., 2022; Farmer et al., 2019). Again, from the set of simulations presented
in the literature, varying these settings has a subdominant impact on the
location of the upper mass gap.

Finally, the most fundamental assumption in this interpretation is that the
turnover in the $M_{\rm CO}-M_{\rm BH}$ mapping can in fact be associated with
the pair-instability process at all. Studies have suggested that associating
the observed peak in the mass distribution with the PPISN pileup is in tension
with known stellar physics and observed supernovae rates (Hendriks et al.,
2023; Woosley & Heger, 2021). Inferring the underlying $M_{\rm CO}-M_{\rm BH}$
mapping that gives rise to the observed BH mass distribution, this work also
provides evidence of such a tension in terms of the underlying physics that
would be necessary to generate a turnover in the $M_{\rm CO}-M_{\rm BH}$ map
at the correct location.

The cause of the peak in the observed mass distribution at
${\sim}35\,M_{\odot}$ therefore requires alternative explanations (Hendriks et
al., 2023). Recent studies have proposed that this overdensity could be a
signature from a subpopulation of binaries which had undergone stable mass
transfer (Briel et al., 2023), BBH systems in globular clusters (Antonini et
al., 2023), and stars which have undergone significant wind-driven mass loss.
The model introduced in this work could be used to describe or inform any
mechanism that generates a peak in the high-mass tail of the 1G mass
distribution via a transition to a nonlinear $M_{\rm CO}-M_{\rm BH}$
relationship.

###  5.4 Model Limitations

Our model has a few additional caveats. For example, our model does not
attempt to fit for the features that we know exist at lower masses beyond a
power law (Edelman et al., 2022; LIGO Scientific Collaboration et al., 2023a;
Farah et al., 2023a). We have confirmed that neglecting this does not bias the
inference in the higher-mass region that we care about here, and this will be
explored further in future work. We also ignore the effect of spins in our
population, but we demonstrate in Appendix B that this does not cause a
notable bias in our results of interest. Given that certain mass-spin
correlations have been found in the BBH population (Callister et al., 2021),
it may be insightful to use the spins to help distinguish the 1G and 2G
subpopulations (see, e.g., Fishbach et al. (2017); Gerosa & Berti (2017); Farr
et al. (2017).

##  6 Conclusions

Characterizing the population of BBH masses with direct phenomenological or
nonparametric fits can provide insight into the shape of the mass
distribution, but does not provide direct constraints on the underlying
physics of BBH masses. With the method we propose here, we can infer the
underlying physics by fitting the implied (derived) astrophysical BBH
distribution to the observed data. We demonstrate the use of this method by
evaluating the role of the PPISN process giving rise to the 1G BH mass
distribution and its structure, including an exceess (bump) in the mass
distribution at the lower edge of the PPISN mass gap. Fitting this model to
the observed data, we find that the necessary physical parameters to explain
the excess of BHs at ${\sim}35\,M_{\odot}$ are unrealistic from a nuclear
physics perspective if we take this PPISN model at face value. We therefore
conclude it is highly unlikely that the feature at ${\sim}35\,M_{\odot}$ is
associated with the PPISN process.

This framework motivates future investigations to better constrain the physics
underlying astrophysical populations in general. Future work using additional
observations and enhanced versions of our model may be able to constrain
proposed astrophysical mechanisms underpinning the BBH mass, spin, and
redshift distributions. This approach may offer fruitful applications such as
calibration of “spectral siren” features for cosmology (Farr et al., 2019),
investigating other proposed interpretations of the bumps in the mass
distribution, and understanding progenitor populations by relating back to
population synthesis configurations (Wong et al., 2023; Andrews et al., 2021;
Zevin et al., 2017).

## Acknowledgements

We thank Mathieu Renzo for helpful discussions about MESA simulations and the
physics of stellar evolution and the PPISN process. We also acknowledge Alan
Weinstein for useful comments and discussions, and Maya Fishbach for helpful
comments on the manuscript. The Flatiron Institute is a division of the Simons
Foundation. The authors are grateful for computational resources provided by
the LIGO Lab and supported by National Science Foundation Grants PHY-0757058
and PHY-0823459. JG is supported by NSF award No. 2207758. This research has
made use of data or software obtained from the Gravitational Wave Open Science
Center (gwosc.org), a service of the LIGO Scientific Collaboration, the Virgo
Collaboration, and KAGRA.

This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation, as well as 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 construction of Advanced LIGO and construction and operation of the
GEO600 detector. Additional support for Advanced LIGO was provided by the
Australian Research Council. Virgo is funded, through the European
Gravitational Observatory (EGO), by the French Centre National de Recherche
Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN)
and the Dutch Nikhef, with contributions by institutions from Belgium,
Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain.
KAGRA is supported by Ministry of Education, Culture, Sports, Science and
Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan;
National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in
Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC)
in Taiwan.

This manuscript carries LIGO Document Number #P2300418.

##  Appendix A Details of the Likelihood and Differential Rate Calculation

Setting $\theta$ to be the set of single-event parameters, we can write the
contribution from i$th$-event to the population likelihood as (Mandel et al.,
2019):

$p(d_{i}|\Lambda)=\frac{\int
d\theta_{i}\,p(d_{i}|\theta_{i})\,p(\theta_{i}|\Lambda)\,p_{\rm
det}(\theta_{i},d_{i})}{\int\int
dd_{i}\,d\theta_{i}\,p(d_{i}|\theta_{i})\,p(\theta_{i}|\Lambda)\,p_{\rm
det}(\theta_{i},d_{i})}$ (A1)

Recalling that the probability density should be normalized over the arguments
on the left side of the bar, the denominator is included to explicitely
normalize the numerator in terms of the data from the i$th$ detection, and is
commonly known as the “selection effects” term. We write the detection
probability as $p_{\rm det}(\theta_{i},d_{i})$ in order to include the general
possiblity of thresholding detection in terms of the event parameters, which
may be implemented when considering, for example, a simulated catalog. For our
purposes, the detection probability depends on the data, as this is the input
to a detection pipeline when assigning a FAR. The normalization in the
denominator also corresponds to the fraction of detectable events expected
from the population given by $\Lambda$ (Farr, 2019). We make the following
definition of the denominator:

$\mu(\Lambda)\equiv\int\int
dd_{i}\,d\theta_{i}\,p(d_{i}|\theta_{i})\,p(\theta_{i}|\Lambda)\,p_{\rm
det}(\theta_{i},d_{i})$ (A2)

The total likelihood comes from considering the probability of the entire
dataset $\\{d_{i}\\}$ of $N_{d}$ detections (where the i$th$ event is detected
if $p_{\rm det}$ is 1), given a population with parameters $\Lambda$ that
predicts $N$ total events, $N\mu\equiv K(\Lambda)$ of which are expected to be
detected. The total likelihood is just the product of the contributions from
all the detected events, and the likelihood of detecting $N_{d}$ events,
considering the realization of $N_{d}$ comes from a Poisson distribution with
expected value $K$:

$p(\\{d\\}|\Lambda,K)=p(N_{d}|K(\Lambda))\prod^{N_{d}}_{i}p(d_{i}|\Lambda)\propto
K(\Lambda)^{N_{d}}e^{-K(\Lambda)}\mu(\Lambda)^{-N_{d}}\prod^{N_{d}}_{i}\int
d\theta_{i}p(d_{i}|\theta_{i})p(\theta_{i}|\Lambda).$ (A3)

If we assume a prior of $\pi(K)\propto{1}/{K}$, we can write the posterior
over $\Lambda$ and analytically integrate out the distribution over $K$:

$\begin{split}p(\Lambda|\\{d\\})&\propto\pi(\Lambda)\int
dK\frac{K(\Lambda)^{N_{d}}e^{-K(\Lambda)}}{K}\mu(\Lambda)^{-N_{d}}\prod^{N_{d}}_{i}\int
d\theta_{i}p(d_{i}|\theta_{i})p(\theta_{i}|\Lambda)\\\
&=\Gamma(N_{d})\pi(\Lambda)\mu(\Lambda)^{-N_{d}}\prod^{N_{d}}_{i}\int
d\theta_{i}p(d_{i}|\theta_{i})p(\theta_{i}|\Lambda)\\\
&\propto\pi(\Lambda)\mu(\Lambda)^{-N_{d}}\prod^{N_{d}}_{i}\int
d\theta_{i}p(d_{i}|\theta_{i})p(\theta_{i}|\Lambda)\end{split}$ (A4)

where $\Gamma(N_{d})$ does not depend on $\Lambda$, so marginalizing over
$K(\Lambda)$ with this choice of $\pi(K)$ allows us to factorize the above
equation, without explicitly considering its dependence on the Poisson term.

In practice, $p(\theta|\Lambda)$ does not need to be normalized, as any
prefactors will divide out in Eq. A1. We therefore only need to calculate
something proportional to $p(\theta|\Lambda)$. For reasons that will become
apparent, we compute $p(\theta|\Lambda)$ in terms of something proportional to
$\frac{dN}{d\theta}(\Lambda)$. We want to define a normalization factor for
the population distribution such that:

$\frac{1}{\alpha(\Lambda)}\left.m_{1}\frac{dN}{dm_{1}dqdVdt_{s}}\right|_{(m_{\rm
ref},q_{\rm ref},z_{\rm ref})}=1$ (A5)

where the differential rate is evaluated at a set of reference parameters.

With the distributions in Sec. 3 defined in terms of $\frac{dN}{dm}$ and
$\frac{dN}{dVdt_{s}}$ (i.e. source frame merger rate density
$\mathcal{R}(z)$), we can compute a normalization factor $\alpha(\Lambda)$:

$\begin{split}\alpha(\Lambda)=\left.m_{1}\frac{dN}{dm_{1}dqdVdt_{s}}(\Lambda)\right|_{(m_{\rm
ref},q_{\rm ref},z_{\rm
ref})}=&\left.m_{1}\frac{dN}{dm_{1}}\frac{dN}{dm_{2}}\frac{dm_{2}}{dq}\frac{dN}{dVdt_{s}}\right|_{(m_{\rm
ref},q_{\rm ref},z_{\rm ref})}\\\
&=\left.m_{1}^{2}\frac{dN}{dm}\frac{dN}{dm}\frac{dN}{dVdt_{s}}\right|_{(m_{\rm
ref},q_{\rm ref},z_{\rm ref})}\end{split}$ (A6)

Technically, we only know the $dN$ distributions up to a constant. As we will
see below, we will only be considering ratios of values that share the same
unknown constant, so we are free to leave it out for now.

Instead of computing $p(\theta|\Lambda)$ as
$p(\theta|\Lambda)=\frac{1}{N}\frac{dN}{d\theta}(\Lambda)$ exactly, we instead
make the following transformation in Equations A1 and A4:

$p(\theta|\Lambda)\rightarrow\frac{1}{\alpha(\Lambda)}\frac{dN}{d\theta}(\Lambda)$
(A7)

which is directly proportional to the differential rate and
$p(\theta|\Lambda)$.

For each draw of $\Lambda$, we have the normalization factor
$\alpha(\Lambda)$, related to the differential rate at our reference
parameters, as defined in Eq. A6. We outline below how we use this re-
expression to construct the rate independent of the likelihood. Note that this
change in Eq. A7 does not affect the likelihood, as it only affects
$\theta$-independent prefactors, which factor out of both the numerator and
denominator in Eq. A1.

Given the values of $\alpha(\Lambda)$ we have calculated, we wish to draw new
samples of $\alpha(\Lambda)$, given that the number of detections is a
Poisson-distributed realization. Recalling $K\equiv N\mu$, Eq. A7 means that
when we compute the denominator of Eq. A1, we are actually calculating the
ratio $\frac{K}{\alpha}$, and not $\mu$. Noting that the K-dependent integrand
of Eq. A4 is a Gamma-distribution for $K$ with shape parameter $N_{d}$ and a
scale parameter of 1, we can make the identification that
$\left<K\right>=N_{d}$ under this distribution. With $K\equiv N\mu$, we can
express the expectation value for $\alpha$ as
$\left<\alpha\right>=\frac{N_{d}}{K/\alpha}$.555Since $\left<K\right>=N_{d}$,
it follows that $\sigma^{2}_{\alpha}=\frac{N_{d}}{(K/\alpha)^{2}}$ As a final
step in post-processing, we can construct the true underlying distribution for
$\alpha$ by drawing samples $\alpha\sim\rm Gamma(\frac{N_{d}}{K/\alpha},1)$.
This gives us a distribution for the predicted merger rate at the reference
coordinates, given the normalization factor $\alpha$ we computed during the
hierarchical inference, assuming this is Poisson-distributed about the true
value and assuming a $\frac{1}{K}$ prior. With the distribution of $\alpha$,
we can scale $\frac{dN}{d\theta}$ to get the differential merger rate at any
set of coordinates $\theta$.

Note that we have written everything in this section in terms of $\theta$ as
if it is always the parameters in the population model, suppressing the fact
that there will be Jacobians in Eq. A1 to transform from these coordinates to
those in the detector-frame (or the priors from the single-event analyses).

##  Appendix B Accounting for Spin Distribution

In the analysis presented in the body of this work, we assume the
(uninformative) parameter estimation priors in the population reweighting.
Based on population-level mass-spin and mass-redshift correlations presented
in the literature (see, e.g., Callister et al. (2021); Biscoveanu et al.
(2022); LIGO Scientific Collaboration et al. (2023a)), we may expect the
assumed spin distribution can have an effect on our results. However, with the
relatively poor spin constraints in the population, we empiricall demonstrate
that this is likely not the case.

In Figure 12, we compare posteriors obtained from our main analysis ignoring
spins, with those obtained by reweighting the posterior samples from each
event and the sensitivty injections to a fiducial spin distribution. For this
fiducial spin distribution, the spin magnitudes are from a half-Gaussian
centered at $a=0$ with a standard deviation of 0.3, meant to model the
preferentially-small spin mangitudes inferred in LIGO Scientific Collaboration
et al. (2023a). For the contribution aligned projection of the spin tilt angle
($\cos\theta$), we use the mixture model introduced in Talbot & Thrane (2017),
with an aligned-spin fraction of $\xi=0.8$ and an aligned-spin spread of
$\sigma_{t}=1.9$, consistent with the results reported in LIGO Scientific
Collaboration et al. (2023a) (see references for definitions of these model
parameters). We find that reweighting to this spin distribution has a
negligible effect on our inferred population.

Figure 12: Selected population-level parameters from the evolving mass
distribution analysis, without reweighting spins (pink, as presented in
Section 4.2) and with reweighting spins to a fiducial population estimate
(blue). The near-identical posteriors show that the spin population
assumptions in this work do not cause a bias.

## 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. (2018) Abbott, B. P., et al. 2018, Living Rev. Rel., 21, 3, doi: 10.1007/s41114-020-00026-9
  * Abbott et al. (2021) Abbott, R., et al. 2021.  https://arxiv.org/abs/2112.06861
  * Acernese et al. (2015) Acernese, F., et al. 2015, Class. Quant. Grav., 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Andrews et al. (2021) Andrews, J. J., Cronin, J., Kalogera, V., Berry, C. P. L., & Zezas, A. 2021, ApJ, 914, L32, doi: 10.3847/2041-8213/ac00a6
  * Antonini et al. (2023) Antonini, F., Gieles, M., Dosopoulou, F., & Chattopadhyay, D. 2023, Monthly Notices of the Royal Astronomical Society, 522, 466–476, doi: 10.1093/mnras/stad972
  * Baxter et al. (2021) Baxter, E. J., Croon, D., McDermott, S. D., & Sakstein, J. 2021, ApJ, 916, L16, doi: 10.3847/2041-8213/ac11fc
  * Belczynski et al. (2016) Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016, Nature, 534, 512
  * Bingham et al. (2019) Bingham, E., Chen, J. P., Jankowiak, M., et al. 2019, J. Mach. Learn. Res., 20, 28:1.  http://jmlr.org/papers/v20/18-403.html
  * Biscoveanu et al. (2022) Biscoveanu, S., Callister, T. A., Haster, C.-J., et al. 2022, ApJ, 932, L19, doi: 10.3847/2041-8213/ac71a8
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs, 0.3.13.  http://github.com/google/jax
  * Briel et al. (2023) Briel, M. M., Stevance, H. F., & Eldridge, J. J. 2023, Monthly Notices of the Royal Astronomical Society, 520, 5724–5745, doi: 10.1093/mnras/stad399
  * Callister et al. (2020) Callister, T., Fishbach, M., Holz, D. E., & Farr, W. M. 2020, ApJ, 896, L32, doi: 10.3847/2041-8213/ab9743
  * Callister & Farr (2023) Callister, T. A., & Farr, W. M. 2023, arXiv e-prints, arXiv:2302.07289, doi: 10.48550/arXiv.2302.07289
  * Callister et al. (2021) Callister, T. A., Haster, C.-J., Ng, K. K. Y., Vitale, S., & Farr, W. M. 2021, ApJ, 922, L5, doi: 10.3847/2041-8213/ac2ccc
  * deBoer et al. (2017) deBoer, R., Görres, J., Wiescher, M., et al. 2017, Reviews of Modern Physics, 89, 035007, doi: 10.1103/RevModPhys.89.035007
  * Doctor et al. (2020) Doctor, Z., Wysocki, D., O’Shaughnessy, R., Holz, D. E., & Farr, B. 2020, The Astrophysical Journal, 893, 35
  * Duane et al. (1987) Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. 1987, Physics Letters B, 195, 216–222, doi: https://doi.org/10.1016/0370-2693(87)91197-X
  * Edelman et al. (2022) Edelman, B., Doctor, Z., Godfrey, J., & Farr, B. 2022, ApJ, 924, 101, doi: 10.3847/1538-4357/ac3667
  * Ezquiaga & Holz (2022) Ezquiaga, J. M., & Holz, D. E. 2022, Phys. Rev. Lett., 129, 061102, doi: 10.1103/PhysRevLett.129.061102
  * Farag et al. (2022) Farag, E., Renzo, M., Farmer, R., Chidester, M. T., & Timmes, F. X. 2022, The Astrophysical Journal, 937, 112, doi: 10.3847/1538-4357/ac8b83
  * Farah et al. (2023a) Farah, A. M., Edelman, B., Zevin, M., et al. 2023a, arXiv e-prints, arXiv:2301.00834, doi: 10.48550/arXiv.2301.00834
  * Farah et al. (2023b) Farah, A. M., Fishbach, M., & Holz, D. E. 2023b, arXiv e-prints, arXiv:2308.05102, doi: 10.48550/arXiv.2308.05102
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S., Fishbach, M., & Justham, S. 2020, The Astrophysical Journal, 902, L36, doi: 10.3847/2041-8213/abbadd
  * Farmer et al. (2019) Farmer, R., Renzo, M., De Mink, S. E., Marchant, P., & Justham, S. 2019, The Astrophysical Journal, 887, 53, doi: 10.3847/1538-4357/ab518b
  * Farr (2019) Farr, W. M. 2019, Research Notes of the American Astronomical Society, 3, 66, doi: 10.3847/2515-5172/ab1d5f
  * Farr et al. (2019) Farr, W. M., Fishbach, M., Ye, J., & Holz, D. E. 2019, ApJ, 883, L42, doi: 10.3847/2041-8213/ab4284
  * Farr et al. (2017) Farr, W. M., Stevenson, S., Miller, M. C., et al. 2017, Nature, 548, 426–429, doi: 10.1038/nature23453
  * Fishbach & Holz (2017) Fishbach, M., & Holz, D. E. 2017, ApJ, 851, L25, doi: 10.3847/2041-8213/aa9bf6
  * Fishbach & Holz (2020) —. 2020, ApJ, 891, L27, doi: 10.3847/2041-8213/ab7247
  * Fishbach et al. (2017) Fishbach, M., Holz, D. E., & Farr, B. 2017, ApJ, 840, L24, doi: 10.3847/2041-8213/aa7045
  * Fishbach et al. (2018) Fishbach, M., Holz, D. E., & Farr, W. M. 2018, The Astrophysical Journal, 863, L41, doi: 10.3847/2041-8213/aad800
  * Fishbach et al. (2022) Fishbach, M., Kimball, C., & Kalogera, V. 2022, ApJ, 935, L26, doi: 10.3847/2041-8213/ac86c4
  * Fishbach & van Son (2023) Fishbach, M., & van Son, L. 2023, ApJ, 957, L31, doi: 10.3847/2041-8213/ad0560
  * Fishbach et al. (2021) Fishbach, M., Doctor, Z., Callister, T., et al. 2021, The Astrophysical Journal, 912, 98, doi: 10.3847/1538-4357/abee11
  * Fowler & Hoyle (1964) Fowler, W. A., & Hoyle, F. 1964, ApJS, 9, 201, doi: 10.1086/190103
  * Gerosa & Berti (2017) Gerosa, D., & Berti, E. 2017, Phys. Rev. D, 95, 124046, doi: 10.1103/PhysRevD.95.124046
  * Gerosa & Fishbach (2021) Gerosa, D., & Fishbach, M. 2021, Nature Astronomy, 5, 749, doi: 10.1038/s41550-021-01398-w
  * Ghirlanda et al. (2016) Ghirlanda, G., Salafia, O. S., Pescalli, A., et al. 2016, A&A, 594, A84, doi: 10.1051/0004-6361/201628993
  * Golomb & Talbot (2023) Golomb, J., & Talbot, C. 2023, Phys. Rev. D, 108, 103009, doi: 10.1103/PhysRevD.108.103009
  * Hendriks et al. (2023) Hendriks, D. D., Van Son, L. A. C., Renzo, M., Izzard, R. G., & Farmer, R. 2023, Monthly Notices of the Royal Astronomical Society, stad2857, doi: 10.1093/mnras/stad2857
  * Karathanasis et al. (2023) Karathanasis, C., Mukherjee, S., & Mastrogiovanni, S. 2023, Monthly Notices of the Royal Astronomical Society, 523, 4539–4555, doi: 10.1093/mnras/stad1373
  * Kimball et al. (2021) Kimball, C., Talbot, C., Berry, C. P. L., et al. 2021, ApJ, 915, L35, doi: 10.3847/2041-8213/ac0aef
  * Kippenhahn & Weigert (1994) Kippenhahn, R., & Weigert, A. 1994, Stellar structure and evolution, 1st edn., Astronomy and astrophysics library (Berlin; New York: Springer-Verlag)
  * Kovetz et al. (2017) Kovetz, E. D., Cholis, I., Breysse, P. C., & Kamionkowski, M. 2017, Phys. Rev. D, 95, 103010, doi: 10.1103/PhysRevD.95.103010
  * Kroupa (2001) Kroupa, P. 2001, Monthly Notices of the Royal Astronomical Society, 322, 231–246, doi: 10.1046/j.1365-8711.2001.04022.x
  * Kroupa & Jerabkova (2019) Kroupa, P., & Jerabkova, T. 2019, Nature Astronomy, 3, 482–484, doi: 10.1038/s41550-019-0793-0
  * Kunz et al. (2002) Kunz, R., Fey, M., Jaeger, M., et al. 2002, ApJ, 567, 643, doi: 10.1086/338384
  * LIGO Scientific Collaboration & Virgo Collaboration (2021) LIGO Scientific Collaboration, & Virgo Collaboration. 2021, ApJ, 913, L7, doi: 10.3847/2041-8213/abe949
  * LIGO Scientific Collaboration et al. (2021) LIGO Scientific Collaboration, Virgo Collaboration, & KAGRA Collaboration. 2021, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run — Parameter estimation data release, Zenodo, doi: 10.5281/ZENODO.5546663
  * LIGO Scientific Collaboration et al. (2023a) LIGO Scientific Collaboration, VIRGO Collaboration, & KAGRA Collaboration. 2023a, Physical Review X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * LIGO Scientific Collaboration et al. (2023b) LIGO Scientific Collaboration, VIRGO Collaboration, & Kagra Collaboration. 2023b, ApJ, 949, 76, doi: 10.3847/1538-4357/ac74bb
  * LIGO Scientific Collaboration et al. (2023c) LIGO Scientific Collaboration, Virgo Collaboration, & KAGRA Collaboration. 2023c, Physical Review X, 13, doi: 10.1103/physrevx.13.011048
  * LIGO Scientific Collaboration et al. (2023d) —. 2023d, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run — O3 search sensitivity estimates, Zenodo, doi: 10.5281/ZENODO.7890437
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, Annual Review of Astronomy and Astrophysics, 52, 415–486, doi: 10.1146/annurev-astro-081811-125615
  * Maiolino et al. (2008) Maiolino, R., Nagao, T., Grazian, A., et al. 2008, Astronomy & Astrophysics, 488, 463
  * Mancarella et al. (2023) Mancarella, M., Iacovelli, F., & Gerosa, D. 2023, Phys. Rev. D, 107, L101302, doi: 10.1103/PhysRevD.107.L101302
  * Mandel et al. (2017) Mandel, I., Farr, W. M., Colonna, A., et al. 2017, MNRAS, 465, 3254, doi: 10.1093/mnras/stw2883
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086, doi: 10.1093/mnras/stz896
  * Marchant et al. (2019) Marchant, P., Renzo, M., Farmer, R., et al. 2019, The Astrophysical Journal, 882, 36, doi: 10.3847/1538-4357/ab3426
  * Mehta et al. (2022) Mehta, A. K., Buonanno, A., Gair, J., et al. 2022, The Astrophysical Journal, 924, 39, doi: 10.3847/1538-4357/ac3130
  * Miller & Hamilton (2002) Miller, M. C., & Hamilton, D. P. 2002, ApJ, 576, 894, doi: 10.1086/341788
  * Ng et al. (2021) Ng, K. K. Y., Vitale, S., Farr, W. M., & Rodriguez, C. L. 2021, ApJ, 913, L5, doi: 10.3847/2041-8213/abf8be
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Payne et al. (2023) Payne, E., Isi, M., Chatziioannou, K., & Farr, W. M. 2023.  https://arxiv.org/abs/2309.04528
  * Phan et al. (2019) Phan, D., Pradhan, N., & Jankowiak, M. 2019, arXiv preprint arXiv:1912.11554
  * Planck Collaboration (2020) Planck Collaboration. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Rahman et al. (2022) Rahman, N., Janka, H.-T., Stockinger, G., & Woosley, S. E. 2022, Monthly Notices of the Royal Astronomical Society, 512, 4503–4540, doi: 10.1093/mnras/stac758
  * Rakavy & Shaviv (1967) Rakavy, G., & Shaviv, G. 1967, ApJ, 148, 803, doi: 10.1086/149204
  * Ray et al. (2023) Ray, A., Magaña Hernandez, I., Mohite, S., Creighton, J., & Kapadia, S. 2023, arXiv e-prints, arXiv:2304.08046, doi: 10.48550/arXiv.2304.08046
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 100, 043027, doi: 10.1103/PhysRevD.100.043027
  * Safarzadeh & Farr (2019) Safarzadeh, M., & Farr, W. M. 2019, The Astrophysical Journal, 883, L24, doi: 10.3847/2041-8213/ab40bd
  * Sallaska et al. (2013) Sallaska, A. L., Iliadis, C., Champange, A. E., et al. 2013, ApJS, 207, 18, doi: 10.1088/0067-0049/207/1/18
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Schneider et al. (2018) Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69, doi: 10.1126/science.aan0106
  * Stevenson et al. (2019) Stevenson, S., Sampson, M., Powell, J., et al. 2019, ApJ, 882, 121, doi: 10.3847/1538-4357/ab3981
  * Talbot et al. (2019) Talbot, C., Smith, R., Thrane, E., & Poole, G. B. 2019, Phys. Rev. D, 100, 043030, doi: 10.1103/PhysRevD.100.043030
  * Talbot & Thrane (2017) Talbot, C., & Thrane, E. 2017, Phys. Rev. D, 96, 023012, doi: 10.1103/PhysRevD.96.023012
  * Talbot & Thrane (2018) Talbot, C., & Thrane, E. 2018, ApJ, 856, 173, doi: 10.3847/1538-4357/aab34c
  * Thrane & Talbot (2019) Thrane, E., & Talbot, C. 2019, PASA, 36, e010, doi: 10.1017/pasa.2019.2
  * Tiwari (2018) Tiwari, V. 2018, Classical and Quantum Gravity, 35, 145009, doi: 10.1088/1361-6382/aac89d
  * Tiwari & Fairhurst (2021) Tiwari, V., & Fairhurst, S. 2021, ApJ, 913, L19, doi: 10.3847/2041-8213/abfbe7
  * Van Son et al. (2022) Van Son, L. A. C., De Mink, S. E., Callister, T., et al. 2022, The Astrophysical Journal, 931, 17, doi: 10.3847/1538-4357/ac64a3
  * Vangioni et al. (2015) Vangioni, E., Olive, K. A., Prestegard, T., et al. 2015, MNRAS, 447, 2575, doi: 10.1093/mnras/stu2600
  * Vitale et al. (2019) Vitale, S., Farr, W. M., Ng, K. K. Y., & Rodriguez, C. L. 2019, ApJ, 886, L1, doi: 10.3847/2041-8213/ab50c0
  * Wong et al. (2023) Wong, K. W. K., Breivik, K., Farr, W. M., & Luger, R. 2023, ApJ, 950, 181, doi: 10.3847/1538-4357/acc863
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244, doi: 10.3847/1538-4357/836/2/244
  * Woosley (2019) Woosley, S. E. 2019, The Astrophysical Journal, 878, 49, doi: 10.3847/1538-4357/ab1b41
  * Woosley & Heger (2021) Woosley, S. E., & Heger, A. 2021, The Astrophysical Journal Letters, 912, L31, doi: 10.3847/2041-8213/abf2c4
  * Zevin et al. (2017) Zevin, M., Pankow, C., Rodriguez, C. L., et al. 2017, ApJ, 846, 82, doi: 10.3847/1538-4357/aa8408