# Cosmology and Astrophysics with Standard Sirens and Galaxy Catalogs in View
of Future Gravitational Wave Observations

Nicola Borghi Dipartimento di Fisica e Astronomia “Augusto Righi”–Università
di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy  INAF -
Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Piero
Gobetti 93/3, I-40129 Bologna, Italy  Michele Mancarella Dipartimento di
Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della
Scienza 3, 20126 Milano, Italy  INFN, Sezione di Milano-Bicocca, Piazza della
Scienza 3, 20126 Milano, Italy  Michele Moresco Dipartimento di Fisica e
Astronomia “Augusto Righi”–Università di Bologna, via Piero Gobetti 93/2,
I-40129 Bologna, Italy  INAF - Osservatorio di Astrofisica e Scienza dello
Spazio di Bologna, via Piero Gobetti 93/3, I-40129 Bologna, Italy  Matteo
Tagliazucchi Dipartimento di Fisica e Astronomia “Augusto Righi”–Università di
Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy  INAF - Osservatorio
di Astrofisica e Scienza dello Spazio di Bologna, via Piero Gobetti 93/3,
I-40129 Bologna, Italy  Francesco Iacovelli Département de Physique Théorique,
Université de Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland
Gravitational Wave Science Center (GWSC), Université de Genève, CH-1211
Genève, Switzerland  Andrea Cimatti Dipartimento di Fisica e Astronomia
“Augusto Righi”–Università di Bologna, via Piero Gobetti 93/2, I-40129
Bologna, Italy  INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5,
I-50125, Firenze, Italy  Michele Maggiore Département de Physique Théorique,
Université de Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland
Gravitational Wave Science Center (GWSC), Université de Genève, CH-1211
Genève, Switzerland Nicola Borghi<EMAIL_ADDRESS>
###### Abstract

With the growing number of gravitational-wave detections and the advent of
large galaxy redshift surveys, a new era in cosmology is unfolding. This study
explores the synergies between gravitational waves and galaxy surveys to
jointly constrain cosmological and gravitational-wave population parameters.
We introduce CHIMERA , a novel code for gravitational-wave cosmology combining
information from the population properties of compact binary mergers and
galaxy catalogs. We study constraints for scenarios representative of LIGO-
Virgo-KAGRA O4 and O5 observing runs, assuming to have a complete catalog of
potential host galaxies with either spectroscopic or photometric redshift
measurements. We find that a percent-level measurement of $H_{0}$ could be
achieved with the best 100 binary black holes in O5 using a spectroscopic
galaxy catalog. In this case, the intrinsic correlation that exists between
$H_{0}$ and the binary black hole population mass scales is broken. Instead,
by using a photometric catalog the accuracy on $H_{0}$ is degraded up to a
factor of $\sim\\!9$, leaving a significant correlation between $H_{0}$ and
the mass scales that must be carefully modeled to avoid bias. Interestingly,
we find that using spectroscopic redshift measurements in the O4 configuration
yields a better constraint on $H_{0}$ compared to the O5 configuration with
photometric measurements. In view of the wealth of gravitational-wave data
that will be available in the future, we argue the importance of obtaining
spectroscopic galaxy catalogs to maximize the scientific return of
gravitational-wave cosmology.

††software: GWFAST (Iacovelli et al., 2022b), MGCosmoPop (Mancarella et al.,
2022), ChainConsumer (Hinton, 2016), numpy (Harris et al., 2020), matplotlib
(Hunter, 2007).

revtex4-1.clsrevtex4-2.cls

##  1 Introduction

Gravitational waves (GWs) emitted by merging compact binaries can serve as
“standard sirens” (Schutz, 1986) as their signal provides a direct measurement
of the luminosity distance. This makes GWs a new and powerful cosmological
probe for tracing the expansion history of the Universe and testing the
validity of General Relativity on cosmological scales. First of all, the
measurement of the expansion history is a longstanding challenge in the field.
The discovery of a discrepancy between the value of the Hubble constant
$H_{0}$ measured from the local distance ladder (e.g. Riess et al., 2022) and
the one inferred from the cosmic microwave background assuming a $\Lambda$CDM
cosmology (Planck Collaboration, 2020) motivated a rising interest within the
community in developing alternative and emerging cosmological probes (see
Moresco et al., 2022, for a recent review). Besides being self-calibrated
distance tracers, GWs are a particularly promising candidate due to the
expected rapid increase in the rate and redshift range of the available
catalogs with current and next-generation detectors (Baibhav et al., 2019;
Iacovelli et al., 2022a; Branchesi et al., 2023). While the measurement of
$H_{0}$ is driven by low-redshift events, those at redshifts $\mathcal{O}(1)$
and beyond open the possibility of testing the distance–redshift relation in
full generality. This allows to constrain the dark energy sector, notably the
phenomenon of “modified GW propagation” (Amendola et al., 2018; Belgacem et
al., 2018a, b; Lagos et al., 2019), which is a general prediction of models
modifying General Relativity at cosmological scales (Belgacem et al., 2019).

The peculiarity of GW cosmology is that determining the redshift with GW data
alone is not possible because of its inherent degeneracy with binary masses.
External information is required to provide cosmological constraints111An
exception is provided by binary neutron star mergers, which carry a direct
redshift signature through the effect of tidal couplings in the inspiral phase
(Messenger & Read, 2012; Messenger et al., 2014). This will be measurable with
next-generation ground-based detectors (Dhani et al., 2022)..

A direct measurement of the redshift is possible in the case of a “bright
siren”, i.e. an event with the coincident detection of an electromagnetic (EM)
counterpart and the consequent identification of the host galaxy which allows
the redshift to be directly measured from spectroscopy (Schutz, 1986; Holz &
Hughes, 2005; Nissanke et al., 2010; Chen et al., 2018). Such events are rare,
as they typically require mergers involving at least one neutron star and a
coincident detection of the EM emission. So far, only the Binary Neutron Star
(BNS) event GW170817 has an associated counterpart, the gamma-ray burst GRB
170817A (Abbott et al., 2017a, b), which led to the identification of the host
galaxy and the first standard siren measurement of $H_{0}$ (Abbott et al.,
2017c), while being a too low redshift to provide a stringent constraint on
modified GW propagation (Belgacem et al., 2018b). Detection rates for the
ongoing and upcoming observing runs of the LIGO-Virgo-Kagra (LVK)
collaboration are quite uncertain but in general not optimistic (Colombo et
al., 2022, 2023).

The largest fraction of events is instead made of Binary Black Holes (BBH)
without EM counterparts, known as “dark sirens”. Using these sources for
cosmology has recently become a concrete option with the latest LVK data
(Abbott et al., 2021b, c). Two techniques have been primarily employed. The
first consists in statistically inferring the redshift from potential host
galaxies within the GW localization volume (Schutz, 1986; Del Pozzo, 2012;
Chen et al., 2018; Gray et al., 2020; Gair et al., 2023)222A more refined
approach involves matching the spatial clustering of GWs and galaxies (Oguri,
2016; Nair et al., 2018; Yu et al., 2020; Vijaykumar et al., 2023b; Bera et
al., 2020; Mukherjee et al., 2021).. This led to measurements of the Hubble
constant from LVK data in conjunction with the GLADE (Dálya et al., 2018) and
GLADE+ (Dálya et al., 2022) galaxy catalogs (Fishbach et al., 2019; Abbott et
al., 2021a; Finke et al., 2021; Gray et al., 2022; Abbott et al., 2023a), as
well as DES (Soares-Santos et al., 2019; Palmese et al., 2020), the DESI
Legacy Survey (Palmese et al., 2023), the DECam Local Volume Exploration
Survey (DELVE) (Alfradique et al., 2023) and DESI (Ballard et al., 2023). This
method also led to the first dark siren bound on modified gravitational wave
propagation (Finke et al., 2021). Forecasts indicate that this technique can
lead to percent-level measurements with future GW experiments such as LISA
(Laghi et al., 2021; Liu et al., 2023), the Einstein Telescope, and Cosmic
Explorer (Muttoni et al., 2022, 2023).

Alternatively, the degeneracy between mass and redshift can be broken by
modeling intrinsic astrophysical properties as the source-frame mass
distribution. In particular, this is possible with “spectral sirens” - sources
whose source-frame mass distribution contains features such as breaks, peaks
or changes in slope (Chernoff & Finn, 1993; Taylor et al., 2012; Farr et al.,
2019; Mastrogiovanni et al., 2021a; Ezquiaga & Holz, 2022). This method
allowed in particular to obtain constraints on $H_{0}$ (Abbott et al., 2023a)
and on modified GW propagation (Ezquiaga, 2021; Mancarella et al., 2022;
Magana Hernandez, 2021; Leyde et al., 2022) from the presence of a mass scale
in the BBH distribution (Abbott et al., 2021d, 2023b).

Multiple pipelines have been publicly released, both for the correlation with
galaxy catalogs (gwcosmo, Gray et al. 2020; DarkSirensStat, Finke et al. 2021;
cosmolisa, Laghi & Del Pozzo 2020) and for spectral sirens (icarogw,
Mastrogiovanni et al. 2021b; MGCosmoPop, Mancarella et al. 2022). However, it
is known that these two techniques are, in fact, special cases of a unique
concept: incorporating prior knowledge about the population in a comprehensive
analysis of the ensemble of GW sources (Mastrogiovanni et al., 2021a; Finke et
al., 2021; Abbott et al., 2023a; Moresco et al., 2022). Since population
properties are known in the source frame (mass and redshift), while GW signals
are only sensitive to detector frame quantities (detector-frame mass and
luminosity distance), the cosmology-dependent mapping between the two allows
to break the mass–redshift degeneracy statistically. In the two aforementioned
cases, the population prior used has been either a redshift distribution
constructed from a catalog of galaxies, or a model of the source-frame mass
distribution. It is natural to seek a unified approach. Moreover, it has been
largely recognized that the correlation with a galaxy catalog actually relies
on assumptions about the population model for GW sources, and in particular of
a source-frame mass, redshift, and spin distribution (Mastrogiovanni et al.,
2021a; Finke et al., 2021; Abbott et al., 2023a). These are needed in
particular to consistently account for selection effects, i.e. the fact that
not all GW sources are equally likely to be observed. Moreover, an implicit
assumption on the mass, redshift, and spin distributions is inherently present
in the form of Bayesian priors adopted in the analysis of individual events,
and some assumption has also to be used when accounting for the
(in-)completeness of the galaxy catalog. These two facts actually make
mandatory a joint analysis of the population properties of compact objects and
of the cosmological parameters. Doing so while also incorporating the galaxy
catalog information is however particularly costly, and until very recently,
the population properties of GW sources (and in particular their source-frame
mass distribution) have been kept fixed when including the galaxy catalog
data.

While this paper was in preparation, two pipelines implementing the joint
population and cosmological inference including a galaxy catalog were
presented in Mastrogiovanni et al. (2023) and Gray et al. (2023), namely
icarogw2.0 and gwcosmo2.0, building on the homonyms aforementioned codes.
Their application to the GWTC-3 catalog (Abbott et al., 2021c) with the GLADE+
galaxy catalog (Dálya et al., 2022) led to updated and more robust
measurements of the Hubble constant, as well as of the parameters describing
modified GW propagation (the latter being further constrained in Chen et al.
2023), and of alternative cosmological models (Raffai et al., 2023). A
machine-learning-based approach has also been proposed by Stachurski et al.
(2023).

In this paper, we present and release CHIMERA333Available at
https://github.com/CosmoStatGW/CHIMERA (Combined Hierarchical Inference Model
for Electromagnetic and gRavitational-wave Analysis), a novel, independent
Python code for the joint analysis of GW transient catalogs and galaxy
catalogs. CHIMERA is built starting from DarkSirensStat (Finke et al., 2021)
and MGCosmoPop (Mancarella et al., 2022) and allows the joint inference of
cosmology and population properties of GW sources. Currently, the information
that can be extracted from data with the combined method just described is
limited by the lack of a sufficient number of GW events with small enough
localization regions and covered by a complete galaxy catalog. In order to
fully appreciate the gain achievable from this method, it is therefore of
interest to study its potential with upcoming GW data combined with a complete
galaxy catalog. This will become a concrete possibility in the coming years,
thanks to galaxy surveys already in a data-release stage such as DESI (DESI
Collaboration et al., 2016), and in data-taking period such as Euclid
(Laureijs et al., 2011). In this context, it is also important to understand
how large the impact of the redshift measurement uncertainties can be and to
assess the relative importance of the information coming from the galaxy
survey and from the population model respectively. This is the second goal of
the paper, where we explore for the first time the application of this
technique in the context of future GW observations and study in detail the
impact of different galaxy catalog properties. We focus here in particular on
the measurement of the Hubble constant. We simulate GW detections with
detector networks and sensitivities corresponding to the ongoing “O4” run of
the LVK collaboration, as well as of the future “O5” observing run. We study
the constraining power of the best 100 events observed in one year,
considering a complete galaxy catalog with both spectroscopic and photometric
uncertainties on the galaxies’ redshifts (for a study of the limited impact of
different assumptions on the shape of galaxy redshift uncertainties with
current data, see Turski et al. 2023). We also compare to the case where no
information from the galaxy catalog is present.

The paper is organized as follows. In Section 2 we present the statistical
framework and its implementation in CHIMERA. In Section 3 we describe the
generation of mock galaxy and GW catalogs. Section 4 provides results for both
the O4-like and O5-like scenarios. Finally, Section 5 concludes the paper and
discusses future applications of this method.

##  2 Statistical framework

We consider a population of GW sources, individually described by source-frame
parameters $\boldsymbol{\theta}$ which globally follow a probability
distribution $p_{\mathrm{pop}}(\boldsymbol{\theta}\,|\,\boldsymbol{\lambda})$
described by hyperparameters $\boldsymbol{\lambda}$. We summarize here the
hierarchical Bayesian inference formalism to infer $\boldsymbol{\lambda}$
given a catalog of GW detections and a catalog of their potential hosts with
measured redshifts444For a comprehensive introduction, see Gair et al.
(2023).. We separate the hyperparameters describing the underlying cosmology
$\boldsymbol{\lambda}_{\mathrm{c}}$ from the ones describing the astrophysical
population of GW sources, i.e. mass distribution
$\boldsymbol{\lambda}_{\mathrm{m}}$ and redshift distribution
$\boldsymbol{\lambda}_{\mathrm{z}}$, so that
$\boldsymbol{\lambda}=\\{\boldsymbol{\lambda}_{\mathrm{c}},\boldsymbol{\lambda}_{\mathrm{m}},\boldsymbol{\lambda}_{\mathrm{z}}\\}$555The
formalism can be trivially extended to include the spin distribution, which we
do not consider in this work.. Given a set $\boldsymbol{d}^{\rm
GW}=\\{\boldsymbol{d}^{\rm GW}_{i}\\}$ ($i=1,...,N_{\rm ev}$) of data from
independent GW events, the population likelihood is (Mandel et al., 2019;
Vitale et al., 2022),

