# Testing the nature of gravitational wave propagation using dark sirens and
galaxy catalogues

Anson Chen  Rachel Gray  Tessa Baker

###### Abstract

The dark sirens method enables us to use gravitational wave events without
electromagnetic counterparts as tools for cosmology and tests of gravity.
Furthermore, the dark sirens analysis code gwcosmo can now robustly account
for information coming from both galaxy catalogues and the compact object mass
distribution. We present here an extension of the gwcosmo code and methodology
to constrain parameterized deviations from General Relativity that affect the
propagation of gravitational waves. We show results of our analysis using data
from the GWTC-3 gravitational wave catalogues, in preparation for application
to the O4 observing run. After testing our pipelines using the First Two Years
mock data set, we reanalyse 46 events from GWTC-3, and combine the posterior
for BBH and NSBH sampling results for the first time. We obtain joint
constraints on $H_{0}$ and parameterized deviations from General Relativity in
the Power Law + Peak BBH population model. With increased galaxy catalogue
support in the future, our work sets the stage for dark sirens to become a
powerful tool for testing gravity.

##  1 Introduction

As the imprint of a propagating gravitational field perturbation,
gravitational waves (GWs) are the most natural tool with which to test the
nature of gravity. The nature of their generation and subsequent propagation
means that they uniquely access both the strong-field regime of gravity near
compact objects, as well as the weak-field regime of gravity on cosmological
scales [1, 2]. Furthermore, unlike many electromagnetic (EM) observables, they
carry direct information about the luminosity distance of their source. When
combined with redshift information, we can therefore use gravitational wave
data to probe the luminosity distance-redshift relation – and hence key
cosmological parameters like the Hubble constant – in a manner distinct to
traditional methods such as the Cosmic Microwave Background or the local
distance ladder [3, 4]. Some departures from General Relativity (GR) will also
impact the luminosity distance-redshift relation, and hence can be tested in
this way.

For both tests of gravity and measurements of $H_{0}$, the presence of an EM
counterpart to a GW source – a “bright siren” event first introduced by [5] –
provides the most direct and (usually) powerful constraints. This was seen
most spectacularly with binary neutron star (BNS) merger GW170817 and its EM
counterpart [6, 7]. Here, the localization of the accompanying kilonova
allowed the galaxy NGC4993 to be uniquely identified as the host galaxy of
GW170817 [8]. With an additional correction for the peculiar velocity of
NGC4993, from this single event alone the Hubble constant was estimated to be
$H_{0}=70.0^{+12.0}_{-8.0}~{}{\rm km}~{}{\rm s}^{-1}{\rm Mpc}^{-1}$ [9]. The
implications of GW170817 for tests of gravity were even more impactful. The
arrival of the GRB counterpart just $\sim 1.7$s after the GW merger signal
constrained the propagation speed of GWs ($c_{T}$) in the frequency band of
terrestrial detectors111A caveat is that an effective field approach to
modified gravity still allows a deviation of $c_{T}$ from $c$ at an energy
scale lower than the terrestrial detection band [10]. This scenario can be
probed with massive BBHs by the Laser Interferometer Space Antenna (LISA)
alone [11], or with stellar origin BBHs by joint detection of LVK and LISA
[12, 13], or possibly by pulsar timing arrays. at the level
$|1-c_{T}/c|\leqslant 3\times 10^{-15}$ [7]. This led to stringent constraints
on various modified gravity models, such as scalar-tensor and vector-tensor
theories[14, 15, 16], Born-Infeld gravity[17] and parity-violating
gravity[18].

However, so far bright sirens have proved relatively elusive. With current
broad uncertainties on the underlying event rates, it remains unclear to what
extent we can rely on them as cosmological tools. This has led to the
development of an alternative set of techniques for GW cosmology that do not
rely on EM counterparts, known as dark siren methods. These techniques instead
employ galaxy catalogues, together with features in the mass spectrum of
compact objects, to supply redshift information about the possible hosts of GW
sources. Under these methods, the redshift of a GW host galaxy cannot be
uniquely identified; instead the galaxy catalogue is used to construct a prior
distribution for the redshift along specific lines of sight (in a pixellated
fashion [19]). This prior is combined with GW data in a hierarchical Bayesian
formalism [20, 21, 22, 23] to infer cosmological parameters. The population
model of compact objects also informs this result because, assuming a known
distribution of compact object masses in the source frame, this is related to
the distribution that GW detectors actually measure (in ‘detector frame’) by
factors of $(1+z)$. Hence for any distinct features in the mass distribution
we can effectively compare their assumed and measured locations to infer
information about source redshifts, again in a statistical sense only.

To date, the predominant focus of dark sirens development has been to
constrain the Hubble constant (but see also [24, 25]). The primary goal of
this paper is to introduce an extension of a key dark sirens tool, described
below, which allows the method to test extensions of GR impacting the
luminosity distance-redshift relation, i.e. one of the most generic effects of
cosmologically-motivated modified gravity theories. With the fourth observing
run of the LIGO-Virgo-KAGRA (LVK) collaboration underway, we ready our tools
for application to the increasing stack of GW detections, to yield (we
anticipate) increasingly precise results.

At present, two sophisticated code packages have been independently developed
to measure cosmological parameters with the bright siren and/or the dark siren
methods: gwcosmo [19, 26] and Icarogw [27, 28]. Both pipelines are now capable
of drawing on information from both the statistical host galaxy distribution
from galaxy catalogues, and also via constraining compact object population
models and merger rate redshift evolution models. The ability to use
information from both galaxy catalogues and the population model
simultaneously is an important recent development, detailed in [26], which
ensures a robust $H_{0}$ result222Previously, the first version of Icarogw
[29] constrained $H_{0}$, population and merger rate models without utilising
galaxy catalogue information; meanwhile gwcosmo v1.0.0 [19] measured $H_{0}$
using galaxy catalogue information, but fixing the population and merger rate
model..

After the O3 observation run [30], 46 dark sirens in GWTC-3 with signal-to-
noise ratios (SNRs) larger than 11 were selected for statistical host galaxy
identification with the GLADE+ galaxy catalogue [19] using gwcosmo v1.0.0
[31]. These yield a measurement of $H_{0}=67^{+13}_{-12}~{}{\rm km}~{}{\rm
s}^{-1}{\rm Mpc}^{-1}$. Combining the dark siren result with the bright siren
GW170817 yields $H_{0}=68^{+8}_{-6}~{}{\rm km}~{}{\rm s}^{-1}{\rm Mpc}^{-1}$.
Meanwhile the analysis with 42 BBHs using Icarogw 2.0 gives
$H_{0}=71^{+35}_{-30}~{}{\rm km}~{}{\rm s}^{-1}{\rm Mpc}^{-1}$ with the
assumption of the compact object binary merger rate proportional to the galaxy
luminosity, and $H_{0}=43^{+48}_{-18}~{}{\rm km}~{}{\rm s}^{-1}{\rm Mpc}^{-1}$
without the assumption [27]. Although these measurements are not more precise
than those from EM probes, they serve as an early step on a pathway towards
and independent and competitive probe of $H_{0}$. We note here that there
exists a third approach to constraining $H_{0}$ with GW data, by fully cross-
correlating GW events and galaxy clustering [32, 33, 34].

We wish to capitalise on these early successes of the dark sirens methodology,
and likewise apply them to test the laws of gravity. The purpose of this paper
is to describe the extension of the gwcosmo software to achieve precisely
this. By reanalysing GWTC-3 events and First Two Years mock data, we
demonstrate the constraints on deviations obtained by gwcosmo and lay the
groundwork for future tests of gravity at a cosmic scale in O4 that have not
been possible in previous observing runs. The inclusion of similar modified
gravity effects in the Icarogw code is described in [24, 27], whilst an
independent code pipeline and results also appears in [25].

The paper is constructed in five sections. In Section 2 we review possible
modifications to the GW luminosity distance outside of GR, and introduce three
classes of parameterization for this phenomenon. In Section 3 we revisit the
theoretical framework of gwcosmo and develop its extension to features beyond
GR. Section 4 and Section 5 show the reanalysis of mock bright siren data and
GWTC-3 events respectively, now including modified gravity parameters with the
latest version of gwcosmo. We conclude and discuss further lines of
development in Section 6.

##  2 GW Propagation Beyond GR

A common feature of gravity theories beyond GR is a modification to the
effective friction term that impacts GW propagation. In theories with this
feature, in the Jordan frame (where matter is minimally coupled to the
metric), the modified GW propagation equation can be written as [35, 36, 37,
38, 39, 25]

$\tilde{h}_{A}^{\prime\prime}+2{\cal
H}[1-\delta(\eta)]\tilde{h}_{A}^{\prime}+c^{2}k^{2}\tilde{h}_{A}=0,$ (2.1)

where $\tilde{h}(\eta,k)$ is the GW amplitude in Fourier space, subscript $A$
denotes the $+$ or $\times$ polarization of GWs, $\eta$ is conformal time, ′
denotes the derivative with respect to $\eta$, and ${\cal H}=a^{\prime}/a$ is
the conformal Hubble parameter. Here $\delta(\eta)$ is the modification to the
friction term introduced by modified gravity theories, where GR is recovered
when $\delta(\eta)=0$. It is shown in Sec. 4.1.4 of [40] that the solution to
Eq. (2.1) in GR gives $\tilde{h}(\eta,k)\propto 1/d_{L}^{\rm EM}(z)$, where
$d_{L}^{\rm EM}(z)$ is the EM luminosity distance defined by

$d_{L}^{\rm
EM}(z)=(1+z)\int_{0}^{z}\frac{c~{}d\tilde{z}}{H_{0}\sqrt{\Omega_{m,0}(1+\tilde{z})^{3}+\Omega_{\Lambda,0}}}$
(2.2)

in the late-time universe. $\Omega_{m,0}$ and $\Omega_{\Lambda,0}$ are the
matter density today and the dark energy density today respectively. When
$\delta(\eta)\neq 0$, it is shown in [35] that one can absorb the changes to
the GW propagation equation by defining a new effective conformal Hubble
factor. To do this, we introduce a new scale factor $\tilde{a}$ such that

$\frac{\tilde{a}^{\prime}}{\tilde{a}}={\cal H}[1-\delta(\eta)],$ (2.3)

When expressed in terms of the new scale factor, the friction term now has the
same form as in GR. From the solution to the transformed equation we obtain a
new quantity called the GW luminosity distance $d_{L}^{\rm GW}(z)$, which is
related to $d_{L}^{\rm EM}(z)$ by

$d_{L}^{\rm GW}(z)=\frac{a(z)}{\tilde{a}(z)}d_{L}^{\rm EM}(z),$ (2.4)

In a non-GR theory, we now have $\tilde{h}(\eta,k)\propto 1/d_{L}^{\rm
GW}(z)$. Eq. (2.3) can be written as

$\left(\log\frac{a}{\tilde{a}}\right)^{\prime}=\delta(\eta){\cal H}(\eta).$
(2.5)

Integrating Eq. (2.5) allows us to re-express the relation between $d_{L}^{\rm
GW}(z)$ and $d_{L}^{\rm EM}(z)$ as:

$d_{L}^{\rm GW}(z)=d_{L}^{\rm
EM}(z)\exp\left\\{-\int^{z}_{0}\frac{dz^{\prime}}{1+z^{\prime}}\,\delta(z^{\prime})\right\\}.$
(2.6)

The exact form of $\delta(z)$ depends on the specific gravity model under
consideration (some subtleties related to this parameterisation are discussed
in [41]).We can also choose to parameterise the entire term in curly brackets
above in some sensible manner, which we demonstrate below. In the next section
we present the analysis of three models/parameterizations implemented in
gwcosmo.

