# Joint cosmological and gravitational-wave population inference using dark
sirens and galaxy catalogues

Rachel Gray  Freija Beirnaert  Christos Karathanasis  Benoît Revenu  Cezary
Turski  Anson Chen  Tessa Baker  Sergio Vallejo  Antonio Enea Romano
Tathagata Ghosh  Archisman Ghosh  Konstantin Leyde  Simone Mastrogiovanni
Surhud More

###### Abstract

In the absence of numerous gravitational-wave detections with confirmed
electromagnetic counterparts, the “dark siren” method has emerged as a leading
technique of gravitational-wave cosmology. The method allows redshift
information of such events to be inferred statistically from a catalogue of
potential host galaxies. Due to selection effects, dark siren analyses
necessarily depend on the mass distribution of compact objects and the
evolution of their merger rate with redshift. Informative priors on these
quantities will impact the inferred posterior constraints on the Hubble
constant ($H_{0}$). It is thus crucial to vary these unknown distributions
during an $H_{0}$ inference. This was not possible in earlier analyses due to
the high computational cost, restricting them to either excluding galaxy
catalogue information, or fixing the gravitational-wave population mass
distribution and risking introducing bias to the $H_{0}$ measurement. This
paper introduces a significantly enhanced version of the Python package
gwcosmo, which allows joint estimation of cosmological and compact binary
population parameters. This thereby ensures the analysis is now robust to a
major source of potential bias. The gravitational-wave events from the Third
Gravitational-Wave Transient Catalogue are reanalysed with the GLADE+ galaxy
catalogue, and an updated, more reliable measurement of $H_{0}=$
$69^{+12}_{-7}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ is found (maximum a
posteriori probability and 68% highest density interval). This improved method
will enable cosmological analyses with future gravitational-wave detections to
make full use of the information available (both from galaxy catalogues and
the compact binary population itself), leading to promising new independent
bounds on the Hubble constant.

##  1 Introduction

The disagreement surrounding the value of the Hubble constant ($H_{0}$) is one
of the dominant unanswered questions in current cosmology. Current
measurements of $H_{0}$ are in tension with each other to a level above
$4\sigma$ (see, _e.g._[1, 2]), with local (low-redshift) measurements
typically preferring a higher value, and early-universe (high-redshift, model-
dependent) measurements preferring a lower one. The cause of this tension is
unknown, but potential explanations include unaccounted-for measurement errors
(which are more likely to impact the local measurements), or a flaw in the
current standard cosmological model, $\Lambda$-cold-dark-matter ($\Lambda$CDM)
(which would impact the early-universe measurements) [3, 4].

The relationship between redshift $z$, luminosity distance $d_{L}$, and
$H_{0}$ can be expressed as

$d_{L}=\dfrac{c(1+z)}{H_{0}}\int^{z}_{0}\dfrac{dz^{\prime}}{\sqrt{(1+z^{\prime})^{3}\Omega_{M}+\Omega_{\Lambda}}},$
(1.1)

for a flat universe. At low redshifts this simplifies to

$d_{L}\approx\dfrac{cz}{H_{0}}.$ (1.2)

Hence, for a local measurement of $H_{0}$ to be made, all that is required is
redshift and luminosity distance information. The amplitude of a gravitational
wave signal from a compact binary coalescence (CBC) directly encodes the
luminosity distance to the source, requiring no external calibration. This
property has led to these types of gravitational wave (GW) events being
commonly termed “standard sirens”. If information about the redshift of the
source can be supplied, the GW data can be used to produce a measurement of
$H_{0}$ which is independent of the cosmic distance ladder, and of all
electromagnetic (EM) measurements of $H_{0}$ to date. Herein lies significant
potential for the field of gravitational wave astronomy to resolve the long-
standing Hubble tension.

If a uniquely-identified EM counterpart can be associated to the GW event
(making it a “bright siren”), the redshift of the GW host galaxy can be
determined via spectroscopy, allowing for a direct measurement of $H_{0}$.
Thus far this has only been possible for the detection of the binary neutron
star GW170817 and its EM counterpart [5, 6]. The small localisation region of
the event (possible because the two LIGO detectors [7] and Virgo detector [8]
were online at the time) meant that multiple telescopes were able to search
the relevant area over a short time-frame, leading to the identification of
the kilonova associated with the GW event, and subsequently the identification
of its host galaxy. This led to the first ever measurement of $H_{0}$ from GWs
[9], a remarkably informative measurement from a single event.

Since then, the third observing run of Advanced LIGO, Virgo and KAGRA [10] has
taken place, with observations lasting for just under a year and increasing
the current catalogue of GW transients to $\sim 90$ [11]. Despite improvements
in the sensitivity of the LIGO and Virgo detectors [12, 13], no EM
counterparts have been uniquely associated with any of the subsequent events.
This distinct lack of EM counterparts has led to the blossoming of methods for
supplying redshift information to these “dark sirens”.

Current dark sirens methods can be broadly split into two further
categories:111Here we make a distinction between dark siren methods and the
cross-correlation method, in which GWs are cross-correlated with galaxy
surveys or other tracers in order to obtain cosmological constraints. For
details on the cross-correlation method see, _e.g._[14, 15, 16, 17, 18, 19,
20].

  * •

The galaxy catalogue method: here the redshift information is provided by a
catalogue of potential host galaxies within the sky localisation region of a
GW event (see [21, 22, 23, 24, 25, 26, 27, 28]). Each galaxy contributes to a
hypothetical measurement of $H_{0}$, such that the galaxy structure within a
GW event’s localisation volume is reflected in the $H_{0}$ posterior it
produces. How informative the individual events are will depend strongly on
their localisation volumes. By combining the contributions of many events, the
true value of $H_{0}$ will be measured as other values will statistically
average out.

  * •

The population method, also known “spectral sirens”: here features in the mass
distribution of the compact object population break the mass-redshift
degeneracy (see [29, 30, 31, 32, 33]). This method can be thought of as
analogous to the redshifting of emission lines in EM spectra: if the location
of a feature is known in the source frame, its measured position in the
detector frame – shifted by a factor of $(1+z)$ – informs us about the
redshift of the source.

These two varieties of dark sirens methods (as presented in [28, 32]) were
applied in the most recent LIGO-Virgo-KAGRA cosmology paper [34], making use
of the GW data from the third Gravitational-Wave Transient Catalogue (GWTC-3)
[11]. Until recently, the two dark sirens methods have always been carried out
separately due to the computational difficulty in incorporating galaxy
catalogue information (which tends to be highly discrete) into a sampling
method, and the computational cost of incorporating a reliable correction for
galaxy catalogue incompleteness. This means that the population method has not
benefited from the additional constraining power of galaxy redshifts.
Conversely, the galaxy catalogue method required fixing the GW population when
computing the selection function, which opens it up to potential bias if the
assumed model (with fixed hyper-parameters) is incorrect. This second point is
clearly more serious, as it could lead to a significant bias in the
measurement of $H_{0}$ if not corrected.

This paper presents a novel method for carrying out cosmological and GW
population parameter estimation jointly, such that the resulting posterior on
$H_{0}$ is free from bias due to GW population assumptions.222Of course this
assumes that the correct population models have been assumed. See, _e.g._[32].
This method combines the best of both dark siren techniques above: for the
large numbers of binary black holes at high redshifts the galaxy catalogues
are relatively uninformative (due to incompleteness), but these events offer
significant constraining power through their mass distribution. Near-by, well-
localised binary neutron stars, neutron star - black hole pairs and BBHs are
relatively lower in numbers and hence do not strongly probe the mass
distribution, but instead can access the highly informative nature of a galaxy
catalogue to provide redshift information. This updated method is implemented
within the Python package gwcosmo,333https://git.ligo.org/lscsoft/gwcosmo and
used to reanalyse the GWTC-3 data in combination with the GLADE+ galaxy
catalogue [35], leading to a final updated value of $H_{0}$ which utilises 47
GW events and incorporates redshift information from the black hole (BH) mass
distribution, a full-sky galaxy catalogue, and the one confirmed EM
counterpart to date.

This version of gwcosmo, which is about 1000 times faster than its
predecessor, is well suited to the challenges that lie ahead in the field of
GW cosmology, in terms of the large numbers of gravitational wave events which
are expected over the coming years, and the increasing depths of galaxy
surveys which means that the sheer quantity of data that will be available
will require efficient handling and analyses which can scale to large
datasets. Looking to the future, beyond the $H_{0}$ tension, this version of
gwcosmo can additionally be used to put constraints on models of modified
gravitational-wave propagation [36].

Section 2 presents the updated methodology in a Bayesian framework, including
introducing the crucial concept of a “ line-of-sight (LOS) redshift prior”.
Section 3 discusses the choice of GW data and galaxy catalogues, including how
the GLADE+ galaxy catalogue is converted into a LOS redshift prior. Section 4
reanalyses the GWTC-3 data with the updated methodology, under the same
assumptions as [34], and quantifies the impact of various assumptions, then
presents a full population + galaxy catalogue reanalysis of GWTC-3, including
an updated measurement of $H_{0}$ which utilises the greatest amount of GW
information to date. Finally, section 5 concludes the paper and discusses the
application of this method to future data.

Finally, a paper presenting an alternative method of joint cosmological and
population inference using the Icarogw package was released while this paper
was in preparation [37]. The methodology is derived using a different approach
(expressing the likelihood in terms of CBCs merger rates) which is consistent
with the gwcosmo methodology where equivalent assumptions are made. The way in
which gravitational-wave data is handled between the two codes differs, as
Icarogw computes numerical integrals involved in the evaluation of the
likelihood by summing over posterior samples, while gwcosmo constructs a
kernel density estimate (KDE) for different lines of sight. The Icarogw
approach enables a faster analysis, at the expense of requiring that pixel
sizes and galaxy redshift distributions are large enough to contain a
representative number of posterior samples. The gwcosmo approach, while
potentially slower, does not suffer from any of these resolution effects due
to the smoothing of the GW posterior.444No formal timing tests have yet been
carried out between the two codes, so a definitive comparison is left for
future work. Differences aside, we find good agreement with their results when
considering the same set of 42 BBHs for a full catalogue+population analysis,
and when considering the highly informative dark siren GW190814 for a
catalogue-only analysis.

##  2 Method

In general, the posterior on a set of hyper-parameters $\Lambda$, given a set
of GW data $\\{x_{\text{GW}}\\}$, corresponding to $N_{\text{det}}$
detections, $\\{D_{\text{GW}}\\}$, is given by [38, 39]:

$\displaystyle p(\Lambda|\\{x_{\text{GW}}\\},\\{D_{\text{GW}}\\})$
$\displaystyle\propto
p(\Lambda)p(N_{\text{det}}|\Lambda)\prod^{N_{\text{det}}}_{i}p(x_{\text{GW}i}|D_{\text{GW}i},\Lambda).$
(2.1)

The parameter $D_{\text{GW}i}$ is a binary parameter (taking the value 0 or 1)
denoting that the GW data of event $i$ passed some detection threshold
statistic in order to be deemed “detected”. The term $\Lambda$ includes the
cosmological parameters of interest (nominally $H_{0}$ for this analysis) as
well as a set of population hyper-parameters which describes the underlying
mass distribution of compact binaries and their distribution in redshift. This
expression expands to:

$\displaystyle p(\Lambda|\\{x_{\text{GW}}\\},\\{D_{\text{GW}}\\},I)$
$\displaystyle\propto
p(\Lambda|I)p(N_{\text{det}}|\Lambda,I)\prod^{N_{\text{det}}}_{i}\dfrac{\int
p(x_{\text{GW}i}|\theta,\Lambda,I)p(\theta|\Lambda,I)d\theta}{\int
p(D_{\text{GW}i}|\theta,\Lambda,I)p(\theta|\Lambda,I)d\theta},$ (2.2)
$\displaystyle\propto p(\Lambda|I)p(N_{\text{det}}|\Lambda,I)\left[\int
p(D_{\text{GW}}|\theta,\Lambda,I)p(\theta|\Lambda,I)d\theta\right]^{-N_{\text{det}}}$
$\displaystyle\hskip 100.0pt\times\prod^{N_{\text{det}}}_{i}\int
p(x_{\text{GW}i}|\theta,\Lambda,I)p(\theta|\Lambda,I)d\theta,$

where it’s been written explicitly that the individual GW likelihood need to
be marginalised over individual event parameters (such as sky location,
masses, inclination, spins, redshift etc.), grouped together and denoted by
$\theta$. GW selection effects are incorporated in the term
$p(D_{\text{GW}}|\Lambda)=\int
p(D_{\text{GW}}|\theta,\Lambda)p(\theta|\Lambda)d\theta$. This is the
“probability of detection” for a CBC drawn from a population described by
parameters $\Lambda$, and because it applies to the whole population of CBCs
(as opposed to being event-specific) it can be taken outside the
product.555The term $p(D_{\text{GW}i}|x_{\text{GW}i},\Lambda,I)$, which
appears when you first expand the expression using Bayes theorem, is equal to
1 as the GW data is, by definition, detected. The term $I$ on the right-hand
side of each expression is a catch-all term which holds any additional
assumptions that have not been explicitly expressed.

