# A Study of Systematics on the Cosmological Inference of the Hubble Constant
from Gravitational Wave Standard Sirens

Grégoire Pierra  0000-0003-3970-7970<EMAIL_ADDRESS>Université Lyon,
Université Claude Bernard Lyon 1, CNRS, IP2I Lyon/IN2P3,UMR 5822, F-69622
Villeurbanne, France  Simone Mastrogiovanni  0000-0003-1606-4183 INFN, Sezione
di Roma, I-00185 Roma, Italy  Stéphane Perriès  0000-0003-2213-3579 Université
Lyon, Université Claude Bernard Lyon 1, CNRS, IP2I Lyon/IN2P3,UMR 5822,
F-69622 Villeurbanne, France  Michela Mapelli  0000-0001-8799-2548 Institut
für Theoretische Astrophysik, ZAH, Universität Heidelberg, 69120 Heidelberg,
Germany

###### Abstract

Gravitational waves (GWs) from compact binary coalescences (CBCs) can
constrain the cosmic expansion of the universe. In the absence of an
associated electromagnetic counterpart, the spectral sirens method exploits
the relation between the detector frame and the source frame masses to jointly
infer the parameters of the mass distribution of black holes (BH) and the
cosmic expansion parameter $H_{0}$. This technique relies on the choice of the
parametrization for the source mass population of BHs observed in binary black
holes merger (BBHs). Using astrophysically motivated BBH populations, we study
the possible systematic effects affecting the inferred value for $H_{0}$ when
using heuristic mass models like a broken power law, a power law plus peak and
a multi-peak distributions. We find that with 2000 detected GW mergers, the
resulting $H_{0}$ obtained with a spectral sirens analysis can be biased up to
$3\sigma$. The main sources of this bias come from the failure of the
heuristic mass models used so far to account for a possible redshift evolution
of the mass distribution and from their inability to model unexpected mass
features. We conclude that future dark siren GW cosmology analyses should make
use of source mass models able to account for redshift evolution and capable
to adjust to unforeseen mass features.

Hubble constant - Dark sirens - Cosmology - Gravitational waves - Binary black
hole - Redshift evolution - Mass spectrum

††preprint: Short-author PP

##  I Introduction

September 2015 marks the date of the first detection of a gravitational-wave
(GW) [1] emitted by the merger of a binary black hole (BBH). Since then, the
LIGO-Virgo-KAGRA scientific collaboration has detected more than 90 GW signals
from compact binary coalescences (CBCs) [2, 3, 4]. These detected GW signals
enable us to directly measure the luminosity distances $d_{L}$ of the sources
without making any assumption on the cosmological model. However, GW signals
do not directly provide the source redshift unless an electromagnetic (EM)
counterpart is observed. The redshift is required to measure the Universe’s
expansion. This is an interesting prospect given the current tension on the
measurements of the Universe’s local expansion rate, the Hubble constant
$H_{0}$ [5]. The multiple GW detections opened the field of GW cosmology. To
exploit the tens of GW sources observed without EM counterpart, additional
methods have been developed to obtain the source redshift. The “Spectral
Siren” [6] method relies on the intrinsic degeneracy between the redshift $z$
of a source and its detector $\rm m^{\rm det}$ and source frame $\rm m^{\rm
s}$ masses, through the relation $\mathrm{m}^{\rm det}=(1+z)\mathrm{m}^{\rm
s}$. The redshift of the sources is implicitly estimated by measuring the
detector masses and jointly fitting the source mass distribution and
cosmological parameters [7, 8, 9]. This method depends on the choice of the
phenomenological models for the BBH source mass distribution and their
capacity to describe the true underlying BBH population [10, 11, 12, 7, 13,
14, 15, 16, 17, 18, 6]. The “spectral sirens” method is also tightly connected
with the “galaxy catalog method” [19, 15, 16, 17], which identifies the
possible galaxy hosts of the GW source by using galaxy surveys. Even with the
extra information from galaxy surveys, BBH source mass models still have to be
used. In fact, in the case that the galaxy survey is highly incomplete, the
source mass models are the ones that provide the implicit redshift of the GW
sources [8].

In this study, we explore how the commonly used phenomenological source mass
models can fit complex and astrophysically motivated distributions of BBH
mergers. The rationale of this study is to understand how a mismatch between
the source mass models and the population of BBHs would translate to an
$H_{0}$ bias. To make this study, we start with a BBH population generated
synthetically by modelling several astrophysical processes. Then we perform a
hierarchical Bayesian inference [9] to measure $H_{0}$ alongside the
phenomenological BBH mass models.

This paper is organized as follows. In Sec. II we describe the framework built
to simulate the GW detections, and introduce the basic principles of
hierarchical Bayesian inference. We also present the phenomenological BBH
population models for the mass spectrum. In Sec. III, we discuss several
sanity checks to test the validity of the framework. We also investigate the
robustness of the spectral sirens method when a wrong mass model is used for
the inference and when the source mass spectrum includes an unmodelled
redshift evolution. In Sec. IV, we consider a spectral sirens analysis using a
realistic catalog of BBH mergers, finding a biased estimate of $H_{0}$ if
simple phenomenological mass models are employed. Sec. IV also argues about
the possible different sources of the bias, specifically the importance of the
redshift evolution of the mass spectrum when estimating the cosmological
parameters. Finally, in Sec. V we draw the conclusions for this study.

##  II Bayesian inference, GW observations and CBC population modelling

In this section, we describe our simulation procedure. We first simulate BBH
observations using a detection criteria based on the signal-to-noise ratio
(SNR) that we describe in Sec. II.1. Then we use the hierarchical Bayesian
inference, summarized in Sec. II.2, to reconstruct cosmological and population
properties.

###  II.1  A fast generator of GW observations

The workflow that we use to generate the GW observations is the following.
From a population of BBH mergers described by the two source masses and the
redshift of the merger, we compute an approximate value of the optimal SNR.
Instead of calculating the SNR from the full waveform and the matched filter,
we use the relation below as a proxy for the optimal SNR [20, 7],

$\rho=\delta\times 9\left[\frac{\mathcal{M}_{c}}{\mathcal{M}_{c,{\rm
9}}}\right]^{\frac{5}{6}}\left[\frac{d_{L,{\rm 9}}}{d_{L}}\right],$ (1)

where $\mathcal{M}_{c}$ is the binary detector chirp mass, $d_{L}$ is the
luminosity distance, $\rho$ is the optimal SNR and $\delta$ a projection
factor. Although Eq. 1 can underestimate the SNR of events with chirp masses
above $30\rm{M}_{\odot}$, it should suffice to provide us with a
representative sample of the BBH population. The luminosity distance and the
binary chirp mass are computed assuming a flat $\Lambda$CDM cosmology with
$H_{0}=67.8$ ${\rm km\,s^{-1}\,Mpc^{-1}}$and $\Omega_{m,0}=0.308$, consistent
with cosmic microwave background measurements from Planck [21]. The parameters
with a superscript of “$\rm 9$” are the detector chirp mass and luminosity
distance values for which a signal with $\delta=1$ would have a SNR of 9. We
set the reference chirp mass to $\mathcal{M}_{c,{\rm 9}}=25$ $\rm{M}_{\odot}$,
the luminosity distance $d_{L,{\rm 9}}=1.5$ $\rm Gpc$ [22]. The projection
factor $\delta$ captures the variations of the optimal SNR when dealing with a
3-detector network and encapsulates the combined impact of the polarization of
the GWs and the sky localization. We draw it according to a cumulative density
function (CDF) introduced in [23].

To mimic the presence of detector noise in the data, we also calculate a
“detected” SNR $\rho^{\rm det}$ from the optimal SNR. The detected SNR is
drawn from a $\chi^{2}$ distribution with $6$ degrees of freedom (2 for each
detector). A GW is detected if its detected SNR exceeds a SNR threshold of
$12$ and the GW frequency corresponding to the innermost stable circular orbit
(ISCO) is higher than 15 Hz. The GW and ISCO frequencies are related by

$f_{\rm GW}=2\,{}f_{\rm ISCO},$ (2)

where $f_{\rm ISCO}$, is approximated as [24]

$f_{\rm
ISCO}=\frac{1}{2\pi}\frac{1}{6^{3/2}}\frac{c^{3}}{G}\left(\frac{1\rm{M}_{\odot}}{M^{d}_{\rm
tot}}\right)10^{3}\hskip 4.26773pt\rm Hz.$ (3)

For this study, we chose to not generate errors on the detector frame masses
and luminosity distances. We assume that we are perfectly able to measure them
from the GW signal. We make this choice for two motivations. The first is that
we want to ease the computational load of the analysis and explore multiple
test cases with thousands of GW detections. The second is that we want to
maximize the potential effect of the systematics arising from the
phenomenological mass models, which could be hidden by the error budget on the
masses.

Figure 1: Simple representation of the three phenomenological mass models, the
Broken Power Law in the left plot, the Power Law plus Peak in the middle plot
and the Multi Peak in the right plot.

###  II.2 The hierarchical likelihood

The detection of GW events can be described by an inhomogeneous Poisson
process in presence of selection biases [25, 9, 8, 7]. For a set of $N_{\rm
GW}$ GW signals detected over an observation time $T_{\rm obs}$, the
probability of having a specific GW data set $\\{\rm x\\}$ given population
hyper-parameters $\Lambda$ can be written as

$\displaystyle\mathcal{L}(\\{x\\}|\Lambda)$ $\displaystyle\propto$
$\displaystyle e^{-N_{\rm exp}(\Lambda)}\prod_{i}^{N_{\rm GW}}T_{\rm
obs}\int{\rm d}\theta{\rm d}z\;\times$ (4)
$\displaystyle\times\mathcal{L}_{\rm
GW}(x_{i}|\theta,z,\Lambda)\frac{1}{1+z}\frac{{\rm d}\rm N_{CBC}}{{\rm
d}\theta{\rm d}z{\rm d}t_{s}}(\Lambda).$