###  2.1 $(\Xi_{0},n)$ parameterization

As mentioned above, there exist many gravity theories that can introduce a
non-standard friction term in the GW propagation equation. Under these
circumstances, rather than work with a large number of individual models, it
can be efficient to instead work with a general parameterized form of the
ratio $d_{L}^{\rm GW}(z)/d_{L}^{\rm EM}(z)$. Constraints on parameters
belonging to specific gravity models are then obtainable via a mapping from
the parameterised constraint.

A widely-adopted and simple parameterization for the ratio of $d_{L}^{\rm
GW}(z)/d_{L}^{\rm EM}(z)$ is given by [36]

$\frac{d_{L}^{\rm GW}(z)}{d_{L}^{\rm
EM}(z)}=\Xi_{0}+\frac{1-\Xi_{0}}{(1+z)^{n}},$ (2.7)

where $\Xi_{0}$ and $n$ are the parameters controlling the variation from GR.
This parameterization was originally obtained by fitting to a class of non-
local gravity models in [36, 42, 43], however, its form can be motivated on
more general grounds as we describe below momentarily. As a result, this
parameterization has been widely adopted in testing GR with dark sirens, for
example in [25, 24, 27, 44, 45], and with other techniques such as using BNS
mass distribution [46] and strongly lensed GW events [47]. Some mappings of
$\Xi_{0}$ and $n$ to parameters in exact gravity models are listed in Table 1
of [38].

The GR limit of Eq.(2.7) corresponds to $\Xi_{0}=1$, so that the GW distance
is the same as the EM luminosity distance. This limit is also recovered for
$z\sim 0$, largely irrespective of $\Xi_{0}$. Conversely, when $z\gg 1$ the
ratio $d_{L}^{\rm GW}(z)/d_{L}^{\rm EM}(z)$ approaches a constant set by the
value of $\Xi_{0}$. This is motivated by the hypothesis that any deviation
from GR predominantly affects the late universe where dark energy dominates,
so $d_{L}^{\rm GW}(z)/d_{L}^{\rm EM}(z)$ deviates from 1 as redshift increases
under effects beyond GR, and approaches constant at high redshift since the
deviation from GR is suppressed in early universe.

The value of $n$ determines how sharply the transition of the ratio from 1 to
$\Xi_{0}$ take place as redshift increases. In the low redshift limit $z\ll
1$, the ratio can be approximated at the linear order in $z$:

$\frac{d_{L}^{\rm GW}(z)}{d_{L}^{\rm EM}(z)}=1-zn(1-\Xi_{0})+O(z^{2}),$ (2.8)

which is often useful in the regime of the local universe. When $n$ is fixed
as a parameter, its fiducial value is often taken to be $n=1.91$, the value
predicted by the RT non-local gravity model when measuring $\Xi_{0}$ [25].
This may seem a somewhat arbitrary choice, but it will not significantly
impact the results of this work (typically we will allow $n$ to vary as a free
parameter with broad priors anyway). Note that GR is also recovered as $n$
approaches 0, meaning we can potentially expect to see some degeneracy between
constraints on $\Xi_{0}$ and $n$, if the data are consistent with GR.

In Fig. 1 we plot the ratio of $d_{L}^{\rm GW}(z)/d_{L}^{\rm EM}(z)$ for some
different values of $\Xi_{0}$. The ratio is greater than 1 when $\Xi_{0}>1$
and lower than 1 when $\Xi_{0}<1$, for $n\neq 0$.

###  2.2 Horndeski class parameterization

In addition to the generic $(\Xi_{0},n)$ parameterization, we also implement a
specific parameterization originating from the Horndeski family of extensions
to GR. The Horndeski action describes the most generic scalar-tensor theory of
gravity with a second-order equation of motion [48, 49]. The observation of
the luminal speed of gravitational waves by GW170817 and its EM counterpart
[6] strongly constrains particular terms in the original Horndeski action,
though see the caveats discussed in [10]. Setting these terms to zero, the
remaining reduced Horndeski action is [14, 50, 51, 52]

$S=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}M^{2}_{\rm
eff}(\phi)R+K(\phi,X)-G_{3}(\phi,X)\Box\phi\right]+S_{m}(g_{\mu\nu},\psi_{m}),$
(2.9)

where $\phi$ is the scalar degree of freedom, $M_{\rm eff}$ is the effective
Planck mass which evolves with $\phi$, $X$ is a kinetic term defined by
$X\equiv-\bigtriangledown^{\mu}\phi\bigtriangledown_{\mu}\phi/2$, $S_{m}$ is
the matter action, and $\psi_{m}$ are matter fields minimally coupled to the
metric $g_{\mu\nu}$. When $M_{\rm eff}$ equals the Planck mass $M_{\rm P}$,
and $K=G_{3}=0$, the action reduces to GR. It has been shown in the literature
(e.g. [53, 54, 39]) that the GW propagation equation derived from the reduced
Horndeski action is given by

$\tilde{h}^{\prime\prime}+{\cal
H}[2+\alpha_{M}(z)]\tilde{h}^{\prime}+c^{2}k^{2}\tilde{h}=0,$ (2.10)

where $\alpha_{M}$ is the running rate of the effective Planck mass defined by

$\alpha_{M}(z)\equiv\frac{d\ln(M_{\rm eff}/M_{\rm P})^{2}}{d\ln a}.$ (2.11)

We see from Eq. (2.1) that the GW friction term $\delta(z)$ is related to
$\alpha_{M}(z)$ by $\delta(z)=-\alpha_{M}(z)/2$. A widely-used and simple
parameterization for $\alpha_{M}(z)$ is that it is proportional to the
fractional dark energy density [55, 56, 57, 58] as

$\alpha_{M}(z)=c_{M}\frac{\Omega_{\Lambda}(z)}{\Omega_{\Lambda,0}}=c_{M}\frac{1}{\Omega_{m,0}(1+z)^{3}+\Omega_{\Lambda,0}},$
(2.12)

where $c_{M}$ is the free parameter to be constrained, and
$\Omega_{\Lambda}(z)$ is the fractional dark energy density. This assumption
associates the deviation from GR to the growth of the dark energy density
parameter. On one hand, this seems a plausible assumption if a cosmological
modified gravity theory fulfils its principle motivation to replace dark
energy. On the other, it has been argued that Eq. (2.12) is not representative
of the true evolution of mainstream gravity theories [59]. We will use Eq.
(2.12) as a general agnostic expression; a different parameterization can be
straightforwardly replaced in our calculations. For instance another way of
parameterization can be found in [60].

With the relation $\delta(z)=-\alpha_{M}(z)/2$, one obtains by substituting
the above parameterization into Eq. (2.6) that

$d_{L}^{\rm GW}(z)=d_{L}^{\rm
EM}(z)\exp\left\\{\frac{c_{M}}{2\Omega_{\Lambda,0}}\ln\frac{1+z}{[\Omega_{m,0}(1+z)^{3}+\Omega_{\Lambda,0}]^{1/3}}\right\\}.$
(2.13)

This parameterization only depends on one free parameter, $c_{M}$. GR is
recovered when $c_{M}=0$. We have $d_{L}^{\rm GW}(z)/d_{L}^{\rm EM}(z)>1$ when
$c_{M}>0$, and $d_{L}^{\rm GW}(z)/d_{L}^{\rm EM}(z)<1$ when $c_{M}<0$, as
shown in Fig. 1.

###  2.3 Extra dimensions

Apart from scalar-tensor theories, additional spacetime dimensions are
introduced in some extended theories of gravity, which can also leave an
imprint on GW propagation. For instance, in Dvali-Gabadadze-Porrati (DGP)
theory [61, 62], gravitational forces propagate in higher dimensions beyond
4-dimension spacetime, while ordinary matter is confined on the 3-dimension
spatial brane. This leads to an effective ‘leakage’ of the perceived
gravitational forces over cosmological distance scales, which transforms into
an extra damping term in the GW propagation equation. It is shown in [63] that
the GW and EM luminosity distances in some theories with extra dimensions are
related by

$d_{L}^{\rm GW}(z)=d_{L}^{\rm EM}(z)\left\\{1+\left[\frac{d_{L}^{\rm
EM}(z)}{(1+z)R_{c}}\right]^{n_{D}}\right\\}^{\frac{D-4}{2n_{D}}},$ (2.14)

where $D$ is the number of spacetime dimensions, $R_{c}$ is the comoving scale
which controls the transitions from the GR 4-dimensional regime at small
scales to the higher-dimensional regime on large scales, and $n_{D}$ controls
the sharpness of the transition around $R_{c}$. On small scales where
$d_{L}^{\rm EM}\ll R_{c}$, GR is recovered and hence $d_{L}^{\rm GW}$ is
indistinguishable from $d_{L}^{\rm EM}$. On large scales where $d_{L}^{\rm
EM}\gg R_{c}$, the ratio of $d_{L}^{\rm GW}/d_{L}^{\rm EM}$ grows with the
comoving distance to the power of $(D-4)/2$. Constraints on $D$ with bright
and dark sirens have been studied in [63, 24].

The comparison between $d_{L}^{\rm GW}$ and $d_{L}^{\rm EM}$ for the three
modified gravity parameterizations with a selection of typical values is shown
in Fig. 1. We can see that in all of the three models, $d_{L}^{\rm
GW}(z)/d_{L}^{\rm EM}(z)$ can either be larger or smaller than 1, depending on
the parameter value. It starts from 1 at low redshift and deviates from 1 as
redshift increases. Therefore the deviation from GR is more significant at
high redshift, so at the simplest level, GW events further away from may be
more helpful in testing these gravity models. As we will see in sections 4 and
5, the reality of this statement will be complicated by the sky localisations
and distance errors of each event.

Figure 1: GW luminosity distances in the three modified gravity models
analysed in this paper compared to GR. We fix $n=1.91$ in the $(\Xi_{0},n)$
parameterization, and $R_{c}=100$ Mpc in the extra-dimension model.

##  3 Modifications to gwcosmo

###  3.1 Bayesian framework of gwcosmo

The detailed Bayesian statistical framework of the gwcosmo code is presented
in [23, 19, 26], which we will briefly review in this section. Then we will
show how the posterior probabilities of modified gravity parameters are built.
They are included in a set of hyper-parameters $\Lambda$ to be measured with
the bright siren and the dark siren methods. These hyper-parameters also
include cosmological parameters like $H_{0}$, and parameters describing both
the compact object population distribution and the merger rate redshift
evolution, which will be discussed in section 3.2 .

Given $N_{\text{det}}$ detected GW events with observation data
$\\{x_{\text{GW}}\\}$, the posterior for a set of hyper-parameters $\Lambda$
is given by [64, 65]

$\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}p(x_{\text{GW}i}|D_{\text{GW}i},\Lambda,I),$
(3.1)

where $D_{\text{GW}}$ is a parameter representing whether a GW event is
detected or not, which takes the value of 1 or 0. $I$ represents any other
information or assumptions not explicitly contained in the parameter set
$\Lambda$. The detection criteria of an event is that its SNR surpasses the
SNR threshold we choose. The selected SNR threshold in our analysis will be
discussed in the next two sections. The prior $p(\Lambda|I)$ is uniform in the
analysis for all parameters in $\Lambda$. $p(N_{\text{det}}|\Lambda,I)$ is the
probability of detecting $N_{\text{det}}$ events. It depends on the intrinsic
astrophysical rate of GW events $R$. As discussed in [66, 64], by choosing a
flat prior on $R$ in a log scale, we obtain the combined prior
$p(\Lambda,R|I)\propto 1/R$, so that the dependence on $R$ is dropped out when
marginalizing $p(N_{\text{det}}|\Lambda,I)$ over $R$. Beyond this
$p(N_{\text{det}}|\Lambda,I)$ has no further dependence on hyper-parameters.
With Bayes’ theorem, the likelihood for obtaining a set of GW observations
$x_{\text{GW}}$ can be expanded into