The method, as written to this point, is the same between the population and
the galaxy catalogue methods of computing $H_{0}$, but at this point it is
worth reiterating the manner in which these analyses will differ: the prior
redshift information. In the population method, the redshift prior would
commonly be taken to be uniform in comoving volume (in order to model the
galaxy distribution with the least informative assumptions when lacking any
galaxy positional data), modified by some additional term for the evolution of
the merger rate of binaries with redshift. In this scenario, the redshift
information necessary for cosmological inference comes mainly from the
distribution of the observed GW masses $p(m^{d}_{1},m^{d}_{2})$, in the
detector frame (denoted by the superscript “$d$”). Gravitational-wave
detectors are sensitive to detector-frame mass, rather than source-frame mass.
This means that the values observed in the detector frame have suffered the
effect of the cosmological expansion and are related to their value in the
source frame (labelled with the superscript “$s$”) by the relation
$m^{s}=m^{d}/(1+z)$. If the mass distribution of compact binaries contains
features this breaks the mass-redshift degeneracy, meaning that the observed
population of compact objects becomes cosmologically informative. If the
attributes of these features are known then they can be strongly informative,
with the caveat that assuming an incorrect mass distribution will have a
strong impact on the inferred cosmology (see _e.g._[34, 32]).

In order to avoid biasing the cosmological inference by assuming a mass
distribution, the uncertainty in the population can be marginalised over, by
carrying out multi-parameter estimation on cosmological and population hyper-
parameters together. Explicitly, the $p(\theta|\Lambda,I)$ terms in Eq. 2.2
can be expanded as
$p(\theta|\Lambda,I)=p(m^{s}_{1},m^{s}_{2}|\Lambda,I)p(\phi|\Lambda,I)$, where
$\theta=\\{m^{s}_{1},m^{s}_{2},\phi\\}$ and $\Lambda$ is the set of
cosmological and population hyper-parameters of interest, including those
which describe source-frame mass distribution prior. For the sake of
compactness, we do not make this expansion for the rest of the derivation, but
it is important to remember that redshift information from the mass
distribution is entering the analysis in this way.

In the galaxy catalogue case, instead of taking the redshift prior to be
uniform in comoving volume, it is built from a galaxy catalogue, and so
instead of being a smooth function of redshift, becomes something highly
discrete due to the point-like nature of galaxies.666Galaxies are, for our
purposes, point-like in right ascension and declination. In terms of redshift,
however, galaxies often come with a non-negligible redshift uncertainty. In
this scenario, information from the structures observed in the mass
distribution still contributes to the cosmological inference as before, though
it can be a subdominant effect. 777Information from the mass distribution is
likely to be subdominant when the GW events used in the analysis are well
localised in regions of the sky which have a high level of galaxy catalogue
completeness. Conversely, where the localisation volumes of GW events are very
large, or when the galaxy catalogue is very incomplete, the redshift
information from the galaxy catalogue can be subdominant to the information
coming from the overall population, and the resulting cosmological inference
will converge to the same result as a population-only analysis.  In addition
to galaxies being discrete, galaxy catalogues are incomplete at the redshifts
to which GWs are currently detectable, meaning that an incompleteness
correction must be applied to account for the fact that the host galaxy of the
GW event may not be contained within the catalogue at all, due to being
fainter than the flux limit of the galaxy survey. To further complicate
matters, the completeness of galaxy catalogues is not uniform across the sky
due to 1) obscuration by the Milky Way and 2) different surveys in different
patches of the sky, leading to varying depths of completeness as a function of
sky location, $\Omega$. The variation of incompleteness across the sky was
addressed extensively in [28], in which the sky was separated into
$N_{\text{pix}}$ equally-sized pixels and the completeness correction was
based on the catalogue’s apparent magnitude threshold within each pixel. Here,
a similar approach is applied, but by deriving the rest of the method in a
different order to the one presented in that paper, the analysis can be made
much more computationally tractable. With this in mind, Eq. 2.2 can be
rewritten to explicitly include the marginalisation over $\Omega$ as a sum
over pixels:

$\displaystyle p(\Lambda|\\{x_{\text{GW}}\\},\\{D_{\text{GW}}\\},I)\propto
p(\Lambda|I)p(N_{\text{det}}|\Lambda,I)$
$\displaystyle\left[\int\sum^{N_{\text{pix}}}_{j}p(D_{\text{GW}}|\Omega_{j},\theta,\Lambda,I)p(\theta|\Omega_{j},\Lambda,I)p(\Omega_{j}|I)d\theta\right]^{-N_{\text{det}}}$
(2.3)
$\displaystyle\times\prod^{N_{\text{det}}}_{i}\left[\int\sum^{N_{\text{pix}}}_{j}p(x_{\text{GW}i}|\Omega_{j},\theta,\Lambda,I)p(\theta|\Omega_{j},\Lambda,I)p(\Omega_{j}|I)d\theta\right],$
$\displaystyle\propto p(\Lambda|I)p(N_{\text{det}}|\Lambda,I)$
$\displaystyle\left[\int\sum^{N_{\text{pix}}}_{j}p(D_{\text{GW}}|\Omega_{j},\theta,\Lambda,I)p(\theta|\Omega_{j},\Lambda,I)d\theta\right]^{-N_{\text{det}}}$
$\displaystyle\times\prod^{N_{\text{det}}}_{i}\left[\int\sum^{N_{\text{pix}}}_{j}p(x_{\text{GW}i}|\Omega_{j},\theta,\Lambda,I)p(\theta|\Omega_{j},\Lambda,I)d\theta\right].$

The $p(\Omega_{j}|I)$ terms are recognised as being identical between same-
sized pixels, so come out the front and cancel in numerator and denominator.

Now it is worth explicitly writing the marginalisation over $z$ (with
$\theta=\\{z,\theta^{\prime}\\}$):

$\displaystyle p(\Lambda|\\{x_{\text{GW}}\\},\\{D_{\text{GW}}\\},I)$ (2.4)
$\displaystyle\propto p(\Lambda|I)p(N_{\text{det}}|\Lambda,I)$
$\displaystyle\left[\iint\sum^{N_{\text{pix}}}_{j}p(D_{\text{GW}}|\Omega_{j},z,\theta^{\prime},\Lambda,I)p(\theta^{\prime}|\Omega_{j},\Lambda,I)p(z|\Omega_{j},\Lambda,I)d\theta^{\prime}dz\right]^{-N_{\text{det}}}$
$\displaystyle\times\prod^{N_{\text{det}}}_{i}\left[\iint\sum^{N_{\text{pix}}}_{j}p(x_{\text{GW}i}|\Omega_{j},z,\theta^{\prime},\Lambda,I)p(\theta^{\prime}|\Omega_{j},\Lambda,I)p(z|\Omega_{j},\Lambda,I)d\theta^{\prime}dz\right]$
$\displaystyle\propto p(\Lambda|I)p(N_{\text{det}}|\Lambda,I)$
$\displaystyle\left[\iint
p(D_{\text{GW}}|z,\theta^{\prime},\Lambda,I)p(\theta^{\prime}|\Lambda,I)\sum^{N_{\text{pix}}}_{j}p(z|\Omega_{j},\Lambda,I)d\theta^{\prime}dz\right]^{-N_{\text{det}}}$
$\displaystyle\times\prod^{N_{\text{det}}}_{i}\left[\iint\sum^{N_{\text{pix}}}_{j}p(x_{\text{GW}i}|\Omega_{j},z,\theta^{\prime},\Lambda,I)p(\theta^{\prime}|\Lambda,I)p(z|\Omega_{j},\Lambda,I)d\theta^{\prime}dz\right].$

Here the assumption has been made that $p(D_{\text{GW}}|\Lambda,I)$ is uniform
across the sky (_i.e._ pixel-independent), and hence can be taken outside the
sum over pixels. As $p(D_{\text{GW}}|\Lambda,I)$ is averaged over an observing
run, the rotation of the Earth blurs out much of the right ascension (RA) and
declination (dec) dependence of the detectors’ antenna patterns, leaving only
a mild declination dependence that for the time being can be safely ignored.
The remaining parameters contained within $\theta$ are independent of $\Omega$
and so lose their dependence. If the prior on the GW merger rate, $R$, is set
to $\propto 1/R$, the term $p(N_{\text{det}}|\Lambda,I)$ loses its dependence
on $\Lambda$ and can be ignored (see [23]), a simplification which is used
throughout the rest of this paper.

###  2.1 The line of sight redshift prior

Now the LOS redshift prior, $p(z|\Omega_{j},\Lambda,I)$, can be considered
more closely. As in the original gwcosmo methodology, the redshift prior must
include the cases where the host galaxy is ($G$), and is not ($\bar{G}$),
inside the galaxy catalogue, which also requires marginalising over apparent
and absolute magnitude ($m$ and $M$) as these determine, to leading order,
which galaxies will end up inside a galaxy catalogue which is formed from a
flux-limited EM survey. An important distinction to make at this point is that
between a _galaxy_ and a _host galaxy_. So far in the derivation, the
assumption that a real GW signal has been emitted has been implicit (_i.e._
$x_{\text{GW}}$ is assumed to correspond to a real GW trigger, not a false
alarm). This can be made explicit by including the parameter $s$ (denoting the
presence of a GW source, which until now has been hidden inside $I$ on the
right-hand side of every expression).

$\displaystyle p(z|\Omega_{i},\Lambda,s,I)$
$\displaystyle=\iint\sum_{g=G,\bar{G}}p(z,M,m,g|\Omega_{i},\Lambda,s,I)\;dMdm,$
(2.5) $\displaystyle=p(G|\Omega_{i},\Lambda,s,I)\iint
p(z,M,m|G,\Omega_{i},\Lambda,s,I)\;dMdm$ $\displaystyle\hskip
100.0pt+p(\bar{G}|\Omega_{i},\Lambda,s,I)\iint
p(z,M,m|\bar{G},\Omega_{i},\Lambda,s,I)\;dMdm,$

As in [28], $p(z,M,m|G,\Omega_{i},\Lambda,s,I)$ is the prior on $z$, $M$ and
$m$ for host galaxies, informed by the galaxy catalogue, within the sky area
covered by pixel $i$. Its counterpart,
$p(z,M,m|\bar{G},\Omega_{i},\Lambda,s,I)$ is the prior on the same parameters
but outside the catalogue, and so will have some associated apparent magnitude
threshold ($m_{\text{th}}$).

####  2.1.1 The in-catalogue part

Taking a closer look at the in-catalogue part,
$p(z,M,m|G,\Omega_{i},\Lambda,s,I)$, galaxies come as (approximately) delta-
function-like points in $\Omega$ and $m$, with some larger uncertainty on $z$.
The absolute magnitude of the galaxy is not provided in the catalogue, and
requires a cosmological assumption to compute. For galaxy $j$,
$M_{j}=f(z_{j},m_{j},\Lambda)$. Applying Bayes’ theorem to the term
$p(z,M,m|G,\Omega_{i},\Lambda,s,I)$ gives:

$\displaystyle p(z,M,m|G,\Omega_{i},\Lambda,s,I)$
$\displaystyle=\dfrac{p(z,M,m|G,\Omega_{i},\Lambda,I)p(s|z,M,m,G,\Omega_{i},\Lambda,I)}{p(s|G,\Omega_{i},\Lambda,I)},$
(2.6)
$\displaystyle=\dfrac{1}{p(s|G,\Omega_{i},\Lambda,I)}p(M|z,m,G,\Omega_{i},\Lambda,I)p(z,m|G,\Omega_{i},\Lambda,I)p(s|z,M,\Lambda,I),$
$\displaystyle=\dfrac{1}{p(s|G,\Omega_{i},\Lambda,I)}\delta(M-M(z,m,\Lambda))p(z,m|G,\Omega_{i},I)p(s|z,M,\Lambda,I).$

The value for $M$ is fixed for given values of $z$, $m$ and a cosmological
model, and so turns into a delta function. The term $p(z,m|G,\Omega_{i},I)$ is
simply the prior on $z$ and $m$ given the galaxies in the catalogue within
pixel $i$. Note that this is now the prior on _all_ galaxies, not just host
galaxies, as the parameter $s$ has been separated out. Because this
information comes purely from the galaxy catalogue, the dependence on
$\Lambda$ is lost. The probability for a galaxy to host a GW merger retains
some dependence on $z$ and $M$ to allow for possible luminosity weighting, or
evolving merger rates with redshift (_i.e._ the $p(s|z,M,\Lambda,I)$ term,
which allows galaxies with certain properties to be considered as more likely
hosts for CBCs). Now integrating $p(z,M,m|G,\Omega_{i},\Lambda,s,I)$ over $m$
and $M$ leads to