The hyper-parameters $\Lambda$ describe the BBH population and also include
cosmological expansion parameters. For a flat $\rm\Lambda CDM$ cosmology, the
expansion parameters are the Hubble constant $H_{0}$ and the fraction of
matter density today $\rm\Omega_{m,0}$. Other population hyper-parameters are
the ones controlling the shape of the BBH mass spectrum. We describe the
phenomenological models for the BBH mass spectrum in Sec. II.3 and in more
detail App. A.1. In Eq. 4 the individual GW likelihood $\mathcal{L}_{\rm
GW}(x_{i}|\theta,\Lambda)$ gauges the errors on the estimation of the
intrinsic parameters $\theta$ for a certain GW event $x_{i}$. For our study,
the intrinsic parameters $\rm\theta=(m_{1}^{s},m_{2}^{s})$ denotes the masses
in the source frame. In Eq. 4, the factor $\frac{1}{1+z}$ accounts for the
time dilatation between the detector and source frame times $t_{\rm d}$ and
$t_{\rm s}$.

$\frac{{\rm d}\rm N_{\rm CBC}}{{\rm d}\theta{\rm d}z{\rm d}t_{s}}(\Lambda)$
(5)

is the CBC merger rate per source time, comoving volume and masses. Finally,
$N_{\rm exp}$ is the expected number of GW detections during a time $T_{\rm
obs}$. This term encapsulates the GW selection biases, and can be written as

$N_{\rm exp}(\Lambda)=T_{\rm obs}\int{\rm d}\theta{\rm d}z\;P_{\rm
det}(z,\theta,\Lambda)\frac{{\rm d}\rm N_{\rm CBC}}{{\rm d}\theta{\rm d}z{\rm
d}t_{s}}(\Lambda)\frac{1}{1+z},$ (6)

where $P_{\rm det}(z,\theta,\Lambda)$ represents the probability of detecting
a certain CBC at a redshift $z$, characterized by some specific GW parameters
$\theta$ and given cosmological hyper-parameters $\Lambda$. The detection
probability is given by the integral of the GW likelihood over all the
detectable data realizations, namely

$\displaystyle P_{\rm det}(z,\theta,\Lambda)=\int_{x\in{\rm Det}}{\rm
d}x\mathcal{L}_{\rm GW}(x|\theta,z,\Lambda).$ (7)

The Spectral Siren CBC rate is then parameterized as

$\displaystyle\frac{{\rm d}\rm N_{\rm CBC}}{{\rm d}\theta{\rm d}z{\rm
d}t_{s}}$ $\displaystyle=$ $\displaystyle\frac{{\rm d}\rm N_{\rm CBC}}{{\rm
d}\overrightarrow{m}_{s}{\rm d}V_{c}{\rm d}t_{s}}\frac{{\rm d}V_{c}}{{\rm
d}z}(\Lambda)$ $\displaystyle=$ $\displaystyle R_{0}\psi(z;\Lambda)p_{\rm
pop}(\overrightarrow{m}_{s}|\Lambda)\times\frac{{\rm d}V_{c}}{{\rm d}z}.$

In the Eq. LABEL:eq:cbc_rate_full, $R_{0}$ is the local ($z=0$) CBC merger
rate per time per comoving volume, $\psi(z;\Lambda)$ is the phenomenological
function for the CBC merger rate that evolves with redshift with
$\psi(z=0;\Lambda)=1$. The $\rm p_{\rm pop}$ term is a probability density
function of the phenomenological models for the source frame masses $\rm
p_{\rm pop}(\vec{m}_{s}|\Lambda)$. The last term in Eq. LABEL:eq:cbc_rate_full
is the differential of the comoving volume per redshift, that for a flat
$\Lambda$CDM cosmology is

$\frac{{\rm d}V_{c}}{{\rm
d}z}=\left[\frac{c}{H_{0}}\right]^{3}\left[\int_{0}^{z}\frac{{\rm
d}z^{\prime}}{E(z^{\prime})}\right]^{2},$ (9)

where the term $E(z)$ in Eq. 9 is the dimensionless Hubble parameter defined
as

$E(z)=\frac{H(z)}{H_{0}}=\sqrt{\Omega_{m}(1+z)^{3}+(1-\Omega_{m})}.$ (10)

An important assumption that we implicitly made in Eq. LABEL:eq:cbc_rate_full
when writing $\rm p_{\rm pop}(\vec{m}_{s}|\Lambda)$ is that the CBC mass
spectrum does not evolve in redshift. This assumption is typically made in
most of the state-of-the-art literature [26, 27, 9, 8]. In the following, we
will discuss how the assumption for the source mass distribution to be
independent on redshift could translate to a $H_{0}$ bias.

Figure 2: Posterior predictive checks for the three catalogs of BBHs simulated
and inferred with the BPL (in blue), the PLP (in light blue) and the MLTP mass
model (in green) with 2000 GW events. The left panels correspond to the source
frame masses and redshift, from top to bottom $\rm m_{1}^{s}$, $\rm m_{2}^{s}$
and $z$. The injected populations are under the form of an histogram, and the
90% CL reconstructed spectrum as coloured contours. The right panels displays
the detector frame masses and luminosity distances, with selection effects are
taken into account, from top to bottom $\rm m_{1}^{det}$, $\rm m_{2}^{det}$,
$\rm d_{L}$.

###  II.3 The phenomenological population models of BBHs

Phenomenological CBC merger rate models of masses and redshift are central to
this study. Currently, spectral sirens analysis as in [26] use simple
parametric models describing the distribution of source frame masses and CBC
merger rate. For the masses, we employ three of these models: the Broken Power
Law (BPL), the Power Law plus Peak (PLP) and the Multipeak (MLTP). The three
mass models are represented in Fig. 1. The CBC merger rate as a function of
redshift is parameterized following the Madau $\&$ Dickinson [28] star
formation rate. In App. A.1 and App. A.2, we provide the detailed probability
density functions for the models and the set of population parameters and
prior ranges used in the rest of the paper. In the following, we provide a
general description of these models.

Figure 3: Marginal posterior of the Hubble constant derived from the spectral
sirens inference, with the PLP, BPL and MLTP mass models and a catalog of 2000
detected GW events. The vertical red line is the injected value in the
simulation, fixed to the Planck15 value of $H_{0}=67.7\hskip 4.26773pt{\rm
km\,s^{-1}\,Mpc^{-1}}$.

The Broken Power Law (BPL) [29] is the simplest model that we employ. It is
composed of a power law, with a smooth tapering at its low mass edge plus a
break in the power law at a mass $\rm m_{break}$. The low-mass smoothing is
introduced to model the effects of the stellar progenitor metallicity [29],
which could create a smooth transition for black hole production. The break in
the power law mimics the left edge of the pair-instability supernovae [30]
gap. The second power law, defined after the break, could explain the
existence of a second population of BBHs inside the pair-instability
supernovae gap due to dynamically formed BBHs.

Figure 4: The left plot presents the marginal distribution of the Hubble
constant inferred with the PLP and the MLTP mass models while generating the
BBH catalog with the BPL model. The middle plot is produced simulating the BBH
events with the PLP and inferring with the BPL and the MLTP. Finally the plot
on the right shows the $H_{0}$ results obtained while generating the GW events
with the MLTP model and inferring with the BPL and PLP mass models. On each of
the plots, the true value of the Hubble constant in each population is
represented by the vertical red line.

The Power Law plus Peak (PLP) [31] is based on a power law with the same low
mass smoothing as the BPL and incorporates a gaussian peak introduced to
account for a potential accumulation of BBH just before the pair-instability
supernovae gap [31]. The PLP mass model is more flexible than the BPL, the
position and the width of the peak can vary and allow the model to potentially
catch high mass events, without changing the lower part of the spectrum. The
fraction of events in the Gaussian peak is given by the parameter
$\lambda_{\mathrm{peak}}$ (see App. A.1).

The Multi Peak (MLTP) mass model is built similarly to the PLP, but with a
second Gaussian peak available at higher masses. The motivations for the
second Gaussian component are that BBH systems could arise from second-
generation black holes produced by hierarchical mergers. The MLTP model is
more flexible and can reduce into a PLP, if the secondary component is not
supported by the Bayesian inference.

Figure 5: Posterior predictive checks for the wrong model analysis. The blue
histograms are the injected population of BBHs, generated with a MLTP mass
model, with the same set of parameters as depicted in Sec. III.1. The two
coloured contours are the 90% CL reconstructed populations, inferred
respectively with the PLP in light blue and the BPL mass model in green. This
posterior predictive check corresponds to the result presented on the right
panel of Fig. 4.

In previous studies [26, 7, 8], a significant correlation between the Hubble
constant and possible local over-densities in the BBH mass spectrum was found.
In particular, the position of the gaussian peak for the PLP model ($\mu_{g}$)
shown a strong correlation with $H_{0}$.

##  III  Application to analytical BBH populations

Here we consider diverse populations of BBHs generated from parametric models
to study how the PLP, BPL and MLTP models infer $H_{0}$ for several scenarios.
In Sec. III.1, we use the same mass models (BPL, PLP and MLTP) for the
generation of the BBHs catalog for the inference. In Sec. III.2, we still
generate BBHs from one of our three models, but we perform the inference using
the other two. In the Sec. III.3 we simulate GW events with an extended
version of the PLP that includes a redshift evolution of the Gaussian peak in
order to test the bias on $H_{0}$ obtained by the redshift-independent models.
We will quantify the $H_{0}$ bias according to what confidence interval (C.I.)
the true values of $H_{0}$ are found from the hierarchical posterior. In
principle, biases should be quantified using parameters-parameters (PP) plots
that indicate what is the fraction of time that the true value of $H_{0}$ is
found in a given C.I. for the posterior. As an example, if no bias is present,
the true value of $H_{0}$ would be found in the 50% of C.I. for 50% of the
cases and so on. However, here we do not calculate PP plots as they are
computationally demanding. We argue for the presence of possible biases when
the true value of $H_{0}$ is excluded at 2$\sigma$ or more C.I. The motivation
to choose this criterion is that the probability of obtaining such discrepancy
due to statistical fluctuations, and not due to biases, is only $5\%$.