$\displaystyle\prod^{N_{\text{det}}}_{i}p(x_{\text{GW}i}|D_{\text{GW}i},\Lambda,I)$
$\displaystyle\propto\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}$ (3.2)
$\displaystyle\propto{\bigg{[}\prod^{N_{\text{det}}}_{i}\int
p(x_{\text{GW}i}|\theta,\Lambda,I)p(\theta|\Lambda,I)d\theta\bigg{]}}{\bigg{[}\int
p(D_{\text{GW}}|\theta,\Lambda,I)p(\theta|\Lambda,I)d\theta\bigg{]}^{-N_{\text{det}}}},$

where $\theta$ represents a parameter set describing a detected GW signal,
which must be marginalised over in order to find posterior for $\Lambda$. The
parameters in $\theta$ of interest to us in this work include the detector
frame masses $m_{1}^{\rm det}$ and $m_{2}^{\rm det}$, redshift $z$, and sky
location $\Omega$. $p(x_{\text{GW}i}|\theta,\Lambda,I)$ in the numerator is
obtained from the posterior samples of GW signal parameters estimated for each
detected event, assuming a cosmological model. Here only cosmological hyper-
parameters are relevant, while population and merger rate hyper-parameters
enter in $p(\theta|\Lambda,I)$. $p(D_{\text{GW}i}|\theta,\Lambda,I)$ in the
denominator is called the detection probability of GW events given a cosmology
model. In gwcosmo it is computed with a large number of injections of GW
events. The injections are simulated detections of GW signals that are
generated with the GW waveform injected with noises. The injections we used
are generated with the IMRPhenomPv2 waveform over the sensitivity curves of
specific detector network configurations such as O2, O3 or O4 set ups.
Provided the selected SNR threshold, injections meeting the SNR threshold are
used for computing the detection probability.

The posterior computed above is affected by the selection effects, which are
the effects caused by the sensitivity of GW detectors, our choices of the
population and cosmology model, prior ranges of $\theta$ and the SNR threshold
in computing detection probability. If the prior range of parameters cannot
cover the space that may be detected from GW events, the posterior will be
biased. However, it is not always more beneficial for the prior range to be
larger, because some parameters may be degenerate, so that one of the
unconstrained parameters may biase constraints on another one (we will see an
example of this in Section 5). In addition, the SNR threshold determines
whether an event is counted as detected. It must be consistent with the SNR
threshold of LVK detection, so it can’t be too low. On the other hand, if the
threshold is too high, we will have too few events to obtain a tighter
constraint. Therefore the prior range and the SNR threshold need to be
selected carefully.

The probabilistic distributions of the redshift and sky location of the GW
event are marginalised over in different ways for the bright siren and the
dark siren methods. For a bright siren, a prior on the redshift and the sky
location can be retrieved from the host galaxy identified from an EM
counterpart with high accuracy. In gwcosmo the redshift prior of a bright
siren $p(z|\Lambda,I)$ is approximated by a Gaussian distribution with the
mean and the standard deviation being the same as the measured redshift of the
host galaxy and its redshift uncertainty. (This in principle could be
generalised to non-Gaussian redshift measurements in future.)

On the other hand, the prior on redshift and the sky location in the dark
siren method is obtained with galaxy catalogue information. Without an EM
counterpart, a GW event is localized within a patch of the sky that is too
large to identify its host galaxy. However, galaxies recorded in catalogues
within this patch of the sky are potential host galaxies that can provide
redshift information of the GW source. The redshift prior of the GW source is
then obtained by combining the redshift distribution of galaxies in the
catalogues within the sky localization area of the event upon the prior of the
comoving redshift volume.

In gwcosmo the redshift prior is computed for each pixel on the sky [19]. The
resolution in dividing the sky is determined by a parameter called nside
introduced by healpy [67, 68]. For larger nside, there are more pixels in the
whole sky and correspondingly fewer galaxies in each pixel; this allows for
more finely-grained features in across different patches in the sky, but takes
takes longer to evaluate numerically. For a typical value of nside=64, the sky
is divided into 49,152 pixels. The pixelation of the line-of-sight (LOS)
redshift prior reduces the bias caused by variations in galaxy catalogue
incompleteness – the completeness is different for sky areas explored by
different galaxy surveys that make up the GLADE galaxy catalogue [69].
Assuming a uniform completeness across the sky will overestimate or
underestimate the completeness of some areas, creating bias in GW host
identification. Therefore using galaxy catalogue with specific completeness in
each pixel, which is defined according to the apparent magnitude threshold of
the survey, will reduce such bias.

The final form of the posterior with the dark siren approach is given by

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

where $p(z|\Omega_{j},\Lambda,I)$ is the LOS redshift prior for each pixel at
a sky location $\Omega_{j}$, which can be pre-computed before estimating the
posterior, of which the details are presented in [26]. In this paper, the
GLADE+ galaxy catalogue in the $K$ band is used for the analysis. Although the
$K$ band is less complete than the $B_{J}$ band, its incompleteness model is
better understood than the $B_{J}$ band, and hence gives more reliable
results. Unlike the $K$ band, the $B_{J}$ band suffers from an issue that its
incompleteness is not well-captured by a simple magnitude threshold.

In the GR version of gwcosmo, $H_{0}$ is the only cosmological parameter in
$\Lambda$. Here we extend gwcosmo to the estimation of sets of modified
gravity parameters such as $\Xi_{0}$, $n$ in Eq. (2.7), $c_{M}$ in Eq. (2.13),
or $D$ in Eq. (2.14). Since these parameters modify $d_{L}^{\rm GW}(z)$, but
leave $d_{L}^{\rm EM}(z)$ unchanged, the LOS redshift prior from the galaxy
catalogue is not affected. However, modification is required when
marginalising over the source frame masses of compact object binaries, because
$z$ inferred from $d_{L}^{\rm GW}(z)$ now differs from that from $d_{L}^{\rm
EM}(z)$ under different values of modified gravity parameters. Inside a pixel,
the numerator integrand of the likelihood is given by

$p(x_{\text{GW}i}|z,\Lambda,I)=\int_{m_{\rm min}}^{m_{\rm max}}\int_{m_{\rm
min}}^{m_{1}^{\rm s}}p(x_{\text{GW}i}|z,m_{1}^{\rm s},m_{2}^{\rm
s},\Lambda,I)p(m_{1}^{\rm s},m_{2}^{\rm s}|\Lambda,I)dm_{2}^{\rm s}dm_{1}^{\rm
s},$ (3.4)

where $m_{1}^{\rm s}$ and $m_{2}^{\rm s}$ are the component masses in the GW
source frame, with the assumption $m_{1}^{\rm s}>m_{2}^{\rm s}$. $m_{\rm min}$
and $m_{\rm max}$ are the minimum and the maximum mass for a source of that
type, given by the population model. The prior on source masses $p(m_{1}^{\rm
s},m_{2}^{\rm s}|I)$ is inferred from the compact object mass distribution
models [70]. However, the posterior samples from GW observation data are
generated in the detector frame. In addition, it is $d_{L}^{\rm GW}(z)$ that
is measured from the GW signal instead of $z$. Therefore we have

$p(x_{\text{GW}i}|z,\Lambda,I)=\int_{m_{\rm min}}^{m_{\rm max}}\int_{m_{\rm
min}}^{m_{1}^{\rm s}}p(x_{\text{GW}i}|d_{L}^{\rm GW}(z,\Lambda),m_{1}^{\rm
det}[m_{1}^{\rm s},z],m_{2}^{\rm det}[m_{2}^{\rm s},z],I)p(m_{1}^{\rm
s},m_{2}^{\rm s}|I)dm_{2}^{\rm s}dm_{1}^{\rm s},$ (3.5)

Applying Bayes’ theorem we obtain

$p(x_{\text{GW}i}|d_{L}^{\rm GW},m_{1}^{\rm det},m_{2}^{\rm
det},I)=\frac{p(d_{L}^{\rm GW},m_{1}^{\rm det},m_{2}^{\rm
det}|x_{\text{GW}i},I)p(x_{\text{GW}i}|I)}{\pi(d_{L}^{\rm GW},m_{1}^{\rm
det},m_{2}^{\rm det}|I)}\propto\frac{p(d_{L}^{\rm GW},m_{1}^{\rm
det},m_{2}^{\rm det}|x_{\text{GW}i},I)}{\pi(d_{L}^{\rm GW},m_{1}^{\rm
det},m_{2}^{\rm det}|I)},$ (3.6)

where $\pi(d_{L}^{\rm GW},m_{1}^{\rm det},m_{2}^{\rm det}|I)$ is the prior
applied to $d_{L}^{\rm GW},m_{1}^{\rm det}$ and $m_{2}^{\rm det}$ when
computing the posterior samples of a GW event. The posterior samples are
effectively the numerator of Eq. (3.6). For the GWTC-3 data set, the posterior
samples are generated with two options: $\pi(d_{L}^{\rm GW})\propto(d_{L}^{\rm
GW})^{2}$ or $\pi(d_{L}^{\rm GW})$ uniform in comoving volume prior. In this
paper we used the data set with $\pi(d_{L}^{\rm GW})\propto(d_{L}^{\rm
GW})^{2}$. Meanwhile $\pi(m_{1}^{\rm det},m_{2}^{\rm det})$ is uniform. Then
we convert the posterior samples of $d_{L}^{\rm GW}$ into those of $z$ by

$p(d_{L}^{\rm GW}(z,\Lambda),m_{1}^{\rm det},m_{2}^{\rm
det}|x_{\text{GW}i},I)=p(z,m_{1}^{\rm det},m_{2}^{\rm
det}|x_{\text{GW}i},\Lambda,I)\bigg{|}\frac{\partial z}{\partial d_{L}^{\rm
GW}}\bigg{|}.$ (3.7)

The conversion from $d_{L}^{\rm GW}$ to $z$, and the Jacobian $|{\partial
z}/{\partial d_{L}^{\rm GW}}|$ for posterior samples needs to be modified
under different parameterizations, which is the main modification compared to
the GR version of gwcosmo. Finally, the transformation from source masses into
detected masses gives

$dm_{1}^{\rm s}dm_{2}^{\rm s}=\bigg{|}\frac{\partial(m_{1}^{\rm s},m_{2}^{\rm
s})}{\partial(m_{1}^{\rm det},m_{2}^{\rm det})}\bigg{|}dm_{1}^{\rm
det}dm_{2}^{\rm det}=\frac{1}{(1+z)^{2}}dm_{1}^{\rm det}dm_{2}^{\rm det}$
(3.8)

Following the steps above, the expression for the integrand in Eq. (3.5)
becomes (remember this is part of the numerator in Eq. (3.3)):

$\displaystyle p(x_{\text{GW}i}|z,\Lambda,I)\propto\int_{m^{\rm det}(m_{\rm
min},z)}^{m^{\rm det}(m_{\rm max},z)}\int_{m^{\rm det}(m_{\rm min},z)}^{m^{\rm
det}(m_{1}^{s},z)}$ $\displaystyle\frac{p(z,m_{1}^{\rm det},m_{2}^{\rm
det}|x_{\text{GW}i},\Lambda,I)}{(d_{L}^{\rm GW})^{2}}\bigg{|}\frac{\partial
z}{\partial d_{L}^{\rm GW}}\bigg{|}$ (3.9)
$\displaystyle\times\frac{p(m_{1}^{\rm s}[m_{1}^{\rm det},z],m_{2}^{\rm
s}[m_{2}^{\rm det},z]|I)}{(1+z)^{2}}dm_{2}^{\rm det}dm_{1}^{\rm det}.$