$\displaystyle\iint p(z,M,m|G,\Omega_{i},\Lambda,s,I)\;dMdm$
$\displaystyle=\dfrac{1}{p(s|G,\Omega_{i},\Lambda,I)}\int
p(z,m|G,\Omega_{i},I)p(s|z,M(z,m,\Lambda),\Lambda,I)\;dm.$ (2.7)

The delta-function relationship between $M$, $z$ and $m$ means that the
integral over $M$ can be carried out by simply replacing all occurrences of
$M$ with $M(z,m,\Lambda)$.

The galaxy catalogue information comes into the term $p(z,m|G,\Omega_{i},I)$,
which can be replaced by a sum over galaxies. The measured apparent magnitudes
(in a given observation band) can be taken as delta functions,888In reality,
these galaxies will have some uncertainty on their measured apparent magnitude
which is expected to be small relative to the measured value, and we ignore
this as a simplifying assumption (as has previously been done in the
literature). An in-depth look into the impact of including this uncertainty is
left for future investigation. while the measured redshift of each galaxy will
come with some non-negligible uncertainty, particularly in the case of
photometric redshifts.

$\displaystyle p(z,m|G,\Omega_{i},I)$
$\displaystyle=\dfrac{1}{N_{\text{gal}}(\Omega_{i})}\sum^{N_{\text{gal}}(\Omega_{i})}_{k}p(z|\hat{z}_{k})\delta(m-\hat{m}_{k}).$
(2.8)

Here the hat notation is used for measured quantities. The term
$p(z|\hat{z}_{k})$ is the _posterior_ on the redshift of galaxy $k$, given the
_measured_ redshift $\hat{z}_{k}$. The rest of this paper will use the
assumption that

$\displaystyle p(z|\hat{z}_{k})=\mathcal{G}(z-\hat{z}_{k};\hat{\sigma}_{k}),$
(2.9)

_i.e._ that a Gaussian centered at $\hat{z}_{k}$ with standard deviation
$\hat{\sigma}_{k}$ is a reasonable approximation for the redshift posterior of
galaxy $k$.999While this implementation was chosen for ease, there is no
strict requirement that the galaxy uncertainties be Gaussian, or even that
they follow the same distribution. Hypothetically this method would very
easily extend to include the full PDF for photometric redshifts, allowing
their often non-Gaussian structure to be captured. Additionally the assumption
is made that the galaxy measurements provided in the catalogue used for the
analyses later in this paper (see Sec. 3.2) are posteriors.101010If the
measured redshift is actually a likelihood, a prior ought to be applied, in
which case a uniform in comoving volume prior might be sensible, to account
for the larger volume at higher redshifts. In the limit that
$\hat{\sigma}_{k}\rightarrow 0$, and $p(z|\hat{z}_{k})$ becomes a delta
function, the choice of prior will become irrelevant. This is yet another
reason to look forward to spectroscopic redshift measurements.

Returning to Eq. 2.7, and substituting the sum over galaxies into the integral
over $m$ returns the following:

$\displaystyle\iint p(z,M,m|G,\Omega_{i},\Lambda,s,I)\;dMdm$
$\displaystyle=\dfrac{1}{p(s|G,\Omega_{i},\Lambda,I)}\dfrac{1}{N_{\text{gal}}(\Omega_{i})}$
(2.10) $\displaystyle\hskip
30.0pt\times\int\sum^{N_{\text{gal}}(\Omega_{i})}_{k}p(z|\hat{z}_{k})\delta(m-\hat{m}_{k})p(s|z,M(z,m,\Lambda),\Lambda,I)\;dm,$
$\displaystyle=\dfrac{1}{p(s|G,\Omega_{i},\Lambda,I)}\dfrac{1}{N_{\text{gal}}(\Omega_{i})}\sum^{N_{\text{gal}}(\Omega_{i})}_{k}p(z|\hat{z}_{k})p(s|z,M(z,\hat{m}_{k},\Lambda),\Lambda,I).$

This in-catalogue part of the LOS redshift prior is essentially a weighted sum
of Gaussians, where the weights are some function of redshift and absolute
magnitude.

####  2.1.2 The out-of-catalogue part

Now turning attention to the out-of-catalogue part,
$p(z,M,m|\bar{G},\Omega_{i},\Lambda,I)$, and expanding it leads to a slightly
different expression:

$\displaystyle p(z,M,m|\bar{G},\Omega_{i},\Lambda,s,I)$
$\displaystyle=\dfrac{p(z,M,m|\bar{G},\Omega_{i},\Lambda,I)p(s|z,M,m,\bar{G},\Omega_{i},\Lambda,I)}{p(s|\bar{G},\Omega_{i},\Lambda,I)},$
(2.11)
$\displaystyle=\dfrac{1}{p(s|\bar{G},\Omega_{i},\Lambda,I)}p(z,M,m|\bar{G},\Omega_{i},\Lambda,I)p(s|z,M,\Lambda,I),$
$\displaystyle=\dfrac{1}{p(s|\bar{G},\Omega_{i},\Lambda,I)}\dfrac{p(\bar{G}|z,M,m,\Omega_{i},\Lambda,I)p(z,M,m|\Omega_{i},\Lambda,I)}{p(\bar{G}|\Omega_{i},\Lambda,I)}p(s|z,M,\Lambda,I).$

Here an additional step has been taken to separate out the term
$p(\bar{G}|z,M,m,\Omega_{i},\Lambda,I)$, which is the probability that a
galaxy with properties $z$, $M$, $m$, etc. is outside the galaxy catalogue,
which requires modeling the EM selection effects of that particular survey.
The main assumption here is that the galaxy catalogue is flux limited, such
that the probability that a galaxy is or is not contained within it depends on
the galaxy’s apparent magnitude, and whether it is greater or less than the
apparent magnitude threshold of the catalogue along the same line of sight,
$m_{\text{th}}(\Omega_{i})$. However, this is not the only reason galaxies may
be missing from a particular catalogue. If the redshift or colour information
of galaxies at high redshifts becomes particularly unreliable then it may be
prudent to discard galaxies above some redshift, $z_{\text{cut}}$, from the
catalogue, in which case their exclusion must be accounted for with the out of
catalogue term. In this case it becomes

$\displaystyle p(\bar{G}|z,M,m,\Omega_{i},H_{0},I)$
$\displaystyle=\Theta[m-m_{\text{th}}(\Omega_{i})]\Theta[z_{\text{cut}}-z]+\Theta[z-z_{\text{cut}}].$
(2.12)

Here, $\Theta$ denotes a Heaviside step function. If no redshift information
is discarded, _i.e._ $z_{\text{cut}}\rightarrow\infty$, then
$\Theta[z_{\text{cut}}-z]=1$ and $\Theta[z-z_{\text{cut}}]=0$, meaning that
the selection effects simplify to only depending on the $m_{\text{th}}$ of the
catalogue.

Another term to examine more closely is $p(z,M,m|\Omega_{i},\Lambda,I)$, which
is the prior on redshift, absolute and apparent magnitude for (all) galaxies
within pixel $i$. Because the parameter $\bar{G}$ is not present here, this
term does not depend on direction-dependent EM selection effects, and so
$\Omega_{i}$ can be dropped from the right-hand side. Due to the relationship
between $z$, $M$ and $m$, priors only need to be defined for two of these
parameters, as the third is specified by its relationship to the other two.
For the out of catalogue case it makes sense to define priors for $z$ and $M$
because these are reasonably well known (galaxies are broadly expected to
follow a uniform in comoving volume distribution in redshift when viewed on a
large enough scale, and galaxy absolute magnitudes are expected to follow a
Schechter function). This means that the term can be rewritten as
$p(z,M,m|\Lambda,I)=\delta(m-m(z,M,\Lambda))p(z,M|\Lambda,I)$. Substituting
this and Eq. 2.12 into Eq. 2.11 gives

$\displaystyle p(z,M,m|\bar{G},\Omega_{i},\Lambda,s,I)$
$\displaystyle=\dfrac{1}{p(s|\bar{G},\Omega_{i},\Lambda,I)p(\bar{G}|\Omega_{i},\Lambda,I)}\left(\Theta[m-m_{\text{th}}(\Omega_{i})]\Theta[z_{\text{cut}}-z]+\Theta[z-z_{\text{cut}}]\right)$
(2.13) $\displaystyle\hskip
100.0pt\times\delta(m-m(z,M,\Lambda))p(z,M|\Lambda,I)p(s|z,M,\Lambda,I).$

Now integrating this term over both $m$ and $M$ leads to

$\displaystyle\iint$ $\displaystyle
p(z,M,m|\bar{G},\Omega_{i},\Lambda,s,I)\;dMdm$ (2.14)
$\displaystyle=\dfrac{1}{p(s|\bar{G},\Omega_{i},\Lambda,I)p(\bar{G}|\Omega_{i},\Lambda,I)}\int\left(\Theta[m(z,M,\Lambda)-m_{\text{th}}(\Omega_{i})]\Theta[z_{\text{cut}}-z]+\Theta[z-z_{\text{cut}}]\right)$
$\displaystyle\hskip 200.0pt\times p(z,M|\Lambda,I)p(s|z,M,\Lambda,I)\;dM,$
$\displaystyle=\dfrac{1}{p(s|\bar{G},\Omega_{i},\Lambda,I)p(\bar{G}|\Omega_{i},\Lambda,I)}\left[\Theta[z_{\text{cut}}-z]\int^{M_{\text{max}}(H_{0})}_{M(z,m_{\text{th}}(\Omega_{i}),\Lambda)}p(z,M|\Lambda,I)p(s|z,M,\Lambda,I)\;dM\right.$
$\displaystyle\hskip
150.0pt+\left.\Theta[z-z_{\text{cut}}]\int^{M_{\text{max}}(H_{0})}_{M_{\text{min}}(H_{0})}p(z,M|\Lambda,I)p(s|z,M,\Lambda,I)\;dM\right].$

The Heaviside step function for the apparent magnitude threshold has been
converted to a redshift-dependent integration limit over $M$ which depends on
the apparent magnitude threshold of pixel $i$. The integration limits on $M$
in general depend on $H_{0}$ because the parameters of the Schechter function
are $H_{0}$-dependent, but the final distribution is insensitive to its value.

####  2.1.3 The full expression for the LOS prior

Both Eq. 2.10 and Eq. 2.14 contain prefactors which are of a format hard to
compute. They take the form $1/p(s|G,\Omega_{i},\Lambda,I)$ for the former,
and $1/(p(s|\bar{G},\Omega_{i},\Lambda,I)p(\bar{G}|\Omega_{i},\Lambda,I))$ for
the latter. However, when substituting Eq. 2.10 and Eq. 2.14 into Eq. 2.5 it
can be seen that these expressions simplify. For ease of notation let the rest
of Eq. 2.10 and Eq. 2.14 be termed $f(\text{IN CAT})$ and $f(\text{OUT OF
CAT})$ respectively. Then the expression becomes

$\displaystyle p(z|\Omega_{i},\Lambda,s,I)$
$\displaystyle=\dfrac{p(G|\Omega_{i},\Lambda,s,I)}{p(s|G,\Omega_{i},\Lambda,I)}f(\text{IN
CAT})+\dfrac{p(\bar{G}|\Omega_{i},\Lambda,s,I)}{p(s|\bar{G},\Omega_{i},\Lambda,I)p(\bar{G}|\Omega_{i},\Lambda,I)}f(\text{OUT
OF CAT}),$ (2.15)
$\displaystyle=\dfrac{p(s|G,\Omega_{i},\Lambda,I)p(G|\Omega_{i},\Lambda,I)}{p(s|G,\Omega_{i},\Lambda,I)p(s|\Omega_{i},\Lambda,I)}f(\text{IN
CAT})$ $\displaystyle\hskip
100.0pt+\dfrac{p(s|\bar{G},\Omega_{i},\Lambda,I)p(\bar{G}|\Omega_{i},\Lambda,I)}{p(s|\Omega_{i},\Lambda,I)p(s|\bar{G},\Omega_{i},\Lambda,I)p(\bar{G}|\Omega_{i},\Lambda,I)}f(\text{OUT
OF CAT}),$
$\displaystyle=\dfrac{1}{p(s|\Omega_{i},\Lambda,I)}\left[p(G|\Omega_{i},\Lambda,I)f(\text{IN
CAT})+f(\text{OUT OF CAT})\right].$

There is a lack of symmetry between the in and out of catalogue terms, as one
might have expected a $p(\bar{G}|\Omega_{i},H_{0},I)$ as a prefactor to the
out of catalogue term, but this term has cancelled with one in the
denominator. With further consideration, this makes sense – the sum over
galaxies of the in-catalogue part by definition integrates to 1 over $z$
(assuming uniform galaxy weights for now), while the integral over the out-of-
catalogue part does not, due to the apparent magnitude threshold appearing in
the integral limits.