$p(\boldsymbol{d}^{\rm
GW}\,|\,\boldsymbol{\lambda})\propto\frac{1}{\xi(\boldsymbol{\lambda})^{N_{\rm
ev}}}\prod_{i=1}^{N_{\rm ev}}\int p(\boldsymbol{d}^{\rm
GW}_{i}\,|\,\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\mathrm{c}})\,p_{\mathrm{pop}}(\boldsymbol{\theta}_{i}\,|\,\boldsymbol{\lambda})\,\mathrm{d}\boldsymbol{\theta}_{i}\,,$
(1)

where $p(\boldsymbol{d}^{\rm
GW}_{i}\,|\,\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\mathrm{c}})$ is the
individual source likelihood and

$\xi(\boldsymbol{\lambda})\equiv\int P_{\rm
det}(\boldsymbol{\theta},\boldsymbol{\lambda}_{\mathrm{c}})\,p_{\mathrm{pop}}(\boldsymbol{\theta}\,|\,\boldsymbol{\lambda})\,\mathrm{d}\boldsymbol{\theta}\,,$
(2)

is the selection function, which corrects for selection effects (Loredo, 2004;
Mandel et al., 2019). When
$p_{\mathrm{pop}}(\boldsymbol{\theta}\,|\,\boldsymbol{\lambda})$ is correctly
normalized to unity, $\xi(\boldsymbol{\lambda})$ measures the overall fraction
of detectable events given $\boldsymbol{\lambda}$. $P_{\rm
det}(\boldsymbol{\theta},\boldsymbol{\lambda}_{\mathrm{c}})$ gives the
probability of detecting an event with parameters $\boldsymbol{\theta}$. We
consider the following set of source parameters:
$\boldsymbol{\theta}=\\{z,\hat{\Omega},m_{1},m_{2}\\}$, where $z$ is the
redshift of the binary, $\hat{\Omega}$ the sky localization, and $m_{1}$,
$m_{2}$ are the primary and secondary masses.

Note that, when the cosmological hyperparameters
$\boldsymbol{\lambda}_{\mathrm{c}}$ are included in the analysis, the single-
event GW likelihood written in source-frame variables in Eq. 1 includes an
explicit dependence on $\boldsymbol{\lambda}_{\mathrm{c}}$. This is a
consequence of the fact that GW observations do not provide information on
source-frame parameters, but on detector-frame quantities
$\boldsymbol{\theta}^{\mathrm{det}}=\\{d_{L},\hat{\Omega},m_{1}^{\mathrm{det}},m_{2}^{\mathrm{det}}\\}$.
Those are then mapped to source frame in a way that depends on the underlying
cosmology. Explicitly, the source-to-detector frame conversion requires
inverting the relations

$\displaystyle d_{L}$ $\displaystyle=d_{L}(z;\boldsymbol{\lambda}_{\rm c})$
(3) $\displaystyle m_{1,2}^{\mathrm{det}}$ $\displaystyle=m_{1,2}\,(1+z)\,,$
(4)

where $d_{L}$ is the luminosity distance, which, in a flat $\Lambda$CDM
cosmology, is given by

$d_{L}(z;\boldsymbol{\lambda}_{\rm
c})=c\,(1+z)\int_{0}^{z}\frac{\mathrm{d}z^{\prime}}{H(z^{\prime};\boldsymbol{\lambda}_{\rm
c})}\,,$ (5)

with $\boldsymbol{\lambda}_{\rm c}\equiv\\{H_{0},\Omega_{\rm m,0}\\}$ and the
Hubble parameter being

$H(z;\boldsymbol{\lambda}_{\rm c})=H_{0}\sqrt{\Omega_{\rm
m,0}(1+z)^{3}+(1-\Omega_{\rm m,0})}\,,$ (6)

where the radiation is not considered as its contribution is negligible in the
late Universe. Alternatively, Eq. 1 could be written in the detector frame to
remove the dependence on $\boldsymbol{\lambda}_{\mathrm{c}}$ (see Finke et al.
(2021) for a detailed discussion). In this work, for the actual implementation
we find more convenient to work in source-frame, as it will be clear below -
essentially, this is due to the fact that the galaxy catalog is naturally
given in redshift space. Note that the single-event GW likelihood is a
probability distribution normalized on the data, not on the parameters, so in
the transformation from $\boldsymbol{\theta}$ to
$\boldsymbol{\theta}^{\mathrm{det}}$ it remains unchanged. However, in
practical applications of hierarchical inference, we are provided with a set
of samples drawn from the posterior distribution
$p(\boldsymbol{\theta}^{\mathrm{det}}_{i}\,|\,\boldsymbol{d}^{\rm GW}_{i})$
for each event in the GW catalog, obtained with priors
$\pi(\boldsymbol{\theta}^{\mathrm{det}}_{i})$. It is thus convenient to re-
express the single-event GW likelihood via Bayes’ theorem as
$p(\boldsymbol{d}^{\rm
GW}_{i}\,|\,\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\mathrm{c}})\propto
p(\boldsymbol{\theta}_{i}\,|\,\boldsymbol{d}^{\rm
GW}_{i},\boldsymbol{\lambda}_{\mathrm{c}})/\pi(\boldsymbol{\theta}_{i})$. The
posterior $p(\boldsymbol{\theta}_{i}\,|\,\boldsymbol{d}^{\rm
GW}_{i},\boldsymbol{\lambda}_{\mathrm{c}})$ can then be obtained from the
(cosmology-dependent) conversion of the posterior samples of
$p(\boldsymbol{\theta}^{\mathrm{det}}_{i}\,|\,\boldsymbol{d}^{\rm GW}_{i})$ to
source-frame, while the prior $\pi(\boldsymbol{\theta}_{i})$ is related to the
one in detector-frame by
$\pi(\boldsymbol{\theta}_{i})=\pi(\boldsymbol{\theta}^{\mathrm{det}}_{i})\times|\mathrm{d}\boldsymbol{\theta}^{\mathrm{det}}_{i}/\mathrm{d}\boldsymbol{\theta}_{i}|$.
Explicitly, for the variables of interest here, the Jacobian factors are given
by

$\frac{\mathrm{d}m^{\mathrm{det}}_{1,2}}{\mathrm{d}m_{1,2}}=1+z\,,\quad\frac{\mathrm{d}d_{L}}{\mathrm{d}z}(z;\boldsymbol{\lambda}_{\rm
c})=\frac{d_{L}}{1+z}+\frac{c\,(1+z)}{H(z;\boldsymbol{\lambda}_{\rm c})}\,.$
(7)

Finally, under the assumption that the mass function does not evolve with
cosmic time666In this paper, we adopt a time-independent mass distribution.
More general models should be envisaged in the future in light of recent
results that support evidence of evolution (Rinaldi et al., 2023)., the
population function can be split as follows:

$p_{\mathrm{pop}}(\boldsymbol{\theta}\,|\,\boldsymbol{\lambda})=p(m_{1},m_{2}\,|\,\boldsymbol{\lambda}_{\rm
m})\,p(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm
c},\boldsymbol{\lambda}_{\rm z})\,,$ (8)

where $\hat{\Omega}$ denotes the direction in the sky. In the subsequent
sections, we will describe in detail the construction of the population prior
taking into account a galaxy catalog and the source-frame mass distribution of
the population of compact binaries.

###  2.1 Redshift prior

The core of the method is the construction of a population prior on the GW
source redshift from a galaxy catalog. This is given by the second term of Eq.
8 that we factorize as follows

$p(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c},\boldsymbol{\lambda}_{\rm
z})=\frac{p_{\rm gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})\,p_{\rm
rate}(z\,|\,\boldsymbol{\lambda}_{\rm z})}{\int p_{\rm
gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})\,p_{\rm
rate}(z\,|\,\boldsymbol{\lambda}_{\rm
z})\,\mathrm{d}z\,\mathrm{d}\hat{\Omega}}$ (9)

where $p_{\rm gal}$ is the probability that there is a galaxy at
($z,\hat{\Omega}$) and $p_{\rm rate}$ the probability of a galaxy at redshift
$z$ to host a GW event. This takes into account that the probability for a
galaxy to host a merger can have a non-trivial redshift dependence. We assume
that this can be parameterized as

$p_{\rm rate}(z\,|\,\boldsymbol{\lambda}_{\rm
z})\propto\frac{\psi(z;\boldsymbol{\lambda}_{\rm z})}{(1+z)}$ (10)

where $\psi(z;\boldsymbol{\lambda}_{\rm z})$ is the merger rate evolution of
compact objects with redshift and the term $(1+z)^{-1}$ takes into account the
conversion between source and detector time. Note that the normalization
integral in Eq. 9, as well as any rate normalization factor in Eq. 10,
simplifies between numerator and denominator in Eq. 1 and does not need to be
computed explicitly.

###  2.2 Galaxy catalog

As a first approximation, we assume to have a “complete” galaxy catalog, i.e.
this contains all the potential host galaxies. In this case, we denote as
$p_{\rm cat}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$ the probability
distribution built from galaxies in the catalog. We can write $p_{\rm
gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})=p_{\rm
cat}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$ and compute the latter
probability as a sum over the contribution of each galaxy. Given a set of
$\boldsymbol{d}^{\rm EM}=\\{\boldsymbol{d}^{\rm EM}_{g}\\}$ ($g=1,...,N_{\rm
gal}$) EM observations of galaxies we have

$p_{\rm cat}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm
c})=\frac{\sum_{g}w_{g}\,p(z\,|\,\boldsymbol{d}^{\rm
EM}_{g},\boldsymbol{\lambda}_{\rm
c})\,\delta(\hat{\Omega}-\hat{\Omega}_{g})}{\sum_{g}w_{g}}$ (11)

where $w_{g}$ weights the probability of each galaxy to host a GW event (e.g.,
by the galaxy luminosity), $p(z\,|\,\boldsymbol{d}^{\rm
EM}_{g},\boldsymbol{\lambda}_{\rm c})$ is the galaxy’s redshift distribution
that we want to use as a prior, and $\delta(\hat{\Omega}-\hat{\Omega}_{g})$ is
a Dirac delta distribution of each galaxy’s sky localization that can be
treated as errorless.

The galaxy catalog contains redshift measurements $\tilde{z}_{g}$ and
associated uncertainties $\tilde{\sigma}_{z,\,g}$ for each observed galaxy,
i.e. $\boldsymbol{d}^{\rm EM}_{g}=\\{\tilde{z}_{g},\tilde{\sigma}_{z,\,g}\\}$
. From these quantities, we construct the likelihoods, which we assume to be
Gaussian,
$p(\tilde{z}_{g}\,|\,z)=\mathcal{N}(z;\tilde{z}_{g},\tilde{\sigma}_{z,\,g}^{2})$.
This is a probability distribution over the observed values $\tilde{z}_{g}$.
To get $p(z\,|\,\boldsymbol{d}^{\rm EM}_{g},\boldsymbol{\lambda}_{\rm c})$ we
need to multiply it by a prior on the redshift distribution, which in the
absence of other information is naturally chosen as uniform in comoving volume
(Gair et al., 2023). Using Bayes’ theorem, we get

$p(z\,|\,\boldsymbol{d}^{\rm EM}_{g},\boldsymbol{\lambda}_{\rm
c})=\frac{\mathcal{N}(z;\tilde{z}_{g},\tilde{\sigma}_{z,\,g}^{2})\,\frac{\mathrm{d}V_{c}}{\mathrm{d}z}(z;\boldsymbol{\lambda}_{\rm
c})}{\int\mathcal{N}(z;\tilde{z}_{g},\tilde{\sigma}_{z,\,g}^{2})\,\frac{\mathrm{d}V_{c}}{\mathrm{d}z}(z;\boldsymbol{\lambda}_{\rm
c})\,\mathrm{d}z}\,,$ (12)

where $\mathrm{d}{V_{c}}/\mathrm{d}{z}$ is the differential comoving volume
element in a flat universe. With this definition, Eq. 11 is normalized so that
if $p(\tilde{z}_{g}\,|\,z)=\mathcal{\delta}(z-\tilde{z}_{g})$ and in case of
uniform weights, we get $n_{\rm
cat}(z)=(1/V_{c})\sum_{g}\delta(z-\tilde{z}_{g})$ which is consistent with eq.
3.39 of Finke et al. (2021), i.e. the comoving density of galaxies $n_{\rm
cat}(z)$ is estimated by counting the objects in the catalog and dividing by
the total volume. If, instead, the likelihood is completely uninformative,
$p(\tilde{z}_{g}\,|\,z)=1$, we correctly get a comoving density constant in
redshift. Note that the assumption of a uniform-in-comoving-volume prior
introduces a dependence on the cosmology in the redshift prior, since
$\mathrm{d}{V_{c}}/\mathrm{d}{z}$ is a function of $\boldsymbol{\lambda}_{\rm
c}$. Once normalizing correctly $p(z\,|\,\boldsymbol{d}^{\rm
EM}_{g},\boldsymbol{\lambda}_{\rm c})$, however, the dependence on $H_{0}$
cancels, as this is an overall multiplicative factor in
$\mathrm{d}{V_{c}}/\mathrm{d}{z}$, but a dependence on $\Omega_{\rm m,0}$
remains and should be accounted for. The same would be true if one wished to
include more general cosmologies, e.g. a Dark Energy equation of state $w_{\rm
DE}\neq-1$. Alternatively, one may adopt a flat prior, which eliminates the
extra dependence on the cosmology, as done e.g. in Gray et al. (2023).

In the real world, observed galaxy catalogs are subject to selection effects
that result in having an incomplete set of potential hosts. To address this
issue, the $p_{\rm gal}$ term in Eq. 9 has to be modified including the
probability of missing galaxies, $p_{\rm
miss}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$, as follows (Chen et
al., 2018; Finke et al., 2021):

$p_{\rm gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm
c})=f_{\mathcal{R}}\,p_{\rm cat}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm
c})+(1-f_{\mathcal{R}})\,p_{\rm
miss}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})\,.$ (13)

The probability $p_{\rm miss}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm
c})$ is constructed with two pieces of information. The first is the amount of
missing galaxies as a function of redshift and sky position, encoded in the
completeness fraction of the catalog $P_{\rm comp}(z,\hat{\Omega})$. This can
be estimated with some assumption on the total comoving density of galaxies or
on the luminosity distribution of a complete sample that is assumed to follow
a Schechter function. We refer to Finke et al. (2021) for a rigorous
definition of this quantity and for the discussion of different ways to
compute it. In CHIMERA, we will use in particular the so-called “mask”
completeness of Finke et al. (2021), which we will describe better in Sec.
2.4. From this, one can compute the weight $f_{\mathcal{R}}$ in Eq. 13 as