The posterior samples on $z$ computed from $d_{L}^{\rm GW}$ given $\Lambda$,
$m_{1}^{\rm det}$ and $m_{2}^{\rm det}$ for a GW event can be approximated as
a sum over delta-functions:

$p(z,m_{1}^{\rm det},m_{2}^{\rm
det}|x_{\text{GW}i},\Lambda,I)\approx\frac{1}{N_{\rm pos}}\sum^{N_{\rm
pos}}_{k=1}\delta(z-z_{k}(d_{L,k}^{\rm GW},\Lambda))\delta(m_{1}^{\rm
det}-m_{1,k}^{\rm det})\delta(m_{2}^{\rm det}-m_{2,k}^{\rm det}),$ (3.10)

where $k$ denotes the $k^{\rm th}$ posterior sample and $N_{\rm pos}$ denotes
the total number of posterior samples. Inserting this back to Eq. (3.9) we
obtain the result as

$p(x_{\text{GW}i}|z,\Lambda,I)\propto\frac{1}{N_{\rm pos}}\sum^{N_{\rm
pos}}_{k=1}p(m_{1}^{\rm s}(m_{1,k}^{\rm det},z),m_{2}^{\rm s}(m_{2,k}^{\rm
det},z)|I)\bigg{|}\frac{\partial z}{\partial d_{L}^{\rm
GW}}\bigg{|}\frac{1}{(d_{L}^{\rm GW})^{2}(1+z)^{2}}\delta(z-z_{k}(d_{L,k}^{\rm
GW},\Lambda)).$ (3.11)

The above quantity is then multiplied by the LOS redshift prior333In reality,
a kernel density estimate (KDE) of $p(x_{\text{GW}i}|z,\Lambda,I)$ is created
first, and multiplied with the redshift prior. $p(z|\Omega_{j},\Lambda,I)$,
and integrated over redshift. Summing the results over all pixels one obtains
the numerator of the likelihood as indicated in Eq. (3.3). The final result of
the numerator is given by

$\displaystyle p(x_{\text{GW}i}|\Lambda,I)\propto\frac{1}{N_{\rm
pos}}\sum^{N_{\rm pos}}_{k=1}p(m_{1}^{\rm s}(m_{1,k}^{\rm
det},z_{k}),m_{2}^{\rm s}(m_{2,k}^{\rm det},z_{k})|I)\bigg{|}\frac{\partial
z}{\partial d_{L}^{\rm GW}}$
$\displaystyle\bigg{|}_{z=z_{k}}\frac{1}{(d_{L,k}^{\rm GW})^{2}(1+z_{k})^{2}}$
$\displaystyle\times\sum^{N_{\text{pix}}}_{j}p(z_{k}|\Omega_{j},\Lambda,I).$
(3.12)

The details of computing the LOS redshift prior
$p(z_{k}|\Omega_{j},\Lambda,I)$ are described in [26]. Similarly, injections
in the detector frame need to be converted into the source frame when
marginalising the detection probability in the denominator of the likelihood
over redshift and component masses. The detection probability is given in the
same form by

$\displaystyle p(D_{\text{GW}}|\Lambda,I)\propto\frac{1}{N_{\rm
inj}}\sum^{N_{\rm inj}}_{k=1}p(m_{1}^{\rm s}(m_{1,{\rm inj},k}^{\rm
det},z_{k}),m_{2}^{\rm s}(m_{2,{\rm inj},k}^{\rm det},z_{k})|I)$
$\displaystyle\bigg{|}\frac{\partial z}{\partial d_{L}^{\rm
GW}}\bigg{|}_{z=z_{{\rm inj},k}}\frac{1}{(d_{L,{\rm inj},k}^{\rm
GW})^{2}(1+z_{{\rm inj},k})^{2}}$
$\displaystyle\times\sum^{N_{\text{pix}}}_{j}p(z_{{\rm
inj},k}|\Omega_{j},\Lambda,I).$ (3.13)

By inserting the results of the numerator and the denominator back into Eq.
(3.3), and multiplying over all the GW events, we can thus obtain the
posterior on the parameter set $\Lambda$.

###  3.2 Prior on source masses and merger rate evolution

The source mass prior $p(m_{1}^{\rm s},m_{2}^{\rm s}|I)$ is obtained from the
mass distribution model. There are several different BBH mass distribution
models adopted in the LVK analysis, such as the Broken Power Law model and the
Power Law + Peak model [31, 70]. In our analysis, we use the Power Law + Peak
model for black hole distribution, because it is preferred over other models
from the population analysis of O3 data [71], and is adopted by a number of
similar works [24, 27, 26]. We use a uniform distribution for neutron stars as
in the LVK GWTC-3 cosmology paper [31]. There are 8 parameters in the
description of the Power Law + Peak model. In our analysis of GWTC-3 events in
Section 5, we will study cases where these parameters are fixed or varied. In
the cases where they are varied, the priors on them are uniformly distributed,
and are listed in Table 1 in Section 5.

The explicit redshift prior actually depends on the presence of GW sources,
which we denote by $s$. By Bayes’ theorem it is given that

$p(z|s,\Omega_{j},\Lambda,I)=\frac{p(s|z,\Omega_{j},\Lambda,I)p(z|\Omega_{j},\Lambda,I)}{p(s|\Omega_{j},\Lambda,I)}.$
(3.14)

The term $p(s|\Omega_{j},\Lambda,I)$ is cancelled for appearing in both the
numerator and the denominator of Eq. (3.3). Assuming that GW sources are
distributed uniformly in the sky on the scale of our pixels, the dependence on
$\Omega_{j}$ can be dropped. Then $p(s|z)$ depends on the redshift evolution
of the merger rate of compact object binaries. Following the LVK cosmology
paper [31], we adopt a parameterization motivated by connection of GW events
to the star formation rate in [72], which has the form

$p(s|z,\gamma,k,z_{p})=[1+(1+z_{p})^{-\gamma-k}]\frac{(1+z)^{\gamma}}{1+[(1+z)/(1+z_{p})]^{\gamma+k}}.$
(3.15)

The free parameters $\gamma$, $k$ and $z_{p}$ control the merger rate redshift
evolution, and are included in the hyper-parameters to be measured. The merger
rate evolution can be scaled by an overall factor $R_{0}$, but it is not
included in the current gwcosmo analysis.

##  4 Tests on Mock Data

###  4.1 Simulation setup

Because of the limited number of detected GW events at the present time, we
first use mock GW events to test the measurement of modified gravity
parameters by gwcosmo. We use the First Two Years of Electromagnetic Follow-Up
with Advanced LIGO and Virgo (F2Y) mock data set [73] to test our bright siren
methodology. This catalogue was generated by injecting around 50,000 BNS
mergers with the LIGO and Virgo sensitivity curves mimicking the O2 run. The
criteria for a detected bright siren is having single detector SNR over 4, and
a network SNR over 12. As a result only approximately 500 BNS events are
marked as detected events, among which 250 of them are randomly selected for
parameter estimation. We used these 250 bright sirens for our test of gravity
with gwcosmo; we note this is the same mock data challenge as employed for
$H_{0}$ in [23].

We generated our own BNS injections for computing the detection probability.
We used the same O2 LIGO and Virgo sensitivity curves as in [73], which are
the “early” and “mid” sensitivity curves for L1 and H1 in LIGO technical
document [74], and the O2 noise curve for V1 from [75]. Each of the three
detectors has a duty factor of 76.8% as indicated in Fig. 9 of [73]. We
generated 50,000 injections with the SNR threshold of 12, GW luminosity
distance between 0 and 500 Mpc, and the BNS component mass between $1.2$ and
$1.6\,M_{\odot}$, assuming $m_{1}>m_{2}$.

A significant feature of the F2Y dataset is that it is generated with a local
universe approximation, given that the redshifts of detectable BNS events are
expected to be very low. The linear cosmology approximation $d_{L}^{\rm
EM}\approx cz/H_{0}$ is applied for the EM counterparts of BNS injections. The
BNS $d_{L}^{\rm GW}$ is then parameterized with the linearized $d_{L}^{\rm
EM}$ in the three modified gravity models introduced in Section 2. We note
that in the case of the Hubble constant analysis, the detection probability
under a low redshift approximation simply grows as $H_{0}^{3}$ [76]; in the
case of modified gravity it acquires additional corrections. On the other
hand, the detected BNS masses have no redshift correction since this is
negligible in the nearby universe. The mass prior we used in the analysis is
uniform in $[1.2,1.6]~{}M_{\odot}$, since the F2Y BNSs were generated with
component masses between 1.2 and 1.6 $M_{\odot}$. Such a mass prior is much
narrower than the range of $[1.0,3.0]~{}M_{\odot}$ used in later LVK analyses:
neutron stars were generally expected to have masses strongly peaked around
$1.4~{}M_{\odot}$ before much heavier neutron stars were detected in events
from the O3 observing run.

###  4.2 Results

We first measure $\Xi_{0}$ in the $(\Xi_{0},n)$ parameterization shown in Eq.
(2.7) with a uniform prior in $[0.05,10]$, while fixing other cosmological
parameters as $H_{0}=70$ ${\rm km}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}$ and
$n=1.91$444Note that other cosmological parameters such as $\Omega_{M0}$ will
not be varied in any of our analyses: they are well-measured by EM probes
already, so we adopt their central values.. The normalized posterior of
$\Xi_{0}$ for each bright siren is plotted in black, and the combined
posterior is plotted in red in Fig. 2. The measured value with $1\sigma$ bound
is $\Xi_{0}=0.96^{+0.19}_{-0.19}$, which is consistent with GR. Our results
show that gwcosmo is capable of constraining $\Xi_{0}$ well in the simplified
case of a one-parameter analysis with a number of bright sirens. The number of
bright sirens needed to obtain this bound (250) may seem quite large; this is
because, as discussed in section 2.1, $\Xi_{0}$ is most relevant to high-
redshift GW events, and the mock bright sirens used here are all at low
redshifts.

Figure 2: The normalized posterior probability distribution for $\Xi_{0}$ for
individual events (black) and the combined one (red). The blue vertical line
labels $\Xi_{0}=1$ recovery of GR.

However, the 1D measurement of $\Xi_{0}$ can be misleading, because it has
assumed $H_{0}$ to be a fixed value. Since $H_{0}$ and $\Xi_{0}$ both
determine the conversion from GW distances to redshifts, they are expected to
be degenerate. Ideally it is more sensible to perform a 2D $H_{0}$-$\Xi_{0}$
joint analysis. However, the number of posterior samples for the F2Y mock BNS
is $\sim 1000$, which is very low compared to the real events that have tens
of thousands of posterior samples. Therefore the estimation of redshift for
F2Y data is less accurate. Our joint analysis shows that the GR value
$(H_{0}=70,\Xi_{0}=1)$ is on the edge of the $3\sigma$ bound of the 2D joint
posterior. Due to the low number of posterior samples for the F2Y data, the 2D
result is unreliable.

We further perform the same analysis on $c_{M}$ in the Horndeski class
parameterization shown in Eq.(2.13) and $D$ in the extra-dimension
parameterization shown in Eq.(2.14). We fix $H_{0}$, mass distribution and
merger rate redshift evolution parameters as in the $\Xi_{0}$ measurement. In
the measurement of $c_{M}$, the values of $\Omega_{m,0}$ and
$\Omega_{\Lambda,0}$ are needed to reconstruct the evolution of the fraction
of dark energy density as described in Eq. (2.12). We adopt
$\Omega_{m,0}=0.3065$ from the Planck result [77] as in the LVK cosmology
paper [31], and $\Omega_{\Lambda,0}=1-\Omega_{m,0}$ as other components of
energy density are negligible today. We apply the uniform prior for $c_{M}$ in
[-5, 10]. The normalized posterior of $c_{M}$ for individual events and the
combined posterior are plotted in the left panel of Fig. 3. The $1\sigma$
bound measured value is $c_{M}=-0.18^{+0.71}_{-0.72}$, which is consistent
with GR where $c_{M}=0$.