The term $1/p(s|\Omega_{i},\Lambda,I)$ can be ignored as a normalisation
constant. It will be the same for all pixels, assuming GW sources are
isotropically distributed, and as such will eventually cancel with an
identical term in denominator once the LOS prior is substituted back into Eq.
2.4.

The one term we have not yet looked at is the probability of a galaxy in pixel
$i$ being inside the galaxy catalogue (notice now that there is no $s$ on the
right-hand side, so here we are considering _all_ galaxies, not just host
ones):

$\displaystyle p(G|\Omega_{i},\Lambda,I)$ $\displaystyle=\iiint
p(G|z,M,m,\Omega_{i},\Lambda,I)p(z,M,m|\Omega_{i},\Lambda,I)\;dzdMdm.$ (2.16)

The probability that a galaxy is inside the catalogue can then be expanded as

$\displaystyle p(G|z,M,m,\Omega_{i},\Lambda,I)$
$\displaystyle=\Theta[m_{\text{th}}(\Omega_{i})-m]\Theta[z_{\text{cut}}-z],$
(2.17)

which leads to

$\displaystyle p(G|\Omega_{i},\Lambda,I)$
$\displaystyle=\int^{z_{\text{cut}}}_{0}\int^{M(z,m_{\text{th}}(\Omega_{i}),\Lambda)}_{M_{\text{min}}(H_{0})}p(z,M|\Lambda,I)\;dMdz.$
(2.18)

The fact that all probability densities must be properly normalised may cause
the reader to wonder whether the choice of redshift range over which this is
done so has an impact on this analysis. After all the term
$p(G|\Omega_{i},\Lambda,I)$ will take a very different value if
$p(z,M|\Lambda,I)$ is normalised over a redshift range which is truncated at
$z_{\text{max}}=2$ and one where $z_{\text{max}}=10$. However, it is useful to
note that as long as the same expression for $p(z,M|\Lambda,I)$ is also be
used in Eq. 2.14, the normalisation will come out the front as a constant. As
such, the choice of what redshift range over which to normalise
$p(z,M|\Lambda,I)$ has no bearing on the result.

Combining Eqs. 2.15, 2.10 and 2.14 leads to a full expression for the redshift
prior of

$\displaystyle
p(z|\Omega_{i},\Lambda,s,I)=\dfrac{1}{p(s|\Omega_{i},\Lambda,I)}\Bigg{[}$
$\displaystyle
p(G|\Omega_{i},\Lambda,I)\dfrac{1}{N_{\text{gal}}(\Omega_{i})}\sum^{N_{\text{gal}}(\Omega_{i})}_{k}p(z|\hat{z}_{k})p(s|z,M(z,\hat{m}_{k},\Lambda),\Lambda,I)$
(2.19)
$\displaystyle+\bigg{(}\Theta[z_{\text{cut}}-z]\int^{M_{\text{max}}(H_{0})}_{M(z,m_{\text{th}}(\Omega_{i}),\Lambda)}p(z,M|\Lambda,I)p(s|z,M,\Lambda,I)\;dM$
$\displaystyle\hskip
10.0pt+\Theta[z-z_{\text{cut}}]\int^{M_{\text{max}}(H_{0})}_{M_{\text{min}}(H_{0})}p(z,M|\Lambda,I)p(s|z,M,\Lambda,I)\;dM\bigg{)}\Bigg{]}.$

Looking ahead, Fig. 8 shows an example of the LOS redshift prior,
demonstrating the discrete galaxies which are visible at low redshift, and how
it converges to a uniform in comoving distribution at higher redshifts as
galaxy catalogue completeness decreases.

####  2.1.4 The separability of GW mass and rate evolution hyper-parameters
from the line-of-sight redshift prior

The LOS redshift prior, as presented in Eq. 2.19, retains clear dependence on
$\Lambda$ throughout, where $\Lambda$ is a set of cosmological and population
hyper-parameters of interest, many of which we would like to jointly estimate.
However, if Eq. 2.19 has to be recomputed for every newly drawn set of
parameter values this will slow things down to the point of impracticality. To
address this, $\Lambda$ can be broken into its constituent parts, namely
$\Lambda=\\{\Lambda_{\text{cosmo}},\Lambda_{\text{mass}},\Lambda_{\text{rate}}\\}$.
$\Lambda_{\text{cosmo}}$ denotes the parameters of the cosmological model,
namely $\\{H_{0},\Omega_{M},\Omega_{\Lambda}\\}$. $\Lambda_{\text{mass}}$
denotes the parameters of the GW mass model, such as minimum and maximum mass,
slope, the position of any features, etc. $\Lambda_{\text{rate}}$ denotes the
parameters of the rate-evolution model which determines how the intrinsic rate
of CBC mergers evolves with redshift (see Sec. 4 for details).

The immediate impact of substituting this definition into Eq. 2.19 is that the
priors on galaxy redshift and absolute magnitude depend only on
$\Lambda_{\text{cosmo}}$. Beyond this, some simplifying assumptions must be
made. We assume that the GW mass model doesn’t evolve with redshift (there is,
as yet, no strong evidence for a mass model which evolves with redshift [40],
though this may change once the data from future observing runs has been
analysed [41, 42]). In this scenario, none of the terms in Eq. 2.19 are
dependent on $\Lambda_{\text{mass}}$. The host galaxy weights term retains
dependence on $\Lambda_{\text{rate}}$ (and, in the case where those galaxies
come from the catalogue, such that their absolute magnitude is computed from
their apparent magnitude and redshift, $\Lambda_{\text{cosmo}}$).

So far in Eq. 2.19 the only assumption about the relationship between GW
mergers and host galaxies is that it will have some (unspecified) dependence
on redshift and absolute magnitude. At this point in time, the true
relationship between galaxy properties and GW merger probability is not well
constrained, but there are a couple of reasons for parameterising the problem
this way. The usual justification is to reason that there is likely a causal
link between galaxy mass and/or star formation rate (SFR) and that galaxy’s
probability of hosting a GW merger. Higher mass means more matter available to
be turned into compact binaries. Higher SFR means more stars forming, which
means more binary systems forming and eventually turning into compact
binaries. Of course there may well be some significant time-delay between SFR
and a binary merging, but this is also currently uncertain. In any case, both
mass and SFR are traced by galaxy luminosity in different bands, and so
estimates for both can be constructed from observable data.111111Of course
other things such as galaxy metallicity and galaxy type may all also play a
role in host galaxy probability, but a full exploration of these parameters
and how they could be incorporated into the current methodology is left for
the future.

As we do not yet know the exact relationship between GW mergers and host
galaxy properties, it is worth remaining flexible and allowing that there may
be some additional redshift-dependence to the merger rate model. As such, for
the remainder of this paper, the following form is taken:

$\displaystyle p(s|z,M,\Lambda_{\text{rate}},I)$ $\displaystyle\propto
p(s|z,\Lambda_{\text{rate}},I)p(s|M,I).$ (2.20)

where

$p(s|M,I)\propto\begin{cases}L(M)&\text{if using luminosity weighting}\\\
\text{const}&\text{if assuming uniform weighting.}\end{cases}$ (2.21)

By keeping the parameter(s) $\Lambda_{\text{rate}}$ flexible, such that
$p(s|z,\Lambda_{\text{rate}},I)$ can return to a uniform distribution, this
allows additional unknown redshift dependence to be captured.

Making this substitution in Eq. 2.19 leads to

$\displaystyle p(z|\Omega_{i},\Lambda,s,I)=$
$\displaystyle\dfrac{p(s|z,\Lambda_{\text{rate}},I)}{p(s|\Omega_{i},\Lambda,I)}\Bigg{[}p(G|\Omega_{i},\Lambda,I)\dfrac{1}{N_{\text{gal}}(\Omega_{i})}\sum^{N_{\text{gal}}(\Omega_{i})}_{k}p(z|\hat{z}_{k})p(s|M(z,\hat{m}_{k},\Lambda_{\text{cosmo}}),I)$
(2.22)
$\displaystyle+p(z|\Lambda_{\text{cosmo}},I)\bigg{(}\Theta[z_{\text{cut}}-z]\int^{M_{\text{max}}(H_{0})}_{M(z,m_{\text{th}}(\Omega_{i}),\Lambda_{\text{cosmo}})}p(M|z,\Lambda_{\text{cosmo}},I)p(s|M,I)\;dM$
$\displaystyle\hskip
100.0pt+\Theta[z-z_{\text{cut}}]\int^{M_{\text{max}}(H_{0})}_{M_{\text{min}}(H_{0})}p(M|z,\Lambda_{\text{cosmo}},I)p(s|M,I)\;dM\bigg{)}\Bigg{]}.$

The term $p(s|z,\Lambda_{\text{rate}},I)$ has come out the front, meaning that
$\Lambda_{\text{rate}}$ parameters are not needed for computing the bulk of
the LOS prior, and the aim becomes to simply precompute the part of the
expression within the square brackets.

Within the square brackets there is still a clear dependence on
$\Lambda_{\text{cosmo}}$ in both the in- and out-of-catalogue parts. The term
$p(z|\Lambda_{\text{cosmo}},I)$ is a uniform in comoving volume distribution
and its dependence on $H_{0}$ drops out when normalised, leaving only the
dependence on $\\{\Omega_{M},\Omega_{\Lambda}\\}$. The galaxy luminosity
function, $p(M|z,\Lambda_{\text{cosmo}},I)$, will depend only on $H_{0}$
(Schechter function parameters are quoted with a dependence on $h\equiv
H_{0}/100$). The relation between redshift, apparent and absolute magnitude,
$M(z,m,\Lambda_{\text{cosmo}})$ will retain dependence on all cosmological
parameters due to the necessity of converting redshift to luminosity distance.

Despite the dependence of individual terms on $H_{0}$, when computed for
different values of $H_{0}$ the LOS prior is independent of its value down to
some normalisation constant. Because the same prior is used to evaluate the GW
likelihood and GW selection effects, this dependence cancels, meaning that the
LOS prior can in fact be computed just once and used in an analysis in which
$\\{H_{0},\Lambda_{\text{mass}},\Lambda_{\text{rate}}\\}$ vary.

####  2.1.5 The resolution of the LOS prior

Looking at Eq. 2.22 it should be noted that pixel area does not explicitly
appear anywhere within the current definition, but there is an implicit
assumption present which is worth taking care with. The quantities computed in
the method thus far are based on assumptions about the distribution of
galaxies which we cannot see (_i.e._ that they are distributed, on average,
uniformly in comoving volume, and follow some kind of luminosity distribution
such as a Schechter function). Consider $p(G|\Omega_{i},\Lambda,I)$ for
example, and how it is expressed in Eq. 2.18. This term is the fraction of
galaxies which within pixel $i$ which are inside the galaxy catalogue. The
“true” value is given by
$N_{\text{gal}}(\Omega_{i})/N_{\text{tot}}(\Omega_{i})$ – that is the ratio of
the number of galaxies in the catalogue (in that pixel) to the total number of
galaxies (including the ones we cannot see) which lie in that pixel. Because
we do not have access to the true total number of galaxies we must use an
approximation based on our chosen priors. This approximation is okay, as long
as the pixel size is large enough that $N_{\text{gal}}(\Omega_{i})\approx
p(G|\Omega_{i},\Lambda,I)\times N_{\text{tot}}\times\Delta\Omega$, where
$N_{\text{tot}}$ is the total number of galaxies in the universe, and
$\Delta\Omega$ is the fractional area covered by one pixel. This will clearly
break down in the limit of infinitesimally small pixels which may not contain
any galaxies, but which would still have a non-zero value of
$p(G|\Omega_{i},\Lambda,I)$. In order to avoid this limit, pixel size must be
chosen to be large enough that the approximation is valid.121212Alternatively,
if pixels with the same level of completeness could be grouped together, such
that $N_{\text{gal}}(\Omega_{i})$ can be estimated over a larger patch of the
sky, this issue with small number statistics may be similarly avoided. To do
so is beyond the scope of this paper and is left for future work.

This fine-resolution limit to the method has a drawback, in that all galaxies
within a pixel are treated as having the same $\Omega$, effectively reducing
the galaxies’ positional (RA and dec) information. If the GW likelihood
changes significantly across the area of one pixel, valuable information will
be lost. In order to carry out the analysis with suitably high resolution, we
make the choice to define two different levels of resolution:
$n_{\text{map}}$, which determines the resolution on which
$N_{\text{gal}}(\Omega_{i})$ and $m_{\text{th}}$ are estimated,131313Here we
continue to rely on estimating $m_{\text{th}}(\Omega_{i})$ using the median
$m$ of the galaxies inside pixel $i$, an assumption which also breaks down in
the high-resolution limit. and $n_{\text{high}}$, a finer resolution which
pixels are subdivided into. The specific choices for these resolutions are
discussed in section 3.2. These parameters, $n_{\text{map}}$ and
$n_{\text{high}}$ are values of nside, a parameter which relates to the number
of equal-sized pixels the surface of a sphere is separated into, using the
formula $n_{\text{pix}}=12\times\text{nside}^{2}$.