$f_{\mathcal{R}}\equiv\frac{1}{V_{c}(\boldsymbol{\lambda}_{\rm c})}\int P_{\rm
comp}(z,\hat{\Omega})~{}\mathrm{d}V_{c}\,,$ (14)

where the integral extends to a region sufficiently large to encompass all the
GW events under consideration and $V_{c}(\boldsymbol{\lambda}_{\rm c})$
denotes the corresponding comoving volume.

The second bit of information amounts to specifying how the missing galaxies
are distributed. Two extreme possibilities are a “homogeneous completion”,
where the distribution is assumed to be uniform in comoving volume and in sky
position, and a “multiplicative completion”, where missing galaxies are
assumed to trace the distribution of those present in the catalog (Finke et
al., 2021). The two correspond to

$\displaystyle p_{\rm miss}^{\rm
HOM}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$
$\displaystyle=\frac{1-P_{\rm
comp}(z,\hat{\Omega})}{(1-f_{\mathcal{R}})\,V_{c}(\boldsymbol{\lambda}_{\rm
c})}\,\frac{\mathrm{d}V_{c}}{\mathrm{d}z}(z;\boldsymbol{\lambda}_{\rm c})\,,$
(15) $\displaystyle p_{\rm miss}^{\rm
MULT}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$
$\displaystyle=\frac{f_{\mathcal{R}}}{1-f_{\mathcal{R}}}\,\frac{1-P_{\rm
comp}(z,\hat{\Omega})}{P_{\rm comp}(z,\hat{\Omega})}\,p_{\rm
cat}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})\,.$ (16)

In practice, as suggested in Finke et al. (2021) and done in the code
DarkSirensStat, it is convenient to interpolate between a multiplicative
completion where the catalog is fairly complete, and a homogeneous completion
where the catalog is largely incomplete. Finally, an even better option,
suggested in Finke et al. (2021) and recently studied in detail in Dalang &
Baker (2023), would be to complete the catalog using the typical correlation
length of galaxies.

###  2.3 Full form of the likelihood

By putting together Eqs. 1–14 we can write the likelihood in its final form.
In the following, we explicitly restrict to the sub-space
$\bar{\boldsymbol{\theta}}=\\{z,\hat{\Omega},m_{1},m_{2}\\}$ and
$\bar{\boldsymbol{\theta}}^{\mathrm{det}}=\\{d_{L},\hat{\Omega},m_{1}^{\mathrm{det}},m_{2}^{\mathrm{det}}\\}$.
This relies on the assumption that the remaining parameters of the GW waveform
have a population distribution that coincides with the prior used in the
analysis of individual events. We have:

$\displaystyle p(\boldsymbol{d}^{\rm GW}\,|\,\boldsymbol{\lambda})$
$\displaystyle\propto\frac{1}{\xi(\boldsymbol{\lambda})^{N_{\rm
ev}}}\prod_{i=1}^{N_{\rm
ev}}\int\mathrm{d}z\,\mathrm{d}\hat{\Omega}\,\mathcal{K}_{\mathrm{gw},i}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\mathrm{c}},\boldsymbol{\lambda}_{\mathrm{m}})\,p_{\rm
gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\mathrm{c}})\,\frac{\psi(z;\boldsymbol{\lambda}_{\mathrm{z}})}{1+z}\,,$
(17)
$\displaystyle\mathcal{K}_{\mathrm{gw},i}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\mathrm{c}},\boldsymbol{\lambda}_{\mathrm{m}})$
$\displaystyle\equiv\int\mathrm{d}m_{1}\mathrm{d}m_{2}\,\frac{p(z,m_{1},m_{2},\hat{\Omega}\,|\,\boldsymbol{d}^{\rm
GW}_{i},\boldsymbol{\lambda}_{\mathrm{c}})}{\pi(d_{L})\pi(m_{1}^{\mathrm{det}})\pi(m_{2}^{\mathrm{det}})}\,\frac{1}{\frac{\mathrm{d}d_{L}}{\mathrm{d}z}(z;\boldsymbol{\lambda}_{\mathrm{c}})\,(1+z)^{2}}\,p(m_{1},m_{2}\,|\,\boldsymbol{\lambda}_{\mathrm{m}})\,,$
(18) $\displaystyle\xi(\boldsymbol{\lambda})$
$\displaystyle=\int\mathrm{d}\boldsymbol{\theta}\,P_{\rm
det}(\boldsymbol{\theta},\boldsymbol{\lambda}_{\mathrm{c}})\,\,p(m_{1},m_{2}\,|\,\boldsymbol{\lambda}_{\mathrm{m}})\,p_{\rm
gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\mathrm{c}})\,\frac{\psi(z;\boldsymbol{\lambda}_{\mathrm{z}})}{1+z}\,.$
(19)

Figure 1: Visual representation of the underlying workflow of CHIMERA. The GW
probability is approximated by a three-dimensional KDE
$\mathcal{K}_{\mathrm{gw}}$, while the galaxy probability $p_{\mathrm{gal}}$
is evaluated by summing the individual galaxies contributions within the pixel
enclosed in the GW sky localization area.

Note that we have omitted the overall normalization of the redshift prior, as
this is independent of $z,\hat{\Omega}$ and is present also in the selection
function, leading to a simplification between numerator and denominator in Eq.
17. A different normalization of this quantity just amounts to a different
interpretation of the selection function, which does not correspond to the
fraction of detectable events, as a consequence of the fact that
$p_{\mathrm{pop}}$ looses its correct normalization to unity. We find the form
of Eqs. 17-19 useful since it makes clear what quantities actually need to be
computed in the numerical evaluation. On the other hand, we note that the
normalization of the term $p_{\rm gal}$ has to be computed carefully, in order
to correctly balance the contributions of the catalog and missing terms in Eq.
13.

Figure 2: Analysis of GW170817 with CHIMERA using GLADE+ galaxies with a
luminosity $L_{K}>0.02\,L_{K}^{\ast}$. Left: Sky distribution of the potential
host galaxies (red crosses) contained in the 90% credible pixels (gray
squares). The true host NGC4993 is identified with a gold circle. Center:
Redshift distribution of the GW kernel (blue) assuming
$H_{0}=70~{}\mathrm{km/s/Mpc}$ and galaxy probability $p_{\mathrm{gal}}$ (red)
in each pixel. The pixel including NGC 4993 is colored in black. Right:
Posterior distributions on $H_{0}$ assuming that the true host is or is not
observed.

###  2.4 Implementation in CHIMERA

We now discuss the implementation of the method in CHIMERA. The code
accurately computes Eq. 17 across different regimes, i.e. when the EM redshift
information in highly constraining (bright sirens), more loosely constraining
(dark sirens), or not present at all (spectral sirens). We also prioritize
computational efficiency to meet the demands of next-generation GW
observatories and upcoming galaxy surveys.

In this section, we resume the key aspects of the code, while a more detailed
description is provided in Appendix A. We will use the real event GW170817 as
an illustrative example, considering it a dark siren in conjunction with the
GLADE+ galaxy catalog. This is useful because GW170817 is the best localized
GW event so far and its host galaxy, NGC 4993, has been uniquely identified
(Abbott et al., 2017a), providing an ideal benchmark for dark siren
techniques. Indeed, this is often used as a test case and allows comparison to
other pipelines (Mastrogiovanni et al., 2023; Gray et al., 2023). In Appendix
B, we provide details about the data and cuts used for this example.

Figure 1 shows a compact representation of the workflow of CHIMERA, in the
case of GW170817. In short, the galaxy catalog term is pre-computed as a
function of redshift (red lines) on a pixel-by-pixel basis inside the GW
localization region (given by the gray squares). The GW kernel in Eq. 18 is
given by a kernel density estimation (KDE) of the GW posterior samples re-
weighted for the population prior on source-frame masses (blue area).

##### Pixelization and pre-computation of the galaxy term