In the measurement of $D$ in the extra-dimension parameterization shown in Eq.
(2.14), we choose to fix $R_{c}=100$ Mpc, which is the same magnitude as that
of the luminosity distances for BNS events, so that the modified gravity
effect can be probed at such distances. This value may not match the one
motivated from fundamental physics, but it allows us to test the constraint on
$D$ at a scale of the BNS mock data we have available. In addition we also fix
$n_{D}=1$. The prior of $D$ is uniform in [3.5, 5.0]. The effect of varying
$D$ with fixed $R_{c}$ and $n_{D}$ is shown in Fig. 1. The normalized
posterior of $D$ for individual events and the combined posterior are plotted
in the right panel of Fig. 3. The $1\sigma$ bound measured value is
$D=4.015^{+0.027}_{-0.027}$, which is consistent with the GR limit of this
class of theories.

Figure 3: The normalized posterior probability distribution for $c_{M}$ (left
panel) and $D$ (right panel) for individual events (dark) and the combined
constraint (red). The blue vertical lines indicate the GR limits of $c_{M}=0$
and $D=4$.

In summary, we obtain good 1D constraints on modified gravity parameters in
all of the three alternative gravity models using 250 mock bright sirens, and
recover consistency with GR in all cases (as of course expected, as the mock
data is generated under GR). However, because of the low number of posterior
samples of the F2Y data, 2D joint constraints of $H_{0}$ and modified gravity
parameters are not reliable. Still, our test shows that gwcosmo is capable of
measuring deviations beyond GR with the bright siren method. It will be used
to test gravity with bright sirens that are expected to be detected in the O4
run.

##  5 Reanalysing the GWTC-3 data

###  5.1 Parameter ranges

In this section we show the results of reanalysing the GWTC-3 events with
gwcosmo, including modified gravity parameters. We use the same 46 dark sirens
from the GWTC-3 catalogues as in the LVK GWTC-3 cosmology paper [31] with net
SNR larger than 11, including 42 BBH mergers, 1 BNS merger GW190425, and 3
NSBH mergers GW190814, GW200105_162426 and GW200115_042309. The posterior
samples of GWTC-2.1 [78] are used in computing the numerator of the likelihood
for events from O1 to O3a, and GWTC-3 samples [79] for events in O3b. The
denominator of the likelihood is computed using $2\times 10^{6}$ GW injections
with the SNR threshold set to 11. The injected detected masses range from 1 to
500 $M_{\odot}$, and the injected GW distances are varied between 0.1 and
$\sim 20,000$ Mpc. Such a high value for the maximum GW distance is needed to
cover the GW distance for extreme values of the modified gravity parameters
for $z\in(0,10)$.

The LOS redshift prior we used is pre-computed using the GLADE+ galaxy
catalogue in the $K$ band, with resolution of nside = 64. We will focus on the
$\Xi_{0}-n$ parameterization for the next few subsections, but the other
modified gravity parameters will be treated analogously, and their final
results will be discussed in Section 5.3. Given the maximum redshift $z=10$ of
the LOS redshift prior, we need to choose joint prior bounds of $H_{0}$ and
$\Xi_{0}$ such that the maximum injected GW distance corresponds to redshift
not larger than 10. As a result, we vary $\Xi_{0}$ in a range of [0.3,10] when
fixing $H_{0}=70~{}{\rm km}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}$. A combination of
the maximum $H_{0}$ and the minimum $\Xi_{0}$ gives the largest redshift for
the maximum injected GW distance, so we need to reduce the higher bound of
$H_{0}$ or raise the lower bound of $\Xi_{0}$ when varying both of them. We
chose a range of [20,140] ${\rm km}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}$ for
$H_{0}$ and [0.35,10] for $\Xi_{0}$ for the joint analysis.

###  5.2 Dark sirens results

####  5.2.1 1D measurement of $\Xi_{0}$

We first measure the posterior of $\Xi_{0}$ in Eq. (2.7) with a flat prior on
[0.3,10], while fixing the other two cosmological parameters to be
$H_{0}=70~{}{\rm km}~{}{\rm s}^{-1}{\rm Mpc}^{-1}$ and $n=1.91$. We also fix
the parameters in the Power Law + Peak mass prior as in the LVK cosmology
paper [31], which are $\alpha=3.78$, $\beta=0.81$, $m_{\rm min}^{\rm
BH}=4.98M_{\odot}$, $m_{\rm max}^{\rm BH}=112.5M_{\odot}$,
$\delta_{m}=4.8M_{\odot}$, $\mu_{g}=32.27M_{\odot}$,
$\sigma_{g}=3.88M_{\odot}$ and $\lambda_{p}=0.03$. A brief description of the
physical significance of these parameters is given in Table 1. In the cases of
NSBH and BNS, we fix $m_{\rm min}^{\rm NS}=1.0M_{\odot}$ and $m_{\rm max}^{\rm
NS}=3.0M_{\odot}$. The merger rate evolution parameters are also fixed with
the same values $\gamma=4.59$, $k=2.86$ and $z_{p}=2.47$. The normalized
posterior of $\Xi_{0}$ for each individual event is plotted in Fig. 4. We can
see that there are three typical patterns in the distribution of $p(\Xi_{0})$
for different events: favoring low $\Xi_{0}$ such as GW190910_112807, favoring
high $\Xi_{0}$ such as GW191216_213338, and peaking at a certain value of
$\Xi_{0}$ such as GW170809_082821. The reasons for these patterns can be found
by looking at the posterior sample distribution in $z$ converted from
$d_{L}^{\rm GW}$ for different $\Xi_{0}$ for each event, which we will
elaborate on in Appendix B. The measured value from the combined posterior
with $1\sigma$ bound is $\Xi_{0}=1.06^{+0.25}_{-0.18}$, which is consistent
with GR.

Figure 4: The one-dimensional likelihood of $\Xi_{0}$ from individual selected
events used in our analysis. Note that all other parameters are held fixed.
These 1D curves demonstrate the range of possible curve shapes obtained when
galaxy catalogue support is moderately weak; see Appendix B for a detailed
explanation.

For most dark siren events in GWTC-3, the redshift support from the galaxy
catalogue does not have a significant effect on improving the precision of the
$\Xi_{0}$ measurement. The event with the most galaxy support is GW190814,
which is shown in Fig. 5; here we compare the analysis with the galaxy
catalogue and an empty catalogue (no galaxies). $p(\Xi_{0})$ with galaxy
catalogue support peaks at lower $\Xi_{0}$ compared to the flat distribution
without the support, which shows that galaxy catalogues have the potential to
increase the precision of modified gravity constraints with dark sirens. The
galaxy catalogue information has only a weak effect on the current dark siren
analysis; this is because the present localization of GW events has a
relatively low accuracy, so that the information from clustering of galaxies
gets washed out as it is averaged over too many pixels. If the localization is
small enough, our method is sensitive to clustering of galaxies so that the
results will favour some values of $H_{0}$ and $\Xi_{0}$. In addition, the
galaxy catalogue is incomplete, especially at high redshift, so it provides
limited support to GW events further away.

In future, with an increase in the sensitivity of the KAGRA detector, and the
addition of LIGO India to the terrestial network, the localization of GW
events will be improved. Moreover, galaxy surveys by new detectors such as the
Dark Energy Spectroscopic Instrument (DESI) [80] and the Euclid space
telescope [81] in the future will provide millions of further spectroscopic
galaxy redshifts, particuarly at higher redshifts where catalogue support is
currently lacking. This will further accelerate the precision of dark sirens
constraints.