The $n_{\text{map}}$ resolution is used to compute HEALPix maps [43, 44] of
two quantities:141414HEALPix is available at http://healpix.sf.net.
$m_{\text{th}}$ and $N_{\text{gal}}(\Omega_{i})$. When computing the full LOS
redshift prior on the $n_{\text{high}}$ resolution, $m_{\text{th}}$ is taken
from the $m_{\text{th}}$ map, and used for the incompleteness corrections.
$N_{\text{gal}}(\Omega_{i})$ is also taken from the coarser map, and divided
through by the number of sub-pixels that the coarse pixel is divided into, and
this value is then used in place of $N_{\text{gal}}(\Omega_{i})$ in Eq. 2.22.
For clarity, the sum is still carried out over the true number of galaxies in
the fine pixel, but the normalising factor of $1/N_{\text{gal}}(\Omega_{i})$
is constructed from the coarse map. This means that you will get the same
results when computing the full LOS prior on a coarse map as when computing it
on a higher resolution, then summing pixels to return to the same resolution
as the coarse map. As such it is robust in the limit of infinitely high
resolutions.

Figure 1: Schematic demonstrating the impacts of pixel resolution. The same
patch of sky containing four galaxies is shown with 3 different pixel
resolutions.

Figure 1 may help to visualise this. First assume you do not pre-compute
$m_{\text{th}}$ or $N_{\text{gal}}(\Omega_{i})$ on a coarser resolution.
Looking at the left-hand panel, if the LOS prior was constructed each galaxy
would have a normalising factor of $1/4$ in front of it. In the middle panel,
however, the pixel has been split into 4 sub-pixels, and the same galaxies
have normalisation factors of 1 or $1/2$, depending on which pixel they fall
in. If you then summed the 4 sub-pixels you would not recover the same result
as the left-hand panel; instead two of your galaxies would be down-weighted
relative to the other two, simply due to small number statistics.
Additionally, in the high-resolution limit (right-hand panel) there will be a
maximum of 1 galaxy per pixel, but there can be an arbitrarily high number of
empty pixels, which means the out-of-catalogue contribution will dominate. If,
however, $N_{\text{gal}}(\Omega_{i})$ is computed for a coarse pixel then
divided through by the number of sub-pixels (leading to an average of 4
galaxies per sub-pixel in the left-hand panel of Fig. 1, 1 for the middle
panel, and 1/4 for the right-hand panel), then the effective
$N_{\text{gal}}(\Omega_{i})$ of each high resolution pixel remains inversely
proportional to the total number of out-of-catalogue contributions (1 per sub-
pixel).

###  2.2 Gravitational-wave selection effects

In order to obtain an unbiased estimate of cosmological (and population)
parameters, it is crucial to accurately model GW selection effects. The
“probability of detection” of a GW event (often shortened to $P_{\text{det}}$)
can be written as $p(D_{\text{GW}}|\Lambda,I)=\int
p(D_{\text{GW}}|\theta,\Lambda,I)p(\theta|\Lambda,I)\mathrm{d}\theta$. In
previous versions of gwcosmo, where the only parameter being constrained was
$H_{0}$ (_i.e._ $\Lambda=\Lambda_{\text{cosmo}}=H_{0}$), and all other mass
model and cosmological parameters were set to fixed values, this was done by
simulating a population of GW events and computing their probability of
detection as a function of redshift and $H_{0}$ (marginalising over masses,
sky position, polarisation, and inclination). Now, with a dozen or more
parameters to jointly constrain, this approach is no longer practical due to
the high computational burden.

In order for Eq. 2.4 to be used by a Markov Chain Monte Carlo (MCMC) or Nested
Sampling method, the probability of detection must be quickly computable for
any values of
$\Lambda=\\{\Lambda_{\text{cosmo}},\Lambda_{\text{mass}},\Lambda_{\text{rate}}\\}$.
Probability of detection is equivalent to the ratio of detectable events
$N_{\text{exp}}(\Lambda)$ over the total number of mergers $N(\Lambda)$ up to
a certain redshift. This integral can be computed numerically using a Monte
Carlo approximation.

$\displaystyle\frac{N_{\text{exp}}}{N}(\Lambda)$ $\displaystyle=\iiint
p(D_{\text{GW}}|m_{1}^{s},m_{2}^{s},z,\Lambda,I)p(m_{1}^{s},m_{2}^{s},z|\Lambda,I)\mathrm{d}m_{1}^{s}\mathrm{d}m_{2}^{s}\mathrm{d}z,$
(2.23)
$\displaystyle\approx\frac{1}{N_{\text{sim}}}\sum_{i=1}^{N_{\text{det}}}\frac{p(m^{s}_{1,i},m^{s}_{2,i},z_{i}|\Lambda,I)}{\pi_{\text{inj}}^{s}(m^{s}_{1,i},m^{s}_{2,i},z_{i}|\Lambda,I)},$

where $\pi_{\text{inj}}^{s}$ is used to denote a (prior) probability density
function which is distinct from another prior on the same parameters (see
later in this section). The discrete summation is done over a large set of
simulated GW signals (typically of the order of
$N_{\text{sim}}=10^{5}-10^{6}$) using a prior on the characteristics of the
simulated CBC (masses, redshift or luminosity distance). These signals are
injected in realistic time sequences of noise following the detector
sensitivities and duty cycles in order to reproduce observational conditions
as closely as possible. Then, we can compute the corresponding signal-to-noise
ratio (SNR) for each simulated signal so that we have in the end a number
$N_{\text{det}}$ of detected events (injections) above the specified network
SNR threshold in the analysis, which, when summed produces the probability of
detection.

The probabilities appearing in the numerator of Eq. 2.23,
$p(m^{s}_{1,i},m^{s}_{2,i},z_{i}|\Lambda,I)$, are computed for each value of
$\Lambda$ drawn during the sampling process, applying the desired priors on
compact binary source-frame masses and redshift (as an example, the prior for
the distribution of the primary mass $m_{1}$ in the source frame can be seen
in Eq. 4.1). In particular, the model used for $p(z|\Lambda,I)$ is the LOS
redshift prior of Eq. 2.22. The probabilities of the denominator,
$\pi_{\text{inj}}^{s}(m^{s}_{1,i},m^{s}_{2,i},z_{i}|\Lambda)$, are the
individual probabilities of the detected injections after taking into account
the current value of $\Lambda$ and must be recomputed accordingly. These
probabilities are therefore closely related to their known values, when the
injections have been drawn from the chosen prior distribution
$\pi_{\text{inj}}$. In the following analyses, we chose to draw injections in
the detector frame,
$\pi_{\text{inj}}(m^{d}_{1,i},m^{d}_{2,i},d_{\mathrm{L},i})$, as detectability
of a GW event is a function of these parameters, and they therefore provide
coverage for a wide range of cosmological values. With this choice, it is
possible to obtain the required values in terms of source frame masses,
$\pi_{\text{inj}}^{s}(m^{s}_{1,i},m^{s}_{2,i},z_{i}|\Lambda)$, using the
relation:

$\pi_{\text{inj}}^{s}(m^{s}_{1,i},m^{s}_{2,i},z_{i}|\Lambda)=\pi_{\text{inj}}(m^{d}_{1,i},m^{d}_{2,i},d_{\mathrm{L},i})(1+z_{i})^{2}\left|\frac{\partial
d_{L}}{\partial z}\right|_{z=z_{i}}$ (2.24)

where the Jacobians for transforming between source frame and detector frame,
and $d_{L}$ and $z$ have been written explicitly. Finally, substituting Eq.
2.24 into the denominator of Eq. 2.23 gives:

$\frac{N_{\text{exp}}}{N}(\Lambda)\approx\frac{1}{N_{\text{sim}}}\sum_{i=1}^{N_{\text{det}}}\frac{p(m^{s}_{1,i},m^{s}_{2,i},z_{i}|\Lambda)}{\pi_{\text{inj}}(m^{d}_{1,i},m^{d}_{2,i},d_{\mathrm{L},i})(1+z_{i})^{2}\left|\frac{\partial
d_{L}}{\partial z}\right|_{z=z_{i}}}.$ (2.25)

The denominator acts as a weight of each detected injection. The sampling
process modifies these weights for each new value of $\Lambda$; this is a fast
operation and we are therefore able to get, at every point of the hyper-
parameters space, the estimate of $N_{\text{exp}}/N$ for any value of
$\Lambda$. Assuming that the redshift prior in the numerator of Eq. 2.25 is
constructed from the LOS redshift prior, this is (to within a normalisation
constant) the quantity between the first set of square brackets appearing in
Eq. 2.4:

$\left[\iint
p(D_{\text{GW}}|z,\theta^{\prime},\Lambda,I)p(\theta^{\prime}|\Lambda,I)\sum^{N_{\text{pix}}}_{j}p(z|\Omega_{j},\Lambda,I)d\theta^{\prime}dz\right].$

##  3 Gravitational-wave data and galaxy catalogues

###  3.1 Gravitational-wave data

In this paper we select all GW events with SNR $>11$ from GWTC-3 (as was done
in [34]). This leads to a selection of 42 BBHs, 2 NSBHs (GW200105 and
GW200115), the asymmetric mass binary GW190814 (which is treated as an NSBH in
all the analyses which follow) and as well as two BNSs (GW170817 and
GW190425). gwcosmo makes use of both GW posterior samples (for computing a KDE
of the GW distance distribution along every line-of-sight) and skymaps (for
weighting the normalised distance KDEs to retain the correct 3-dimensional
distribution of GW probability in RA, dec and luminosity distance).151515The
data, in the form of posterior samples and skymaps, for O1, O2 and O3a [45]
can be found at https://zenodo.org/record/6513631 , and for O3b [11] at
https://zenodo.org/record/5546663.

A set of injections (as described in section 2.2) was generated using gwcosmo
to characterise the detectability of GW events assuming a network SNR
threshold of 11.

###  3.2 Galaxy catalogues

For this analysis we use the GLADE+ galaxy catalogue, which is a composite
catalogue made by cross-matching the Gravitational Wave Galaxy Catalogue [46],
HyperLEDA [47], the 2 Micron All-Sky Survey Extended Source Catalog [48, 49],
the 2MASS Photometric Redshift Catalog [50], the WISExSCOS Photometric
Redshift Catalogue [51], and the Sloan Digital Sky Survey quasar catalogue
from the 16th data release [52]. We choose galaxies in the $K$-band in order
for our results to be comparable to [34]. We assume that the galaxy luminosity
function follows Schechter function that does not evolve with
redshift,161616This is a reasonable assumption at low redshifts, but will not
remain true at higher redshifts ($z\sim 1$). The flexibility in the model
which evolves GW merger rate with redshift in gwcosmo will compensate for some
mis-match between the chosen Schechter function here and the truth, but this
is a point which it will be worth investigating thoroughly in the future. with
parameters $M^{*}=-23.39+5\log(h)$ (where $h\equiv H_{0}/100$) and
$\alpha=-1.09$, following [53]. To curtail the faint-end of the Schechter
function we choose a limit of $M_{\text{max}}=-19.0+5\log(h)$. The default
assumption for host galaxy weighting is that host probability is proportional
to $K$-band luminosity, which traces galaxy stellar mass. As mentioned in
section 2.1.4, an assumption of CBC merger rate evolution with redshift is not
required at this stage as it is incorporated at a later stage in the analysis.
$K$-corrections are also applied when converting galaxy apparent magnitudes to
rest-frame absolute magnitudes, following the prescription in [53].

With the galaxy catalogue selected, and the observation band chosen, the LOS
redshift prior can now be pre-computed. The first stage of this is to pre-
compute maps describing the $m_{\text{th}}$ of the catalogue across the sky,
as well as the corresponding numbers of galaxies. In order to assess which
coarse resolution satisfies the conditions for $n_{\text{map}}$ listed in Sec.
2.1.5, we have computed maps for two different $n_{\rm map}$ values: 32 and
64. These maps are shown in Figs. 2 and 3. For the $m_{\text{th}}$ maps, the
minimum number of galaxies a pixel must contain in order to not be considered
empty is 10. For $n_{\rm map}$=32 this results in a map for which the only
empty pixels are those for which the Milky Way band obscures the sky. For
$n_{\rm map}$=64, due to the limited number of galaxies with available
$K$-band luminosities and the small pixel sizes, many of the pixels evaluate
as empty. For this reason, $n_{\text{map}}=32$ is chosen to be the fiducial
value (though the impact of changing to $n_{\text{map}}=64$ is investigated in
Sec. 4.2), and we do not investigate coarser resolutions, as they will not be
able to accurately capture the varying galaxy catalog completeness across the
sky.

(a) $m_{th}$ map for $n_{\rm map}$ = 32. The sky is divided into 12,288
pixels, each covering an area of 3.36 $\text{deg}^{2}$. (b) $m_{th}$ map for
$n_{\rm map}$ = 64. The sky is divided into 49,152 pixels, each covering an
area of 0.84 $\text{deg}^{2}$.