The quantity $p_{\rm gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$
can be cosmology-dependent, but only when adopting a cosmology-dependent prior
on the galaxies’ redshifts, and in any case never on $H_{0}$, as explained in
Sec. 2. If adopting a flat prior on the galaxies’ redshifts, or a narrow prior
on the other cosmological parameters (e.g. a Planck prior on the matter
density), the quantity $p_{\rm gal}$ can be computed once for all, given a
specific galaxy catalog. For each event in the GW catalog, CHIMERA selects
galaxies enclosed in the GW localization volume, with a redshift range
computed taking into account all the prior range for $H_{0}$. Then, the sky
area of the GW event is divided into equal-area pixels in the pixelization
scheme of Healpix (Górski et al., 2005; Zonca et al., 2019). CHIMERA includes
an adaptive pixelization procedure that allows to set the pixel size for each
GW event depending on its localization area (see App. A for details). This
step corresponds to the left panel of Fig. 2, where we show the galaxies
enclosed in the $90\%$ localization volume of GW170817 as red crosses,
overlapped with the pixels. The quantity $p_{\rm
gal}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$ is then approximated as
a redshift-dependent function of the form given in Eqs. 11-12 for each pixel,
by summing the contribution of all the galaxies enclosed within it. The
weights in Eq. 11 can be chosen to be either uniform or corresponding to the
galaxies’ luminosity in the catalog. The redshift distributions in each pixel
are shown in red in the central panel of Fig. 2 for the example of GW170817.
This method allows to compress the information coming from all galaxies in
each pixel in a single function of redshift. In the case of GW170817, we find
88 galaxies in the 90% localization volume, while using 20 pixels, which
corresponds to a factor of $\sim 5$ improvement in the computational cost.
This procedure is even more effective in the case of events with thousands of
galaxies in the localization region. We note that, in general, the GW skymap
changes on much larger angular scales than the galaxy distribution. This means
that the compression introduced here (which basically amounts at approximating
the galaxies’ angular positions with the center of each pixel) does not lead
to a loss of information as far as the pixel size is chosen appropriately.

##### Catalog completeness

The completeness fraction $P_{comp}(z,\hat{\Omega})$ is computed following the
mask method presented in Finke et al. (2021). First, similar pixels are
grouped together in masks by applying an agglomerative clustering algorithm.
The chosen feature for this purpose is the number of galaxies within each
pixel. The number of masks is arbitrarily chosen depending on the survey
properties. For each mask, the completeness fraction is computed by comparing
the luminosity function at each redshift to a reference Schechter function. We
refer to Finke et al. (2021) for a more detailed description of this
procedure. In particular, this allows to easily cut regions particularly
under-covered by the galaxy survey, for example the region of the sky covered
by the Milky Way. It also makes straightforward the use of galaxy catalogs
with partial sky coverage. In the example of GW170817, it is complete.

##### 3D GW kernel

We now turn to the term $\mathcal{K}_{\mathrm{gw}}$ in Eq. 18. For each GW
event, we have a set of posterior samples of the detector-frame quantities
$\bar{\boldsymbol{\theta}}^{\mathrm{det}}=\\{d_{L},\hat{\Omega},m_{1}^{\mathrm{det}},m_{2}^{\mathrm{det}}\\}$
approximating the posterior probability
$p(m_{1}^{\mathrm{det}},m_{2}^{\mathrm{det}},d_{L},\hat{\Omega}\,|\,\boldsymbol{d}^{\rm
GW}_{i})$. CHIMERA assumes that the prior on detector-frame masses used for
the individual event analysis is flat,
$\pi(m_{1}^{\mathrm{det}})=\pi(m_{2}^{\mathrm{det}})=1$, while the prior on
the luminosity distance is $\pi(d_{L})\propto d_{L}^{2}$. Then, for each value
of $\boldsymbol{\lambda}$, the samples are converted to source-frame
quantities by inverting Eqs. 3-5. Each sample, labeled by $j$, is assigned a
weight equal to

$w_{\mathcal{K},j}(\boldsymbol{\lambda})=\frac{p(m_{1,j},m_{2,j}\,|\,\boldsymbol{\lambda}_{\rm
m})}{d_{L}^{2}\,\frac{\mathrm{d}d_{L}}{\mathrm{d}z}(z_{j};\boldsymbol{\lambda}_{\mathrm{c}})\,\left(1+z_{j}\right)^{2}}\,,$
(20)

and $\mathcal{K}_{\mathrm{gw}}$ is obtained by a 3D weighted KDE in the space
$(z,\mathrm{RA},\mathrm{Dec})$ with weights
$w_{\mathcal{K},j}(\boldsymbol{\lambda})$. Using the re-weighted samples in
this subspace corresponds to performing the integral in Eq. 18 with the
correct population prior $p(m_{1},m_{2}\,|\,\boldsymbol{\lambda}_{\rm m})$ on
the _source-frame_ masses. This will encode a dependence on the cosmological
parameters, which enters in the conversion of the posterior samples from
detector-frame to source-frame, resulting in a re-shaping of the kernel
$\mathcal{K}_{\mathrm{gw}}$. A non-flat source-frame mass spectrum will add
redshift information at the statistical level trough this term, even in
absence of a galaxy catalog, in which case we recover the spectral siren case.
For GW170817, the GW kernel is plotted in blue in Fig. 2b, assuming a value of
$H_{0}=70~{}\mathrm{km/s/Mpc}$ and a flat mass distribution between $1$ and
$3$ solar masses.

##### Integration

The double integral in Eq. 17 is performed by integrating numerically
$\mathrm{d}z$ in each pixel first. Then, the integral in
$\mathrm{d}\hat{\Omega}$ can be approximated as a discrete sum of the results
in each pixel, multiplied by the pixel area. The latter, for each given GW
event, is equal for all pixels and inversely proportional to the number of
Healpix pixels of the map. So, in the end, the angular integral in Eq. 17 is
obtained as an average over pixels.

##### Selection effects

The selection bias term $\xi(\boldsymbol{\lambda})$ from Eq. 2 is computed by
using the standard re-weighted Monte Carlo method (Tiwari, 2018; Farr, 2019).
In this case, it is crucial to use _detector-frame_ quantities, since the
detectability in a GW experiment is a function of those only, i.e. $P_{\rm
det}=P_{\rm det}(\boldsymbol{\theta}^{\rm det})$. This allows to compute this
function only once. The integral in Eq. 19 can then be computed with the
change of variables $\boldsymbol{\theta}\mapsto\boldsymbol{\theta}^{\rm
det}(\boldsymbol{\theta},\boldsymbol{\lambda}_{\mathrm{c}})$, where the
Jacobian factors are given in Eq. 7. CHIMERA takes as input a set of simulated
events, each with detector-frame parameters $\boldsymbol{\theta}_{i}^{\rm
det}$, drawn from a reference distribution with probability
$p_{\mathrm{draw}}(\boldsymbol{\theta}_{i}^{\rm det})$. These must have been
previously injected in the same detection pipeline used to obtain the GW
catalog, storing those that pass the detection threshold.

Then, the selection function in Eq. 19 is computed with Monte Carlo
integration as

$\begin{split}\xi(\boldsymbol{\lambda})=&\frac{1}{N_{\mathrm{inj}}}\sum_{i=1}^{N_{\mathrm{det}}}\frac{1}{p_{\mathrm{draw}}(\boldsymbol{\theta}_{i}^{\rm
det})}\frac{1}{\frac{\mathrm{d}d_{L}}{\mathrm{d}z}(z_{i},\boldsymbol{\lambda}_{\mathrm{c}})(1+z_{i})^{2}}\\\
&\times p(m_{1,i},m_{2,i}\,|\,\boldsymbol{\lambda}_{\mathrm{m}})\,p_{\rm
gal}(z_{i},\hat{\Omega}_{i}\,|\,\boldsymbol{\lambda}_{\rm
c})\,\frac{\psi(z_{i};\boldsymbol{\lambda}_{\rm z})}{1+z_{i}}\end{split}$ (21)

where $N_{\mathrm{inj}}$ is the total number of injections, $N_{\mathrm{det}}$
the number of detected ones, and the source-frame quantities
$z_{i},m_{1,i},m_{2,i}$ are understood to be computed from the (fixed)
detector-frame ones via inversion of Eqs. 3-5, and are thus functions of the
cosmological parameters $\boldsymbol{\lambda}_{\mathrm{c}}$. Finally, when
computing the selection bias term, CHIMERA employs a galaxy catalog
interpolant $p_{\rm gal}$ constructed on the full sky.

##### Accuracy of the MC integrals

Finally, when computing Eq. 21, we ensure that the effective number of
independent draws after re-weighting, $N_{\mathrm{eff}}$, for a given
population is sufficiently high ($N_{\mathrm{eff}}>5N_{\mathrm{det}}$, see
Farr 2019). This criterion is also applied to the GW kernel weights to ensure
numerical stability.

In conclusion, the right panel of Fig. 2 shows the results for the posterior
probability on $H_{0}$ (solid line). We obtain a value of
$H_{0}=73^{+58}_{-22}~{}\mathrm{km/s/Mpc}$ in the dark siren case, in good
agreement with Fishbach et al. (2019). As can be seen from Fig. 2 (central
panel), the primary contribution to the posterior arises from galaxies at
approximately $z\sim 0.01$. This value combined with the GW measurement of
$d_{L}\sim 40$ Mpc, implies $H_{0}\sim 70~{}\mathrm{km/s/Mpc}$. The galaxy
“groups” at $z\sim 0.02~{}(0.04)$ provide a shallower contribution to
$H_{0}\sim 150~{}(300)~{}\mathrm{km/s/Mpc}$ due to the presence of selection
effects which disfavor high values of $H_{0}$ [$\xi(H_{0})\sim H_{0}^{3}$]. In
this case, the assumption of a population model for GW170817 provides no
significant improvements to $H_{0}$.

For comparison, and to test the performance of the code in a the bright siren
regime, we repeat the analysis assuming the identification of the host galaxy
with $z=0.0100\pm 0.0005$ (Abbott et al., 2017a). In this case, we assign a
weight $w=1$ to NGC 4993 and $w=0$ to all the other galaxies (see Eq. 11). The
resulting posterior on $H_{0}$ is shown in Fig. 2c (dashed line). In the
bright siren case, we obtain a value of
$H_{0}=69^{+15}_{-8}~{}\mathrm{km/s/Mpc}$, in very good agreement with Abbott
et al. (2023a).

Figure 3: Main properties of the simulated O4-like (blue) and O5-like (red) GW
catalogs. Upper panels: GW sky localization areas at $1$ and $2\,\sigma$
overlaid to their potential host galaxies (gray) extracted from MICEv2. Middle
panels: distribution of the relative uncertainty on the luminosity distance
and sky localization area as a function of network SNR. Lower panels:
distribution of the GW events as a function of (left to right): redshift,
primary mass, and number of galaxies contained in the GW localization volume.

##  3 Mock Catalogs

In this section, we describe the procedure to build realistic mock galaxy and
GW catalogs. These will be used to both validate the code and provide updated
forecasts for O4-like and O5-like detector configurations.

###  3.1 Parent Galaxy Sample

We generate our mock galaxy catalog (hereafter, parent sample) from the MICE
Grand Challenge light-cone simulation (v2), which populates one octant of the
sky (close to 5157 deg2) designed to mimic a complete DES-like (Dark Energy
Survey) survey up to an observed magnitude $i<24$ at redshift $z<1.4$ (Fosalba
et al., 2015; Crocce et al., 2015; Fosalba et al., 2015; Carretero et al.,
2015; Hoffmann et al., 2015). MICE assumes a flat $\Lambda$CDM cosmology with
$H_{0}=70~{}\mathrm{km/s/Mpc}$, $\Omega_{m,0}=0.25$,
$\Omega_{\Lambda,0}=0.75$.

While we ideally require a complete catalog with high number density, as a
simplifying assumption in this paper we consider only galaxies with stellar
masses $\log M_{\star}/\mathrm{M}_{\odot}>10.5$. This cut is consistent with
the idea that the binary merger rate is traced by stellar mass, as also
adopted in current standard sirens analysis via absolute magnitude cuts and
luminosity weighting (Fishbach et al., 2019; Finke et al., 2021; Gray et al.,
2022; Abbott et al., 2023a; Mastrogiovanni et al., 2023; Gray et al., 2023). A
similar cut in mass is also considered in (Muttoni et al., 2023) in the
context of simulations for the Einstein Telescope.

We subsample the MICEv2 catalog to reproduce the density for the cut described
above extracting the galaxies to get a uniform in comoving volume
distribution. In the end, we obtain a parent sample of about 1.6 million
massive galaxies.

For the redshift uncertainties, we consider two cases. First of all, we
explore the possibility of maximizing the galaxy catalog information by having
a spectroscopic catalog. This could be done by expanding the currently
available catalogs (GLADE+ Dálya et al., 2022) in the future by exploiting the
information provided by the next large spectroscopic surveys. As an example,
the ESA mission Euclid (Laureijs et al., 2011) will provide an all-sky map of
spectroscopic redshift in the range $0.9<z<1.8$, with an accuracy of
$\sigma_{z}/(1+z)\lesssim 0.001$, and the Dark Energy Spectroscopic Instrument
(DESI DESI Collaboration et al., 2016) is planned to observe $\sim 14000$ deg2
covering the redshift range $0.4<z<2.1$. Second, we study how the information
extracted changes when using photometric redshift, assuming an uncertainty
$\sigma_{z}/(1+z)=0.05$. This is easily accessible with current ongoing
surveys like the Dark Energy Survey (e.g., DES has reached $\sigma_{z}\sim
0.01$, Myles et al. 2021, and this limit can be pushed to $\sigma_{z}\sim
0.007$ with improved techniques, e.g. Buchs et al. 2019) on a smaller area,
and in future surveys like Euclid and Rubin Observatory (Ivezić et al., 2019)
are planned to extend it to the entire sky to a depth of Euclid H magnitude
$H_{E}\sim 24$ ($H_{E}\sim 26$ in the Deep Survey) and an expected uncertainty
$\sigma_{z}/(1+z)\lesssim 0.05$ (Euclid Collaboration: Desprez et al., 2020;
Euclid Collaboration: Schirmer et al., 2022).

Therefore, in this paper we consider two regimes of photometric (hereafter,
“phot”) and spectroscopic (hereafter, “spec”) redshift uncertainties:

$\sigma_{z}=\begin{cases}0.001\,(1+z)&(z_{\rm spec})\\\ 0.05\,(1+z)&(z_{\rm
phot})\end{cases}$ (22)

Table 1: Summary of the population hyperparameters and priors adopted. $\mathcal{U(\cdot)}$ denotes a uniform distribution. Parameter | Description | Fiducial Value | Prior
---|---|---|---
| Cosmology (flat $\Lambda$CDM) |  |
$H_{0}$ | Hubble constant [km/s/Mpc] | 70.0 | $\mathcal{U}({10.0},{200.0})$
$\Omega_{\rm m,0}$ | Matter energy density | 0.25 | Fixed
| Rate evolution (Madau-like) |  |
$\gamma$ | slope at $z<z_{p}$ | 2.7 | $\mathcal{U}({0.0},{12.0})$
$\kappa$ | slope at $z>z_{p}$ | 3 | $\mathcal{U}({0.0},{6.0})$
$z_{\rm p}$ | peak redshift | 2 | $\mathcal{U}({0.0},{4.0})$
| Mass distribution (PowerLaw+Peak) |  |
$\alpha$ | (primary) slope of the power law | 3.4 | $\mathcal{U}({1.5},{12.0})$
$\beta$ | (secondary) slope of the power law | 1.1 | $\mathcal{U}({-4.0},{12.0})$
$\delta_{m}$ | (primary) smoothing parameter [$\mathrm{M}_{\odot}$] | 4.8 | $\mathcal{U}({0.01},{10.0})$
$m_{\rm low}$ | lower value $[\mathrm{M}_{\odot}]$ | 5.1 | $\mathcal{U}({2.0},{50.0})$
$m_{\rm high}$ | upper value $[\mathrm{M}_{\odot}]$ | 87.0 | $\mathcal{U}({50.0},{200.0})$
$\mu_{\rm g}$ | (primary): mean of the Gaussian component $[\mathrm{M}_{\odot}]$ | 34.0 | $\mathcal{U}({2.0},{50.0})$
$\sigma_{\rm g}$ | (primary): standard deviation of the Gaussian component $[\mathrm{M}_{\odot}]$ | 3.6 | $\mathcal{U}({0.4},{10.0})$
$\lambda_{\rm g}$ | (primary): fraction of the Gaussian component | 0.039 | $\mathcal{U}({0.01},{0.99})$

###  3.2 Sample of GW events

We generate mock GW events from the parent sample by fixing cosmological
hyperparameters $\boldsymbol{\lambda}_{\rm c}$ and astrophysical population
hyperparameters $\boldsymbol{\lambda}_{\rm z}$ and $\boldsymbol{\lambda}_{\rm
m}$. We describe the redshift and mass distributions in turn.

####  3.2.1 Rate evolution

The source-frame merger rate is parameterized as follows (see Madau &
Dickinson, 2014):

$\psi(z;\boldsymbol{\lambda}_{\rm
z})=\frac{(1+z)^{\gamma}}{1+\left(\frac{1+z}{1+z_{\mathrm{p}}}\right)^{\gamma+\kappa}},$
(23)

where $\psi(z)\propto(1+z)^{\gamma}$ at low $z$, then reaches its peak near
$z_{\mathrm{p}}$, and subsequently declines as
$\psi(z)\propto(1+z)^{-\kappa}$. To build the catalog, we use $\gamma=2.7$
consistent with the LVK GWTC-3 results (Abbott et al., 2023b). The limited
detection range of current gravitational wave detectors restricts our ability
to determine the merger rate at higher redshifts, therefore we adopt
$z_{\mathrm{p}}=2$, $\kappa=3$, consistent with the idea that
$\psi(z;\boldsymbol{\lambda}_{\rm z})$ follows the galaxy’s star formation
rate density with parameters from Madau & Dickinson (2014). The catalog of
potential sources is then obtained by sampling the parent sample using a
weight proportional to the detector-frame merger rate
$\psi(z;\boldsymbol{\lambda}_{\rm z})/(1+z)$.

####  3.2.2 Mass distribution

For the mass distribution, we adopt the phenomenological “PowerLaw+Peak” (PLP)
model following the LVK GWTC-3 results (Abbott et al., 2023b). The mass term
of Eq. 8 is factorized as follows

$p(m_{1},m_{2}\,|\,\boldsymbol{\lambda}_{\rm
m})=p(m_{1}\,|\,\boldsymbol{\lambda}_{\rm
m})\,p(m_{2}\,|\,m_{1},\boldsymbol{\lambda}_{\rm m})\,.$ (24)

The probability of the primary BH mass is given by

$p(m_{1}\,|\,\boldsymbol{\lambda}_{\rm m})\propto\left[(1-\lambda_{\rm
g})\,\mathcal{P}(m_{1})+\lambda_{\rm
g}\,\mathcal{G}(m_{1})\right]\,\mathcal{S}(m_{1})\,,$ (25)

where $\mathcal{P}(m)\propto m^{-\alpha}$ is a power-law truncated in the
domain $m\in[m_{\rm low},m_{\rm high}]$,
$\mathcal{G}(m)\propto\mathcal{N}(\mu_{\mathrm{g}};\sigma_{\mathrm{g}}^{2})$
is a Gaussian component whose contribution relative to $\mathcal{P}$ is
regulated by $\lambda_{\rm g}$, and $\mathcal{S}(m)\in[0,1]$ is a smoothing
piece-wise function with a tapering parameter $\delta_{m}$ fully described in
Abbott et al. (2021d) appendix B. The secondary BH mass is modeled by a power-
law with an index $\beta$ in the domain $m\in[m_{\rm low},m_{1}]$.

####  3.2.3 Summary

With the above assumptions, the cosmological and astrophysical hyperparameters
to be studied are:

$\begin{split}\boldsymbol{\lambda}_{\rm c}&=\\{H_{0},\Omega_{\rm m,0}\\}\\\
\boldsymbol{\lambda}_{\rm z}&=\\{\gamma,k,z_{\rm p}\\}\\\
\boldsymbol{\lambda}_{\rm m}&=\\{\alpha,\beta,\delta_{m},m_{\rm low},m_{\rm
high},\mu_{\rm g},\sigma_{\rm g},\lambda_{\rm g}\\}.\end{split}$ (26)

The fiducial values and the prior ranges chosen for this work are reported in
Table 1. For the cosmological parameters, we will consider the value of the
matter density to be fixed to its fiducial value.

###  3.3 GW data generation

For the simulation of GW detections we use the simulation pipeline
GWFAST777Available at https://gwfast.readthedocs.io/en/latest/ (Iacovelli et
al., 2022a, b).

We assume quasi-circular non-precessing BBH systems. Their waveform is
characterized by the detector-frame parameters

$\boldsymbol{\theta}^{\mathrm{det}}=\left\\{\mathcal{M}_{c},\eta,\chi_{1,z},\chi_{2,z},d_{L},\theta,\phi,\iota,\psi,t_{c},\Phi_{c}\right\\}\,,$
(27)

where $\mathcal{M}_{c}$ is the detector-frame chirp mass, $\eta$ is the
symmetric mass ratio, $\chi_{1/2,z}$ are the adimensional spin parameters
along the direction of the orbital angular momentum, $d_{L}$ is the luminosity
distance, $\theta=\pi/2-\mathrm{Dec}$ and $\phi=\mathrm{RA}$ are the sky
position angles, $\iota$ refers to the inclination angle of the binary’s
orbital angular momentum with respect to the line of sight, $\psi$ is the
polarization angle, $t_{c}$ is the coalescence time, and $\Phi_{c}$ is the
phase at coalescence.