Figure 5: The 1D posterior of $\Xi_{0}$ for GW190814, with both the GLADE+
galaxy catalogue and an empty catalogue analysis (not galaxies. We see that
for this event the clustering in the galaxy catalogue (dis)favours certain
values of $\Xi_{0}$, creating a bumpy structure in the likelihood. Comparable
features can be seen in the $H_{0}$ likelihood shown in Fig. 9 of [26].

####  5.2.2 BBH joint analysis

Next we perform a joint posterior measurement for all hyper-parameters under
the $(\Xi_{0},n)$ MG parameterization with nested sampling in gwcosmo. We
adopt uniform priors for all hyper-parameters including the mass distribution
model and the merger rate redshift evolution model as listed in Table 1. We
present the corner plot of the joint posteriors for 42 BBH events from nested
samplings for 6 of the most interesting parameters in Fig. 6, which are
$\gamma$, $m_{\rm max}^{\rm BH}$, $\mu_{g}$, $H_{0}$, $\Xi_{0}$ and $n$. The
full corner plot with all parameters is shown in Fig. 11 in Appendix A. We can
see that some BBH mass distribution parameters are moderately well-
constrained; However, $k$ and $z_{p}$ in the merger rate redshift evolution
model are poorly constrained, which is similar to results in the LVK cosmology
paper [31]. In addition, the estimated value of $H_{0}$ is
$55.93^{+38.35}_{-26.29}~{}{\rm km}~{}{\rm s}^{-1}{\rm Mpc}^{-1}$, again an
uncertainty similar to previous results. The estimated value for $\Xi_{0}$ is
$1.29^{+1.22}_{-0.67}$, which is consistent with GR; as expected, the result
is less constrained than that in the 1D measurement because of co-variance
with other parameters. Similar results have been obtained using the Icarogw
software [24, 27]. Meanwhile $n$ is measured to be $0.79^{+3.86}_{-0.65}$,
which is lower than the fiducial value of 1.91, but with quite broad
uncertainties. The reason that $p(n)$ favors low values of $n$ is likely that
the $\Xi(z)$ parameterization reduces to GR when $n=0$, so a shift towards low
values of $n$ is an alternative fit to the data (instead of $\Xi_{0}=1$), if
the data are consistent with GR. Moreover, as shown by the 2D $\Xi_{0}$-$n$
posterior contour, $\Xi_{0}$ is better constrained when $n$ is larger. This is
because the redshift-dependent term with the inverse power of $n$ in Eq. (2.7)
becomes less significant when $n$ is large, leaving the constant $\Xi_{0}$ as
the more dominant term. Therefore the constraint is put more directly on
$\Xi_{0}$ when $n$ is large.

Parameter | Definition | Prior
---|---|---
$H_{0}$ | Hubble constant | ${\cal U}(20.0,140.0)$
$\Xi_{0}$ | Modified gravity parameter controlling high-z limit of | ${\cal U}(0.35,10.0)$
| distance ratio in Eq. (2.7) |
$n$ | Modified gravity parameter controlling steepness of | ${\cal U}(0.1,10.0)$
| distance ratio in Eq. (2.7) |
$\alpha$ | The power of the power law component in the primary | ${\cal U}(1.5,8.0)$
| mass distribution |
$\beta$ | The power of the power law component in the mass | ${\cal U}(-4.0,6.0)$
| ratio distribution |
$m_{\rm min}^{\rm BH}~{}[M_{\odot}]$ | The minimum mass of the mass distribution | ${\cal U}(2.0,10.0)$
$m_{\rm max}^{\rm BH}~{}[M_{\odot}]$ | The maximum mass of the mass distribution | ${\cal U}(50.0,200.0)$
$\lambda_{p}$ | Fraction of the model in the Gaussian component | ${\cal U}(0.0,0.5)$
$\mu_{g}$ | Mean of the Gaussian component in the primary mass | ${\cal U}(10.0,50.0)$
| distribution |
$\sigma_{g}$ | Width of the Gaussian component in the primary mass | ${\cal U}(0.1,20.0)$
| distribution |
$\delta_{m}$ | Range of mass tapering at the lower end of the mass | ${\cal U}(0.0,15.0)$
| distribution |
$\gamma$ | The power of the power law distribution of the rate | ${\cal U}(0.0,12.0)$
| evolution before redshift $z_{p}$ |
$k$ | The power of the power law distribution of the rate | ${\cal U}(0.0,10.0)$
| evolution after redshift $z_{p}$ |
$z_{p}$ | The redshift turning point between two power law | ${\cal U}(0.0,10.0)$
| distributions |
Table 1: Prior of the parameters in the Power Law + Peak BBH mass
distribution. Figure 6: Selected corner plot of joint posteriors for 42 BBH
events from the GWTC-3 catalogue computed by the gwcosmo sampling run. The
contours show 1, 2 and 3$\sigma$ bound for the 2D posteriors. The vertical
lines in the 1D histogram plot correspond to the peak values and 1$\sigma$
bound of the parameters in the title. The full corner plot including all
parameters is shown in Fig. 11. A brief description of the mass and merger
rate hyper-parameters is given in Table 1.

From the contour plot we can see that $H_{0}$ and $\Xi_{0}$ are partially
degenerate with each other. This is because $d_{L}^{\rm
GW}\propto\Xi_{0}/H_{0}$ as shown in Eq. (2.2) and (2.7), so that in the
conversion of posterior samples from $d_{L}^{\rm GW}$ to $z$ in computing the
likelihood, varying $H_{0}$ and $\Xi_{0}$ while keeping the ratio
$\Xi_{0}/H_{0}$ unchanged gives the same redshift. So the result that the
estimated value of $H_{0}$ is higher than that in the LVK cosmology paper [31]
is likely caused by $\Xi_{0}$ favoring a slightly higher value than GR.
Similar degeneracy can also be found in Fig. 13 of [24].

Apart from $H_{0}$, we can also see that $\Xi_{0}$ is strongly degenerate with
$\gamma$. In the computation of the likelihood, the LOS redshift prior is
multiplied by the merger rate redshift evolution, which is approximately
proportional to $(1+z)^{\gamma}$ when $z\ll z_{p}$. Given that all GW events
in GWTC-3 have redshift less than 1, and from star formation arguments $z_{p}$
is typically considered to have a value around 2, the power of $\gamma$
dominates the redshift evolution of the merger rate. This can also explain why
$k$, which is the power of the redshift evolution when $z>z_{p}$, and $z_{p}$,
are poorly constrained with GWTC-3 events. So the numerator of the likelihood
is proportional to $\Xi_{0}$ and $(1+z)^{\gamma}$, leading to a strong
degeneracy between $\Xi_{0}$ and $\gamma$. A similar feature is also found in
Fig. 12 of the Icarogw results [24], where the distribution of $p(\gamma)$
with the Power Law + Peak BBH population model is similar to ours. Given that
the distribution of $p(\gamma)$ favors a higher value, if we allow the prior
bound of $\gamma$ to be larger, the distribution of $p(\Xi_{0})$ will migrate
to higher values because of the degeneracy. The prior ranges for $\gamma$, $k$
and $z_{p}$ are chosen to be sufficiently wide in consideration of the effect
of a possible time delay between the formation and the merger of the binary.

####  5.2.3 NSBH joint analysis

We then perform joint analysis on all hyper-parameters with the three NSBH
mergers by nested samplings implemented in gwcosmo. The NSBH events selected
are GW190814, GW200105_162426 and GW200115_042309. Note that GW190814 is not a
confirmed NSBH, but we choose to treat it as one for this analysis. In
addition to all of the parameters in the analysis with BBHs, parameters
$m_{\rm min}^{\rm NS}$ and $m_{\rm max}^{\rm NS}$ for the uniform neutron star
population model are included in the estimation. We use the uniform prior of
$m_{\rm min}^{\rm NS}\in[0.5,1.5]M_{\odot}$ and $m_{\rm max}^{\rm
NS}\in[1.501,5]M_{\odot}$. Due to the small number of events, most of the
parameters in the BBH population model and the merger rate evolution model are
poorly constrained. Furthermore, the marginalized posteriors for some of the
parameters show misleading results, as shown in the selected contours in Fig.
7. For example, the posteriors for $\mu_{g}$ and $\sigma_{g}$ indicate a
narrow Gaussian peak at $\sim 24M_{\odot}$ for the black hole population
model. This is caused alone by the black hole with mass $\sim 23M_{\odot}$ in
GW190814. The same appearance is probably also shown in the constraint of
$m_{\rm max}^{\rm NS}$, which is driven by the assumption that GW190814’s
secondary mass is a neutron star with mass $\sim 2.6M_{\odot}$. Meanwhile the
constraints on cosmological parameters are generally weaker than those with
the BBHs due to the lower number of NSBH events.

Figure 7: Selected corner plot of joint posteriors for 3 NSBH events from the
GWTC-3 catalogue computed by the gwcosmo sampling run. The contours show 1, 2
and 3$\sigma$ bound for the 2D posteriors. The vertical lines in the 1D
histogram plot correspond to the peak values and 1$\sigma$ bound of the
parameters in the title. A brief description of the mass and merger rate
hyper-parameters is given in Table 1.

###  5.3 Combined BBH and NSBH joint analysis

We can improve the constraints on cosmological parameters by combining the
results of BBHs and NSBHs. We believe this is the first time a combined dark
sirens analysis of BBHs and NSBHs has been applied to tests of gravity. Since
the sampling analysis for NSBHs includes two more parameters, $m_{\rm
max}^{\rm NS}$ and $m_{\rm min}^{\rm NS}$, the sampling data from BBHs and
NSBHs cannot be combined simply. We first compute a kernel density estimation
(KDE) function for the multi-dimensional joint posterior of the cosmological
parameters with the sampling data. Then we compute the joint posterior for
each gridded value of the cosmological parameter space with the KDE for BBHs
and NSBHs, and multiply the two posteriors over the space. Finally we obtain
the posterior of each parameter by marginalizing the joint multi-dimensional
posterior over other parameters. The combined posteriors for $H_{0}$,
$\Xi_{0}$ and $n$ are shown in Fig. 8. The estimated values of $H_{0}$,
$\Xi_{0}$ and $n$ for the combined posterior are between the BBH posterior and
the NSBH posterior, with a slightly smaller $1\sigma$ uncertainty than the
individual BBH and NSBH posteriors, which is as expected. Again, we see the
partial degeneracies between $\Xi_{0}$ and $H_{0}$, and $\Xi_{0}$ and $n$ are
maintained.

Figure 8: Posteriors of $H_{0}$, $\Xi_{0}$ and $n$ by combining results of
BBHs and NSBHs from the gwcosmo sampling run.

The results for the other two modified gravity models introduced in this work
for the combined posterior of BBHs and NSBHs are shown in Fig. 9 and 10. For
the Horndeski class model, we apply a uniform prior for $c_{M}$ in [-4,10].
The lower prior bound is restricted by the maximum injected GW distance. The
combined posterior gives $c_{M}=1.5_{-2.1}^{+2.2}$, which is consistent with
the GR prediction of $c_{M}=0$ within a $1\sigma$ bound. Principally due to
the higher redshifts of the dark siren events, this constraint is much
narrower than that obtained by the bright siren GW170817 [82], and is
competitive with that from spectral sirens [83]. We can also see slight
degeneracy between $H_{0}$ and $c_{M}$ in the contour.

On the other hand, for the extra-dimension model, we use a uniform prior for
$D\in[3.7,6]$, $\log(R_{c}/{\rm Mpc})\in[0.5,6]$ and $n_{D}\in[0.3,10]$. We
obtain $D=4.07_{-0.23}^{+1.01}$, which is consistent with 4-dimensional
spacetime within 1$\sigma$ bound, but the distribution of $p(D)$ has a long
tail at higher values. This constraint is slightly better than that obtained
with the spectral siren method [84]. Meanwhile the constraint on $R_{c}$ is
$\log(R_{c}/{\rm Mpc})=3.93_{-0.55}^{+1.55}$, corresponding to $R_{c}\sim 10$
Gpc, which is around the same scale as the Hubble radius. The drop of the
posterior of $\log(R_{c}/{\rm Mpc})$ at the upper bound is artificial due to
the over-smoothing of KDE; the distribution of samples in $R_{c}$ is actually
flat at the upper limit, as GR is recovered when $R_{c}$ approaches infinity.
In Fig. 10 we trim the range of $R_{c}$ where the over-smoothing takes place.
We have also included $n_{D}$ in the joint posterior, but $n_{D}$ is
essentially unconstrained. Compared with the other two modified gravity
models, there is little degeneracy among $H_{0}$, $D$ and $\log R_{c}$.

Figure 9: Posteriors of $H_{0}$ and $c_{M}$ by combining results of BBHs and
NSBHs from the gwcosmo sampling run.  Figure 10: Posteriors of $H_{0}$, $D$
and $\log R_{c}$ by combining results of BBHs and NSBHs from the gwcosmo
sampling run.

##  6 Conclusions

The detection of gravitational waves offers us a powerful tool for testing
theories of gravity beyond GR. The detection of bright siren GW170817 kicked
off these investigations by providing a strong constraint on a particular
aspect of modified GW propagation, namely the propagation speed. However, due
to its relative proximity, the bounds obtained on other aspects of modified GW
propagation remain relatively weak. We stress that an event ideal for
constraining the Hubble parameter (e.g. one at low redshift) is not
necessarily an ideal system for tests of gravity and fundamental physics.
Therefore, high redshift events should not be always discarded in favour of
apparently simpler analyses at low redshift.

Although more bright siren detections are expected going forward, their rate
and redshift distribution remains currently rather unknown. The dark sirens
technique, which statistically identifies host galaxies of GW events with
galaxy catalogues, enables us to avoid pinning our hopes for GW cosmology and
tests of fundamental physics on the uncertain bright siren landscape. Although
dark sirens are individually less constraining than bright sirens as
cosmological probes, their advantage lies in both their number and security
(meaning that, they rely only on regular BBH and NSBH mergers, which we are
confident will continue to be detected in increasing numbers). Any moderately
well-localised GW event can be used for a dark sirens analysis, provided it
meets a modest SNR threshold and has satisfactory parameter estimation
results. As such, at least tens (or possibly $\sim 100$) of suitable events
are expected during the recently-commenced O4 observation run. We can expect
to see corresponding improvements in the parameter constraints presented here
on the timescale of a year or two.

In this work we described the modification to the gwcosmo software pipeline
for constraining modified gravity parameters that affect GW luminosity
distances, with both bright sirens and dark sirens. A crucial step forward is
that these additional MG parameters can now be constrained concurrently with
hyper-parameters describing the mass distribution and redshift evolution of
compact objects, as well as the Hubble constant.

To validate our pipeline and explore the constraining power of future data
sets, we first performed our tests of gravity on 250 mock bright sirens from
the F2Y mock data catalogue. We computed the 1D posteriors on $\Xi_{0}$,
$c_{M}$ and $D$ in the three modified gravity parameterizations of Section 2
respectively, recovering the expected consistency with GR. The joint posterior
on $H_{0}$ and $\Xi_{0}$ with the F2Y mock bright sirens remained weak: it is
difficult to constrain modified propagation effects with bright sirens located
at low redshifts, where changes to the GW luminosity distances are small.

Turning our attention to the current data, we then computed the posterior
probability on the modified gravity parameters in addition to mass
distribution and rate evolution hyper-parameters with 46 dark sirens from the
GWTC-3 catalogue. The constraint on $\Xi_{0}$ from the marginalized posterior
with 42 BBHs is $1.29_{-0.67}^{+1.22}$, which is consistent with GR. We
observed high degeneracy between $\Xi_{0}$ and $\gamma$ (a parameter
controlling the redshift evolution of the merger rate), which makes the choice
of the prior bound of $\gamma$ important for tests of gravity. This is a good
example of how constraints on cosmology, astrophysics and fundamental physics
are entangled, at least for the present.

Finally, we combined for the first time the results of 42 BBHs and 3 NSBHs,
which yields $\Xi_{0}=1.67_{-0.94}^{+0.93}$. Furthermore, the combined
posterior of BBHs and NSBHs gives the constraint of $c_{M}=1.5_{-2.1}^{+2.2}$
in the Horndeski class model, and $D=4.07_{-0.23}^{+1.01}$ in the extra-
dimension model, all of which are consistent with GR.

At present, the effects of adding galaxy catalogue information are most
pronounced for event GW190814, due to its good localisation and moderately low
redshift; they are not highly significant for most other events. This is
primarily due to the relatively low redshift range for which the GLADE+
catalogue employed our analysis has high completeness. Therefore, for most of
the GW events to date, the bulk of their GW luminosity distance distribution
lies beyond the catalogue. The good news is that other spectroscopic galaxy
catalogues with higher completeness and redshift range are forthcoming, such
as those from the DESI [80] and Euclid [81] experiments (see [85] for an
initial dark siren analysis with the DESI catalogue). In some sense, this is
the much better issue to have: we know that galaxy catalogue completeness will
improve dramatically on a timescale of five years or less. We do not currently
know what the corresponding improvements for bright siren counts will be. On
longer timescales, the measurement of cosmological and modified gravity
parameters with direct cross-correlation between GW events and galaxy surveys
has been forecast in [86, 87, 88, 89, 90] (see also [91, 92]). This can serve
as an independent probe in the 3G era.

A secondary effect in the current constraints is localization of GW events,
which was particularly large for some events in the first two LVK observing
runs. When ‘averaged’ over a large sky area we expect the galaxy clustering
signal – which is key for (dis)favouring some parameter values – to be washed
out. Given the incompleteness issue above, this is not a major stumbling block
at present. With developments at the KAGRA and LIGO-India sites progressing,
observations by three widely-separated detectors in future will lead to
corresponding improvements in event localization, and should prevent this
effect from ever becoming a limiting factor.

In this analysis we selected three specific parameterizations of modified GW
luminosity distances to study. Whilst these functional forms are commonly used
and motivated by general properties of current MG models, they are intended as
representative examples only. The analysis pipelines developed in this work
can be easily adapted to use alternative parameterizations or specific gravity
models, if desired. In future work we intend to broaden the range of options
in gwcosmo, including effects beyond modified GW luminosity distances.

As the fourth LVK observing run progresses, there may well be significant
advances or new discoveries in some of the ingredients that make up the dark
sirens recipe, e.g. the mass distribution of compact objects, its possible
evolution with redshift, or the merger rates of different source types (as
analysed with GWTC-3 [71]). With our current framework in hand, these are all
straightforwards to accommodate. Such improvements will help to pin down the
astrophysical parts of the analysis, lifting degeneracies with the
cosmological and beyond-GR parameters. We have set the stage for improved dark
siren tests of gravity, and eagerly anticipate the forthcoming GW and galaxy
data.

## Acknowledgments

The authors would like to thank Benoit Revenu and Christos Karathanasis for
assistance with GW injection computations. A.C. is supported by a PhD grant
from the Chinese Scholarship Council (grant no.202008060014), and the STSM
GravNet grant during his visit to the NBI. The research of R.G. was supported
by ERC Starting Grant SHADE 949572 and STFC grant ST/V005634/1. 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). 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.