In the rest of the paper, we will use the posterior predictive checks (PPCs)
to assess how well the PLP, BPL and MLTP can fit the mass spectrum. The
astrophysical posterior predictive distribution, given the observed GW events,
for redshift and source masses is defined as

$p_{\rm pop}(z,\mathrm{m}_{1}^{\rm s},\mathrm{m}_{2}^{\rm s}|\\{x\\})=\int{\rm
d}\Lambda\,p_{\rm pop}(z,\mathrm{m}_{1}^{\rm s},\mathrm{m}_{2}^{\rm
s}|\Lambda)p(\Lambda|\\{x\\}).$ (11)

The posterior predictive distribution for the luminosity distance and detector
frame masses is defined as

$\displaystyle p_{\rm obs}(d_{L},\mathrm{m}_{1}^{{\rm d}},\mathrm{m}_{2}^{{\rm
d}}|\\{x\\})=$ $\displaystyle\int{\rm d}\Lambda\,p_{\rm
pop}(d_{L},\mathrm{m}_{1}^{{\rm d}},\mathrm{m}_{2}^{{\rm d}}|\Lambda)\times$
(12) $\displaystyle p(\Lambda|\\{x\\})P_{\rm det}(d_{L},\mathrm{m}_{1}^{{\rm
d}},\mathrm{m}_{2}^{{\rm d}}).$

We will discuss how the mapping between source and detector frame, and the
mismatch of the BBH mass spectrum will translate to an $H_{0}$ bias.

###  III.1 Using the same mass models for simulation and inference

We simulate three test populations of BBH mergers, with the three source mass
distributions and the CBC merger rate evolution described before. The
parameters of the population models used for the simulation of the catalogs
are summarized in App. B.2. For the three populations, we generate 2000
detected GW events, that will be used in the spectral sirens analysis. We
choose to work with 2000 GW events since for this number of detections it is
typically possible to have a data-informed posterior bounded in the prior
range $H_{0}\in[20,140]{\rm km\,s^{-1}\,Mpc^{-1}}$.

Figure 6: Evolution of the $H_{0}$ value inferred with respect to the final
position of the PLP peak at $z=1$. The blue points correspond to the median of
the inferred posterior of the Hubble constant, for each of the analysis. The
point of reference (no redshift evolution) is for $\mu_{g}=30$ at $z=1$. The
error bars are the $1\sigma$ intervals of the posteriors and the green line is
the best linear fit passing by all points.

We jointly estimate the cosmological parameter $H_{0}$, as well as all the
mass and the CBC merger rate parameters using hierarchical Bayesian inference.

Fig. 3 presents the marginal posterior distributions for $H_{0}$ obtained for
each source mass model. For all the simulations, the true value of the Hubble
constant is estimated within the $90\%$ C.I. We verified that the true values
of all the other population parameters describing the CBC merger rate and mass
spectrum are contained in the 90% C.I. of their marginal posterior. The
inferred $H_{0}$ is more precise for the PLP and the MLTP mass models rather
than the BPL. The more precise constraint on $H_{0}$ is given by the fact that
the PLP and MLTP contain sharper features in the source mass spectrum.

Fig. 2 shows the PPCs for the BPL, PLP and MLTP mass models. The figure
demonstrates that all the structures in the mass and rate spectrum are well
reconstructed by the analysis. In particular, the inference can correctly
reconstruct all the features present in the mass spectra. In conclusion, when
using the correct mass model for inference, both the correct population and
cosmological expansion parameters are recovered. Additionally, this test
demonstrates the importance of having features in the BBH mass spectrum to
improve the precision for cosmological expansion parameters.

Figure 7: Posterior predictive check for the simulated populations,
incorporating an additional linear evolution of the gaussian component
$\mu_{g}(z)$ with respect of the redshift, between $z=0$ and $z=1$. The
coloured contour are the 90%CL estimated spectrum obtained with the spectral
sirens analysis. The coloured histograms display the injected populations in
each run, they are denoted by the position of their gaussian peak at $z=1$.
Only three out of the seven populations are shown: The vanilla analysis in
purple (reference run), the $\mu_{g}^{1}=35\rm{M}_{\odot}$ in orange and the
$\mu_{g}^{1}=25$ in blue.

###  III.2 Using different mass models for simulation and inference

We generate 3 populations of 2000 detected CBCs from the BPL, PLP and MLTP
mass models and then we perform the estimation of populations and cosmological
parameters with the other two mass models.

Fig. 4 shows the marginal posterior distributions of the Hubble constant
inferred. We find that, when the detections are simulated with the BPL and the
inference performed with PLP and MLTP, the true value of $\rm H_{0}$ is
contained in the $68\%$ C.I. despite the wrong mass model was used for the
inference. In contrast, when simulating with the PLP and MLTP and inferring
with the BPL, we find the true value of $H_{0}$ is excluded at a 98.7% C.I.
level. This bias is linked to the difficulties of the BPL model in catching
all the sharp mass features contained in the PLP and MLTP mass spectra. The
mismatching of the mass features result in a systematically mismatched
redshift for the sources and then on a bias on $H_{0}$.

Similarly to Sec. III.1, we look at the PPC to understand the cosmological
results. Fig. 5 shows the PPCs for the case in which the simulated population
was generated from a MLTP model and the Bayesian inference performed with the
BPL and PLP models. From the source frame mass spectrum of $\rm m_{1}$ and
$\rm m_{2}$, we see that the PLP inference can reconstruct a peak located
around $50M_{\odot}$. While the reconstructed mass spectrum from the BPL
misses the high mass feature by underestimating it. The resulting bias in
$H_{0}$ can be understood as, for lower values of $H_{0}$, GW events are
placed at lower redshifts hence with higher source frame masses. By lowering
the $H_{0}$ value, the BPL is trying to fit the over-density of high-mass
events produced by the MLTP. Concerning the simulations with the most complex
of the models, the MLTP, we find that both the PLP and BPL are not able to
correctly reconstruct the lowest part of the mass spectrum. However, the
mismatch is not sufficient to introduce a $H_{0}$ bias for the reconstruction
with the PLP. The CBC merger rate as a function of redshift, in all cases, is
correctly reconstructed. In conclusion, the sharper mass features in the BBH
mass spectrum can source a significant bias on $H_{0}$.

###  III.3 Mass models response to a redshift evolution of the mass spectrum

In this section, we perform a test to examine the robustness of redshift-
independent mass models to a possible evolution of the mass spectrum in
redshift. In doing so, we will assume a linear dependence between the features
of the mass spectrum and the redshift. This simple model has already been
explored in [6] for the evolution of the low and high edges of the mass
spectrum, finding a resulting mild $H_{0}$ bias when the low edge of the mass
spectrum is not modelled correctly. In this paper, we modify the PLP mass
model to include a peak position linearly evolving with redshift. The redshift
evolution of the mass spectrum for the simulated BBH population is
incorporated in the parameter that governs the position of the Gaussian
component, while the edges are kept fixed in redshift. This type of evolution
could be the result of an interplay of the Pair Instability Supernova mass
scale and time-delays between binary formation and merger [14, 32] that highly
suppress the mass spectrum above the PISN scale (the gaussian peak). The
impact of such models has been explored in [33] on GWTC-3 finding no
statistically significant preference and a consistent $H_{0}$ posterior
against the non-evolving PLP.

Here, we do not suppress the mass spectrum above the PLP peak but just assume
a linear dependency:

$\mu_{g}(z)=\mu_{g}^{0}+z\left(\mu_{g}^{1}-\mu_{g}^{0}\right),$ (13)

where $\mu_{g}^{0}=\mu_{g}(z=0)$ and $\mu_{g}^{1}=\mu_{g}(z=1)$ are the
positions of the gaussian peak at $z=0$ and $z=1$. We fix
$\mu_{g}^{0}=30\rm{M}_{\odot}$ and we considered 7 values for $\mu_{g}^{1}$
from $25\rm{M}_{\odot}$ to $35\rm{M}_{\odot}$.

For each simulation, we run a full hierarchical Bayesian inference using the
redshift-independent PLP to reconstruct $H_{0}$ as well as the other
population parameters from 2000 GW detections. We report the $H_{0}$
estimation with 68.3% C.I. for these analyses in Fig. 6, showing the
variations of the inferred $H_{0}$ value as a function of the position of
$\mu_{g}^{1}$. For instance, if the Gaussian peak shifts of about 5 solar
masses from redshift 0 to redshift 1, $H_{0}$ is found outside the 99.7% C.I.
Fig. 6 shows that the redshift evolution of the structures in the mass
spectrum can bias the inferred value of the Hubble constant. Even when
considering a mild evolution of the mass spectrum in redshift, the bias on the
inferred $H_{0}$ is significant. We also notice that the systematic bias
obtained on $H_{0}$ is linearly proportional to the magnitude of the redshift
evolution of the mass spectrum feature. As we will describe in the following
paragraph, the bias is introduced by the fact that the source-frame mass
spectrum is not reconstructed correctly.