First of all, we generate a population of GW events following the
prescriptions given in Section 3.2, so that each source is characterized by a
set of source-frame parameters $\boldsymbol{\theta}$. For the parameters in
Eq. 27 that are not explicitly mentioned in Sec. 3.2, we assume the following
distributions: the sky position angles are uniform on the octant covered by
the parent sample, the inclination angles have a uniform distribution in
$\cos\iota$ in the range $[0,\,\pi]$, the polarization angle and coalescence
phase have a uniform distribution in $[0,\,\pi]$ and $[0,\,2\pi]$,
respectively, and the time of coalescence, given as Greenwich mean sidereal
time and expressed in units of fraction of a day, is uniform in the interval
$[0,\,1]$. For the spin parameters instead, we adopt a flat distribution
between $[-1,\,1]$ for the components aligned with the orbital angular
momentum. The source-frame parameters are then converted to detector frame
using the fiducial cosmology (see Table 1).

For each source, we simulate GW emission and compute the network matched-
filtered signal-to-noise ratio ($\mathrm{SNR}$) using the IMRPhenomHM (London
et al., 2018; Kalaghatgi et al., 2020) waveform approximant, which includes
the contribution of sub-dominant modes to the signal, that are of great
relevance to describe in particular the merger phase of BBH systems. We
consider two network configurations. The first, denoted as O4-like, is
composed of a network of the two LIGO interferometers at Hanford and
Livingston, USA, the Virgo interferometer in Cascina, Italy, and the KAGRA
interferometer in Japan. We assume sensitivity curves representative of the O4
run of the LVK collaboration, with public sensitivity curves from Abbott et
al. (2016). For the second configuration, denoted as O5-like the network
includes the two LIGO, Virgo, and KAGRA instruments, as well as a LIGO
detector located in India. For O4, we use sensitivity curves representative of
the O5 run of the LVK collaboration888The amplitude spectral densities can be
found at https://dcc.ligo.org/LIGO-T2000012/public. For O4, we use kagra_3Mpc
for KAGRA. For O5, we use AplusDesign for the three LIGO detectors,
avirgo_O5low_NEW for Virgo, and kagra_80Mpc for KAGRA. (Abbott et al., 2016).
We assume a $100\%$ duty cycle in all cases.

Then, we select samples of GW events as follows,

  * •

O4-like: 100 events with a network $\mathrm{SNR}>12$;

  * •

O5-like: 100 events with a network $\mathrm{SNR}>25$.

These cuts are designed to yield the 100 best events for each configuration
over approximately one year of observation. We determined these numbers
simulating with GWFAST a one-year observing run for each of the two scenarios,
with a population of BBHs with overall merger rate calibrated on the latest
LVK constraints (see e.g. Iacovelli et al., 2022a). We note that the fact that
the simulation is performed on one octant of the sky is irrelevant in our
case, as we do not constrain direction-dependent hyperparameters.

We approximate the GW posterior probability with the Fisher Information Matrix
(FIM) approximation, where the GW likelihood $p(\boldsymbol{d}^{\rm
GW}_{i}\,|\,\boldsymbol{\theta}_{i}^{\rm det})$ is assumed to be given by a
multivariate Gaussian distribution. This approximation is valid for the high
SNR events that we are considering in this work (see Iacovelli et al. (2022a)
for details). We compute the FIM with GWFAST. For each event, we then draw
$5000$ samples from a multivariate Gaussian distribution with covariance equal
to the inverse of the FIM, further imposing priors that are chosen to be flat
in detector-frame masses with the condition $m_{2}^{\rm det}<m_{1}^{\rm det}$,
$\propto d_{L}^{2}$ for the luminosity distance, and equal to the
distributions used to sample the events for the remaining parameters. We then
use the samples of $m_{1}^{\rm det},m_{2}^{\rm
det},d_{L},\mathrm{RA},\mathrm{Dec}$ to compute the GW likelihood as described
in Sec. 2.4.

The main properties of the galaxy and GW catalogs are summarized in Fig. 3.
The top panels show the GW skymaps of the events detected in the two
configurations overlaid to the galaxy sky distribution. The central panels
present the scatter plots of the sky localization area versus the error on
$d_{L}$ (with a color scale giving information on the SNR). Finally, the
bottom panels show the redshift and mass distribution of the detected GW
events, as well as the distribution of the number of galaxies found within
their localization volumes.

The O4-like events are at redshift $z\lesssim 0.9$ with sky areas between a
few and a few hundreds $\rm{deg}^{2}$. This typically results in more than
$\sim 500$ potential host galaxies for the 90% of the events (Fig. 3, bottom-
right panel). A particularly lucky event is present with small localization
area ($\sim 2~{}\mathrm{deg^{2}}$) and high SNR ($\sim 55$, see the mid-left
panel of Fig. 3). While this event represents an outlier of the distribution
that can be ascribed to a statistical fluctuation, it still represents a
possibility with this configuration. For the O5-like events, the redshift
distribution remains limited to $z\lesssim 1$ as a consequence of the higher
SNR cut, while the larger and more sensitive detector network substantially
improves the localization capabilities. This typically results in more than
$\sim 50$ potential host galaxies for the 90% of the events (Fig. 3, bottom-
right panel). In this configuration, there is a significant tail of events
with just few tens of galaxies, corresponding to sky localization regions of
less than one degree with SNR that can exceed 100. Overall, while the number
of galaxies in the localization volume depends on the assumption on the galaxy
catalog employed, it is interesting to observe the $10\times$ reduction in the
number of potential host galaxies between the O4 and O5 networks. This
improvement, following from the much smaller localization volumes, is a key
factor to obtain improved dark sirens constraints as discussed in Section 4.

Ultimately, for the computation of the selection bias, we generate injection
sets for both the O4-like and the O5-like scenarios with GWFAST, adopting the
same selection cuts as for the GW catalogs. The injections cover the same sky
area as the catalogs and a distance range that arrives up to the detector
horizon for the given $\mathrm{SNR}$ cuts. The injections set is made of
$N_{\rm inj}=2\times 10^{7}$ and $4\times 10^{7}$ events, resulting in about
$1.5\times 10^{6}$ and $1\times 10^{6}$ detected events in the O4-like and
O5-like scenario, respectively. These are then used to estimate the selection
bias as in Eq. 21.

Table 2: Median and 68% confidence level interval of the hyperparameters
resulting from spectral only and full analyses for the O4- and O5-like
configurations and for spectroscopic and photometric errors on the galaxy
redshifts. The relative uncertainties are given in parentheses. The rate
parameters $k$ and $z_{p}$ are not included as they remain practically
unconstrained.

O4-like (100 events with SNR $>$ 12) O5-like (100 events with SNR $>$ 25)
Spectral Full ($z_{\rm phot}$) Full ($z_{\rm spec}$) Spectral Full ($z_{\rm
phot}$) Full ($z_{\rm spec}$) $H_{0}$ $64^{+32}_{-23}\,(43\,\%)$
$76^{+16}_{-12}\,(18\,\%)$ $75.3^{+5.2}_{-4.9}\,(7\,\%)$
$55^{+20}_{-16}\,(32\,\%)$ $73^{+7}_{-6}\,(9\,\%)$ $70.2\pm 0.8\,(1\,\%)$
$\gamma$ $1.6^{+1.3}_{-1.0}\,(67\,\%)$ $1.8^{+1.1}_{-0.9}\,(54\,\%)$
$1.8^{+1.4}_{-0.9}\,(62\,\%)$ $2.2^{+1.9}_{-1.1}\,(68\,\%)$
$2.3^{+1.6}_{-0.9}\,(56\,\%)$ $2.3^{+1.3}_{-1.0}\,(50\,\%)$ $\alpha$ $3.6\pm
0.4\,(10\,\%)$ $3.7\pm 0.4\,(10\,\%)$ $3.7\pm 0.3\,(9\,\%)$ $3.2\pm
0.3\,(8\,\%)$ $3.3\pm 0.3\,(8\,\%)$ $3.3\pm 0.3\,(8\,\%)$ $\beta$
$2.5^{+0.8}_{-0.7}\,(29\,\%)$ $2.4^{+0.9}_{-0.7}\,(32\,\%)$ $2.6\pm
0.7\,(29\,\%)$ $1.9^{+0.6}_{-0.5}\,(28\,\%)$ $2.0\pm 0.6\,(29\,\%)$ $1.8\pm
0.6\,(30\,\%)$ $\delta_{m}$ $7.1^{+1.9}_{-3.2}\,(34\,\%)$
$7.0^{+2.1}_{-3.1}\,(37\,\%)$ $7.3^{+1.8}_{-2.5}\,(29\,\%)$
$3.1^{+2.5}_{-1.8}\,(69\,\%)$ $2.8^{+2.3}_{-1.6}\,(70\,\%)$
$3.5^{+2.3}_{-1.7}\,(57\,\%)$ $m_{\rm low}$ $4.4^{+0.8}_{-0.7}\,(16\,\%)$
$4.3^{+0.7}_{-0.8}\,(16\,\%)$ $4.3^{+0.6}_{-0.5}\,(13\,\%)$
$6.0^{+0.7}_{-1.1}\,(15\,\%)$ $5.8^{+0.7}_{-1.1}\,(15\,\%)$
$5.6^{+0.7}_{-1.0}\,(15\,\%)$ $m_{\rm high}$ $102^{+52}_{-21}\,(35\,\%)$
$95^{+47}_{-17}\,(33\,\%)$ $94^{+59}_{-15}\,(38\,\%)$
$102^{+23}_{-13}\,(18\,\%)$ $91.5^{+27.8}_{-7.2}\,(19\,\%)$
$93.2^{+35.3}_{-7.1}\,(23\,\%)$ $\mu_{\rm g}$ $34.2^{+3.5}_{-4.8}\,(12\,\%)$
$32.9^{+2.4}_{-3.5}\,(9\,\%)$ $32.8^{+1.7}_{-2.7}\,(7\,\%)$ $36.5\pm
2.7\,(7\,\%)$ $34.1^{+1.3}_{-1.4}\,(4\,\%)$ $34.40^{+0.99}_{-1.05}\,(3\,\%)$
$\sigma_{\rm g}$ $5.3^{+2.6}_{-1.8}\,(42\,\%)$ $5.0^{+2.0}_{-1.6}\,(35\,\%)$
$4.9^{+2.2}_{-1.4}\,(38\,\%)$ $4.47^{+1.21}_{-0.84}\,(23\,\%)$
$4.09^{+0.90}_{-0.65}\,(19\,\%)$ $4.24^{+0.82}_{-0.71}\,(18\,\%)$
$\lambda_{\rm g}$ $0.03^{+0.03}_{-0.01}\,(74\,\%)$
$0.03^{+0.03}_{-0.01}\,(82\,\%)$ $0.03^{+0.02}_{-0.01}\,(62\,\%)$
$0.07^{+0.04}_{-0.03}\,(50\,\%)$ $0.07^{+0.04}_{-0.03}\,(46\,\%)$
$0.07^{+0.04}_{-0.03}\,(49\,\%)$

Figure 4: Marginalized posterior distributions of a selection of
hyperparameters in the O4-like (upper panels) and O5-like (lower panels)
detector networks. The analysis is performed with CHIMERA in three setups,
namely: spectral-only (Spect.) and spectral with the inclusion of a
photometric ($z_{\rm phot}$) or spectroscopic ($z_{\rm spec}$) galaxy catalog.

##  4 Results

In this Section, we report the results of the analyses of the O4- and O5-like
configurations. For both detector networks, we consider three distinct
analysis setups:

  * •

Full ($z_{\rm spec}$): full analysis using the parent catalog of 1.6 million
galaxies derived from MICEv2 adopting spectroscopic redshift uncertainties as
in Eq. 22.

  * •

Full ($z_{\rm phot}$): full analysis as above, adopting broader (photometric)
redshift uncertainties as in Eq. 22.

  * •

Spectral: no galaxy catalog information is incorporated and instead galaxies
are assumed to be uniformly distributed in comoving volume. The constraints
are based only on the assumption of the source-frame mass distribution.

In this way, we obtain a total of six configurations. We adopt wide uniform
priors for all the hyperparameters (Table 1). The posterior distribution is
sampled with affine-invariant Markov Chain Monte Carlo sampler emcee (Foreman-
Mackey et al., 2013) and the convergence is assessed ensuring that the number
of samples is at least 50 larger than the integrated autocorrelation time for
all the parameters. Results are reported using the median statistic with
symmetric 68.3% credible levels.

The final constraints on all individual parameters are reported in Table 2,
while the marginalized posterior distribution for a selection of cosmological
($H_{0}$), mass ($\alpha$, $\mu_{\rm g}$, $\sigma_{\rm g}$), and rate
($\gamma$) parameters is shown in Fig. 4. In Fig. 5 we show the relative
uncertainty on $H_{0}$ in all configurations. Finally, in Appendix D we show
an example of the full corner plot. For all the six configurations, with
CHIMERA we recover the fiducial values with a typical deviation of
$0.2~{}\sigma$, with fluctuations that can be ascribed to the specific
realizations of the data sets.

Figure 5: Relative uncertainty on $H_{0}$ obtained from spectral-only and full
standard sirens analysis of 100 BBHs in the O4- and O5-like network
configurations, including a complete galaxy catalog with photometric or
spectroscopic redshift uncertainties. Figure 6: Joint cosmological and
astrophysical constraints from the full standard sirens analysis of 100 BBHs
in the O4-like (blue) and O5-like (red) configurations. Here we show, for
illustrative purposes, the most relevant parameters, namely
$\boldsymbol{\lambda}_{\rm c}=\\{H_{0}\\}$, the BBH mass function parameters
$\boldsymbol{\lambda}_{\rm m}=\\{\mu_{\rm g},\sigma_{\rm g}\\}$, and the rate
evolution parameter $\boldsymbol{\lambda}_{\rm z}=\\{\gamma\\}$, while in App.
D we show for completeness the entire distribution of parameters. The contours
represent the $1$ and $2\,\sigma$ confidence levels. The dotted lines indicate
the fiducial values adopted. Figure 7: Posterior on the Power Law + Peak BBH
primary mass function from the full standard analysis of 100 BBH in the
O4-like (blue) and O5-like (red) configurations. The colored lines and bands
represent the median and 90% credible level. The black dashed line is the mass
function evaluated at the fiducial hyperparameters.

###  4.1 Full O4- and O5-like scenarios

We start by discussing the results for the best-case scenario, consisting of a
complete galaxy catalog with spectroscopic redshift measurements. The marginal
posterior distributions for the selection of parameters
$\\{H_{0},\alpha,\mu_{\rm g},\sigma_{\rm g},\gamma\\}$ are shown in Fig. 6.