##  Appendix A Corner plots for GWTC-3 reanalysis

We present the full corner plot for the joint posteriors obtained from the
gwcosmo sampling run with 42 BBHs in GWTC-3 in Fig. 11. $\gamma$, $\kappa$ and
$z_{p}$ in the merger rate evolution model are essentially unconstrained. On
the other hand, parameters for BBH population model and cosmological
parameters are constrained at a better level. There exists notable degeneracy
between $H_{0}$-$\Xi_{0}$, $\gamma$-$\Xi_{0}$, $\mu_{g}$-$\sigma_{g}$ and
$m_{\rm min}^{\rm BH}$-$\sigma_{m}$ (see discussion in the main text).

Figure 11: Full corner plot of joint posteriors for 42 BBH events from the
GWTC-3 catalogue. The contours show 1, 2 and 3$\sigma$ bound for the 2D
posteriors. The vertical lines in the 1D histogram plot correspond to the peak
values and 1$\sigma$ bound of parameters in the title.

##  Appendix B Posterior samples for GWTC-3 reanalysis

Here we explain the three varieties of $\Xi_{0}$ posteriors observed from
individual events in Fig 4. In Fig. 12 we plot histograms of the redshift
posterior samples for the 46 events selected from GWTC-2.1 and GWTC-3. The
redshifts of posterior samples are converted from $d_{L}^{\rm GW}$ with two
values, $\Xi_{0}=1$ (GR) and $\Xi_{0}=2$ with $H_{0}=70~{}{\rm km}~{}{\rm
s}^{-1}{\rm Mpc}^{-1}$ and $n=1.91$. The grey dashed vertical line marks the
estimated redshift in GR from GWTC-3. We can see that the posterior samples
shift to lower redshift as $\Xi_{0}$ increases. For the events in which the
peak of the samples is located at higher redshift than the estimated redshift,
the peak crosses through to the other side of the estimated redshift when
$\Xi_{0}$ increases. As a result, there is a peak in the likelihood of
$\Xi_{0}$ over the prior range of $\Xi_{0}$, see for example GW170809_082821
in Fig. 4. For the events in which the peaks of the samples for different
$\Xi_{0}$ values are on the same side of the estimated redshift, the
likelihood of $\Xi_{0}$ continuously decreases instead of having a peak, for
example GW190910_112807 in Fig. 4. For the events in which the posterior
samples shift little when $\Xi_{0}$ increases, the likelihood is dominated by
the denominator that is continuously decreasing, so the likelihood is
continuously increasing, for example GW191216_213338 in Fig. 4.

Figure 12: Posterior samples converted from GW luminosity distance to redshift
for the 46 events from the GWTC-3 catalogue with different values of
$\Xi_{0}$. The grey dashed vertical line marks the estimated redshift in GR
from GWTC-3.