Figure 8: Kernel density estimation of the primary masses $\rm m_{1}$ in
source frame from the A03 BBH catalog. Each plot corresponds to a slice in
redshift, going from $z=0$ up to $z=7$. This figure highlights the overall
mass spectrum evolution with respect to the redshift, especially the
appearance of peaks and their shift. Figure 9: Kernel density estimation of
the primary masses $\rm m_{1}$ in source frame from the A03 BBH catalog. Each
plot corresponds to a slice in redshift, going from $z=0$ up to $z=1.75$. This
figure highlights the overall mass spectrum evolution with respect to the
redshift, especially the appearance of peaks and their shift.

Fig. 7 summarizes the PPC for this test case. The right panel shows the
detected mass distributions for three of the seven simulations. In the
detector frame, the 3 simulations display similar distributions for the masses
and luminosity distance of detected GW events. The distributions are the ones
that are fit by the non-evolving PLP models. However, in the source frame, the
three simulations have a significantly different mass spectrum (due to
redshift evolution). The PLP model is not able to catch this evolution and
produces a wrong reconstruction of the mass spectrum. In particular, the non-
evolving PLP model always reconstructs a peak located at 30 $\rm{M}_{\odot}$.
This is the consequence of the fact that about 80% of the GW detections are
located at low redshift ($z<0.3$), where even for redshift evolving models,
the peak of the Gaussian component is located around $\mu_{\rm
g}^{0}=30\rm{M}_{\odot}$.

The introduction of a $H_{0}$ bias is due to the misplacing of the Gaussian
component for events at higher redshifts. When $\mu_{\rm g}^{1}>\mu_{\rm
g}^{0}$ and the PLP is not able to fit it, we obtain a $H_{0}$ biased to
higher values. The motivation for this bias is that, for high $H_{0}$ values,
GW events are placed at higher redshifts and hence at lower source masses. By
placing the events at lower source masses, the PLP model is able to
accommodate them in the peak that fits around $30\rm{M}_{\odot}$ in the source
frame. When $\mu_{\rm g}^{1}<\mu_{\rm g}^{0}$, we observe a bias to lower
values of $H_{0}$. Similarly to the previous case, this is due to the fact
that GW events are placed at lower redshifts, and therefore they will have a
higher source mass that can be place in the peak at $30\rm{M}_{\odot}$. This
simple test case shows that the evolution in redshift of a feature in the mass
spectrum could potentially bias the GW-based $H_{0}$ estimation.

Figure 10: Marginal posterior distributions of the Hubble constant $H_{0}$
obtained with 2000 GW events, using the PLP in light blue, the BPL in blue and
the MLTP mass model in green. The vertical red line represent the true value
of the Hubble constant intrinsic to the A03 catalog and equal to $H_{0}=67.7$
${\rm km\,s^{-1}\,Mpc^{-1}}$.

###  III.4 Summary of the tests

To summarize, the $H_{0}$ bias is closely linked to the ability of the mass
model to catch and reconstruct the features of the CBC mass spectrum. A bias
on the $H_{0}$ value can be introduced from several factors.

We have shown that in the case where the true underlying population of CBCs
presents sharp features, and the reconstruction model is not able to include
them, a significant bias on $H_{0}$ can be obtained. We have also shown that
if any of the mass features of the spectrum evolve in redshift, a $H_{0}$ bias
can be introduced even if the evolution is of the order of a few solar masses.

##  IV Application to a complex BBH population

In this section, we study how the BPL, PLP and MLTP models respond to a
population of BBHs generated from synthetic astrophysical simulations. We use
a catalogue, referred to as A03, of BBH mergers formed by various formation
channels from [34, 35]. We then present the results of the hierarchical
Bayesian analysis obtained with the BPL, PLP and MLTP mass models in Sec.
IV.2. Finally, in Sec. IV.3, we investigate the possible sources of bias for
the Hubble constant.

###  IV.1  General description of the A03 catalog

For this study we employ the ”A03” catalog from [35]. This model consists of
four different channels: isolated binary evolution, dynamical assembly of BBHs
in young, globular, and nuclear star clusters. They are mixed based on their
corresponding redshift-dependent rates as described in [35]. The dynamical
formation channels assume an initial black hole mass function obtained with
the MOBSE population synthesis code [36]. In model A03, we assume the delayed
model by [37] for core-collapse supernovae, the pair instability supernova
formalism by [38], and we take into account the dependence of the BH mass
function on the metallicity of the progenitor stars, as discussed in [36]. The
BHs then pair up dynamically with other BHs and harden because of dynamical
interactions in their host star clusters [34]; we also allow for hierarchical
mergers of BBHs, as described in [35]. The isolated binary evolution channel
consists of systems evolved with the MOBSE binary population synthesis code
accounting for binary star evolution [36]. For these runs, we include a
treatment for mass transfer, common envelope (with efficiency parameter
$\alpha=1$), tidal evolution, natal kicks, and gravitational-wave decay as
described in [36]. A further discussion of the astrophysical processes
included in model A03 and the corresponding uncertainties are beyond the scope
of this study. Here, we explore the possibility that non-trivial structures of
the mass distribution for BBHs can introduce a bias on the estimation of
$H_{0}$.

Figure 11: Posterior predictive check for the spectral sirens analysis with
2000 GW observations, generated from the A03 BBH catalog. The plain black
histograms are the injected population, in source frame on the left panels and
in detector frame on the right panels. The three coloured contours depict the
90%CL for the three inferences (BPL, PLP and MLTP).

Fig. 8 shows the distribution of the primary mass $\rm m_{1}^{s}$ for
different redshift bins. The BBH mass spectrum presents several features.
Moving from lower to higher masses, a feature at $10\rm{M}_{\odot}$, a
“valley” between $13-15\rm{M}_{\odot}$ and a second large peak at
$16\rm{M}_{\odot}$. In the high mass region of the spectrum, a local over-
density of BBHs is visible around $25M_{\odot}$, plus several less significant
peaks up to $90\rm{M}_{\odot}$. These structures are also visible on the $\rm
m_{2}^{s}$ distribution. Another property of the A03 catalog is the evolution
in redshift of the primary and secondary source frame masses. From Fig. 8, we
notice that the structures in the mass spectrum are evolving as the redshift
increases. The feature at $20\rm{M}_{\odot}$ is nonexistent for mergers below
$z=2$, and emerges while shifting to the lower masses in the higher redshift
bins.

Fig. 9 shows the BBH mass spectrum in a redshift range of $0<z<1.75$. This
range corresponds to the typical range in which our simulation is able to
detect BBH mergers. The evolution of the BBH mass spectrum at low redshifts is
fainter than the one at high redshift, nevertheless peaks around
$11\rm{M}_{\odot}$ appear while shifting to higher masses and the structure
around $25\rm{M}_{\odot}$ is drifting to lower masses. Since real GW events at
high redshift ($z>1$) are hardly detectable in the O4 like network, this
effect is thought to be subdominant for population and cosmological inferences
[39]. If these structures are evolving with respect to the redshift, as
discussed in [40, 41, 42, 43, 44], the mass models could mismatch the features
and produce a biased value of $H_{0}$ as shown in Sec. II.3.

Figure 12: Posterior predictive check of the spectral sirens analysis of 2000
GW observations generated from the A03 BBH catalog. The plain black histograms
are the injected population, in source frame on the left panels and in
detector frame on the right panels. Each coloured histogram corresponds to a
bin in redshift from $z=0$ to $z=2$. The three coloured contours depict the
90%CL computed with the BPL, PLP and MLTP mass models.

###  IV.2 Vanilla analysis of the A03 catalog

Using the A03 BBH catalog, and with the framework described in Sec. II.1, we
generate GW catalogs made of $2000$ detected GW events. For each catalog, we
perform a full hierarchical Bayesian inference reconstructing posteriors on
$H_{0}$ as well as the other population parameters of the redshift rate and
the three redshift-independent BBH mass models. The prior ranges used for the
Bayesian inference are reported in the Tab. 1, Tab. 2 and Tab. 3.

Figure 13: Marginal posterior distributions of the Hubble constant $H_{0}$
obtained with 2000 GW events, using the PLP and the MLTP mass models. The
vertical purple line represent the true value of the Hubble constant intrinsic
to the A03 catalog and equal to $H_{0}=67.7$ ${\rm km\,s^{-1}\,Mpc^{-1}}$.
This distribution is obtained using the A03 catalog, where the redshift
evolution of the mass spectrum has been removed.

Fig. 10 shows the marginal posterior distributions of the Hubble constant
estimated. When using the PLP and MLTP mass models, we find an $H_{0}$ bias
towards higher values. In particular, we estimate for the PLP
$H_{0}=103^{+11}_{-10}{\rm km\,s^{-1}\,Mpc^{-1}}$ and $H_{0}=87^{+9}_{-8}{\rm
km\,s^{-1}\,Mpc^{-1}}$ at 68.3 % C.I. The true value of $H_{0}$ is excluded
from the 99.7% C.I. for the PLP and the 95% C.I. for the MLTP. From Fig. 10,
it is possible to note that the BPL mass model obtains an uninformative
$H_{0}$ posterior. We obtain a value of $H_{0}=82^{+20}_{-19}{\rm
km\,s^{-1}\,Mpc^{-1}}$ which includes the injected $H_{0}$ within the $68\%$
C.I. As we will argue later, the BPL recovers a less informative posterior on
$H_{0}$ as it does not contain any strong mass feature. We note that, even if
the BPL does not display a $H_{0}$ bias for 2000 events, it might still
display it as more and more GW detections are collected. The bias is not
introduced by any of the other population parameters “railing” over the prior
ranges. In fact, we verified that all the parameters are well-constrained in
their prior ranges.

Figure 14: Comparison of posterior predictive checks for the spectral sirens
analysis of the A03, with and without redshift evolution of the mass spectrum,
using 2000 GW events. The blue contour and histograms are the vanilla
population, and are the same results shown in Fig. 11. The purple histograms
are the new A03 population with no redshift evolution of the mass spectrum and
the purple contours are the 90%CL inferred spectrum with the population,
denoted ”A03 no redshift evolution”.