We find that in one year of observations in the O4-like configuration, the LVK
data are able to constrain $H_{0}$ with $7\%$ uncertainty (at the $1\sigma$
level) from BBH in a combined cosmological and astrophysical population
inference. This is a remarkable improvement with respect to the current
constraints, as the analysis of the 42 BBHs at $\mathrm{SNR}>12$ observed so
far yields a $46\%$ measurement of $H_{0}$ (Mastrogiovanni et al., 2023). To
arbitrate the Hubble tension, percent-level measurements are required. If we
assume that the uncertainty on $H_{0}$ scales as $1/\sqrt{N}$, in the O4
configuration it is not possible to reach $1\%$ in the planned schedule of
about two years.

In contrast, with the O5-like configuration it is possible to reach $1\%$
uncertainty in about one year of operation. We stress that this is a best-case
scenario, relying on having a complete catalog of potential hosts. In general,
the actual completeness can vary among different galaxy surveys and galaxy
types. For example, with Euclid it would change between the photometric or the
spectroscopic survey mode and the north or south direction, requiring a more
detailed assessment in a future study. Moreover, even if we marginalise on the
astrophysical parameters, the results still rely on an assumption on the
functional form of the BBH mass distribution. This dependence should be
investigated in more detail in the future. On the other side, this analysis is
based only on the BBH population, and further improvements can be obtained
including NSBH and BNS events and their potential electromagnetic
counterparts.

We now move to the population hyperparameters. Currently, in population
studies the cosmology is typically fixed (e.g. Abbott et al., 2023b). Here we
study how well the fiducial models are recovered when taking into account
potential correlations with cosmological hyperparameters. In Fig. 7 we show
the reconstructed primary mass distribution. For both O4 and O5 scenarios, the
Gaussian peak at $34~{}\mathrm{M}_{\odot}$ is clearly visible and its mean
value $\mu_{\mathrm{g}}$ is recovered with a precision of $7\%$ and $3\%$,
respectively. The second best-constrained mass parameter is the slope $\alpha$
of the primary BBH mass distribution that is recovered with fractional
uncertainties of $15\%$ and $13\%$, respectively. Overall, for the mass
function parameters these are small improvements with respect to the full
cosmological and astrophysical inference with the current GWTC-3 catalog,
which yields $\sigma_{\mu_{\mathrm{g}}}/\mu_{\mathrm{g}}\sim 13\%$ and
$\sigma_{\alpha}/\alpha\sim 11\%$ (see Mastrogiovanni et al., 2023). The same
is true for the population analysis at fixed cosmology, which gives
$\sigma_{\mu_{\mathrm{g}}}/\mu_{\mathrm{g}}\sim 9\%$ and
$\sigma_{\alpha}/\alpha\sim 11\%$ (Abbott et al., 2023b). Similarly, the
constraints on the rate parameter $\gamma$ remain essentially unchanged with
respect to the current knowledge, even when considering the transition from O4
to O5. This can be explained by the higher SNR threshold adopted in our O5
catalog, resulting in the same number of GW events that map a redshift range
comparable to that of O4.

In conclusion, we recall that our results are based on the full astrophysical
and cosmological analysis of the best 100 GW events detectable in one year for
each configuration. Population studies typically benefit from the inclusion of
all confident events detected and are carried out by fixing the cosmological
parameters (e.g., Abbott et al., 2023b). In this sense, our population
constraints should not be taken as representative of the overall performance
of O4 and O5.

###  4.2 Spectroscopic vs. photometric galaxy catalog

While obtaining a complete spectroscopic catalog poses challenges and awaits
future facilities, ongoing surveys such as DES and Euclid are already building
extensive photometric galaxy catalogs. The full $z_{\rm phot}$ and $z_{\rm
spec}$ configurations for both the O4-like and O5-like configurations are
compared in Fig. 4.

With photometric redshifts the constraint on $H_{0}$ for O4 is notably less
accurate, with a measurement uncertainty that is three times greater
($\sigma_{H_{0}}/H_{0}\approx 18\%$) compared to the spectroscopic approach.
In the case of O5 this factor increases to nine ($\sigma_{H_{0}}/H_{0}\approx
9\%$). Interestingly, this shows that from 100 O5 events at $\mathrm{SNR}>25$
it is not possible to achieve percent-level precision on $H_{0}$ using a
photometric catalog, even under the assumption of completeness. Of course, one
may consider lowering the SNR threshold to include more events; however, such
events would also have much larger localization volumes, thus a large number
of potential hosts. This would likely limit their additional constraining
power. Some information can be still retrieved for the mass distribution,
whose reconstruction benefits from a larger sample; however, in the next
section we will show that the constraints obtained in the absence of a galaxy
catalog with our sample are only marginally worse than what can be obtained
from a sample with a much lower SNR cut. This seems to suggest that it would
be difficult to reach percent-level accuracy.

Another interesting result concerns the comparison of the Full ($z_{\rm
phot}$) O5-like and Full ($z_{\rm spec}$) O4-like configurations (see Fig. 5).
We find that considering a galaxy catalog with spectroscopic redshift
uncertainties in the O4-like scenario, we are able to achieve a more precise
constraint on $H_{0}$ compared to having a larger LVK detector network at O5
design sensitivity with a photometric galaxy catalog. This occurs despite the
factor 10 improvement in the GW localization volume (see Fig. 3). Overall,
these findings underline the importance of mapping the GW localization volume
- at least for well-localized events - with dedicated spectroscopic surveys.

###  4.3 Spectral-only analysis

In the absence of the galaxy catalog information (spectral sirens case), the
cosmological constraints are determined by the capability of the source-frame
mass function to break the mass-redshift degeneracy.

In our spectral siren analysis, the Hubble parameter is recovered at
$\sigma_{H_{0}}/H_{0}\approx 43\%$ in O4-like and $\sigma_{H_{0}}/H_{0}\approx
32\%$ in O5-like configuration. Our results are in good agreement with those
of Mancarella et al. (2022) and Leyde et al. (2022). In particular, Leyde et
al. (2022) find $38\%$ ($24\%$) uncertainty on $H_{0}$ using $\mathrm{SNR}>12$
spectral siren events for O4 (O5) obtaining a total of 87 (423) events. When
comparing with our results for O5, we must consider that we applied a higher
$\mathrm{SNR}>25$ threshold, resulting in a factor of $4$ difference in the
number of detected events. It is interesting to note that this large
difference results only in a marginal improvement in the constraining power,
suggesting that the constraints are mainly driven by the events with higher
SNR.

In general, we conclude that the analysis of 100 well-localized GW events with
a complete galaxy catalog with spectroscopic measurements, would provide
better constraints with respect to a pure spectral siren analysis under the
assumed mass functions. This result holds true also when comparing with the
five-year results of O5 spectral sirens by Farr et al. (2019), providing a
relative error on $H(z)$ of $\approx 3\%$ at a pivot redshift of $z\simeq
0.8$. We note however that the authors assumed a much more optimistic event
rate and a mass function with sharper edges, in line with the knowledge of the
BBH population at the time.

Figure 8: Constraints on the $H_{0}-\mu_{\rm g}$ plane for the O5-like
scenario in case of spectral, Full ($z_{\rm phot}$), Full ($z_{\rm spec}$)
analyses with CHIMERA.

###  4.4 The $H_{0}-\mu_{\rm g}$ correlation

In the absence of a galaxy catalog, there is a strong correlation between
$H_{0}$ and $\mu_{\rm g}$, which drives the constraint on $H_{0}$. In this
section, we comment on this correlation and compare the three scenarios
considered in the paper. Figure 8 shows the constraints in the $H_{0}-\mu_{\rm
g}$ plane in the O5-like configuration.

A relatively strong anti-correlation is present between $\mu_{\rm g}$ and
$H_{0}$ in both the Spectral and Full ($z_{\rm phot}$) analyses. In terms the
Spearman’s rank correlation coefficient computed from the MCMC chains, these
correspond to $\rho\sim-0.8$ and $\rho\sim-0.6$, respectively.

This is a consequence of the fact that, in the absence of a precise redshift
measurement, a higher value of $H_{0}$ would move the events with a given
luminosity distance at higher $z$, requiring smaller values of the source-
frame masses to reproduce the observed data (see Eq. 4). This leads to a
smaller value of $\mu_{\rm g}$. This anti-correlation is observed also in the
latest GWTC-3 analyses (Abbott et al., 2023a). Having a complete spectroscopic
galaxy catalog allows instead to break the $H_{0}-\mu_{\rm g}$ degeneracy,
constraining $H_{0}$ with higher precision. This is visible in Fig. 8 where no
significant correlation is present in the Full ($z_{\rm spec}$) case.
Interestingly, in the case where only photometric redshifts are available,
this result shows that the information coming from the mass distribution still
plays a role in the measurement of $H_{0}$ and has therefore to be accurately
modeled. In particular, Fig. 8 shows that a wrong assumption on the value of
the mass scale would lead to a bias on $H_{0}$ even in the presence of a
complete photometric galaxy survey.

##  5 Conclusions

In this work, we study the scientific potential of dark sirens in the context
of a joint astrophysical and cosmological parameter inference including galaxy
catalogs. Here we summarize the main results of the paper.

  1. 1.

We present CHIMERA999Available at https://github.com/CosmoStatGW/CHIMERA, a
novel Python code for gravitational-wave cosmology combining information from
the population properties of compact binary mergers and galaxy catalogs. The
code is designed to be accurate in different scenarios, from the spectral
siren (population inference), to the dark siren (galaxy catalog) and the
bright siren (redshift precisely measured after the identification of an
electromagnetic counterpart) methods, and efficient for large galaxy catalogs.

  2. 2.

We generate two catalogs of BBH events for detector networks representative of
the LIGO-Virgo-KAGRA O4 and O5 observing runs. We start from a parent catalog
of galaxies extracted from the MICEv2 mock galaxy catalog, that we select with
a cut in stellar mass $\log M_{\star}/\mathrm{M}_{\odot}>10.5$. From this
catalog, we generate BBH events with the GWFAST code, considering 100 best
events that can be detected over approximately one year of observation, namely
$\mathrm{SNR}>12$ BBHs for the O4-like configuration and 100 $\mathrm{SNR}>25$
BBHs for the O5-like configuration. We then associate the galaxy redshift
catalog with two different uncertainties on the redshift, representative of
the case of having a catalog of photometric and spectroscopic redshifts.

  3. 3.

We find that in the best-case scenario of a complete spectroscopic galaxy
catalog, the Hubble constant $H_{0}$ can be constrained to 7% with O4-like and
1% with O5-like in about a year of observations, demonstrating the strong
potential of joint standard sirens and galaxy catalog analyses in the context
of the Hubble tension using novel independent probes.

  4. 4.

With a photometric galaxy catalog, the constraints on $H_{0}$ are notably
weaker, i.e. by a factor of $\sim 3$ for O4-like and $\sim 9$ for O5-like. In
this case, the correlation between $H_{0}$ and the mass scales in the BBH
population - which is intrinsic in the absence of a galaxy catalog - is not
fully broken. It is therefore crucial to accurately model the mass
distribution to avoid biases in the cosmological parameters.

  5. 5.

Interestingly, we find that the O4-like configuration with a spectroscopic
catalog provides a more precise measurement of $H_{0}$ than the O5-like
configuration with photometric redshifts. This points to the importance of
having spectroscopic redshift measurements for GW cosmology, in the absence of
which the potential of a factor $\sim 10$ improvement in the localization
capabilities between O4-like and O5-like GW detector networks could be
completely lost.

The results obtained in this work are based on some simplifying assumptions.
First of all, the catalog of potential hosts that we employ includes only the
more massive galaxies, in line with the idea that the binary merger rate is
traced by stellar mass. Second, the catalog is assumed to be complete. Despite
these limitations, the crucial role that the spectroscopic surveys of galaxies
can have in gravitational wave cosmology is already evident. Moreover, while
obtaining a complete catalog of less massive galaxies is challenging at
present, this will become a concrete possibility in the coming years thanks to
ongoing (e.g., DESI, DESI Collaboration et al. 2016; and Euclid, Laureijs et
al. 2011), future (e.g., WFIRST, Akeson et al. 2019; LSST, LSST Science
Collaboration et al. 2009) and proposed (e.g., WST, Ellis & Dawson 2019) large
galaxy surveys. In this context, further investigations are needed to further
study the impact of incompleteness in the observed catalog.

Potential improvements may come from two directions. Spectroscopic targeted
searches may be carried out, also a posteriori, in the localization volumes of
best-localized events. To study the feasibility of this approach a
comprehensive study should be carried out, considering both the observational
requirements and the optimal GW localization volumes within which this could
be done. Even if this approach is not feasible, the additional information
from observational and physical galaxy properties in combination with binary
population studies (e.g. Santoliquido et al., 2022; Vijaykumar et al., 2023a)
can be included in the inference, potentially improving the constraints.

From the GW data side, a second assumption, which mostly impacts the
photometric and spectral-only analyses, is the specific distributions for the
binary population. In fact, even if we marginalize on the astrophysical
parameters, the results still rely on the functional form of the BBH mass
distribution and the presence of features that can be used as standard rulers
to obtain cosmological results. This dependence should be investigated in more
detail in the future. We also emphasize that the GW catalogs analyzed in this
work are built assuming a 100% duty-cycle, though, in reality, detectors may
undergo temporary maintenance and upgrade phases, requiring the modeling of
the specific duty cycle for each detector. Moreover, we recall that these
results are obtained assuming at least three GW detectors working at nominal
sensitivity, which is crucial in order to achieve a sufficiently precise
localization.

In the future, it will be interesting to expand this analysis to 3G networks,
like the Einstein Telescope (Punturo et al., 2010; Branchesi et al., 2023),
LISA (Amaro-Seoane et al., 2017), and Cosmic Explorer (Reitze et al., 2019),
exploring deeper galaxy catalogs needed to properly cover the detectability
range of those detectors. We will address the case of future GW experiments in
a separate study.

To conclude, our results highlight that the synergy between future GW
detectors and large galaxy surveys is crucial to fully attain the potential of
dark sirens analyses. In particular, they suggest the importance of targeted
spectroscopic campaigns to map the galaxy distribution within the GW
localization volumes, at least for the best-localized events.