Figure 2: $m_{\text{th}}$ maps for different $n_{\rm map}$ values. The minimum
number of galaxies a pixel must contain in order to not be considered empty is
10. Empty pixels are shown in white.

(a) $m_{th}$ map for $n_{\rm map}$ = 32. (b) $m_{th}$ map for $n_{\rm map}$ =
64.

Figure 3: Effective $N_{\text{gal}}$ maps for different $n_{\rm map}$ values.

With the $m_{\text{th}}$ and $N_{\text{gal}}$ maps precomputed, the full LOS
redshift prior can now be computed at a higher resolution. The LOS redshift
prior (the part within square brackets of Eq. 2.22) is computed for every
pixel assuming $n_{\rm high}$ = 32, 64 and 128. The impact of changing the
choice of resolution is investigated in section 4.2, and $n_{\text{high}}=128$
is chosen to be the fiducial value for the analyses moving forwards, as it
allows the maximum information from galaxy RAs and decs to be utilised. Higher
resolutions are not investigated in this paper due to 1) the length of time it
takes to compute the LOS redshift prior (though this is parallelisable) and 2)
the resulting file sizes: a LOS prior with $n_{\text{high}}=128$ is about 7GB,
and each subsequent increase in resolution divides each pixel into 4, and so
quadruples the resulting file size. Usefully, however, this file size is
completely independent of the number of galaxies which were originally in the
galaxy catalogue.171717Additionally, even though the file size is large,
because only the pixels which cover the _e.g._ 99.9% sky area of the GW event
are required, the full gwcosmo analysis does not require holding all 7GB in
memory – the data can be swiftly reduced to a manageable quantity. The second
criterion which determines the overall file size, other than pixel resolution,
is the redshift resolution.

Every pixel in the LOS redshift prior is computed for the same redshift array,
in order to facilitate fast operations summing pixels at later stages. For
this reason, the redshift array resolution must be high enough to capture the
redshift profile of every galaxy in the catalogue while maintaining an output
file with manageable size. Looking at Fig. 4, which shows a scatter plot of
absolute redshift uncertainty against redshift for all the galaxies in GLADE+,
it is clear that there are several distinct regions. These correspond to the
different galaxy surveys which make up GLADE+, for which different uncertainty
models were used.181818Some galaxies can be seen to have $\sigma_{z}>z$; these
galaxies make up a minority of the sample ($\sim 0.7\%$) and are an artefact
of the photometric redshift models which were used. We identified 4 redshift
ranges in GLADE+ where the galaxies with the smallest $z$ uncertainties
($\sigma_{z}$) show differing dependence on $z$ and we construct the grid
accordingly, shown in orange. The grid starts at $z=10^{-6}$ and ends at
$z=10$, switching between linear and logarithmic spacing in order to closely
match the data. By setting a requirement of having at least 1 data point in
the range $[z,z+\sigma_{z}]$ for each galaxy, the resolution is chosen in such
a way that is guaranteed to capture all of the galaxy structure in the
catalogue. This does mean that in the future, as redshift uncertainties for
galaxy surveys decrease, the resolution will need to increase in order to
ensure this information is captured.

Figure 4: Scatter plot of galaxy redshift absolute uncertainty vs redshift for
$K$-band galaxies in the GLADE+ catalogue, where the colour shows the number
of galaxies. The spacing of the corresponding redshift array against redshift
is shown in orange, demonstrating the changes in spacing required to capture
all the GLADE+ galaxies with the minimum number of array elements.

The choice to use the GLADE+ catalogue, rather than a deeper galaxy survey, is
motivated in part so that the results which follow are directly comparable to
[34], and partially because, by the method’s construction, the same LOS prior
should be used for every GW event. Using GLADE+ prioritises the catalogue’s
broad sky coverage, at the expense of greater depth in other parts of the sky.
It is known that the Dark Energy Survey (DES) [54] covers the footprints of
GW190814 and GW170814 (two well-localised binary mergers) to a much higher
level of completion, and will thus produce a more informative result for those
events than the GLADE+ catalogue [24, 26]. However, the inclusion of these
deeper surveys is left to future work.

##  4 Results

We choose to use the same fiducial population models as [34]. For mass models,
the Power Law + Peak mass distribution is used for BBH. The distribution of
the primary mass, $m_{1}$, is given by

$p(m_{1}|M_{\text{min}}^{\text{BH}},M_{\text{max}}^{\text{BH}},\alpha,\lambda_{g},\mu_{g},\sigma_{g})=[(1-\lambda_{g})\mathcal{P}(m_{1}|M_{\text{min}}^{\text{BH}},M_{\text{max}}^{\text{BH}},-\alpha)+\lambda_{g}\mathcal{G}(m_{1}|\mu_{g},\sigma_{g})],$
(4.1)

while the secondary mass, $m_{2}$ is described by a truncated power-law with
slope $\beta$ between $M_{\text{min}}^{\text{BH}}$ and $m_{1}$. The NSBH
events also use the Power Law + Peak mass distribution for the primary mass,
but with the secondary mass using a uniform prior between
$M_{\text{min}}^{\text{NS}}$ and $\min(M_{\text{max}}^{\text{NS}},m_{1})$. The
BNSs use a uniform prior between $M_{\text{min}}^{\text{NS}}$ and
$M_{\text{max}}^{\text{NS}}$. Table 1 describes the parameters used in the
mass distributions and gives the prior ranges which are used for the analyses
in Secs. 4.1 and 4.3.

POWER LAW + PEAK
Parameter Description Prior $M_{\text{min}}^{\text{BH}}$ Minimum mass of the
PL component of the black hole mass distribution.
$\mathcal{U}(2.0$M${}_{\odot},10.0$M${}_{\odot})$ $M_{\text{max}}^{\text{BH}}$
Maximum mass of the PL component of the black hole mass distribution.
$\mathcal{U}(50.0$M${}_{\odot},200.0$M${}_{\odot})$ $\alpha$ Spectral index
for the PL of the primary mass distribution. $\mathcal{U}(1.5,12.0)$ $\mu_{g}$
Mean of the Gaussian component in the primary mass distribution.
$\mathcal{U}(20.0$M${}_{\odot},50.0$M${}_{\odot})$ $\sigma_{g}$ Width of the
Gaussian component in the primary mass distribution.
$\mathcal{U}(0.4$M${}_{\odot},10.0$M${}_{\odot})$ $\lambda_{g}$ Fraction of
the model in the Gaussian component. $\mathcal{U}(0.0,1.0)$ $\delta_{m}$ Range
of mass tapering on the lower end of the mass distribution.
$\mathcal{U}(0.0$M${}_{\odot},10.0$M${}_{\odot})$ $\beta$ Spectral index for
the PL of the secondary mass distribution. $\mathcal{U}(-4.0,12.0)$
$M_{\text{min}}^{\text{NS}}$ Minimum mass of the uniform neutron star mass
distribution. $\mathcal{U}(1.0$M${}_{\odot},1.5$M${}_{\odot})$
$M_{\text{max}}^{\text{NS}}$ Maximum mass of the uniform neutron star mass
distribution. $\mathcal{U}(2.0$M${}_{\odot},5.0$M${}_{\odot})$

Table 1: Summary of mass distribution parameters with the corresponding prior
ranges.

The CBC merger rate evolution model is assumed to follow a Madau-Dickinson-
like distribution [55] which is given by:

$R(z)=R_{0}(1+z)^{\gamma}\,\frac{1+(1+z_{p})^{-(\gamma+\kappa)}}{1+\left(\frac{1+z}{1+z_{p}}\right)^{\gamma+\kappa}}.$
(4.2)

Definitions of the parameters appearing in Eq. 4.2 are given in Tab. 2, as
well as the prior ranges which are used for the analyses in Secs. 4.1 and 4.3.

MERGER RATE SHAPE PARAMETERS  Parameter Description Prior $R_{0}$ Local merger
rate. $1/R_{0}$ (implicit) $\gamma$ Power-law index describing the merger rate
at low redshift. $\mathcal{U}(0,12.0)$ $\kappa$ Power-law index describing the
merger rate at high redshift. $\mathcal{U}(0,6.0)$ $z_{p}$ The redshift where
the slope of the merger rate changes. $\mathcal{U}(0,4.0)$

Table 2: Summary of merger rate shape parameters with the corresponding
priors.

Using the selection of events described in 3.1, three separate analyses are
carried out: a pure population analysis using 42 BBHs (see Sec. 4.1), a galaxy
catalogue analysis with fixed population assumptions using all the BBHs, NSBHs
and BNSs (see Sec. 4.2) and, finally, a population + galaxy catalogue analysis
which uses all GW events and galaxy catalogue information, and jointly
estimates cosmological and population parameters (Sec. 4.3). For the multi-
parameter estimation analyses, gwcosmo utilises nested sampling with the
sampler dynesty [56],191919https://zenodo.org/record/7995596 while for the
$H_{0}$-only analysis a gridded method is used.

###  4.1 The population analysis of GWTC-3 BBHs

In this section the 42 BBHs are reanalysed using gwcosmo, allowing
cosmological and population parameters to vary. Fig. 5 shows the result of the
analysis for $H_{0}$ and the parameters which correlate most strongly with it,
namely $\mu_{g},M_{\text{min}}^{\text{BH}},M_{\text{max}}^{\text{BH}}$ and
$\sigma_{g}$ from the mass model (see Table 1) and $\gamma$ from the merger
rate model (see Table 2). The full corner plot showing all parameters can be
found in appendix A. Marginalising over population parameters results in a
value of $H_{0}=$$33^{+51}_{-18}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ (maximum
a posteriori probability (MAP) and 68% highest density interval (HDI)). All
results from this point on will be quoted using the same, unless specified
otherwise. The result is fully compatible with the Icarogw result from [34],
with minor deviations due to the versions of GW parameter estimation which
were used in each analysis ([34] used parameter estimation from the GWTC-2
release, whereas we make use of the updated GWTC-2.1 release), and the fact
that gwcosmo implicitly marginalises over $R_{0}$ with a $1/R_{0}$ prior,
while in [34] Icarogw used a uniform prior on $R_{0}$.

Figure 5: Posteriors on parameters $H_{0}$, $\gamma$, $\mu_{g}$,
$M_{\text{max}}^{\text{BH}}$, $M_{\text{min}}^{\text{BH}}$ and $\sigma_{g}$
using the 42 BBH events of the GWTC-3 catalog, obtained with a pure population
analysis with gwcosmo.

###  4.2 The galaxy catalogue analysis of GWTC-3 with fixed population
assumptions

Figure 6: $H_{0}$ posteriors (uniform prior) for the 46 dark sirens from
GWTC-3 with SNR $>11$. The dotted line shows the contribution coming purely
from the GW population assumptions and no galaxy catalogue information (the
“empty catalogue” scenario). The dashed orange line includes galaxy catalogue
information with uniform weighting of host galaxies, while the blue dot-dashed
line includes galaxy catalogue information and applies weighting proportional
to the galaxy’s luminosity.

In order to best assess the individual contributions of GW events to the
$H_{0}$ measurement, as well as determine to what extent a result is driven by
population information versus galaxy catalogue information, it is useful to
run an analysis with fixed population assumptions. This analysis is carried
out for the 46 dark siren events from GWTC-3 with SNR $>11$. For the Power Law
+ Peak mass model we use $\alpha=3.78$, $\beta=0.81$,
$M_{\text{max}}^{\text{BH}}=112.5M_{\odot}$,
$M_{\text{min}}^{\text{BH}}=4.98M_{\odot}$ , $\delta_{m}=4.8M_{\odot}$
$\mu_{g}=32.27M_{\odot}$ , $\sigma_{g}=3.88M_{\odot}$ and $\lambda_{g}=0.03$.
For neutron stars we fix $M_{\text{min}}^{\text{NS}}=1M_{\odot}$and
$M_{\text{max}}^{\text{NS}}=3M_{\odot}$. For the evolving merger rate model we
use $\gamma=4.59$, $k=2.86$ and $z_{p}=2.47$. These values match those used
for the gwcosmo analysis in [34], using the previous version of gwcosmo,
making the results directly comparable.

The posteriors for individual events (assuming uniform priors on $H_{0}$) can
be seen in Fig. 6, for three different analysis assumptions: luminosity
weighting of host galaxies in the $K$-band, uniform weighting of host
galaxies, and no galaxy catalogue (also known as the “empty catalogue”
analysis). The empty catalogue analysis demonstrates the information which is
driven by the chosen GW population, and as such it is the difference between
this line and the others which should be used to assess how informative the
galaxy catalogue method is for any given event. For the majority of events,
lying at distances higher than that to which the GLADE+ catalogue is complete
in the $K$-band, the results of all 3 analyses are the same. Some nearer-by
events show deviations at the low-$H_{0}$ end, which corresponds to low
redshift (for a fixed luminosity distance) where the galaxy catalogue is more
complete. In general, the luminosity weighted results are the most
informative, as the majority of the GW host probability gets condensed into a
smaller number of galaxies, while in the unweighted case the presence of many
low-luminosity galaxies means that the final LOS redshift prior will appear
much closer to uniform in comoving volume when averaged over the volume of an
event.