## References

  * [1] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Tests of general relativity with binary black holes from the second ligo-virgo gravitational-wave transient catalog_ , _Phys. Rev. D_ 103 (2021) 122002.
  * [2] LIGO Scientific, VIRGO, KAGRA collaboration, _Tests of General Relativity with GWTC-3_ , 2112.06861.
  * [3] L. Verde, T. Treu and A.G. Riess, _Tensions between the early and late universe_ , _Nature Astronomy_ 3 (2019) 891.
  * [4] E. Di Valentino et al., _Snowmass2021 - letter of interest cosmology intertwined ii: The hubble constant tension_ , _Astroparticle Physics_ 131 (2021) 102605.
  * [5] B. Schutz, _Determining the hubble constant from gravitational wave observations_ , _Nature_ 323 (1986) 310.
  * [6] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Gw170817: Observation of gravitational waves from a binary neutron star inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101.
  * [7] B.P. Abbott et al., _Gravitational waves and gamma-rays from a binary neutron star merger: Gw170817 and grb 170817a_ , _The Astrophysical Journal Letters_ 848 (2017) L13.
  * [8] LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT collaboration, _Multi-messenger Observations of a Binary Neutron Star Merger_ , _Astrophys. J. Lett._ 848 (2017) L12 [1710.05833].
  * [9] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _A gravitational-wave standard siren measurement of the hubble constant_ , _Nature_ 551 (2017) 85.
  * [10] C. de Rham and S. Melville, _Gravitational Rainbows: LIGO and Dark Energy at its Cutoff_ , _Phys. Rev. Lett._ 121 (2018) 221101 [1806.09417].
  * [11] LISA Cosmology Working Group collaboration, _Measuring the propagation speed of gravitational waves with LISA_ , _JCAP_ 08 (2022) 031 [2203.00566].
  * [12] T. Baker, E. Barausse, A. Chen, C. de Rham, M. Pieroni and G. Tasinato, _Testing gravitational wave propagation with multiband detections_ , _JCAP_ 03 (2023) 044 [2209.14398].
  * [13] I. Harry and J. Noller, _Probing the speed of gravity with LVK, LISA, and joint observations_ , _Gen. Rel. Grav._ 54 (2022) 133 [2207.10096].
  * [14] T. Baker, E. Bellini, P.G. Ferreira, M. Lagos, J. Noller and I. Sawicki, _Strong constraints on cosmological gravity from GW170817 and GRB 170817A_ , _Phys. Rev. Lett._ 119 (2017) 251301 [1710.06394].
  * [15] J.M. Ezquiaga and M. Zumalacárregui, _Dark energy after gw170817: Dead ends and the road ahead_ , _Phys. Rev. Lett._ 119 (2017) 251304.
  * [16] P. Creminelli and F. Vernizzi, _Dark energy after gw170817 and grb170817a_ , _Phys. Rev. Lett._ 119 (2017) 251302.
  * [17] S. Jana, G.K. Chakravarty and S. Mohanty, _Constraints on born-infeld gravity from the speed of gravitational waves after gw170817 and grb 170817a_ , _Phys. Rev. D_ 97 (2018) 084011.
  * [18] A. Nishizawa and T. Kobayashi, _Parity-violating gravity and gw170817_ , _Phys. Rev. D_ 98 (2018) 124018.
  * [19] 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].
  * [20] W. Del Pozzo, _Inference of cosmological parameters from gravitational waves: Applications to second generation interferometers_ , _Phys. Rev. D_ 86 (2012) 043011.
  * [21] H.-Y. Chen, M. Fishbach and D.E. Holz, _A two per cent hubble constant measurement from standard sirens within five years_ , _Nature_ 562 (2018) 545.
  * [22] M. Fishbach et al., _A standard siren measurement of the hubble constant from gw170817 without the electromagnetic counterpart_ , _The Astrophysical Journal Letters_ 871 (2019) L13.
  * [23] R. Gray, I.M.n. Hernandez, H. Qi, A. Sur, P.R. Brady, H.-Y. Chen et al., _Cosmological inference using gravitational wave standard sirens: A mock data analysis_ , _Phys. Rev. D_ 101 (2020) 122001.
  * [24] K. Leyde, S. Mastrogiovanni, D. Steer, E. Chassande-Mottin and C. Karathanasis, _Current and future constraints on cosmology and modified gravitational wave friction from binary black holes_ , _Journal of Cosmology and Astroparticle Physics_ 2022 (2022) 012.
  * [25] 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].
  * [26] R. Gray et al., _Joint cosmological and gravitational-wave population inference using dark sirens and galaxy catalogues_ , 2308.02281.
  * [27] S. Mastrogiovanni, D. Laghi, R. Gray, G.C. 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 preprint_ (2023) [2305.10488].
  * [28] S. Mastrogiovanni, G. Pierra, S. Perriès, D. Laghi, G.C. Santoro, A. Ghosh et al., _ICAROGW: A python package for inference of astrophysical population properties of noisy, heterogeneous and incomplete observations_ , 2305.17973.
  * [29] 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.
  * [30] LIGO Scientific, VIRGO, KAGRA collaboration, _GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run_ , 2111.03606.
  * [31] R. Abbott et al., _Constraints on the cosmic expansion history from GWTC-3_ , _arXiv e-prints_ (2021) arXiv:2111.03604 [2111.03604].
  * [32] 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_, 2203.03643.
  * [33] 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].
  * [34] M. Oguri, _Measuring the distance-redshift relation with the cross-correlation of gravitational wave standard sirens and galaxies_ , _Phys. Rev. D_ 93 (2016) 083511.
  * [35] E. Belgacem, Y. Dirian, S. Foffa and M. Maggiore, _Gravitational-wave luminosity distance in modified gravity theories_ , _Phys. Rev. D_ 97 (2018) 104066.
  * [36] E. Belgacem, Y. Dirian, S. Foffa and M. Maggiore, _Modified gravitational-wave propagation and standard sirens_ , _Phys. Rev. D_ 98 (2018) 023510.
  * [37] A. Nishizawa, _Generalized framework for testing gravity with gravitational-wave propagation. i. formulation_ , _Phys. Rev. D_ 97 (2018) 104037.
  * [38] E. Belgacem, G. Calcagni, M. Crisostomi, C. Dalang, Y. Dirian, J.M. Ezquiaga et al., _Testing modified gravity at cosmological distances with lisa standard sirens_ , _Journal of Cosmology and Astroparticle Physics_ 2019 (2019) 024.
  * [39] T. Baker and I. Harrison, _Constraining scalar-tensor modified gravity with gravitational waves and large scale structure surveys_ , _Journal of Cosmology and Astroparticle Physics_ 2021 (2021) 068.
  * [40] M. Maggiore, _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford University Press (2007), 10.1093/acprof:oso/9780198570745.001.0001.
  * [41] A.E. Romano and M. Sakellariadou, _Mirage of luminal modified gravitational-wave propagation_ , _Phys. Rev. Lett._ 130 (2023) 231401.
  * [42] E. Belgacem, Y. Dirian, S. Foffa and M. Maggiore, _Gravitational-wave luminosity distance in modified gravity theories_ , _Phys. Rev. D_ 97 (2018) 104066 [1712.08108].
  * [43] E. Belgacem, Y. Dirian, S. Foffa and M. Maggiore, _Nonlocal gravity. Conceptual aspects and cosmological predictions_ , _JCAP_ 03 (2018) 002 [1712.07066].
  * [44] M. Mancarella, E. Genoud-Prachex and M. Maggiore, _Cosmology and modified gravitational wave propagation from binary black hole population models_ , _Phys. Rev. D_ 105 (2022) 064030 [2112.05728].
  * [45] M. Mancarella, N. Borghi, S. Foffa, E. Genoud-Prachex, F. Iacovelli, M. Maggiore et al., _Gravitational-wave cosmology with dark sirens: state of the art and perspectives for 3G detectors_ , _PoS_ ICHEP2022 (2022) 127 [2211.15512].
  * [46] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore and M. Mancarella, _Modified gravitational wave propagation and the binary neutron star mass function_ , _Phys. Dark Univ._ 36 (2022) 100994 [2108.04065].
  * [47] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore and M. Mancarella, _Probing modified gravitational wave propagation with strongly lensed coalescing binaries_ , _Phys. Rev. D_ 104 (2021) 084057 [2107.05046].
  * [48] G.W. Horndeski, _Second-order scalar-tensor field equations in a four-dimensional space_ , _Int. J. Theor. Phys._ 10 (1974) 363.
  * [49] C. Deffayet, X. Gao, D.A. Steer and G. Zahariade, _From k-essence to generalised Galileons_ , _Phys. Rev. D_ 84 (2011) 064039 [1103.3260].
  * [50] J.M. Ezquiaga and M. Zumalacárregui, _Dark Energy After GW170817: Dead Ends and the Road Ahead_ , _Phys. Rev. Lett._ 119 (2017) 251304 [1710.05901].
  * [51] P. Creminelli and F. Vernizzi, _Dark Energy after GW170817 and GRB170817A_ , _Phys. Rev. Lett._ 119 (2017) 251302 [1710.05877].
  * [52] J. Sakstein and B. Jain, _Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories_ , _Phys. Rev. Lett._ 119 (2017) 251303 [1710.05893].
  * [53] E. Bellini and I. Sawicki, _Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity_ , _Journal of Cosmology and Astroparticle Physics_ 2014 (2014) 050.
  * [54] M. Lagos, E. Bellini, J. Noller, P.G. Ferreira and T. Baker, _A general theory of linear cosmological perturbations: stability conditions, the quasistatic limit and dynamics_ , _Journal of Cosmology and Astroparticle Physics_ 2018 (2018) 021.
  * [55] E. Bellini, A.J. Cuesta, R. Jimenez and L. Verde, _Constraints on deviations from $\lambda$cdm within horndeski gravity_, _Journal of Cosmology and Astroparticle Physics_ 2016 (2016) 053.
  * [56] D. Alonso, E. Bellini, P.G. Ferreira and M. Zumalacárregui, _Observational future of cosmological scalar-tensor theories_ , _Phys. Rev. D_ 95 (2017) 063502.
  * [57] M. Lagos, M. Fishbach, P. Landry and D.E. Holz, _Standard sirens with a running planck mass_ , _Phys. Rev. D_ 99 (2019) 083504.
  * [58] S. Mastrogiovanni, D.A. Steer and M. Barsuglia, _Probing modified gravity theories and cosmology using gravitational-waves and associated electromagnetic counterparts_ , _Phys. Rev. D_ 102 (2020) 044009.
  * [59] E.V. Linder, _Challenges in connecting modified gravity theory and observations_ , _Phys. Rev. D_ 95 (2017) 023518 [1607.03113].
  * [60] R. D’Agostino and R.C. Nunes, _Probing observational bounds on scalar-tensor theories from standard sirens_ , _Phys. Rev. D_ 100 (2019) 044041 [1907.05516].
  * [61] G. Dvali, G. Gabadadze and M. Porrati, _4d gravity on a brane in 5d minkowski space_ , _Physics Letters B_ 485 (2000) 208.
  * [62] A. Padilla, _Ghost-free braneworld bigravity_ , _Classical and Quantum Gravity_ 21 (2004) 2899.
  * [63] M. Corman, C. Escamilla-Rivera and M. Hendry, _Constraining extra dimensions on cosmological scales with lisa future gravitational wave siren data_ , _Journal of Cosmology and Astroparticle Physics_ 2021 (2021) 005.
  * [64] I. Mandel, W.M. Farr and J.R. Gair, _Extracting distribution parameters from multiple uncertain observations with selection biases_ , _Monthly Notices of the Royal Astronomical Society_ 486 (2019) 1086 [https://academic.oup.com/mnras/article-pdf/486/1/1086/28390969/stz896.pdf].
  * [65] 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_ , in _Handbook of Gravitational Wave Astronomy_ , C. Bambi, S. Katsanevas and K.D. Kokkotas, eds., (Singapore), pp. 1–60, Springer Singapore (2020), DOI.
  * [66] 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].
  * [67] 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.
  * [68] 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_ , _The Astrophysical Journal_ 622 (2005) 759.
  * [69] 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].
  * [70] R. Abbott et al., _Population properties of compact objects from the second ligo–virgo gravitational-wave transient catalog_ , _The Astrophysical Journal Letters_ 913 (2021) L7.
  * [71] KAGRA, VIRGO, LIGO Scientific collaboration, _Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3_ , _Phys. Rev. X_ 13 (2023) 011048 [2111.03634].
  * [72] P. Madau and M. Dickinson, _Cosmic star-formation history_ , _Annual Review of Astronomy and Astrophysics_ 52 (2014) 415–486.
  * [73] L.P. Singer, L.R. Price, B. Farr, A.L. Urban, C. Pankow, S. Vitale et al., _The first two years of electromagnetic follow-up with advanced ligo and virgo_ , _The Astrophysical Journal_ 795 (2014) 105.
  * [74] L. Barsotti and P. Fritschel, _Early aligo configurations: Example scenarios toward design sensitivity_ , _Tech. Rep. LIGO-T1200307-v4_ (2012) .
  * [75] KAGRA, LIGO Scientific, Virgo, VIRGO collaboration, _Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA_ , _Living Rev. Rel._ 21 (2018) 3 [1304.0670].
  * [76] H.-Y. Chen, M. Fishbach and D.E. Holz, _A two per cent Hubble constant measurement from standard sirens within five years_ , _Nature_ 562 (2018) 545 [1712.06531].
  * [77] Planck Collaboration collaboration, _Planck 2015 results. XIV. Dark energy and modified gravity_ , _aap_ 594 (2016) A14 [1502.01590].
  * [78] The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration collaboration, _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 preprint_ (2021) [2108.01045].
  * [79] 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) [2111.03606].
  * [80] A. Dey et al., _Overview of the desi legacy imaging surveys_ , _The Astronomical Journal_ 157 (2019) 168.
  * [81] G.D. Racca et al., _The Euclid mission design_ , in _Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave_ , H.A. MacEwen, G.G. Fazio, M. Lystrup, N. Batalha, N. Siegler and E.C. Tong, eds., vol. 9904, p. 99040O, International Society for Optics and Photonics, SPIE, 2016, DOI.
  * [82] M. Lagos, M. Fishbach, P. Landry and D.E. Holz, _Standard sirens with a running Planck mass_ , _Phys. Rev. D_ 99 (2019) 083504 [1901.03321].
  * [83] J.M. Ezquiaga, _Hearing gravity from the cosmos: GWTC-2 probes general relativity at cosmological scales_ , _Phys. Lett. B_ 822 (2021) 136665 [2104.05139].
  * [84] I. Magana Hernandez, _Constraining the number of spacetime dimensions from GWTC-3 binary black hole mergers_ , _Phys. Rev. D_ 107 (2023) 084033 [2112.07650].
  * [85] A. Palmese, C.R. Bom, S. Mucesh and W.G. Hartley, _A Standard Siren Measurement of the Hubble Constant Using Gravitational-wave Events from the First Three LIGO/Virgo Observing Runs and the DESI Legacy Survey_ , _Astrophys. J._ 943 (2023) 56 [2111.06445].
  * [86] S. Mukherjee and B.D. Wandelt, _Beyond the classical distance-redshift test: cross-correlating redshift-free standard candles and sirens with redshift surveys_ , 1808.06615.
  * [87] S. Mukherjee, B.D. Wandelt and J. Silk, _Probing the theory of gravity with gravitational lensing of gravitational waves and galaxy surveys_ , _Mon. Not. Roy. Astron. Soc._ 494 (2020) 1956 [1908.08951].
  * [88] S. Mukherjee, B.D. Wandelt and J. Silk, _Testing the general theory of relativity using gravitational wave propagation from dark standard sirens_ , _Mon. Not. Roy. Astron. Soc._ 502 (2021) 1136 [2012.15316].
  * [89] S. Mukherjee, B.D. Wandelt, S.M. Nissanke and A. Silvestri, _Accurate precision Cosmology with redshift unknown gravitational wave sources_ , _Phys. Rev. D_ 103 (2021) 043520 [2007.02943].
  * [90] C.C. Diaz and S. Mukherjee, _Mapping the cosmic expansion history from LIGO-Virgo-KAGRA in synergy with DESI and SPHEREx_ , _Mon. Not. Roy. Astron. Soc._ 511 (2022) 2782 [2107.12787].
  * [91] G. Cañas Herrera, O. Contigiani and V. Vardanyan, _Learning How to Surf: Reconstructing the Propagation and Origin of Gravitational Waves with Gaussian Processes_ , _Astrophys. J._ 918 (2021) 20 [2105.04262].
  * [92] J. Fonseca, S. Zazzera, T. Baker and C. Clarkson, _The observed number counts in luminosity distance space_ , 2304.14253.