We thank Simone Mastrogiovanni, Sofia Contarini, Massimiliano Romanello,
Giorgio Lesci, Alberto Traina, and Laura Leuzzi for useful discussions. We
acknowledge the use of computational resources from the parallel computing
cluster of the Open Physics Hub (https://site.unibo.it/openphysicshub/en) at
the Physics and Astronomy Department in Bologna. N.B. and M.Mor. acknowledge
support from MIUR, PRIN 2017 (grant 20179ZF5KS). M.Mor. and A.C. acknowledge
the grants ASI n.I/023/12/0 and ASI n.2018-23-HH.0. A.C. acknowledges the
support from grant PRIN MIUR 2017 - 20173ML3WW_001. M.Mor. acknowledges the
support from grant PRIN MIUR 2022 (grant 2022NY2ZRS_001) M.T. acknowledges the
funding by the European Union - NextGenerationEU, in the framework of the HPC
project – “National Centre for HPC, Big Data and Quantum Computing” (PNRR -
M4C2 - I1.4 - CN00000013 – CUP J33C22001170001)). The work of M.Manc. is
supported by European Union’s H2020 ERC Starting Grant No. 945155–GWmining,
Cariplo Foundation Grant No. 2021-0555, the ICSC National Research Centre
funded by NextGenerationEU, and MIUR PRIN Grant No. 2022-Z9X4XS. The work of
F.I. and M.Mag. is supported by the Swiss National Science Foundation, grant
200020$\\_$191957, and by the SwissMap National Center for Competence in
Research. The views and opinions expressed are those of the authors only and
do not necessarily reflect those of the European Union or the European
Commission. Neither the European Union nor the European Commission can be held
responsible for them. This work is (partially) supported by ICSC – Centro
Nazionale di Ricerca in High Performance Computing, Big Data and Quantum
Computing, funded by European Union – NextGenerationEU.

##  Appendix A CHIMERA: code implementation

CHIMERA (Combined Hierarchical Inference Model for Electromagnetic and
gRavitational wave Analysis) is a new, flexible Python code that enables
jointly fitting the cosmological and astrophysical population parameters by
using information from galaxy catalogs. The code is designed to be accurate
for different scenarios, encompassing bright, dark, and spectral sirens
methods, and computationally efficient in view of next-generation GW
observatories and upcoming galaxy surveys. The ultimate goal of CHIMERA is to
compute the result of Eq. 17 in a computationally efficient way in view of
next-generation GW observatories and upcoming galaxy surveys. For this reason,
we use the LAX-backend implementation and Just in time computation
capabilities of JAX (Bradbury et al., 2018). The main functions are developed
starting from DarkSirensStat (Finke et al., 2021) and MGCosmoPop (Mancarella
et al., 2022). The workflow of CHIMERA is illustrated in Fig. 1.

Figure 9: The flowchart of CHIMERA. Main modules include functions related to
the full likelihood computation (violet), computation of probabilities (red),
data file I/O (yellow), and model functions (gray).

In this appendix, we provide an overview of the code. As standard units for
CHIMERA, we adopt radians for the angles and Gpc for the distances. Data sets
are given in the form of dictionaries of numpy arrays. All model functions
accept one argument in the form of a dictionary containing the related
parameters.

###  A.1 Full likelihood

The core modules of CHIMERA are Likelihood.py and Bias.py. In particular, the
file Likelihood.py contains classes to compute the product of the integrals in
Eq. 17, while Bias.py allows for the computation of the selection bias term
$\xi(\boldsymbol{\lambda})$. The computation of the likelihood uses functions
to analyze the gravitational wave (from GW.py) and electromagnetic (from
EM.py) information. These modules contain specific methods to perform all the
preliminary computations that do not change during the likelihood evaluation
(a more extended discussion can be found in the next paragraph). Data are
loaded with specific classes present in DataGW.py (e.g., DataGWMock and
DataLVK) and DataEM.py (e.g., MockGalaxiesMICEv2 and GLADEPlus), respectively.
The algorithm to compute the Likelihood proceeds as follows:

  1. 1.

First of all, the class MockLike stores all the population models (model_*),
GW data (data_GW*), galaxy data (data_GAL*), pixelization parameters
(nside_list, npix_event, sky_conf), and integration parameters
(z_int_H0_prior, z_int_sigma, z_int_res).

  2. 2.

The GW class is initialized and the pixelization and redshift grids are pre-
computed. To optimize the computation in case of large galaxy catalog
analysis, the code not only restricts the sky localization, but also the
redshift integration grid. The first task is done starting from posterior
distributions in RA and Dec. The user can choose the approximate number of
pixels desired for each event (npix_event) to be found within a confidence
level ellipse (sky_conf), given a list of possible pixelizations (nside_list).
CHIMERA optimizes the pixelization of each event (in particular, the nside
parameter of HealPix) to obtain the number of pixels closest to npix_event.
The second task must take into account that by varying the cosmology, also
$\mathcal{K}_{\mathrm{gw}}$ varies, so to avoid biases, the integration grid
must encompass all the redshift range explored during the inference. This is
obtained starting from the GW posteriors on $d_{L}$ and defining the range
inside z_int_sigma standard deviations at a resolution of z_int_res spanning
all the range101010In practice the grid starts at
$z(d_{L}-N\,\sigma_{d_{L}}\,|\,H_{0,\mathrm{min}})$ and ends at
$z(d_{L}+N\,\sigma_{d_{L}}\,|\,H_{0,\mathrm{max}})$, where $N$ is the number
of standard deviations (z_int_sigma) and
$[H_{0,\mathrm{min}},H_{0,\mathrm{max}}]$ the interval of $H_{0}$ spanned in
the inference. of $H_{0}$ explored.

  3. 3.

The EM class is initialized, and the quantity $p_{\rm
cat}(z,\hat{\Omega}\,|\,\boldsymbol{\lambda}_{\rm c})$ (Eqs. 11 and 12) is
pre-computed pixel-by-pixel on the redshift grids. At this point, based on the
included catalog (e.g., GLADE+, MICEv2), it is possible to activate catalog-
related tasks (e.g., luminosity cut, or associate user-defined redshift
uncertainties). Similarly, if defined, the completeness is computed, and the
pixelized completeness function $P_{\rm compl}(z)$ (see the discussion in the
above paragraph) is stored in the class to be accessible during the inference.

The pixelization approach and the pre-computation of $p_{\mathrm{gal}}$ are
essential for next-gen galaxy surveys. The first allows us to reduce the
dimension, combining the probability of $10-10^{3}$ galaxies in one single
pixel ensuring granularity in the sky analysis. The second allows us to avoid
the computation of $p_{\mathrm{gal}}$ at each step of the MCMC inference.
While this last approximation is not a problem for $H_{0}$ inference since
both the numerator and denominator of Eq. 12 have a $H_{0}^{3}$ dependence,
the impact of a possible bias on $\Omega_{m}$ should be carefully assessed
when future analyses using events at higher $z$ will be carried out. The
algorithm to compute the selection Bias proceeds as follows:

  1. 1.

First of all, the class Bias stores all the population models (model_*), the
directory of the GW injection catalog data (file_inj), and the SNR threshold
to be applied. Optional arguments also include the catalog interpolant
function. If not given, the bias is evaluated for a uniform in comoving volume
galaxy distribution (i.e., $\propto\mathrm{d}{V_{C}}/\mathrm{d}{z}$).

  2. 2.

The injection catalog is loaded by applying the chosen SNR cut. It is
important to ensure that this cut is equivalent to the one adopted when
creating the catalog of GW events to analyze.

Thus, following Eq. 1, the full likelihood is conveniently derived by first
calculating the log-likelihood for all the events and then subtracting the
log-bias term, which is computed only once and multiplied by the number of
events. The computation of the likelihood and selection bias is performed by
calling the .compute() methods or, in logarithmic form, the .compute_ln()
methods. In the latter case, the models should also be given in logarithmic
form.

###  A.2 Included models

Currently, CHIMERA includes the following population models:

  * •

Mass distributions (mass.py): logpdf_TPL (Truncated Power Law), logpdf_BPL
(Broken Power Law), logpdf_PLP (Power Law + Peak), logpdf_PL2P (Power Law + 2
Peaks).

  * •

Rate evolutions (rate.py): logphi_PL (Power Law), logphi_MD (Madau-like).

  * •

Spin distributions (spin.py): logpdf_G (Gaussian), logpdf_U (Uniform).

  * •

Cosmological models (cosmo.py): fLCDM, fLCDM modified gravity

All the above functions accept parameters in the form of dictionaries, e.g.
lambda_cosmo = {'H0':70.0, 'Om0':0.3}.

##  Appendix B Choices for the analysis of GW170817

In this appendix we provide details about the analysis setup for the example
of GW170817 in Sec. 2.4.

We use the GLADE+ galaxy catalog (Dálya et al., 2022), which includes data
from multiple galaxy surveys (GWGC, 2MPZ, 2MASS XSC, HyperLEDA, WISExSCOSPZ,
and SDSS-DR16Q quasars) and is largely employed in standard sirens analyses.
As done in previous analyses, we select galaxies with a luminosity above a
certain threshold to obtain a more complete catalog of potential hosts. The
implicit assumption is that the true host lies above the chosen threshold.
This is well motivated by the fact that the luminosity is a reasonable proxy
for the stellar mass (in K-band) and the star formation rate (in B-band),
therefore more luminous galaxies have a higher chance to host mergers. In this
case, we use K-band luminosities with a threshold $L_{K}>0.02\,L_{K}^{\ast}$,
where $L_{K}^{\ast}\simeq 1.1\times 10^{11}~{}\mathrm{M}_{\odot}$ (Kochanek et
al., 2001) denotes the characteristic K-band luminosity. In this analysis we
weight each galaxy by its luminosity, $w_{g}\propto L_{k}$ (see Eq. 11). The
completeness is computed with the mask method, as outlined in Section 2.4. The
sky is pixelized with a size of 0.83 deg2 and pixels are organized into nine
distinct masks based on the galaxy counts in each pixel. For each mask, the
completeness fraction $P_{c}(z,\hat{\Omega})$ is computed by comparing the
luminosity density to a reference Schechter function with parameters from
Kochanek et al. (2001). We adopt the following population model: flat mass
distribution between $1$ and $3~{}\mathrm{M}_{\odot}$ and power-law rate
evolution with slope $\gamma=2.7$.Then, we fix the population parameters and
study the posterior on $H_{0}$. For the inference, we consider each galaxy
within the $90\%$ GW sky localization area (about 28 deg2) as a potential host
(Fig. 2, left).

##  Appendix C Comparison of CHIMERA with MGCosmoPop

We compare CHIMERA in the spectral-only case (population analysis) with
MGCosmoPop111111Available at https://github.com/CosmoStatGW/MGCosmoPop. The
code is presented and used in Mancarella et al. (2022) to provide joint
cosmological and astrophysical constraints, including deviation from general-
relativity, in the empty catalog case. The codes feature similar model
functions implementations, except for the cosmology class, which has been
rewritten in CHIMERA to improve computational efficiency. We compare these
functions across a wide range of parameters, finding machine-level
differences.

On the contrary, the two codes feature distinct implementations of the
likelihood evaluation. In MGCosmoPop the integral in Eq. 1 is performed with a
Monte Carlo approach, while in CHIMERA we pixelize the sky, evaluate the
redshift prior inside each pixel, and perform the integral on the redshift
grid (see Fig. 1). While this approach has the advantage of improving the
computational efficiency with large galaxy catalogs (Gray et al., 2022) and
the evaluation of their completeness (Finke et al., 2021), it is crucial to
perform a comparative analysis to better assess potential biases.

Figure 10: Posterior distributions obtained in the case of the Spectral Siren
analysis with CHIMERA and MGCosmoPop using the O5-like catalog.

We use the set of 100 BBH events at $\mathrm{SNR}>25$ presented in Section
3.2. We assume a flat $\Lambda$CDM cosmology, PLP mass distribution, and
Madau-like rate evolution, with parameters and priors reported in Table 1. In
CHIMERA, we set npix_event to 15 and nside to $[2^{n}]$ with $n=(3,...,9)$ so
that the pixelization of each event is automatically adjusted to have
approximately 15 pixels in the 90% credible sky area. Then, we perform full
MCMC analyses with CHIMERA and MGCosmoPop. The results are shown in Fig. 10.
Overall, the posteriors obtained with the two codes are in excellent
agreement.

##  Appendix D Full corner plot

In Fig. 11 we show the full corner plot with posteriors on the 12 parameters
$\boldsymbol{\lambda}$ in case of galaxy catalog with spectroscopic redshift
measurements. In this scenarios, as in the other four configurations (see
Table 2), we are able to recover the fiducial values for almost all of the
hyperparameters within the 68% credible levels.