Figure 7: $H_{0}$ posteriors for GW190814 using different values for $n_{\rm
high}$ and $n_{\rm map}$. Figure 8: Ratio of the galaxy catalogue and empty
catalogue LOS redshift prior, weighted by the sky area of GW190814 with
$n_{\rm high}=128$. The blue curve corresponds to the $n_{\rm map}$=32 case
and the orange curve corresponds to $n_{\rm map}$=64. The horizontal dashed
black line at 1.0 corresponds to the scenario when matter is distributed truly
uniformly in comoving volume, so oscillations above and below this correspond
to over- and under-densities in matter. At high redshift, where the catalogue
contains no information in the $K$-band, the curves tend to 1.0 as the
incompleteness correction dominates.

The event of most interest, from a galaxy catalogue perspective, is GW190814.
Due to its nearby distance (around 240 Mpc) and small sky localisation (18.5
$\text{deg}^{2}$ for the 90% HDI) [57] this event is an ideal dark siren. In
order to assess the impact of analysis assumptions relating to the choice of
pixel resolution, GW190814 is reanalysed under a variety of different
assumptions (see Fig. 7). Two different resolutions of $m_{\text{th}}$ and
$N_{\text{gal}}$ maps are used: $n_{\rm map}$=32 and 64 (see Figs. 2 and 3),
and for each of these, different resolutions up to $n_{\rm high}$=128 are used
for the final LOS prior. Interestingly, the choice of map resolution seems to
have a larger impact on the shape of the posterior than that of
$n_{\text{high}}$. Looking at Fig. 8, it can be seen that while the same
galaxies are present in the LOS redshift prior for the $n_{\rm map}$=32 and 64
cases, the relative contributions of these galaxies can differ significantly.
The driving reason for this is the small number statistics discussed in detail
in Sec. 2.1.5: in the $n_{\rm map}$=64 scenario, a lot of pixels within
GW190814’s sky area only contain about a dozen galaxies, meaning that the
relative contributions from adjacent pixels can be significantly impacted. The
larger pixel sizes in the $n_{\rm map}$=32 scenario reduce this impact, and
the results are more robust because of it.

The value of $n_{\text{high}}$ also impacts the final shape of the posterior.
This is most noticeable for the $n_{\text{high}}=32$ result, compared with the
$n_{\text{high}}=64$ and 128 results using $n_{\rm map}$=32 in all cases. A
lower redshift galaxy which contributes a large peak for $n_{\text{high}}=32$
is subsequently suppressed at higher resolutions, because the pixel it lies
within contains significant variation in GW support: a detail which is only
captured at higher resolutions. The $n_{\text{high}}=64$ and 128 results are
very similar, showing a reasonable level of convergence has been reached. As
GW190814 is the most well-localised dark siren in GWTC-3, we conclude that
this resolution is sufficient to capture all the necessary GW information for
this analysis. It is worth remembering, however, that the choice of resolution
which is “good enough” will depend on the localisation area of the GW events
it is used with – a future GW event with a smaller localisation region that
GW190814 will need a higher resolution LOS prior to do it justice.

###  4.3 Combined cosmological and population inference using GWTC-3

With the robustness of this version of gwcosmo established, the final step is
to do a full population + galaxy catalogue analysis of GWTC-3. The BBHs are
analysed as in Sec. 4.1, but this time utilising a LOS redshift prior
constructed from the GLADE+ $K$ band, using a resolutions of
$n_{\text{map}}=32$ and $n_{\text{high}}=128$. The two NSBHs and GW190814 are
also analysed using the full population + galaxy catalogue method. The maximum
and minimum mass of the neutron star (NS) component is allowed to vary in this
analysis, as GW190814’s secondary mass is very heavy for a NS and we want to
avoid a cut off in the mass distribution which could incorporate additional
cosmological information into the analysis (see appendix A for the full
results of this analysis). When considering these 3 events alone, the majority
of population parameters are poorly constrained (as expected, see also the
NSBH analysis in [36]). Because GW190425 is the only BNS without an EM
counterpart we do not carry out a population analysis for it. Given the lack
of features in the assumed BNS mass distribution (assuming the event is not so
heavy or light as to touch the BNS prior bounds), we simply fix the population
parameters to fiducial values and carry out a 1-dimensional analysis over
$H_{0}$.

Figure 9: Posterior on $H_{0}$ produced using 42 BBHs, 2 NSBHs, GW190814, and
GW190425 using the GLADE+ galaxy catalogue and GW170817 with its EM
counterpart (black line). The separate contributions from the BBHs (blue) the
NSBHs and GW190814 (orange) and GW170817 (green) are also shown (normalised).
The $1\sigma$ Planck [1] and SH0ES [2] results are shown as shaded pink and
brown bands respectively.

The result which combines the inference from all 3 GW populations is presented
in Fig. 9. The final constraint on $H_{0}$ is $69^{+12}_{-7}$ $\mbox{km
s}^{-1}\,\mbox{Mpc}^{-1}$ which, as anticipated, is less informative than the
previous result produced by gwcosmo in [34] of $H_{0}=$ $68^{+8}_{-6}$
$\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$, due to the marginalisation over the
uncertainty of the population hyper-parameters. While less informative, this
result is far more robust, and no longer susceptible to bias due to fixed
assumptions about the GW population.

##  5 Conclusions

This paper presents a novel method for cosmological analyses using GW data and
galaxy catalogues, incorporating information from the GW mass distribution in
a way which allows the joint estimation of cosmological and GW population
parameters, a crucial improvement while the true mass distribution of compact
mergers remains uncertain. A key part of this new method is the construction
of a LOS redshift prior which incorporates both galaxy catalogue data and an
incompleteness correction, the pre-computation of which produces a speed-up in
the analysis of roughly 3 orders of magnitude over the previous version of
gwcosmo. Comparing the population-only cosmological results using 42 BBHs from
GWTC-3 to the Icarogw results in [34] shows excellent agreement between the
methods. Similarly, the fixed population result using the GLADE+ galaxy
catalogue agrees very well with the previous version of gwcosmo in [34]. The
ability to jointly estimate, or marginalise over, the population of compact
mergers means that this analysis is largely robust to the current
uncertainties in the population. The LVK’s fourth observing run (O4) is
currently underway and is expected to increase the number of GW detections to
several hundreds. This improved version of gwcosmo will allow these new events
to be used in an informative cosmological analysis which incorporates data
from current galaxy surveys, and also utilises information from the GW
population itself. Nearby and well localised sources will make informative
contributions due to the known galaxy structure within their localisation
volumes. The numerous GW events at high distances, where catalogue support is
low will still contribute cosmologically through their redshifted masses. By
joining these two sources of redshift information the best of both is
utilised, enabling reliable and informative future measurements of $H_{0}$
using dark sirens.

However, there are still areas in which further improvements must be made.
Here, for example, in the final result we treat population hyper-parameters as
uncorrelated between the different source populations. Ideally, parameters
which are common between the populations (_e.g._ the hyper-parameters
describing the black hole as the primary mass in both the BBH and NSBH cases)
would be estimated jointly. However for certain parameters, such as those
parameterising the merger rate evolution model, one would not expect to be the
same between the different GW populations. There, and in general perhaps, the
solution could be to infer hyper-parameters for the entire CBC population,
rather than breaking it into sub categories. Given the current difficulties in
firmly placing current GW detections into these categories (GW190814 being an
excellent example, not to mention the current uncertainty surrounding the
existence – or not – of a mass gap between the NS and BH population [40]) this
could simplify the process noticeably.

While the gwcosmo method is now able to jointly constrain cosmological and
population hyper-parameters, it shouldn’t be forgotten that the choice of mass
model can have an impact on the inference, and ensuring that a model is used
which adequately captures the structure in the compact binary mass
distribution will be an important aspect of this analysis moving forwards.
Similarly, while current evidence does not confirm whether the mass
distribution of compact binaries evolves with redshift or remains constant,
this would be an important feature to include should we reach the point where
the data shows a strong preference one way or the other.

From the galaxy catalogue side of things, the current gwcosmo methodology is
much improved to handle future data releases, as the speed of the analysis is
not tied to the number of galaxies contained within the catalogue. It can also
be adapted with much more ease in future to non-gaussian redshift
uncertainties, which may be of high importance in the regime where the galaxy
catalogue information becomes a more dominant contributor to the measurement
of $H_{0}$. With current galaxy surveys now publicly available which are
complete to higher redshifts than the GLADE+ catalogue (such as the DESI
Legacy Imaging Surveys [58] and DES [54]), and upcoming surveys promising ever
greater depths and redshift precision (such as the DESI survey [59], LSST [60]
and Euclid [61]) the era of galaxy catalogue driven dark siren analyses is
just about to begin.

The main results of this paper, $H_{0}=$ $69^{+12}_{-7}$ $\mbox{km
s}^{-1}\,\mbox{Mpc}^{-1}$, provides a robust measurement of the Hubble
constant which combines information from GW events in combination with EM
counterparts, galaxy catalogues, and redshifted mass information, paving the
way forward for future GW cosmology analyses which can utilise information
from an array of different sources. With the number of GW events set to vastly
increase over the coming years, and the galaxy survey information becoming
available, this is a truly exciting time for the field of gravitational-wave
cosmology.

## Acknowledgments

The authors would like to thank Jonathan Gair and Sylvain Marsat for their
contributions to reviewing the gwcosmo code, as well as Danny Laghi and
Suvodip Mukherjee for their useful comments on the paper draft.

This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation. The authors
are grateful for computational resources provided by the LIGO Laboratory and
supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.
The research of RG was suported by ERC starting grant SHADE 949572 and STFC
grant ST/V005634/1. FB, CT and AG are supported by the Ghent University
Special Research Funds (BOF) project BOF/STA/202009/040 and the Fonds
Wetenschappelijk Onderzoek (FWO) iBOF project BOF20/IBF/124. A.C. is supported
by a PhD grant from the Chinese Scholarship Council (grant no.202008060014).
T.B. is supported by ERC Starting Grant SHADE (grant no. StG 949572) and a
Royal Society University Research Fellowship (grant no.
URF$\backslash$R1$\backslash$180009). A.E.R. was supported by the UDEA
Dedicacion exclusiva and Sostenibilidad programs and the CODI projects
2019-28270 and 2021-44670 of UDEA.

##  Appendix A Full population results

Here we present the full parameter constraints resulting from the analyses
described in section 4.

Figure 10: Posteriors on 12 cosmological and population parameters using the
42 BBH events of the GWTC-3 catalog, obtained with a pure population analysis
with gwcosmo. Figure 11: Posteriors on 14 cosmological and population
parameters using the 2 NSBH events (GW200105 and GW200115) and GW190814,
including information from the GLADE+ galaxy catalogue.

