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

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 H0 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.


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 strongfield 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 km s −1 Mpc −1 [9].The implications of GW170817 for tests of gravity were even more impactful.The arrival of the GRB counterpart just ∼ 1.7s after the GW merger signal constrained the propagation speed of GWs (c T ) in the frequency band of terrestrial detectors 1 at the level |1 − c T /c| ⩽ 3 × 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 result2 .
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 km s −1 Mpc −1 .Combining the dark siren result with the bright siren GW170817 yields H 0 = 68 +8 −6 km s −1 Mpc −1 .Meanwhile the analysis with 42 BBHs using Icarogw 2.0 gives H 0 = 71 +35 −30 km s −1 Mpc −1 with the assumption of the compact object binary merger rate proportional to the galaxy luminosity, and H 0 = 43 +48 −18 km s −1 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.

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 [25,[35][36][37][38][39] h′′ where h(η, k) is the GW amplitude in Fourier space, subscript A denotes the + or × polarization of GWs, η is conformal time, ′ denotes the derivative with respect to η, and H = a ′ /a is the conformal Hubble parameter.Here δ(η) is the modification to the friction term introduced by modified gravity theories, where GR is recovered when δ(η) = 0.It is shown in Sec.4.1.4 of [40] that the solution to Eq. (2.1) in GR gives in the late-time universe.Ω m,0 and Ω Λ,0 are the matter density today and the dark energy density today respectively.When δ(η) ̸ = 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 ã such that 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 GW L (z), which is related to d EM L (z) by In a non-GR theory, we now have h(η, k) ∝ 1/d GW L (z).Eq. ( 2.3) can be written as Integrating Eq. (2.5) allows us to re-express the relation between d GW L (z) and d EM L (z) as: (2.6) The exact form of δ(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.

(Ξ 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 GW L (z)/d EM L (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 GW L (z)/d EM L (z) is given by [36] d GW L (z) where Ξ 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 [24,25,27,44,45], and with other techniques such as using BNS mass distribution [46] and strongly lensed GW events [47].Some mappings of Ξ 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 Ξ 0 = 1, so that the GW distance is the same as the EM luminosity distance.This limit is also recovered for z ∼ 0, largely irrespective of Ξ 0 .Conversely, when z ≫ 1 the ratio d GW L (z)/d EM L (z) approaches a constant set by the value of Ξ 0 .This is motivated by the hypothesis that any deviation from GR predominantly affects the late universe where dark energy dominates, so d GW L (z)/d EM L (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 Ξ 0 take place as redshift increases.In the low redshift limit z ≪ 1, the ratio can be approximated at the linear order in z: 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 Ξ 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 Ξ 0 and n, if the data are consistent with GR.
In Fig. 1 we plot the ratio of d GW L (z)/d EM L (z) for some different values of Ξ 0 .The ratio is greater than 1 when Ξ 0 > 1 and lower than 1 when Ξ 0 < 1, for n ̸ = 0.

Horndeski class parameterization
In addition to the generic (Ξ 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]] where ϕ is the scalar degree of freedom, M eff is the effective Planck mass which evolves with ϕ, X is a kinetic term defined by X ≡ − ▽ µ ϕ ▽ µ ϕ/2, S m is the matter action, and ψ m are matter fields minimally coupled to the metric g µν .When M eff equals the Planck mass M P , and K = G 3 = 0, the action reduces to GR.It has been shown in the literature (e.g.[39,53,54]) that the GW propagation equation derived from the reduced Horndeski action is given by where α M is the running rate of the effective Planck mass defined by

.11)
We see from Eq. (2.1) that the GW friction term δ(z) is related to α M (z) by δ(z) = −α M (z)/2.A widely-used and simple parameterization for α M (z) is that it is proportional to the fractional dark energy density [55][56][57][58] as where c M is the free parameter to be constrained, and Ω Λ (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 δ(z) = −α M (z)/2, one obtains by substituting the above parameterization into Eq.(2.6) that (2.13) This parameterization only depends on one free parameter, c M .GR is recovered when c M = 0.
We have d GW L (z)/d EM L (z) > 1 when c M > 0, and d GW L (z)/d EM L (z) < 1 when c M < 0, as shown in Fig. 1.

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 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 grows with the comoving distance to the power of (D − 4)/2.Constraints on D with bright and dark sirens have been studied in [24,63].
The comparison between d GW L and d EM L 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 GW L (z)/d EM L (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.

Bayesian framework of gwcosmo
The detailed Bayesian statistical framework of the gwcosmo code is presented in [19,23,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 Λ 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 det detected GW events with observation data {x GW }, the posterior for a set of hyper-parameters Λ is given by [64,65] p(Λ|{x GW }, {D GW }, I) ∝ p(Λ|I)p(N det |Λ, I) where D 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 Λ.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(Λ|I) is uniform in the analysis for all parameters in Λ. p(N det |Λ, I) is the probability of detecting N det events.It depends on the intrinsic astrophysical rate of GW events R. As discussed in [64,66], by choosing a flat prior on R in a log scale, we obtain the combined prior p(Λ, R|I) ∝ 1/R, so that the dependence on R is dropped out when marginalizing p(N det |Λ, I) over R. Beyond this p(N det |Λ, I) has no further dependence on hyper-parameters.With Bayes' theorem, the likelihood for obtaining a set of GW observations x GW can be expanded into ) where θ represents a parameter set describing a detected GW signal, which must be marginalised over in order to find posterior for Λ.The parameters in θ of interest to us in this work include the detector frame masses m det 1 and m det 2 , redshift z, and sky location Ω. p(x GWi |θ, Λ, 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 hyperparameters are relevant, while population and merger rate hyper-parameters enter in p(θ|Λ, I).p(D GWi |θ, Λ, 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 IMRPhe-nomPv2 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 θ 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|Λ, 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 ) where p(z|Ω j , Λ, I) is the LOS redshift prior for each pixel at a sky location Ω 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 Λ.Here we extend gwcosmo to the estimation of sets of modified gravity parameters such as Ξ 0 , n in Eq. (2.7), c M in Eq. (2.13), or D in Eq. (2.14).Since these parameters modify d GW L (z), but leave d EM L (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 GW L (z) now differs from that from d EM L (z) under different values of modified gravity parameters.Inside a pixel, the numerator integrand of the likelihood is given by where m s 1 and m s 2 are the component masses in the GW source frame, with the assumption m s 1 > m s 2 .m min and m 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 s 1 , m s 2 |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 GW L (z) that is measured from the GW signal instead of z.Therefore we have ) is the prior applied to d GW L , m det 1 and m det 2 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: uniform in comoving volume prior.In this paper we used the data set with π( (3.7) The conversion from d GW L to z, and the Jacobian |∂z/∂d GW L | 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 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)): dm det 2 dm det 1 .
(3.9)The posterior samples on z computed from d GW L given Λ, m det 1 and m det 2 for a GW event can be approximated as a sum over delta-functions: where k denotes the k th posterior sample and N pos denotes the total number of posterior samples.Inserting this back to Eq. (3.9) we obtain the result as (3.11) The above quantity is then multiplied by the LOS redshift prior3 p(z|Ω j , Λ, 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 The details of computing the LOS redshift prior p(z k |Ω j , Λ, 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 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 Λ.

Prior on source masses and merger rate evolution
The source mass prior p(m s 1 , m s 2 |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,26,27].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, Ω j , Λ, I) = p(s|z, Ω j , Λ, I)p(z|Ω j , Λ, I) p(s|Ω j , Λ, I) . (3.14) The term p(s|Ω j , Λ, 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 Ω 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 (3.15) The free parameters γ, 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

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 ⊙ , 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 EM L ≈ cz/H 0 is applied for the EM counterparts of BNS injections.The BNS d GW L is then parameterized with the linearized d EM L 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 3 0 [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 ⊙ , since the F2Y BNSs were generated with component masses between 1.2 and 1.6 M ⊙ .Such a mass prior is much narrower than the range of [1.0, 3.0] M ⊙ used in later LVK analyses: neutron stars were generally expected to have masses strongly peaked around 1.4 M ⊙ before much heavier neutron stars were detected in events from the O3 observing run.

Results
We first measure Ξ 0 in the (Ξ 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 km s −1 Mpc −1 and n = 1.91 4 .The normalized posterior of Ξ 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σ bound is Ξ 0 = 0.96 +0. 19 −0.19 , which is consistent with GR.Our results show that gwcosmo is capable of constraining Ξ 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, Ξ 0 is most relevant to high-redshift GW events, and the mock bright sirens used here are all at low redshifts.However, the 1D measurement of Ξ 0 can be misleading, because it has assumed H 0 to be a fixed value.Since H 0 and Ξ 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 -Ξ 0 joint analysis.However, the number of posterior samples for the F2Y mock BNS is ∼ 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, Ξ 0 = 1) is on the edge of the 3σ 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 Ξ 0 measurement.In the measurement of c M , the values of Ω m,0 and Ω Λ,0 are needed to reconstruct the evolution of the fraction of dark energy density as described in Eq. (2.12).We adopt Ω m,0 = 0.3065 from the Planck result [77] as in the LVK cosmology paper [31], and Ω Λ,0 = 1 − Ω 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σ 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σ bound measured value is D = 4.015 +0.027 −0.027 , which is consistent with the GR limit of this class of theories.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.

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 × 10 6 GW injections with the SNR threshold set to 11.The injected detected masses range from 1 to 500 M ⊙ , and the injected GW distances are varied between 0.1 and ∼ 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 ∈ (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 Ξ 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 Ξ 0 such that the maximum injected GW distance corresponds to redshift not larger than 10.As a result, we vary Ξ 0 in a range of [0.3,10] when fixing H 0 = 70 km s −1 Mpc −1 .A combination of the maximum H 0 and the minimum Ξ 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 Ξ 0 when varying both of them.We chose a range of [20,140] km s −1 Mpc −1 for H 0 and [0.35,10] for Ξ 0 for the joint analysis.

1D measurement of Ξ 0
We first measure the posterior of Ξ 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 km s −1 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 α = 3.78, β = 0.81, m BH min = 4.98M ⊙ , m BH max = 112.5M⊙ , δ m = 4.8M ⊙ , µ g = 32.27M⊙ , σ g = 3.88M ⊙ and λ 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 NS min = 1.0M ⊙ and m NS max = 3.0M ⊙ .The merger rate evolution parameters are also fixed with the same values γ = 4.59, k = 2.86 and z p = 2.47.The normalized posterior of Ξ 0 for each individual event is plotted in Fig. 4. We can see that there are three typical patterns in the distribution of p(Ξ 0 ) for different events: favoring low Ξ 0 such as GW190910_112807, favoring high Ξ 0 such as GW191216_213338, and peaking at a certain value of Ξ 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 GW L for different Ξ 0 for each event, which we will elaborate on in Appendix B. The measured value from the combined posterior with 1σ bound is Ξ 0 = 1.06 +0.25 −0.18 , which is consistent with GR.
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 Ξ 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(Ξ 0 ) with galaxy catalogue support peaks at lower Ξ 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 Ξ 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.The one-dimensional likelihood of Ξ 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.The 1D posterior of Ξ 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 Ξ 0 , creating a bumpy structure in the likelihood.Comparable features can be seen in the H 0 likelihood shown in Fig. 9 of [26].

BBH joint analysis
Next we perform a joint posterior measurement for all hyper-parameters under the (Ξ 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 γ, m BH max , µ g , H 0 , Ξ 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 km s −1 Mpc −1 , again an uncertainty similar to previous results.The estimated value for Ξ 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 Ξ(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 Ξ 0 = 1), if the data are consistent with GR.Moreover, as shown by the 2D Ξ 0 -n posterior contour, Ξ 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 Ξ 0 as the more dominant term.Therefore the constraint is put more directly on Ξ 0 when n is large.
From the contour plot we can see that H 0 and Ξ 0 are partially degenerate with each other.This is because d GW L ∝ Ξ 0 /H 0 as shown in Eq. (2.2) and (2.7), so that in the conversion of posterior samples from d GW L to z in computing the likelihood, varying H 0 and Ξ 0 while keeping the ratio Ξ 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 Ξ 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 Ξ 0 is strongly degenerate with γ.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) γ when z ≪ 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 γ 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 Ξ 0 and (1 + z) γ , leading to a strong degeneracy between Ξ 0 and γ.A similar feature is also found in Fig. 12 of the Icarogw results [24], where the distribution of p(γ) with the Power Law + Peak BBH population model is similar to ours.Given that the distribution of p(γ) favors a higher value, if we allow the prior bound of γ to be larger, the distribution of p(Ξ 0 ) will migrate to higher values because of the degeneracy.The prior ranges for γ, 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.

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, The power of the power law distribution of the rate U(0.0, 12.0) evolution before redshift z p k The power of the power law distribution of the rate U(0.0, 10.0) evolution after redshift z p z p The redshift turning point between two power law U(0.0,10.0) distributions Table 1: Prior of the parameters in the Power Law + Peak BBH mass distribution.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 NS min and m NS max for the uniform neutron star population model are included in the estimation.We use the uniform prior of m NS min ∈ [0.5, 1.5]M ⊙ and m NS max ∈ [1.501, 5]M ⊙ .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 µ g and σ g indicate a narrow Gaussian peak at ∼ 24M ⊙ for the black hole population model.This is caused alone by the black hole with mass ∼ 23M ⊙ in GW190814.The same appearance is probably also shown in the constraint of m NS max , which is driven by the assumption that GW190814's secondary mass is a neutron star with mass ∼ 2.6M ⊙ .Meanwhile the constraints on cosmological parameters are generally weaker than those with the BBHs due to the lower number of NSBH events.

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 NS max and m NS min , the sampling data from BBHs and NSBHs cannot be combined simply.We first compute a kernel density estimation (KDE) function for the multi-    1.
and the NSBH posterior, with a slightly smaller 1σ uncertainty than the individual BBH and NSBH posteriors, which is as expected.Again, we see the partial degeneracies between Ξ 0 and H 0 , and Ξ 0 and n are maintained.
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    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 ∼ 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 Ξ 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 Ξ 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 hyperparameters with 46 dark sirens from the GWTC-3 catalogue.The constraint on Ξ 0 from the marginalized posterior with 42 BBHs is 1.29 +1.22 −0.67 , which is consistent with GR.We observed high degeneracy between Ξ 0 and γ (a parameter controlling the redshift evolution of the merger rate), which makes the choice of the prior bound of γ 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 Ξ 0 = 1.67 +0.93  −0.94 .Furthermore, the combined posterior of BBHs and NSBHs gives the constraint of c M = 1.5At 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.the samples is located at higher redshift than the estimated redshift, the peak crosses through to the other side of the estimated redshift when Ξ 0 increases.As a result, there is a peak in the likelihood of Ξ 0 over the prior range of Ξ 0 , see for example GW170809_082821 in Fig. 4. For the events in which the peaks of the samples for different Ξ 0 values are on the same side of the estimated redshift, the likelihood of Ξ 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 Ξ 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 1 :
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 (Ξ 0 , n) parameterization, and R c = 100 Mpc in the extra-dimension model.

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

Figure 3 :
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.

Figure 4 :
Figure4: The one-dimensional likelihood of Ξ 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.

Figure 5 :
Figure5: The 1D posterior of Ξ 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 Ξ 0 , creating a bumpy structure in the likelihood.Comparable features can be seen in the H 0 likelihood shown in Fig.9of[26].
Modified gravity parameter controlling high-z limit of U(0.35, 10.0) distance ratio in Eq. (2.7) n Modified gravity parameter controlling steepness of U(0.1, 10.0) distance ratio in Eq. (2.7) α The power of the power law component in the primary U(1.5, 8.0) mass distribution β The power of the power law component in the mass U(−4.0, 6.0) ratio distribution m BH min [M ⊙ ] The minimum mass of the mass distribution U(2.0, 10.0) m BH max [M ⊙ ] The maximum mass of the mass distribution U(50.0, 200.0) λ p Fraction of the model in the Gaussian component U(0.0, 0.5) µ g Mean of the Gaussian component in the primary mass U(10.0, 50.0) distribution σ g Width of the Gaussian component in the primary mass U(0.1, 20.0) distribution δ m Range of mass tapering at the lower end of the mass U(0.0, 15.0) distribution γ

Figure 6 :
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σ bound for the 2D posteriors.The vertical lines in the 1D histogram plot correspond to the peak values and 1σ bound of the parameters in the title.The full corner plot including all parameters is in Fig.11.A brief description of the mass and merger rate hyper-parameters is given in Table1.

Figure 7 :
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σ bound for the 2D posteriors.The vertical lines in the 1D histogram plot correspond to the peak values and 1σ bound of the parameters in the title.A brief description of the mass and merger rate hyper-parameters is given in Table1.
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.2 −2.1 , which

Figure 9 :
Figure 9: Posteriors of H 0 and c M by combining results of BBHs and NSBHs from the gwcosmo sampling run.

Figure 10 :
Figure 10: Posteriors of H 0 , D and log R c by combining results of BBHs and NSBHs from the gwcosmo sampling run.

Figure 12 :
Figure 12: Posterior samples converted from GW luminosity distance to redshift for the 46 events from the GWTC-3 catalogue with different values of Ξ 0 .The grey dashed vertical line marks the estimated redshift in GR from GWTC-3.
+2.2 −2.1 in the Horndeski class model, and D = 4.07 +1.01 −0.23 in the extradimension model, all of which are consistent with GR.