Figure 11: Joint cosmological and astrophysical constraints from the full
standard sirens analysis of 100 BBHs in the O4 (blue) and O5 (red)
configurations for the full parameter list described in Tab. 1. The contours
represent the 68% and 95% confidence levels. The dotted lines indicate the
fiducial values adopted.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Living Reviews in Relativity, 19, 1, doi: 10.1007/lrr-2016-1
  * Abbott et al. (2017a) —. 2017a, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2017c) —. 2017c, Nature, 551, 85, doi: 10.1038/nature24471
  * Abbott et al. (2021a) —. 2021a, ApJ, 909, 218, doi: 10.3847/1538-4357/abdcb7
  * Abbott et al. (2021b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, Physical Review X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2021c) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021c, arXiv e-prints, arXiv:2111.03606, doi: 10.48550/ARXIV.2111.03606
  * Abbott et al. (2021d) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021d, ApJ, 913, L7, doi: 10.3847/2041-8213/abe949
  * Abbott et al. (2023a) Abbott, R., Abe, H., Acernese, F., et al. 2023a, ApJ, 949, 76, doi: 10.3847/1538-4357/ac74bb
  * Abbott et al. (2023b) Abbott, R., Abbott, T. D., Acernese, F., et al. 2023b, Physical Review X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * Akeson et al. (2019) Akeson, R., Armus, L., Bachelet, E., et al. 2019, arXiv e-prints, arXiv:1902.05569, doi: 10.48550/arXiv.1902.05569
  * Alfradique et al. (2023) Alfradique, V., Bom, C. R., Palmese, A., et al. 2023, arXiv e-prints, arXiv:2310.13695.  https://arxiv.org/abs/2310.13695
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786, doi: 10.48550/arXiv.1702.00786
  * Amendola et al. (2018) Amendola, L., Sawicki, I., Kunz, M., & Saltas, I. D. 2018, J. Cosmology Astropart. Phys, 2018, 030, doi: 10.1088/1475-7516/2018/08/030
  * Baibhav et al. (2019) Baibhav, V., Berti, E., Gerosa, D., et al. 2019, Phys. Rev. D, 100, 064060, doi: 10.1103/PhysRevD.100.064060
  * Ballard et al. (2023) Ballard, W., Palmese, A., Magaña Hernandez, I., et al. 2023, arXiv e-prints, arXiv:2311.13062, doi: 10.48550/arXiv.2311.13062
  * Belgacem et al. (2018a) Belgacem, E., Dirian, Y., Foffa, S., & Maggiore, M. 2018a, Phys. Rev. D, 97, 104066, doi: 10.1103/PhysRevD.97.104066
  * Belgacem et al. (2018b) —. 2018b, Phys. Rev. D, 98, 023510, doi: 10.1103/PhysRevD.98.023510
  * Belgacem et al. (2019) Belgacem, E., Calcagni, G., Crisostomi, M., et al. 2019, J. Cosmology Astropart. Phys, 2019, 024, doi: 10.1088/1475-7516/2019/07/024
  * Bera et al. (2020) Bera, S., Rana, D., More, S., & Bose, S. 2020, ApJ, 902, 79, doi: 10.3847/1538-4357/abb4e0
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs.  http://github.com/google/jax
  * Branchesi et al. (2023) Branchesi, M., Maggiore, M., Alonso, D., et al. 2023, J. Cosmology Astropart. Phys, 2023, 068, doi: 10.1088/1475-7516/2023/07/068
  * Buchs et al. (2019) Buchs, R., Davis, C., Gruen, D., et al. 2019, MNRAS, 489, 820, doi: 10.1093/mnras/stz2162
  * Carretero et al. (2015) Carretero, J., Castander, F. J., Gaztañaga, E., Crocce, M., & Fosalba, P. 2015, MNRAS, 447, 646, doi: 10.1093/mnras/stu2402
  * Chen et al. (2023) Chen, A., Gray, R., & Baker, T. 2023, arXiv e-prints, arXiv:2309.03833, doi: 10.48550/arXiv.2309.03833
  * Chen et al. (2018) Chen, H.-Y., Fishbach, M., & Holz, D. E. 2018, Nature, 562, 545, doi: 10.1038/s41586-018-0606-0
  * Chernoff & Finn (1993) Chernoff, D. F., & Finn, L. S. 1993, ApJ, 411, L5, doi: 10.1086/186898
  * Colombo et al. (2022) Colombo, A., Salafia, O. S., Gabrielli, F., et al. 2022, ApJ, 937, 79, doi: 10.3847/1538-4357/ac8d00
  * Colombo et al. (2023) Colombo, A., Duqué, R., Sharan Salafia, O., et al. 2023, arXiv e-prints, arXiv:2310.16894, doi: 10.48550/arXiv.2310.16894
  * Crocce et al. (2015) Crocce, M., Castander, F. J., Gaztañaga, E., Fosalba, P., & Carretero, J. 2015, MNRAS, 453, 1513, doi: 10.1093/mnras/stv1708
  * Dalang & Baker (2023) Dalang, C., & Baker, T. 2023, arXiv e-prints, arXiv:2310.08991, doi: 10.48550/arXiv.2310.08991
  * Dálya et al. (2018) Dálya, G., Galgóczi, G., Dobos, L., et al. 2018, MNRAS, 479, 2374, doi: 10.1093/mnras/sty1703
  * Dálya et al. (2022) Dálya, G., Díaz, R., Bouchet, F. R., et al. 2022, MNRAS, 514, 1403, doi: 10.1093/mnras/stac1443
  * Del Pozzo (2012) Del Pozzo, W. 2012, Phys. Rev. D, 86, 043011, doi: 10.1103/PhysRevD.86.043011
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv e-prints, arXiv:1611.00036, doi: 10.48550/arXiv.1611.00036
  * Dhani et al. (2022) Dhani, A., Borhanian, S., Gupta, A., & Sathyaprakash, B. 2022, arXiv e-prints, arXiv:2212.13183, doi: 10.48550/arXiv.2212.13183
  * Ellis & Dawson (2019) Ellis, R., & Dawson, K. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 45, doi: 10.48550/arXiv.1907.06797
  * Euclid Collaboration: Desprez et al. (2020) Euclid Collaboration: Desprez, G., Paltani, S., Coupon, J., et al. 2020, A&A, 644, A31, doi: 10.1051/0004-6361/202039403
  * Euclid Collaboration: Schirmer et al. (2022) Euclid Collaboration: Schirmer, M., Jahnke, K., Seidel, G., et al. 2022, A&A, 662, A92, doi: 10.1051/0004-6361/202142897
  * Ezquiaga (2021) Ezquiaga, J. M. 2021, Physics Letters B, 822, 136665, doi: 10.1016/j.physletb.2021.136665
  * Ezquiaga & Holz (2022) Ezquiaga, J. M., & Holz, D. E. 2022, Phys. Rev. Lett., 129, 061102, doi: 10.1103/PhysRevLett.129.061102
  * 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
  * Finke et al. (2021) Finke, A., Foffa, S., Iacovelli, F., Maggiore, M., & Mancarella, M. 2021, J. Cosmology Astropart. Phys, 2021, 026, doi: 10.1088/1475-7516/2021/08/026
  * Fishbach et al. (2019) Fishbach, M., Gray, R., Magaña Hernandez, I., et al. 2019, ApJ, 871, L13, doi: 10.3847/2041-8213/aaf96e
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fosalba et al. (2015) Fosalba, P., Crocce, M., Gaztañaga, E., & Castander, F. J. 2015, MNRAS, 448, 2987, doi: 10.1093/mnras/stv138
  * Fosalba et al. (2015) Fosalba, P., Gaztañaga, E., Castander, F. J., & Crocce, M. 2015, MNRAS, 447, 1319, doi: 10.1093/mnras/stu2464
  * Gair et al. (2023) Gair, J. R., Ghosh, A., Gray, R., et al. 2023, The Astronomical Journal, 166, 22, doi: 10.3847/1538-3881/acca78
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Gray et al. (2022) Gray, R., Messenger, C., & Veitch, J. 2022, Monthly Notices of the Royal Astronomical Society, 512, 1127, doi: 10.1093/mnras/stac366
  * Gray et al. (2020) Gray, R., Hernandez, I. M., Qi, H., et al. 2020, Phys. Rev. D, 101, 122001, doi: 10.1103/PhysRevD.101.122001
  * Gray et al. (2023) Gray, R., Beirnaert, F., Karathanasis, C., et al. 2023, arXiv e-prints, arXiv:2308.02281, doi: 10.48550/arXiv.2308.02281
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hinton (2016) Hinton, S. R. 2016, The Journal of Open Source Software, 1, 00045, doi: 10.21105/joss.00045
  * Hoffmann et al. (2015) Hoffmann, K., Bel, J., Gaztañaga, E., et al. 2015, MNRAS, 447, 1724, doi: 10.1093/mnras/stu2492
  * Holz & Hughes (2005) Holz, D. E., & Hughes, S. A. 2005, ApJ, 629, 15, doi: 10.1086/431341
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iacovelli et al. (2022a) Iacovelli, F., Mancarella, M., Foffa, S., & Maggiore, M. 2022a, ApJ, 941, 208, doi: 10.3847/1538-4357/ac9cd4
  * Iacovelli et al. (2022b) —. 2022b, ApJS, 263, 2, doi: 10.3847/1538-4365/ac9129
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Kalaghatgi et al. (2020) Kalaghatgi, C., Hannam, M., & Raymond, V. 2020, Phys. Rev. D, 101, 103004, doi: 10.1103/PhysRevD.101.103004
  * Kochanek et al. (2001) Kochanek, C. S., Pahre, M. A., Falco, E. E., et al. 2001, ApJ, 560, 566, doi: 10.1086/322488
  * Laghi & Del Pozzo (2020) Laghi, D., & Del Pozzo, W. 2020, cosmolisa.  https://github.com/wdpozzo/cosmolisa
  * Laghi et al. (2021) Laghi, D., Tamanini, N., Del Pozzo, W., et al. 2021, MNRAS, 508, 4512, doi: 10.1093/mnras/stab2741
  * Lagos et al. (2019) Lagos, M., Fishbach, M., Landry, P., & Holz, D. E. 2019, Phys. Rev. D, 99, 083504, doi: 10.1103/PhysRevD.99.083504
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193, doi: 10.48550/arXiv.1110.3193
  * Leyde et al. (2022) Leyde, K., Mastrogiovanni, S., Steer, D. A., Chassande-Mottin, E., & Karathanasis, C. 2022, arXiv e-prints, arXiv:2202.00025, doi: 10.48550/arXiv.2202.00025
  * Liu et al. (2023) Liu, C., Laghi, D., & Tamanini, N. 2023, arXiv e-prints, arXiv:2310.12813, doi: 10.48550/arXiv.2310.12813
  * London et al. (2018) London, L., Khan, S., Fauchon-Jones, E., et al. 2018, Phys. Rev. Lett., 120, 161102, doi: 10.1103/PhysRevLett.120.161102
  * Loredo (2004) Loredo, T. J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 195–206, doi: 10.1063/1.1835214
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201, doi: 10.48550/arXiv.0912.0201
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Magana Hernandez (2021) Magana Hernandez, I. 2021, arXiv e-prints, arXiv:2112.07650, doi: 10.48550/arXiv.2112.07650
  * Mancarella et al. (2022) Mancarella, M., Genoud-Prachex, E., & Maggiore, M. 2022, Phys. Rev. D, 105, 064030, doi: 10.1103/PhysRevD.105.064030
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086, doi: 10.1093/mnras/stz896
  * Mastrogiovanni et al. (2021a) Mastrogiovanni, S., Haegel, L., Karathanasis, C., Magaña Hernandez, I., & Steer, D. A. 2021a, J. Cosmology Astropart. Phys, 2021, 043, doi: 10.1088/1475-7516/2021/02/043
  * Mastrogiovanni et al. (2021b) Mastrogiovanni, S., Leyde, K., Karathanasis, C., et al. 2021b, Phys. Rev. D, 104, 062009, doi: 10.1103/PhysRevD.104.062009
  * Mastrogiovanni et al. (2023) Mastrogiovanni, S., Laghi, D., Gray, R., et al. 2023, Phys. Rev. D, 108, 042002, doi: 10.1103/PhysRevD.108.042002
  * Messenger & Read (2012) Messenger, C., & Read, J. 2012, Phys. Rev. Lett., 108, 091101, doi: 10.1103/PhysRevLett.108.091101
  * Messenger et al. (2014) Messenger, C., Takami, K., Gossan, S., Rezzolla, L., & Sathyaprakash, B. S. 2014, Physical Review X, 4, 041004, doi: 10.1103/PhysRevX.4.041004
  * Moresco et al. (2022) Moresco, M., Amati, L., Amendola, L., et al. 2022, Living Reviews in Relativity, 25, 6, doi: 10.1007/s41114-022-00040-z
  * Mukherjee et al. (2021) Mukherjee, S., Wandelt, B. D., Nissanke, S. M., & Silvestri, A. 2021, Phys. Rev. D, 103, 043520, doi: 10.1103/PhysRevD.103.043520
  * Muttoni et al. (2023) Muttoni, N., Laghi, D., Tamanini, N., Marsat, S., & Izquierdo-Villalba, D. 2023, Phys. Rev. D, 108, 043543, doi: 10.1103/PhysRevD.108.043543
  * Muttoni et al. (2022) Muttoni, N., Mangiagli, A., Sesana, A., et al. 2022, Phys. Rev. D, 105, 043509, doi: 10.1103/PhysRevD.105.043509
  * Myles et al. (2021) Myles, J., Alarcon, A., Amon, A., et al. 2021, MNRAS, 505, 4249, doi: 10.1093/mnras/stab1515
  * Nair et al. (2018) Nair, R., Bose, S., & Saini, T. D. 2018, Phys. Rev. D, 98, 023502, doi: 10.1103/PhysRevD.98.023502
  * Nissanke et al. (2010) Nissanke, S., Holz, D. E., Hughes, S. A., Dalal, N., & Sievers, J. L. 2010, ApJ, 725, 496, doi: 10.1088/0004-637X/725/1/496
  * Oguri (2016) Oguri, M. 2016, Phys. Rev. D, 93, 083511, doi: 10.1103/PhysRevD.93.083511
  * Palmese et al. (2023) Palmese, A., Bom, C. R., Mucesh, S., & Hartley, W. G. 2023, ApJ, 943, 56, doi: 10.3847/1538-4357/aca6e3
  * Palmese et al. (2020) Palmese, A., deVicente, J., Pereira, M. E. S., et al. 2020, ApJ, 900, L33, doi: 10.3847/2041-8213/abaeff
  * Planck Collaboration (2020) Planck Collaboration. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Punturo et al. (2010) Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Classical and Quantum Gravity, 27, 194002, doi: 10.1088/0264-9381/27/19/194002
  * Raffai et al. (2023) Raffai, P., Pálfi, M., Dálya, G., & Gray, R. 2023, arXiv e-prints, arXiv:2310.16556, doi: 10.48550/arXiv.2310.16556
  * Reitze et al. (2019) Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 35, doi: 10.48550/arXiv.1907.04833
  * Riess et al. (2022) Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7, doi: 10.3847/2041-8213/ac5c5b
  * Rinaldi et al. (2023) Rinaldi, S., Del Pozzo, W., Mapelli, M., Medina, A. L., & Dent, T. 2023, arXiv e-prints, arXiv:2310.03074, doi: 10.48550/arXiv.2310.03074
  * Santoliquido et al. (2022) Santoliquido, F., Mapelli, M., Artale, M. C., & Boco, L. 2022, MNRAS, 516, 3297, doi: 10.1093/mnras/stac2384
  * Schutz (1986) Schutz, B. F. 1986, Nature, 323, 310, doi: 10.1038/323310a0
  * Soares-Santos et al. (2019) Soares-Santos, M., Palmese, A., Hartley, W., et al. 2019, The Astrophysical Journal, 876, L7, doi: 10.3847/2041-8213/ab14f1
  * Stachurski et al. (2023) Stachurski, F., Messenger, C., & Hendry, M. 2023, arXiv e-prints, arXiv:2310.13405.  https://arxiv.org/abs/2310.13405
  * Taylor et al. (2012) Taylor, S. R., Gair, J. R., & Mandel, I. 2012, Phys. Rev. D, 85, 023535, doi: 10.1103/PhysRevD.85.023535
  * Tiwari (2018) Tiwari, V. 2018, Classical and Quantum Gravity, 35, 145009, doi: 10.1088/1361-6382/aac89d
  * Turski et al. (2023) Turski, C., Bilicki, M., Dálya, G., Gray, R., & Ghosh, A. 2023, MNRAS, 526, 6224, doi: 10.1093/mnras/stad3110
  * Vijaykumar et al. (2023a) Vijaykumar, A., Fishbach, M., Adhikari, S., & Holz, D. E. 2023a, arXiv e-prints, doi: 10.48550/ARXIV.2312.03316
  * Vijaykumar et al. (2023b) Vijaykumar, A., Saketh, M. V. S., Kumar, S., Ajith, P., & Choudhury, T. R. 2023b, Phys. Rev. D, 108, 103017, doi: 10.1103/PhysRevD.108.103017
  * Vitale et al. (2022) Vitale, S., Gerosa, D., Farr, W. M., & Taylor, S. R. 2022, in Handbook of Gravitational Wave Astronomy (Springer), 45, doi: 10.1007/978-981-15-4702-7_45-1
  * Yu et al. (2020) Yu, J., Wang, Y., Zhao, W., & Lu, Y. 2020, MNRAS, 498, 1786, doi: 10.1093/mnras/staa2465
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, The Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298