The PPC plot in Fig. 11 gives insights into the origin of the $H_{0}$ bias.
While the reconstructed distributions of detector frame masses and luminosity
distance do not show any particular deviation between the injected population
and the recovered one, the reconstructed source population display larger
differences. In the source frame, the CBC merger rate is correctly
reconstructed by the BPL, whereas the PLP and MLTP significantly underestimate
the presence of BBHs between $40\rm{M}_{\odot}$ and $80\rm{M}_{\odot}$.
Therefore, from observed data, the PLP and MLTP, predict a significantly
larger number of BBHs at lower source masses. This is achieved by pushing the
value of the Hubble constant at higher values so that events are found at
higher redshifts and have on average a lower value of source masses. Another
indicator in support of the finding that the PLP and MLTP models place events
at higher redshifts is the reconstruction of the CBC merger rate in redshift.
In fact, from the bottom left plot of Fig. 11, we can see that the PLP and
MLTP reconstruct a biased merger rate in redshift that prefers higher values
of redshifts for the BBH mergers. Instead, the BPL inference manages to fit
the overall trend of the A03 source frame population. The BPL reconstructs
better the source frame masses between $40\rm{M}_{\odot}$ and
$80\rm{M}_{\odot}$.

Even though the PLP and MLTP have enough degrees of freedom to approximate a
smooth distribution of masses (as the ones present in the A03 catalog), they
still fail to reconstruct the mass spectrum. The motivation for which the PLP
and MLTP are not able to correctly reconstruct the mass spectrum of the A03
catalog and $H_{0}$ is not trivial. A first hypothesis for this result was
that the PLP and MLTP were catching local peaks in the mass spectrum (at a
fixed redshift) that were evolving in redshift. In fact, from the posterior
distributions, we observed that the PLP and MLTP were finding local peaks at
low masses ($\approx 17M_{\odot}$) with a large standard deviation and a
higher peak at $\approx 30M_{\odot}$. To test this hypothesis, we performed a
PPC in Fig. 12 dividing the true population in redshift bins. Unfortunately,
the PPC check does not clearly show a discrepancy between the population
reconstruction and the any of the mass spectra for the redshift bins. As this
test is inconclusive on the origin of the $H_{0}$ bias, we performed
additional studies that we discuss below.

###  IV.3  Investigating the sources of the $H_{0}$ bias: blinding the mass-
redshift relation

In order to test further the possibility that the redshift evolution of the
mass spectrum is introducing a $H_{0}$ bias, we performed a simulation where
we blind the A03’s mass spectra to the redshift evolution. To remove the
redshift dependency of the mass spectrum, we randomly shuffle the pairs of BBH
merger redshifts and masses. This procedure artificially removes the redshift
dependence of the mass spectrum while conserving its non-trivial shape. For
the Bayesian inference, we only use the PLP and MLTP mass models since they
were displaying a non-negligible bias of $H_{0}$.

The result of the joint inference with the PLP mass model is presented in Fig.
13. For the PLP, we obtain $H_{0}=67.0^{+9.7}_{-8.9}{\rm
km\,s^{-1}\,Mpc^{-1}}$ and for the MLTP we obtain
$H_{0}=75.13^{+7.8}_{-7.9}{\rm km\,s^{-1}\,Mpc^{-1}}$, the true value of
$H_{0}$ is contained within the 68.3% C.I. of the posterior. Fig. 14 shows the
PPC plot for this test. The figure overlaps the reconstruction for the
original A03 catalog and the one after the redshift-blinding procedure. From
the detector frame point of view, we see that the two reconstructed
distributions do not significantly differ. However, from the source frame
point of view, the PLP model, for the redshift-blinded A03 catalog, is able to
reconstruct better the BBH masses present in the range
$30\rm{M}_{\odot}-60\rm{M}_{\odot}$. Although the reconstruction of the mass
spectrum above $60\rm{M}_{\odot}$ is still poor, the PLP for the non-redshift
evolving A03 is able to include the true value of $H_{0}$. As a consequence,
also the reconstruction of the BBH merger rate as a function of redshift
slightly improves, resulting in an underestimation of about 20%-30% of the
true CBC merger rate. We have also found similar results for the MLTP model.
We verified by running the simulations with 20 different population
realizations that the value of $H_{0}$ obtainable by the PLP and MLTP from the
redshift-blinded A03 is always unbiased.

It is important to notice that, even after blinding the redshift evolution of
the A03 mass spectrum, the PLP and MLTP are not yet able to completely
reconstruct the source frame distributions. As such, we still expect a
systematic bias on $H_{0}$ to be present and hidden by the large statistical
uncertainties that can be seen in Fig. 13. Identifying the source of these
systematic biases is not trivial, but from the test that we performed, we are
sure that the redshift evolution of the mass spectrum is expected to be the
dominant source of the $H_{0}$ bias.

##  V Conclusions

In this work we have extensively studied the interplay between systematic
biases on $H_{0}$ and the reconstruction of the BBH mass spectrum for Spectral
Siren Cosmology. We performed several simulations using BBH merger catalogs
generated with 3 of the current mass models used in literature and we also
considered a catalog of synthetically simulated BBHs.

In Sec. III, we showed that simple phenomenological models that do not contain
sharp features could introduce a bias on $H_{0}$ if the true population of
BBHs includes a local over-density of sources. Moreover, we showed that the
redshift evolution of source frame mass features can induce a $H_{0}$ bias if
not accounted for in the mass models. In particular, the Hubble constant would
act as a free parameter and shift to higher value if the true source mass
spectrum is underestimated, hence placing more events at higher redshift, and
equivalently decreasing the global source mass spectrum. This results is
respectively true if the source mass spectrum is overestimated, creating a
shift of $H_{0}$ to lower values.

In Sec. IV, using the A03 BBH catalog, we showed that even if the detector
frame is well reconstructed by the inference, we might still get a biased
value of the Hubble constant. Deviations up to $2\sigma$ and $3\sigma$ from
the injected value were observed for the MLTP and PLP inferences,
respectively. We observed that while each model yielded to different
populations in the source frame, they produce remarkably similar populations
in the detector frame. This discrepancy in the source frame mapping appeared
to be related to the response of the models to the GW events, as each model
attempted to accommodate to the true population in its own unique manner.
Consequently, the Hubble constant is used as a compensatory factor, when
necessary, to align the inferred populations with the observed data. Finally,
we find that the $H_{0}$ bias obtained with the PLP and MLTP mass model
disappears when the redshift evolution of the source mass spectrum is removed
from the population.

In light of these findings, we demonstrate that the most commonly used mass
models, namely PLP, BPL, and MLTP, combined with the spectral sirens analysis
for GW cosmology, are not suitable choices if the BBH population exhibits even
a slight redshift evolution of its mass spectrum. In addition, even if the
mass spectrum is not evolving with redshift, the mismatched of sharp
structures in the BBH mass distribution can also be responsible of a biased
estimation of the Hubble constant. We showed that this kind of assumption can
be dangerous for the spectral sirens method and more generally for methods
that make use of non evolving parametric models of the mass spectrum since the
constraining power on the Hubble constant relies on the ability of the
phenomenological mass models to fit the position of structures in the mass
spectrum. Given that the true distribution of BBHs in our universe remains
unknown and may differ significantly from the current state-of-the-art models,
future GW cosmology studies based on the spectral sirens method should
incorporate source mass models capable of absorbing effects such as the
redshift evolution of the mass spectrum and possessing sufficient degrees of
freedom to accommodate unforeseen mass features. Failure to do so may
introduce biases in future Hubble constant measurements, especially when
dealing with a large number of detected GW events.

###### Acknowledgements.

The authors are grateful for computational resources provided by the LIGO
Laboratory and supported by the National Science Foundation Grants PHY-0757058
and PHY-0823459. We thank Viola Sordini for the helpful comments and
discussions on the manuscript. MM acknowledges financial support from the
European Research Council for the ERC Consolidator grant DEMOBLACK, under
contract no. 770017, and from the German Excellence Strategy via the
Heidelberg Cluster of Excellence (EXC 2181 - 390900948) STRUCTURES.

##  Appendix A Mass and redshift population models

###  A.1 Source mass models

In this appendix, we present in more details the three phenomenological source
mass models used in this paper for the spectral sirens inference. These models
are the current mass models used by the LVK scientific collaboration for
cosmology results [3, 4, 2]. We give the mathematical expressions of priors
and probability density functions for the : Power Law plus Peak, the Broken
Power Law and the Multi Peak model. Note that they are the same models
described in [9]. All three models are made of combinations of Truncated Power
Law and Truncated Gaussian distributions.

The usual truncated Power Law distribution is given by :

$\rm\mathcal{P}(x|x_{min},x_{max},\alpha)=\begin{cases}\rm\frac{1}{N}x^{\alpha},&\left(\rm
x_{\rm min}\leqslant x\rm\leqslant x_{\rm max}\right),\\\ \rm
0,&\mathrm{otherwise},\end{cases}$ (14)

where $\rm N$ is the normalization factor defined as :

$\rm N=\begin{cases}\rm ln(\frac{x_{\rm max}}{x_{\rm min}}),&\left(\rm
if\hskip 4.26773pt\alpha=-1\right),\\\ \rm\frac{1}{\alpha+1}(x_{\rm
max}^{\alpha+1}-x_{\rm min}^{\alpha+1}),&\mathrm{otherwise}.\end{cases}$ (15)

The Truncated Gaussian distribution with mean $\mu$ and standard deviation
$\sigma$:

$\mathcal{G}_{[a,b]}(x|\mu,\sigma)=\begin{cases}\frac{1}{N}\frac{1}{\sigma\sqrt{2\pi}}\exp{\left[-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right]},&a\leqslant
x\rm\leqslant b,\\\ \rm 0,&\mathrm{otherwise},\end{cases}$ (16)

where $\rm N$ is also the normalization factor expressed as :

$\rm N=\int^{\rm a}_{\rm
b}\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}dx.$ (17)

The Broken Power Law (BPL): This model introduced in [29], is made of two
truncated Power Law attached at a breaking point b,

$\rm b=m_{min}+(m_{max}-m_{min})f,$ (18)

where $f$ belongs to $[0,1]$. When $f=0$, the breaking point b is equal to
$\rm m_{min}$. The probability density functions for $\rm m_{1}^{s}$ and $\rm
m_{2}^{s}$ is then defined as

$\rm\pi(m_{1}^{s}|m_{min},m_{max},\alpha)=\frac{1}{N}\left[\mathcal{P}(m_{1}^{s}|m_{min},b,-\alpha_{1})+\frac{\mathcal{P}(b|m_{min},b,-\alpha_{1})}{\mathcal{P}(b|b,m_{max},-\alpha_{2})}\mathcal{P}(m_{1}^{s}|m_{max},b,-\alpha_{2})\right],$
(19)
$\rm\pi(m_{2}^{s}|m_{min},m_{1}^{s},\beta)=\mathcal{P}(m_{2}^{s}|m_{min},\beta,m_{1}^{s}).$
(20)

N is the normalization factor defined as :

$\rm
N=1+\frac{\mathcal{P}(b|m_{min},b,-\alpha_{1})}{\mathcal{P}(b|b,m_{max},-\alpha_{2})}.$
(21)

Tab. 1 presents the parameters that control the BPL distribution, as well as
the prior distributions used to perform the joint inference of the cosmology
and the BBH population presented in Sec. IV.

Parameter |  Description |  |  Prior
---|---|---|---
$\alpha_{1}$ |  Power Law index number 1 primary mass. |  |  $\mathcal{U}$($-4$, $10$)
$\alpha_{2}$ |  Power Law index number 2 primary mass. |  |  $\mathcal{U}$($-4$, $10$)
$\beta$ |  Power Law index secondary mass. |  |  $\mathcal{U}$($-4,10$)
$\rm m_{min}$ |  Minimum value of the source mass $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($1\rm{M}_{\odot}$, $10\rm{M}_{\odot}$)
$\rm m_{max}$ |  Maximum value of the source mass $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($100\rm{M}_{\odot}$, $200\rm{M}_{\odot}$)
$\delta_{\rm m}$ |  Smoothing parameter $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($0\rm{M}_{\odot}$, $10\rm{M}_{\odot}$)
$\rm b$ |  Breaking point $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($0\rm{M}_{\odot}$, $1\rm{M}_{\odot}$)
Table 1: Summary of the prior ranges used with the BPL mass models for the
spectral sirens analysis.

The Power Law plus Peak (PLP): This model composed of a Truncated Power Law
and a Gaussian component was first introduced by [45]. The probability density
functions of the primary masses are given by :

$\rm\pi(m_{1}^{s}|m_{min},m_{max},\alpha)=(1-\lambda)\mathcal{P}(m_{1}^{s}|m_{min},m_{max},-\alpha)+\lambda\mathcal{G}(m_{1}^{s}|\mu_{g},\sigma_{g}),(0\leqslant\lambda\rm\leqslant
1),$ (22)
$\rm\pi(m_{2}^{s}|m_{min},m_{max},\alpha)=\mathcal{P}(m_{2}^{s}|m_{1}^{s},m_{min},\beta).$
(23)

Tab. 2 presents the parameters that control the PLP distribution, as well as
the prior distributions used to perform the joint inference of the cosmology
and the BBH population presented in Sec. IV.

Parameter |  Description |  |  Prior
---|---|---|---
$\alpha$ |  Spectral index for the PL of the primary mass distribution. |  |  $\mathcal{U}$($1$, $10$)
$\beta$ |  Spectral index for the PL of the mass ratio distribution. |  |  $\mathcal{U}$($-4$, $10$)
$\rm m_{min}$ |  Minimum mass of the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($1\rm{M}_{\odot}$, $12\rm{M}_{\odot}$)
$\rm m_{max}$ |  Maximum mass of the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($50\rm{M}_{\odot}$, $200\rm{M}_{\odot}$)
$\lambda_{\rm g}$ |  Fraction of the model in the Gaussian component. |  |  $\mathcal{U}$($0$, $1$)
$\mu_{\rm g}$ |  Mean of the Gaussian in the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($10\rm{M}_{\odot}$, $40\rm{M}_{\odot}$)
$\sigma_{\rm g}$ |  Width of the Gaussian in the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($6\rm{M}_{\odot}$, $17\rm{M}_{\odot}$)
$\delta_{m}$ |  Range of mass tapering at the lower end of the mass distribution $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($0\rm{M}_{\odot}$, $12\rm{M}_{\odot}$)
Table 2: Summary of the prior ranges used with the PLP mass models for the
spectral sirens analysis.

The Multi Peak (MLTP): This model is an extension of the PLP model, it
consists of a Truncated Power Law and two Gaussian components. This
parameterization was first applied in [29] and the primary masses PDF are
given by :

$\rm\pi(m_{1}^{s}|m_{min},m_{max},\alpha)=\left[(1-\lambda)\mathcal{P}(m_{1}^{s}|m_{min},m_{max},-\alpha)+\lambda\lambda_{low}\mathcal{G}(m_{1}^{s}|\mu_{g}^{low},\sigma_{g}^{low})+\lambda(1-\lambda_{low})\mathcal{G}(m_{1}^{s}|\mu_{g}^{high},\sigma_{g}^{high})\right],$
(24)
$\rm\pi(m_{2}^{s}|m_{min},m_{max},\beta)=\mathcal{P}(m_{2}^{s}|m_{1}^{s},m_{min},\beta).$
(25)

Tab. 3 presents the parameters that control the MLTP distribution, as well as
the prior distributions used to perform the joint inference of the cosmology
and the BBH population presented in Sec. IV.

Parameter |  Description |  |  Prior
---|---|---|---
$\alpha$ |  Power Law index primary mass. |  |  $\mathcal{U}$($1$, $10$)
$\beta$ |  Power Law index secondary mass. |  |  $\mathcal{U}$($-4$, $10$)
$\rm m_{min}$ |  Minimum value of the source mass $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($1\rm{M}_{\odot}$, $10\rm{M}_{\odot}$)
$\rm m_{max}$ |  Maximum value of the source mass $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($50\rm{M}_{\odot}$, $200\rm{M}_{\odot}$)
$\delta_{\rm m}$ |  Smoothing parameter $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($0\rm{M}_{\odot}$, $10\rm{M}_{\odot}$)
$\rm\mu_{g}^{low}$ |  Mean of the lower gaussian peak $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($11\rm{M}_{\odot}$, $30\rm{M}_{\odot}$)
$\rm\mu_{g}^{high}$ |  Mean of the higher gaussian peak $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($40\rm{M}_{\odot}$, $80\rm{M}_{\odot}$)
$\rm\sigma_{g}^{low}$ |  S.t.d of the higher gaussian peak $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($6\rm{M}_{\odot}$, $17\rm{M}_{\odot}$)
$\rm\sigma_{g}^{high}$ |  S.t.d of the higher gaussian peak $[\rm{M}_{\odot}]$. |  |  $\mathcal{U}$($6\rm{M}_{\odot}$, $17\rm{M}_{\odot}$)
$\rm\lambda_{g}$ |  Proportion of events in the peaks. |  |  $\mathcal{U}$($0$, $1$)
$\rm\lambda_{g}^{low}$ |  Proportion of events in the lower peak. |  |  $\mathcal{U}$($0,1$)
Table 3: Summary of the prior ranges used with the MLTP mass models for the
spectral sirens analysis.

All of the three mass models described above are normalized functions, since
they are all made of normalized PDF.

###  A.2 Merger rate evolution model

In this paper, we use a specific parameterization of the CBC merger rate
$\mathcal{R}(z)$, also denoted $p_{pop}(z|\Lambda)$ in Sec. II as a function
of the redshift. This model is constructed after the Madau $\&$ Dickinson star
formation rate [28]. Using the same notation as in Eq. LABEL:eq:cbc_rate_full,
the CBC merger rate is expressed as :

$\rm
p_{pop}(z|\Lambda)=\left[1+(1+z_{p})^{-\gamma-k}\right]\frac{(1+z)^{\gamma}}{1+\left(\frac{1+z}{1+z_{p}}\right)^{\gamma+k}}.$
(26)

Tab. 1 presents the parameters that control the CBC merger rate distribution,
as well as the prior distributions used to perform the joint inference of the
cosmology and the BBH population presented in Sec. IV.

Parameter |  Description |  |  Prior
---|---|---|---
$\gamma$ |  Slope of the power law regime for the rate evolution before the point $z_{p}$ |  |  $\mathcal{U}$($0$, $10$)
$k$ |  Slope of the power law regime for the rate evolution after the point $z_{\rm p}$ |  |  $\mathcal{U}$($0$, $6$)
$z_{p}$ |  Redshift turning point between the power law regimes with $\gamma$ and $k$ |  |  $\mathcal{U}$($0$, $4$)
Table 4:  Summary of the prior ranges used with the CBC merger rate for the
spectral sirens analysis.

##  Appendix B Injected BBHs populations

This appendix presents the values for each of the mass and CBC merger rate
parameters, chosen to perform the spectral sirens analysis of Sec. III, for
each of the following mass models : BPL, PLP and MLTP.

###  B.1 Mass parameters simulated

####  B.1.1 Broken power law parameters