## References

  * [1] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [2] A.G. Riess, W. Yuan, L.M. Macri, D. Scolnic, D. Brout, S. Casertano et al., _A comprehensive measurement of the local value of the hubble constant with 1 km s-1 mpc-1 uncertainty from the hubble space telescope and the sh0es team_ , _The Astrophysical Journal Letters_ 934 (2022) L7.
  * [3] L. Verde, T. Treu and A.G. Riess, _Tensions between the early and late Universe_ , _Nature Astronomy_ 3 (2019) 891 [1907.10625].
  * [4] E.D. Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri et al., _In the realm of the hubble tension—a review of solutions*_ , _Classical and Quantum Gravity_ 38 (2021) 153001.
  * [5] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _GW170817: Observation of gravitational waves from a binary neutron star inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101.
  * [6] B.P. Abbott, R. Abbott, T.D. Abbott, F. Acernese, K. Ackley, C. Adams et al., _Multi-messenger Observations of a Binary Neutron Star Merger_ , _Astrophys. J._ 848 (2017) L12.
  * [7] J. Aasi et al., _Advanced ligo_ , _Classical and Quantum Gravity_ 32 (2015) 074001.
  * [8] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca et al., _Advanced virgo: a second-generation interferometric gravitational wave detector_ , _Classical and Quantum Gravity_ 32 (2014) 024001.
  * [9] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _A gravitational-wave standard siren measurement of the Hubble constant_ , _Nature_ 551 (2017) 85 [1710.05835].
  * [10] T. Akutsu, M. Ando, K. Arai, Y. Arai, S. Araki, A. Araya et al., _Overview of KAGRA: Detector design and construction history_ , _Progress of Theoretical and Experimental Physics_ 2021 (2020) 05A101 [https://academic.oup.com/ptep/article-pdf/2021/5/05A101/37974994/ptaa125.pdf].
  * [11] The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration collaboration, _GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run_ , _arXiv preprint_ (2021) .
  * [12] A. Buikema, C. Cahillane, G.L. Mansell, C.D. Blair, R. Abbott, C. Adams et al., _Sensitivity and performance of the advanced ligo detectors in the third observing run_ , _Phys. Rev. D_ 102 (2020) 062003.
  * [13] Virgo Collaboration collaboration, _Increasing the astrophysical reach of the advanced virgo detector via the application of squeezed vacuum states of light_ , _Phys. Rev. Lett._ 123 (2019) 231108.
  * [14] M. Oguri, _Measuring the distance-redshift relation with the cross-correlation of gravitational wave standard sirens and galaxies_ , _Phys. Rev. D_ 93 (2016) 083511.
  * [15] S. Mukherjee and B.D. Wandelt, _Beyond the classical distance-redshift test: cross-correlating redshift-free standard candles and sirens with redshift surveys_ , _arXiv e-prints_ (2018) arXiv:1808.06615 [1808.06615].
  * [16] S. Mukherjee, B.D. Wandelt and J. Silk, _Probing the theory of gravity with gravitational lensing of gravitational waves and galaxy surveys_ , _Mon. Not. R. Astron. Soc._ 494 (2020) 1956 [1908.08951].
  * [17] S. Mukherjee, B.D. Wandelt, S.M. Nissanke and A. Silvestri, _Accurate precision cosmology with redshift unknown gravitational wave sources_ , _prd_ 103 (2021) 043520 [2007.02943].
  * [18] C. Cigarrán Díaz and S. Mukherjee, _Mapping the cosmic expansion history from LIGO-Virgo-KAGRA in synergy with DESI and SPHEREx_ , _Mon. Not. R. Astron. Soc._ 511 (2022) 2782 [2107.12787].
  * [19] S. Bera, D. Rana, S. More and S. Bose, _Incompleteness Matters Not: Inference of H 0 from Binary Black Hole-Galaxy Cross-correlations_, _Astrophys. J._ 902 (2020) 79 [2007.04271].
  * [20] S. Mukherjee, A. Krolewski, B.D. Wandelt and J. Silk, _Cross-correlating dark sirens and galaxies: measurement of $H_{0}$ from GWTC-3 of LIGO-Virgo-KAGRA_, _arXiv e-prints_ (2022) arXiv:2203.03643 [2203.03643].
  * [21] B.F. Schutz, _Determining the Hubble constant from gravitational wave observations_ , _Nature_ 323 (1986) 310.
  * [22] W. Del Pozzo, _Inference of cosmological parameters from gravitational waves: Applications to second generation interferometers_ , _Phys. Rev. D_ 86 (2012) 043011 [1108.1317].
  * [23] M. Fishbach, D.E. Holz and W.M. Farr, _Does the Black Hole Merger Rate Evolve with Redshift?_ , _Astrophys. J._ 863 (2018) L41 [1805.10270].
  * [24] DES, LIGO Scientific, Virgo collaboration, _First measurement of the Hubble constant from a dark standard siren using the Dark Energy Survey galaxies and the LIGO/Virgo binary-black-hole merger GW170814_ , _Astrophys. J._ 876 (2019) L7 [1901.01540].
  * [25] R. Gray et al., _Cosmological inference using gravitational wave standard sirens: A mock data analysis_ , _Phys. Rev. D_ 101 (2020) 122001 [1908.06050].
  * [26] A. Palmese, J. deVicente, M.E.S. Pereira, J. Annis, W. Hartley, K. Herner et al., _A statistical standard siren measurement of the hubble constant from the LIGO/virgo gravitational wave compact object merger GW190814 and dark energy survey galaxies_ , _Astrophys. J. Lett._ 900 (2020) L33.
  * [27] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore and M. Mancarella, _Cosmology with LIGO/Virgo dark sirens: Hubble parameter and modified gravitational wave propagation_ , _JCAP_ 08 (2021) 026 [2101.12660].
  * [28] R. Gray, C. Messenger and J. Veitch, _A pixelated approach to galaxy catalogue incompleteness: improving the dark siren measurement of the Hubble constant_ , _Monthly Notices of the Royal Astronomical Society_ 512 (2022) 1127 [https://academic.oup.com/mnras/article-pdf/512/1/1127/45303118/stac366.pdf].
  * [29] D.F. Chernoff and L.S. Finn, _Gravitational Radiation, Inspiraling Binaries, and Cosmology_ , _Astrophys. J. Lett._ 411 (1993) L5 [gr-qc/9304020].
  * [30] S.R. Taylor, J.R. Gair and I. Mandel, _Hubble without the Hubble: Cosmology using advanced gravitational-wave detectors alone_ , _Phys. Rev._ D85 (2012) 023535 [1108.5161].
  * [31] W.M. Farr, M. Fishbach, J. Ye and D.E. Holz, _A Future Percent-level Measurement of the Hubble Expansion at Redshift 0.8 with Advanced LIGO_ , _Astrophys. J. Lett._ 883 (2019) L42 [1908.09084].
  * [32] S. Mastrogiovanni, K. Leyde, C. Karathanasis, E. Chassande-Mottin, D.A. Steer, J. Gair et al., _On the importance of source population models for gravitational-wave cosmology_ , _Phys. Rev. D_ 104 (2021) 062009.
  * [33] J.M. Ezquiaga and D.E. Holz, _Spectral sirens: Cosmology from the full mass distribution of compact binaries_ , _Phys. Rev. Lett._ 129 (2022) 061102.
  * [34] R. Abbott et al., _Constraints on the cosmic expansion history from GWTC-3_ , _arXiv e-prints_ (2021) arXiv:2111.03604 [2111.03604].
  * [35] G. Dálya, R. Díaz, F.R. Bouchet, Z. Frei, J. Jasche, G. Lavaux et al., _GLADE+ : an extended galaxy catalogue for multimessenger searches with advanced gravitational-wave detectors_ , _Monthly Notices of the Royal Astronomical Society_ 514 (2022) 1403 [https://academic.oup.com/mnras/article-pdf/514/1/1403/44001190/stac1443.pdf].
  * [36] A. Chen, R. Gray and T. Baker, _Testing the nature of gravitational wave propagation using dark sirens and galaxy catalogues_ , _arXiv e-prints_ (2023) arXiv:2309.03833 [2309.03833].
  * [37] S. Mastrogiovanni, D. Laghi, R. Gray, G. Caneva Santoro, A. Ghosh, C. Karathanasis et al., _A novel approach to infer population and cosmological properties with gravitational waves standard sirens and galaxy surveys_ , _arXiv e-prints_ (2023) arXiv:2305.10488 [2305.10488].
  * [38] I. Mandel, W.M. Farr and J.R. Gair, _Extracting distribution parameters from multiple uncertain observations with selection biases_ , _Mon. Not. R. Astron. Soc._ 486 (2019) 1086 [1809.02063].
  * [39] S. Vitale, D. Gerosa, W.M. Farr and S.R. Taylor, _Inferring the properties of a population of compact binaries in presence of selection effects_ , 2007.05579.
  * [40] LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration collaboration, _Population of merging compact binaries inferred using gravitational waves through gwtc-3_ , _Phys. Rev. X_ 13 (2023) 011048.
  * [41] S. Mukherjee, _The redshift dependence of black hole mass distribution: is it reliable for standard sirens cosmology?_ , _Mon. Not. R. Astron. Soc._ 515 (2022) 5495 [2112.10256].
  * [42] C. Karathanasis, S. Mukherjee and S. Mastrogiovanni, _Binary black holes population and cosmology in new lights: signature of PISN mass and formation channel in GWTC-3_ , _Mon. Not. R. Astron. Soc._ 523 (2023) 4539 [2204.13495].
  * [43] K.M. Górski, E. Hivon, A.J. Banday, B.D. Wandelt, F.K. Hansen, M. Reinecke et al., _HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere_ , _Astrophys. J._ 622 (2005) 759 [arXiv:astro-ph/0409513].
  * [44] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon et al., _healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in python_ , _Journal of Open Source Software_ 4 (2019) 1298.
  * [45] The LIGO Scientific Collaboration, the Virgo Collaboration, R. Abbott et al., _GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run_ , _arXiv e-prints_ (2021) arXiv:2108.01045 [2108.01045].
  * [46] D.J. White, E.J. Daw and V.S. Dhillon, _A List of Galaxies for Gravitational Wave Searches_ , _Class. Quant. Grav._ 28 (2011) 085016 [1103.0695].
  * [47] D. Makarov, P. Prugniel, N. Terekhova, H. Courtois and I. Vauglin, _HyperLEDA. III. The catalogue of extragalactic distances_ , _Astron. Astrophys._ 570 (2014) A13 [1408.3476].
  * [48] T.H. Jarrett, T. Chester, R. Cutri, S. Schneider, M. Skrutskie and J.P. Huchra, _2MASS Extended Source Catalog: Overview and Algorithms_ , _aj_ 119 (2000) 2498 [astro-ph/0004318].
  * [49] M.F. Skrutskie et al., _The Two Micron All Sky Survey (2MASS)_ , _Astron. J._ 131 (2006) 1163.
  * [50] M. Bilicki, T.H. Jarrett, J.A. Peacock, M.E. Cluver and L. Steward, _2MASS Photometric Redshift catalog: a comprehensive three-dimensional census of the whole sky_ , _Astrophys. J. Suppl._ 210 (2014) 9 [1311.5246].
  * [51] M. Bilicki, J.A. Peacock, T.H. Jarrett, M.E. Cluver, N. Maddox, M.J.I. Brown et al., _Wise × supercosmos photometric redshift catalog: 20 million galaxies over 3 $\pi$ steradians_, _The Astrophysical Journal Supplement Series_ 225 (2016) 5.
  * [52] B.W. Lyke, A.N. Higley, J.N. McLane, D.P. Schurhammer, A.D. Myers, A.J. Ross et al., _The sloan digital sky survey quasar catalog: Sixteenth data release_ , _The Astrophysical Journal Supplement Series_ 250 (2020) 8.
  * [53] C.S. Kochanek, M.A. Pahre, E.E. Falco, J.P. Huchra, J. Mader, T.H. Jarrett et al., _The k-band galaxy luminosity function* **_ , _The Astrophysical Journal_ 560 (2001) 566.
  * [54] The Dark Energy Survey Collaboration, _The Dark Energy Survey_ , _arXiv e-prints_ (2005) astro [astro-ph/0510346].
  * [55] P. Madau and M. Dickinson, _Cosmic star-formation history_ , _Annual Review of Astronomy and Astrophysics_ 52 (2014) 415–486.
  * [56] J.S. Speagle, _DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences_ , _Mon. Not. R. Astron. Soc._ 493 (2020) 3132 [1904.02180].
  * [57] R. Abbott et al., _GW190814: Gravitational waves from the coalescence of a 23 solar mass black hole with a 2.6 solar mass compact object_ , _Astrophys. J. Lett._ 896 (2020) L44.
  * [58] A. Dey, D.J. Schlegel, D. Lang, R. Blum, K. Burleigh, X. Fan et al., _Overview of the DESI Legacy Imaging Surveys_ , _Astron. J._ 157 (2019) 168 [1804.08657].
  * [59] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L.E. Allen et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , _arXiv e-prints_ (2016) arXiv:1611.00036 [1611.00036].
  * [60] Ž. Ivezić, S.M. Kahn, J.A. Tyson, B. Abel, E. Acosta, R. Allsman et al., _LSST: From Science Drivers to Reference Design and Anticipated Data Products_ , _Astrophys. J._ 873 (2019) 111 [0805.2366].
  * [61] Euclid Collaboration, R. Scaramella, J. Amiaux, Y. Mellier, C. Burigana, C.S. Carvalho et al., _Euclid preparation. I. The Euclid Wide Survey_ , _Astron. Astrophys._ 662 (2022) A112 [2108.01201].

  *[$H_{0}$]: the Hubble constant
  *[$\Lambda$CDM]: Λ-cold-dark-matter
  *[CBC]: compact binary coalescence
  *[GW]: gravitational wave
  *[EM]: electromagnetic
  *[GWs]: gravitational wave
  *[GWTC-3]: the third Gravitational-Wave Transient Catalogue
  *[BBHs]: binary black hole
  *[BH]: black hole
  *[LOS]: line-of-sight
  *[CBCs]: compact binary coalescence
  *[KDE]: kernel density estimate
  *[RA]: right ascension
  *[dec]: declination
  *[$m_{\text{th}}$]: apparent magnitude threshold
  *[SFR]: star formation rate
  *[SNR]: signal-to-noise ratio
  *[NSBHs]: neutron star - black hole pair
  *[NSBH]: neutron star - black hole pair
  *[BNSs]: binary neutron star
  *[KDEs]: kernel density estimate
  *[RAs]: right ascension
  *[decs]: declination
  *[DES]: Dark Energy Survey
  *[BBH]: binary black hole
  *[NS]: neutron star
  *[BNS]: binary neutron star