Parameter |  Description |  |  Injected value
---|---|---|---
$\alpha_{1}$ |  Power Law index number 1 primary mass. |  |  1.5
$\alpha_{2}$ |  Power Law index number 2 primary mass. |  |  5.5
$\beta$ |  Power Law index secondary mass. |  |  1.4
$\rm m_{min}$ |  Minimum value of the source mass $[\rm{M}_{\odot}]$. |  |  $5\rm{M}_{\odot}$
$\rm m_{max}$ |  Maximum value of the source mass $[\rm{M}_{\odot}]$. |  |  $100\rm{M}_{\odot}$
$\delta_{\rm m}$ |  Smoothing parameter $[\rm{M}_{\odot}]$. |  |  $5\rm{M}_{\odot}$
$\rm b$ |  Breaking point $[\rm{M}_{\odot}]$. |  |  $0.4\rm{M}_{\odot}$
Table 5: Summary of the prior ranges used with the BPL mass models for the
spectral sirens analysis.

####  B.1.2 Power Law plus peak parameters

Parameter |  Description |  |  Injected value
---|---|---|---
$\alpha$ |  Spectral index for the PL of the primary mass distribution. |  |  2
$\beta$ |  Spectral index for the PL of the mass ratio distribution. |  |  1
$\rm m_{min}$ |  Minimum mass of the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $5\rm{M}_{\odot}$
$\rm m_{max}$ |  Maximum mass of the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $100\rm{M}_{\odot}$
$\lambda_{\rm g}$ |  Fraction of the model in the Gaussian component. |  |  0.1
$\mu_{\rm g}$ |  Mean of the Gaussian in the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $35\rm{M}_{\odot}$
$\sigma_{\rm g}$ |  Width of the Gaussian in the primary mass distribution $[\rm{M}_{\odot}]$. |  |  $5\rm{M}_{\odot}$
$\delta_{m}$ |  Range of mass tapering at the lower end of the mass distribution $[\rm{M}_{\odot}]$. |  |  $5\rm{M}_{\odot}$
Table 6: Summary of the values injected to construct the PLP mass distribution
of the BBH populations analysed in Sec. III.

####  B.1.3 Multi Peak parameters

Parameter |  Description |  |  Injected value
---|---|---|---
$\alpha$ |  Power Law index primary mass. |  |  2
$\beta$ |  Power Law index secondary mass. |  |  1
$\rm m_{min}$ |  Minimum value of the source mass $[\rm{M}_{\odot}]$. |  |  $5\rm{M}_{\odot}$
$\rm m_{max}$ |  Maximum value of the source mass $[\rm{M}_{\odot}]$. |  |  $100\rm{M}_{\odot}$
$\delta_{\rm m}$ |  Smoothing parameter $[\rm{M}_{\odot}]$. |  |  $5\rm{M}_{\odot}$
$\rm\mu_{g}^{low}$ |  Mean of the lower gaussian peak $[\rm{M}_{\odot}]$. |  |  $16\rm{M}_{\odot}$
$\rm\mu_{g}^{high}$ |  Mean of the higher gaussian peak $[\rm{M}_{\odot}]$. |  |  $50\rm{M}_{\odot}$
$\rm\sigma_{g}^{low}$ |  S.t.d of the higher gaussian peak $[\rm{M}_{\odot}]$. |  |  $8\rm{M}_{\odot}$
$\rm\sigma_{g}^{high}$ |  S.t.d of the higher gaussian peak $[\rm{M}_{\odot}]$. |  |  $8\rm{M}_{\odot}$
$\rm\lambda_{g}^{high}$ |  Proportion of events in the peaks. |  |  0.7
$\rm\lambda_{g}^{low}$ |  Proportion of events in the lower peak. |  |  0.8
Table 7: Summary of the values injected to construct the MLTP mass
distribution of the BBH populations analysed in Sec. III.

###  B.2 Rate parameters simulated

Parameter |  Description |  |  Injected value
---|---|---|---
$\gamma$ |  Slope of the power law regime for the rate evolution before the point $z_{p}$ |  |  3
$k$ |  Slope of the power law regime for the rate evolution after the point $z_{\rm p}$ |  |  3
$z_{p}$ |  Redshift turning point between the power law regimes with $\gamma$ and $k$ |  |  2
Table 8:  Summary of the values injected to construct the CBC merger rate of
the BBH populations analysed in Sec. III.

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2023a] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run, Phys. Rev. X 13, 041039 (2023a), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO), GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, arXiv (2021b), arXiv:2108.01045 [gr-qc] .
  * Freedman [2021] W. L. Freedman, Measurements of the Hubble Constant: Tensions in Perspective, Astrophys. J. 919, 16 (2021), arXiv:2106.15656 [astro-ph.CO] .
  * Ezquiaga and Holz [2022] J. M. Ezquiaga and D. E. Holz, Spectral Sirens: Cosmology from the Full Mass Distribution of Compact Binaries, Phys. Rev. Lett. 129, 061102 (2022), arXiv:2202.08240 [astro-ph.CO] .
  * Mastrogiovanni  _et al._ [2021] S. Mastrogiovanni, K. Leyde, C. Karathanasis, E. Chassande-Mottin, D. A. Steer, J. Gair, A. Ghosh, R. Gray, S. Mukherjee, and S. Rinaldi, On the importance of source population models for gravitational-wave cosmology, Phys. Rev. D 104, 062009 (2021), arXiv:2103.14663 [gr-qc] .
  * Mastrogiovanni  _et al._ [2023a] S. Mastrogiovanni, D. Laghi, R. Gray, G. C. Santoro, A. Ghosh, C. Karathanasis, K. Leyde, D. A. Steer, S. Perries, and G. Pierra, A novel approach to infer population and cosmological properties with gravitational waves standard sirens and galaxy surveys, arXiv (2023a), arXiv:2305.10488 [astro-ph.CO] .
  * Mastrogiovanni  _et al._ [2023b] S. Mastrogiovanni, G. Pierra, S. Perriès, D. Laghi, G. C. Santoro, A. Ghosh, R. Gray, C. Karathanasis, and K. Leyde, ICAROGW: A python package for inference of astrophysical population properties of noisy, heterogeneous and incomplete observations, arXiv (2023b), arXiv:2305.17973 [astro-ph.CO] .
  * Taylor  _et al._ [2012] S. R. Taylor, J. R. Gair, and I. Mandel, Hubble without the Hubble: Cosmology using advanced gravitational-wave detectors alone, Phys. Rev. D 85, 023535 (2012), arXiv:1108.5161 [gr-qc] .
  * Wysocki  _et al._ [2019] D. Wysocki, J. Lange, and R. O’Shaughnessy, Reconstructing phenomenological distributions of compact binaries via gravitational wave observations, Phys. Rev. D 100, 043012 (2019), arXiv:1805.06442 [gr-qc] .
  * Farr  _et al._ [2019] W. M. Farr, M. Fishbach, J. Ye, and D. Holz, A Future Percent-Level Measurement of the Hubble Expansion at Redshift 0.8 With Advanced LIGO, Astrophys. J. Lett. 883, L42 (2019), arXiv:1908.09084 [astro-ph.CO] .
  * Mancarella  _et al._ [2022] M. Mancarella, E. Genoud-Prachex, and M. Maggiore, Cosmology and modified gravitational wave propagation from binary black hole population models, Phys. Rev. D 105, 064030 (2022), arXiv:2112.05728 [gr-qc] .
  * Mukherjee [2022] S. Mukherjee, The redshift dependence of black hole mass distribution: is it reliable for standard sirens cosmology?, Mon. Not. Roy. Astron. Soc. 515, 5495 (2022), arXiv:2112.10256 [astro-ph.CO] .
  * Gray [2021] R. Gray,  _Gravitational wave cosmology : measuring the Hubble constant with dark standard sirens_ , Ph.D. thesis, University of Glasgow (2021).
  * Gray  _et al._ [2022] R. Gray, C. Messenger, and J. Veitch, A pixelated approach to galaxy catalogue incompleteness: improving the dark siren measurement of the Hubble constant, Mon. Not. Roy. Astron. Soc. 512, 1127 (2022), arXiv:2111.04629 [astro-ph.CO] .
  * Gair  _et al._ [2022] J. R. Gair _et al._ , The Hitchhiker’s guide to the galaxy catalog approach for gravitational wave cosmology, arXiv (2022), arXiv:2212.08694 [gr-qc] .
  * Leyde  _et al._ [2022] K. Leyde, S. Mastrogiovanni, D. A. Steer, E. Chassande-Mottin, and C. Karathanasis, Current and future constraints on cosmology and modified gravitational wave friction from binary black holes, JCAP 09, 012, arXiv:2202.00025 [gr-qc] .
  * Schutz [1986] B. F. Schutz, Determining the Hubble Constant from Gravitational Wave Observations, Nature 323, 310 (1986).
  * Fishbach  _et al._ [2018] M. Fishbach, D. E. Holz, and W. M. Farr, Does the Black Hole Merger Rate Evolve with Redshift?, Astrophys. J. Lett. 863, L41 (2018), arXiv:1805.10270 [astro-ph.HE] .
  * Ade  _et al._ [2016] P. A. R. Ade _et al._ (Planck), Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (KAGRA, LIGO Scientific, Virgo, VIRGO), Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Rev. Rel. 21, 3 (2018), arXiv:1304.0670 [gr-qc] .
  * Dominik  _et al._ [2015] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik, and F. Pannarale, Double Compact Objects III: Gravitational Wave Detection Rates, Astrophys. J. 806, 263 (2015), arXiv:1405.7016 [astro-ph.HE] .
  * Moore  _et al._ [2016] B. Moore, M. Favata, K. G. Arun, and C. K. Mishra, Gravitational-wave phasing for low-eccentricity inspiralling compact binaries to 3PN order, Phys. Rev. D 93, 124061 (2016), arXiv:1605.00304 [gr-qc] .
  * Vitale  _et al._ [2020] S. Vitale, D. Gerosa, W. M. Farr, and S. R. Taylor, Inferring the properties of a population of compact binaries in presence of selection effects, arXiv 10.1007/978-981-15-4702-7_45-1 (2020), arXiv:2007.05579 [astro-ph.IM] .
  * Abbott  _et al._ [2023b] R. Abbott _et al._ (LIGO Scientific, Virgo,, KAGRA, VIRGO), Constraints on the Cosmic Expansion History from GWTC–3, Astrophys. J. 949, 76 (2023b), arXiv:2111.03604 [astro-ph.CO] .
  * Gray  _et al._ [2023] R. Gray _et al._ , Joint cosmological and gravitational-wave population inference using dark sirens and galaxy catalogues, arXiv (2023), arXiv:2308.02281 [astro-ph.CO] .
  * Madau and Dickinson [2014] P. Madau and M. Dickinson, Cosmic Star Formation History, Ann. Rev. Astron. Astrophys. 52, 415 (2014), arXiv:1403.0007 [astro-ph.CO] .
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, Virgo), Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog, Astrophys. J. Lett. 913, L7 (2021c), arXiv:2010.14533 [astro-ph.HE] .
  * Farmer  _et al._ [2019] R. Farmer, M. Renzo, S. E. de Mink, P. Marchant, and S. Justham, Mind the gap: The location of the lower edge of the pair instability supernovae black hole mass gap, arXiv 10.3847/1538-4357/ab518b (2019), arXiv:1910.12874 [astro-ph.SR] .
  * Talbot and Thrane [2018] C. Talbot and E. Thrane, Measuring the binary black hole mass spectrum with an astrophysically motivated parameterization, Astrophys. J. 856, 173 (2018), arXiv:1801.02699 [astro-ph.HE] .
  * Karathanasis  _et al._ [2023a] C. Karathanasis, B. Revenu, S. Mukherjee, and F. Stachurski, GWSim: Python package for creating mock GW samples for different astrophysical populations and cosmological models of binary black holes, A&A 677, A124 (2023a), arXiv:2210.05724 [astro-ph.CO] .
  * Karathanasis  _et al._ [2023b] C. Karathanasis, S. Mukherjee, and S. Mastrogiovanni, Binary black holes population and cosmology in new lights: signature of PISN mass and formation channel in GWTC-3, MNRAS 523, 4539 (2023b), arXiv:2204.13495 [astro-ph.CO] .
  * Mapelli  _et al._ [2021] M. Mapelli _et al._ , Hierarchical black hole mergers in young, globular and nuclear star clusters: the effect of metallicity, spin and cluster properties, Mon. Not. Roy. Astron. Soc. 505, 339 (2021), arXiv:2103.05016 [astro-ph.HE] .
  * Mapelli  _et al._ [2022] M. Mapelli, Y. Bouffanais, F. Santoliquido, M. A. Sedda, and M. C. Artale, The cosmic evolution of binary black holes in young, globular, and nuclear star clusters: rates, masses, spins, and mixing fractions, Mon. Not. Roy. Astron. Soc. 511, 5797 (2022), arXiv:2109.06222 [astro-ph.HE] .
  * Giacobbo  _et al._ [2018] N. Giacobbo, M. Mapelli, and M. Spera, Merging black hole binaries: the effects of progenitor’s metallicity, mass-loss rate and Eddington factor, Mon. Not. Roy. Astron. Soc. 474, 2959 (2018), arXiv:1711.03556 [astro-ph.SR] .
  * Fryer  _et al._ [2012] C. L. Fryer, K. Belczynski, G. Wiktorowicz, M. Dominik, V. Kalogera, and D. E. Holz, Compact Remnant Mass Function: Dependence on the Explosion Mechanism and Metallicity, Astrophys. J. 749, 91 (2012), arXiv:1110.1726 [astro-ph.SR] .
  * Mapelli  _et al._ [2020] M. Mapelli, M. Spera, E. Montanari, M. Limongi, A. Chieffi, N. Giacobbo, A. Bressan, and Y. Bouffanais, Impact of the Rotation and Compactness of Progenitors on the Mass of Black Holes, Astrophys. J. 888, 76 (2020), arXiv:1909.01371 [astro-ph.HE] .
  * Abbott  _et al._ [2023c] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Phys. Rev. X 13, 011048 (2023c), arXiv:2111.03634 [astro-ph.HE] .
  * Belczynski  _et al._ [2001] K. Belczynski, V. Kalogera, and T. Bulik, A Comprehensive study of binary compact objects as gravitational wave sources: Evolutionary channels, rates, and physical properties, Astrophys. J. 572, 407 (2001), arXiv:astro-ph/0111452 .
  * Spera and Mapelli [2017] M. Spera and M. Mapelli, Very massive stars, pair-instability supernovae and intermediate-mass black holes with the SEVN code, Mon. Not. Roy. Astron. Soc. 470, 4739 (2017), arXiv:1706.06109 [astro-ph.SR] .
  * Mapelli  _et al._ [2017] M. Mapelli, N. Giacobbo, E. Ripamonti, and M. Spera, The cosmic merger rate of stellar black hole binaries from the Illustris simulation, Mon. Not. Roy. Astron. Soc. 472, 2422 (2017), arXiv:1708.05722 [astro-ph.GA] .
  * Vitale  _et al._ [2019] S. Vitale, W. M. Farr, K. Ng, and C. L. Rodriguez, Measuring the star formation rate with gravitational waves from binary black holes, Astrophys. J. Lett. 886, L1 (2019), arXiv:1808.00901 [astro-ph.HE] .
  * Renzo  _et al._ [2020] M. Renzo, R. J. Farmer, S. Justham, S. E. de Mink, Y. Götberg, and P. Marchant, Sensitivity of the lower-edge of the pair instability black hole mass gap to the treatment of time dependent convection, Mon. Not. Roy. Astron. Soc. 493, 4333 (2020), arXiv:2002.08200 [astro-ph.SR] .
  * Talbot  _et al._ [2019] C. Talbot, R. Smith, E. Thrane, and G. B. Poole, Parallelized Inference for Gravitational-Wave Astronomy, Phys. Rev. D 100, 043030 (2019), arXiv:1904.02863 [astro-ph.IM] .
  * Turski  _et al._ [2023] C. Turski, M. Bilicki, G. Dálya, R. Gray, and A. Ghosh, Impact of modelling galaxy redshift uncertainties on the gravitational-wave dark standard siren measurement of the Hubble constant, arXiv (2023), arXiv:2302.12037 [gr-qc] .
  * Abbott  _et al._ [2021d] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, arXiv (2021d), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Multi-messenger Observations of a Binary Neutron Star Merger, Astrophys. J. Lett. 848, L12 (2017a), arXiv:1710.05833 [astro-ph.HE] .
  * Coulter  _et al._ [2017] D. A. Coulter _et al._ , Swope Supernova Survey 2017a (SSS17a), the Optical Counterpart to a Gravitational Wave Source, Science 358, 1556 (2017), arXiv:1710.05452 [astro-ph.HE] .
  * Mapelli [2021] M. Mapelli, Formation Channels of Single and Binary Stellar-Mass Black Holes (2021) arXiv:2106.00699 [astro-ph.HE] .
  * Veitch  _et al._ [2015] J. Veitch _et al._ , Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library, Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Bethe and Brown [1998] H. A. Bethe and G. E. Brown, Evolution of binary compact objects which merge, Astrophys. J. 506, 780 (1998), arXiv:astro-ph/9802084 .
  * Mehta  _et al._ [2022] A. K. Mehta, A. Buonanno, J. Gair, M. C. Miller, E. Farag, R. J. deBoer, M. Wiescher, and F. X. Timmes, Observing Intermediate-mass Black Holes and the Upper Stellar-mass gap with LIGO and Virgo, Astrophys. J. 924, 39 (2022), arXiv:2105.06366 [gr-qc] .
  * Babak  _et al._ [2013] S. Babak _et al._ , Searching for gravitational waves from binary coalescence, Phys. Rev. D 87, 024033 (2013), arXiv:1208.3491 [gr-qc] .
  * Fishbach and Holz [2020] M. Fishbach and D. E. Holz, Picky Partners: The Pairing of Component Masses in Binary Black Hole Mergers, Astrophys. J. Lett. 891, L27 (2020), arXiv:1905.12669 [astro-ph.HE] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), A gravitational-wave standard siren measurement of the Hubble constant, Nature 551, 85 (2017b), arXiv:1710.05835 [astro-ph.CO] .
  * Chatterjee  _et al._ [2021] D. Chatterjee, A. Hegade K R, G. Holder, D. E. Holz, S. Perkins, K. Yagi, and N. Yunes, Cosmology with Love: Measuring the Hubble constant using neutron star universal relations, Phys. Rev. D 104, 083528 (2021), arXiv:2106.06589 [gr-qc] .
  * Messenger and Read [2012] C. Messenger and J. Read, Measuring a cosmological distance-redshift relationship using only gravitational wave observations of binary neutron star coalescences, Phys. Rev. Lett. 108, 091101 (2012), arXiv:1107.5725 [gr-qc] .
  * Woosley [2017] S. E. Woosley, Pulsational Pair-Instability Supernovae, Astrophys. J. 836, 244 (2017), arXiv:1608.08939 [astro-ph.HE] .
  * Marchant  _et al._ [2021] P. Marchant, K. M. W. Pappas, M. Gallegos-Garcia, C. P. L. Berry, R. E. Taam, V. Kalogera, and P. Podsiadlowski, The role of mass transfer and common envelope evolution in the formation of merging binary black holes, Astron. Astrophys. 650, A107 (2021), arXiv:2103.09243 [astro-ph.SR] .