A Model-Independent Precision Test of General Relativity using Bright Standard Sirens from ongoing and upcoming detectors

Gravitational waves (GWs) provide a new avenue to test Einstein's General Relativity (GR) using the ongoing and upcoming GW detectors by measuring the redshift evolution of the effective Planck mass proposed by several modified theories of gravity. We propose a model-independent, data-driven approach to measure any deviation from GR in the GW propagation effect by combining multi-messenger observations of GW sources accompanied by EM counterparts, commonly known as bright sirens (Binary Neutron Star (BNS) and Neutron Star Black Hole systems (NSBH)). We show that by combining the GW luminosity distance measurements from bright sirens with the Baryon Acoustic Oscillation (BAO) measurements derived from galaxy clustering, and the sound horizon measurements from the Cosmic Microwave Background (CMB), we can make a data-driven reconstruction of deviation of the variation of the effective Planck mass (jointly with the Hubble constant) as a function of cosmic redshift. Using this technique, we achieve a precise measurement of GR with redshift (z) with a precision of approximately $7.9\%$ for BNSs at redshift z=0.075 and $10\%$ for NSBHs at redshift z=0.225 with 5 years of observation from LIGO-Virgo-KAGRA network of detectors. Employing Cosmic Explorer and Einstein Telescope for just 1 year yields the best precision of about $1.62\%$ for BNSs and $2\%$ for NSBHs at redshift z=0.5 on the evolution of the frictional term, and a similar precision up to z=1. This measurement can discover potential deviation from any kind of model that impacts GW propagation with ongoing and upcoming observations.


INTRODUCTION
The General Theory of Relativity (GR) has stood as a cornerstone of our understanding of gravity for over a century.It elegantly explains how massive objects warp the fabric of spacetime, causing other objects to move along curved trajectories.However, in extreme astrophysical environments such as black holes and the early universe, the GR faces challenges in its applicability (Thorne 1995;Sathyaprakash & Schutz 2009).Recently, a groundbreaking development occurred with the direct detection of gravitational waves (GWs) by the LIGO-Virgo collaboration (Abbott et al. 2016b,c;Scientific et al. 2017;Abbott et al. 2017a;Goldstein et al. 2017).This achievement has opened a new avenue for testing GR, particularly in the context of compact binary systems existing in relativistic regimes.These systems include Binary Neutron Stars (BNS), Neutron Star-Black Hole pairs (NSBH), and Binary Black Holes (BBH).These compact binaries offer unprecedented opportunities to delve into the fundamental nature of gravity across a vast range of mass scales and cosmological distances.Ground-based detectors like LIGO (Aasi et al. 2015), Virgo (Acernese et al. 2014) and KAGRA (Akutsu et al. 2021) have played pivotal roles in this endeavor.Furthermore, the prospect of future space-based detectors such as LISA (Amaro-Seoane et al. 2017) and ground-based detectors such as LIGO-Aundha (LIGO-India) ★ samsuzzaman.afroz@tifr.res.in† suvodip.mukherjee@tifr.res.in(Saleem et al. 2021) 1 , Cosmic Explorer(CE) (Reitze et al. 2019), and the Einstein Telescope(ET) (Punturo et al. 2010) holds even greater promise for advancing our understanding of gravity and rigorously testing GR's predictions.Through these combined efforts, we aim to refine and expand our comprehension of the intricate interplay between gravity and the cosmos.While various models proposing modifications to the theory of gravity have been put forth to address different scenarios, this study is dedicated to testing GR in the propagation of GW through a model-independent, data-driven approach.
GWs and EM signals' propagation speed and luminosity distance offer insight into modified gravity theories.Factors like graviton mass, frictional term(due to running effect of Planck mass) and anisotropic source term contribute to deviations (Deffayet & Menou 2007;Saltas et al. 2014;Nishizawa 2018;Belgacem et al. 2018a,b;Lombriser & Taylor 2016;Lombriser & Lima 2017).Precise measurements enable rigorous testing of these alternative theories.Employing bright sirens that emit both GW and EM radiation provides a robust avenue for conducting meticulous assessments of GW propagation speed, graviton mass, and the frictional term((z), see in Equation .2).The prompt measurement of the electromagnetic counterpart, occurring approximately 1.7 seconds after the BNS event GW170817, has facilitated the imposition of exceedingly precise constraints on the speed of gravitational wave propagation and the mass of the graviton (Abbott et al. 2017a(Abbott et al. ,b, 2019a)).For LVK, the BNS case considers the total events corresponding to a local merger rate  0 = 500 Gpc −3 yr −1 , while for NSBH,  0 = 20Gpc −3 yr −1 .For CE & ET, the BNS scenario assumes  0 = 100 Gpc −3 yr −1 , and for NSBH,  0 = 20 Gpc −3 yr −1 are used.The plot illustrates the enhanced accuracy in measuring the function F () using future detectors compared to our current instruments is possible in the redshift range 0.2 to 1 which is an interesting cosmic epoch when the acceleration of the Universe becomes important.
Various model-dependent parametrizations of () exist in the literature, and one such example is the (Ξ 0 , ) parametrization (Belgacem et al. 2018b) and its detectability is shown for LVK (Leyde et al. 2022) and LISA (Baker & Harrison 2021).
This parameterization depends on two positive parameters, denoted as Ξ 0 and n.When Ξ 0 is set to 1, it corresponds to GR.There exist several models of modified gravity, each providing predictions for the frictional term in various scenarios.One prominent category is Scalar-Tensor theories.These theories encompass a scalar degree of freedom, the dynamics of which play a pivotal role in the evolution of the universe.Among the most well-known scalar-tensor theories of gravity are the Brans-Dicke theory (Brans & Dicke 1961), f(R) gravity (Hu & Sawicki 2007) and covariant Galileon models (Chow & Khoury 2009).These belong to the broader class known as Horndeski theories (Arai & Nishizawa 2018).Beyond Horndeski theories, we have generalized galileon models (Gleyzes et al. 2014a), Degenerate Higher Order Scalar-Tensor (DHOST) theories (Frusciante et al. 2019;Crisostomi et al. 2016;Achour et al. 2016).These have been constructed and represent, thus far, the most general scalar-tensor theories that propagate a single scalar degree of freedom alongside the helicity-2 mode of a massless graviton.Additionally, there are targeted scalar-tensor theories specifically designed to remain compatible with observed phenomena (Creminelli & Vernizzi 2017;Gleyzes et al. 2014b).There are several other well-known modified gravity exist in the literature such as f(Q) gravity, f(T) gravity (Cai et al. 2016), bigravity De Felice et al. (2021), gravity in some extra dimensions (Dvali et al. 2000;Yamashita & Tanaka 2014;Corman et al. 2021) and so on.In addition, there exist additional parameterizations (Belgacem et al. 2018b;Romano 2022), with the most significant being the polynomial-exponential and exponential parameterizations.
These models exhibit variations from one another (Belgacem et al. 2019) from the power-law form which cannot be captured by Ξ 0 and .It is worth mentioning that for this particular parametrization, the fitting formula for (z) becomes less accurate in both low and high redshifts in comparison to some of the models (see for example (Deffayet & Menou 2007)).Additionally, the intent behind parameterizing them is to facilitate a seamless transition from a unitary value at low redshifts, where the cumulative impacts of modified gravity wave propagation haven't had adequate time to deviations from GR, to a consistent value represented by Ξ 0 at higher redshifts.If there is a loss of gravitons into the extra dimensions, it would result in a reduction in the apparent brightness of a GW source.Consequently, we would observe that the luminosity distance for gravitational waves exceeds that of electromagnetic waves (Dvali et al. 2000;Yamashita & Tanaka 2014;Corman et al. 2021).In intermediate redshifts, the integrated effect of () may diverge from a simple power-law form in reality.Consequently, it becomes crucial to test GR in a model-independent way and measure the integrated effect of () as a function of the data.To address these issues, we introduce a novel model-independent function, denoted as F (), defined in Equation 4. This function F () is specifically designed to account for any deviation from GR.It serves as a metric for assessing the disparity between the distances denoted D GW l (z) and D EM l (z) so in standard GR this F ()=1.The model-independent reconstruction technique offers significant advantages by capturing both parametric and non-parametric deviations across any redshift range.This approach efficiently overcomes the constraints encountered in specific models throughout varying redshifts.
Figure 1, highlighting the enhanced precision in measuring the redshift-dependent variation of the Planck mass (F ()) for NSBH and BNS systems.Notably, the precision achieved is approximately 7.9% for BNS systems at redshift 0.075 and 10% for NSBH systems at redshift 0.225 over a 5-year observation period from the LVK network of detectors.Furthermore, employing CE & ET for just 1 year yields remarkable precision, reaching about 1.62% for BNS systems and 2% for NSBH systems at redshift 0.5, assuming a fixed Hubble constant.However, when considering a varying Hubble constant, the error bar increases to approximately 2 to 2.2 times that of a fixed Hubble constant.Figure 1 also emphasizes that the measurement of F () is particularly improved in the redshift range of 0.2 to 1.0.This improvement is attributed to two main factors: signal-to-noise ratio (SNR) and the number of events.For a more detailed analysis, refer to Section 6.In the (Ξ 0 , n) frictional term parametrizations, Ξ 0 is set to 1 as the fiducial value.Notably, parameter which controls the redshift evolution denoted by n lacks definition at this point.Measuring n becomes challenging as Ξ 0 nears 1 (Belgacem et al. 2018b).Model-independent parametrization enables redshift scaling for F(z), even for the GR value, in contrast to parametric forms.
The paper is structured as follows: In Section 2, we delve into the propagation of GW beyond GR.Subsequently, we present our proposed model-independent, data-driven reconstruction of the variation of the effective Planck mass(frictional term).In Section 3, we shift our attention to the astrophysical population of GW sources and explore their detection using both current and future detectors.Section 4 is dedicated to a comprehensive discussion on the modelindependent measurement of the BAO scale from the galaxy power spectrum.In Section 5, we provide a detailed exposition of the formalism underpinning our work.Moving forward, Section 6 encapsulates the main findings of our methodology, while Section 7 serves as the culmination, summarizing our key discoveries and engaging in a discourse on future prospects.

NON-PARAMETRIC RECONSTRUCTION OF DEVIATION FROM GENERAL RELATIVITY: FORMALISM
The propagation of GW in the fabric of spacetime with speed , as described by GR, can be mathematically expressed as follows: where ℎ (+,×) represents the GW strain for plus (+) and cross (×) polarizations, the prime symbol denotes the derivative with respect to the conformal time (), and H represents the Hubble parameter in comoving coordinates.This equation serves as a foundational framework for examining the propagation of GWs and evaluating the validity of GR.
In the formulation of modified theories of gravity, several additional parameters come into play.Firstly, there's the frictional term denoted by (), representing the influence of friction in the theory.Additionally, parameters such as the speed of GW propagation, the graviton mass, and the anisotropic stress term are considered.Within the framework of GR, all these additional parameters typically assume a fiducial value of zero, except for the speed of gravity, which corresponds to the speed of light ().Therefore, in the context of testing GR through propagation, it's essential to scrutinize these additional parameters to ascertain whether they deviate from the values predicted by GR.
Recent measurement from GW170817 event, have provided stringent constraints indicating that the graviton mass is zero, and thus the speed of GW propagation is equivalent to the speed of light () (Abbott et al. 2017b(Abbott et al. , 2019a)).Modifications to the anisotropic stress term affect the phase of binary waveforms.Recent observations of BBH mergers, such as GW150914 and GW151226, have set some limits on such modifications, although they are not currently very stringent Abbott et al. (2016a).
As a consequence, () becomes the primary parameter of interest for our study.Consequently, the propagation of GW within the framework of the modified theory of gravity of our interest is described by the equation The presence of the extra frictional term (z) introduces a notable alteration in the amplitude of the GW signal received from a source situated at cosmological distances.This distinction is particularly intriguing as it implies that the luminosity distance measured with standard sirens deviates in principle from that obtained using standard candles or other electromagnetic probes such as CMB or BAO (Amendola et al. 2018;Matos et al. 2023).This distinction arises due to the incorporation of the frictional term within the framework of modified gravity theory.In contrast, standard electromagnetic theory assumes no such frictional term, ensuring that the luminosity distance measured with standard sirens remains unaffected as redshift changes.This discrepancy presents a compelling opportunity for the detection of modified gravity, potentially serving as a "smoking gun" signature of such modifications.So as a consequence the GW luminosity distance D GW l is related to the electromagnetic wave-based luminosity distance D EM l by an exponential factor, where the exponent involves integrating (z ′ )/(1 + z ′ ) with respect to  ′ (Belgacem et al. 2018a).Consequently, the expression becomes In our investigation, we focus on tensor perturbations, where the impact is encoded in the non-trivial function D GW l (z)/D EM l (z).In this study, we propose a non-parametric reconstruction of the deviation from GR as The term F () can capture any deviation from GR as a function of redshift.The reconstruction of this quantity from observation of D l GW (z) and D l EM (z) can capture any modified gravity models which predict a modification in the GW propagation (Belgacem et al. 2019).Alternatively, a parametric form using Ξ 0 and 1+z) n is also used to capture all those models which predicts a power-law modification with redshift of the GW propagation effect (Belgacem et al. 2018b).
In the realm of binary systems, GW offers a unique method for determining luminosity distance, denoted as D GW l .This distance measurement is derived from GW events.While the luminosity distance can also be determined using the standard ΛCDM cosmology, such an approach is inherently model-dependent.To obtain a modelindependent data-driven inference of luminosity distance, we utilize EM luminosity distance calculations derived from various length scales.
Data-driven inference of EM luminosity distance: A data-driven approach to infer the EM luminosity distance can be made from the angular diameter distance  EM  () measured using BAO angular peak position ( BAO ()) from large scale structure observation at different redshifts (Peebles 1973), and using the distance duality relation which connects the angular diameter distance with the EM luminosity distance by  EM  () =  EM  ()/(1 + ) 2 .The position of the BAO peak in the galaxy two-point correlation function can be ).These measurements are obtained from different independent observables such as CMB observations, galaxy correlation functions, and gravitational wave sources.The combination of the sound horizon and the BAO scale gives the EM angular diameter distance  EM  , and by using the distance duality relation we can make a data-driven test of GR as a function of redshift using Equation 6.
inferred directly from observation.The angular scales are related to the comoving sound horizon until the redshift of drag epoch (  ∼ 1020) as By replacing  EM  in Equation ( 4) with the  BAO () and by using the distance duality relation for EM probes, one can write the above equation as (Mukherjee et al. 2021a) The comoving sound horizon   at the drag epoch can be related to the comoving sound horizon  * at the redshift of recombination ( * = 1090) by   ≈ 1.02 * .The quantity  * is inferred from CMB observation, using the position of the first peak in the CMB temperature power spectrum  * =  * /(1 +  * )    ( * ) (Spergel et al. 2007;Komatsu et al. 2011;Hinshaw et al. 2013;Aghanim et al. 2020a).
In the above equation,    () is measured from GW sources,  BAO is measured from galaxy two-point correlation function, and   is inferred from the CMB temperature fluctuations.The redshift to the GW source can be inferred from the EM counterpart for a bright stand siren.As a result, the product of all these observable quantities will lead to an identity value at all redshifts if GR and the fiducial ΛCDM model of cosmology is the correct theory.However, if there is any departure from these models, then we can measure a deviation from one with redshift.It is important to note that the value of   depends on the redshift above   = 1020.Though it depends on the model of cosmology, it is a redshift-independent number of a constant value for the low redshifts  < 2. So, even if there is an inaccuracy in the inference of true   , the value of F () will vary by an overall normalization.But the reconstructed redshift evolution of F () will remain the same.In one of the later sections, we will show results for both cases, (i) only F () and (ii) F () and  0 together.The second case will make it possible to marginalize over any inaccurate inference of sound horizon   due to incorrect inference of the Hubble constant.
Several alternative models of GR predict deviation in different natures of the acceleration of the Universe that are predicted from the ΛCDM.As a result, this data-driven approach can make it possible to explore both deviation from GR and  = −1 equation of state of dark energy.The use of different observational probes was done previously (Mukherjee et al. 2021a) for a parametric form of deviation from GR (using Ξ 0 and ) valid for a class of model.However, the nonparametric form using F (z) proposed in this analysis, can make a redshift-dependent reconstruction of any deviation from GR using the multi-messenger observations.The redshift can be deduced through the EM counterparts associated with these bright sirens.Ongoing and upcoming missions dedicated to identifying EM counterparts include DESI (Adame et al. 2023), Vera Rubin (Abell et al. 2009), Fermi (Bissaldi et al. 2009), Swift (Burrows et al. 2005), HST (Scoville et al. 2007), Roman space telescope (Eifler et al. 2021), Chandra (Vikhlinin et al. 2009), VLA (Benz & Guedel 1989), ZTF (Bellm 2014), Astrosat (Agrawal 2006).A comprehensive article on the required multi-messenger observations can be found here (Ronchini et al. 2022).For the measurement of BAO scale, ongoing and upcoming missions are committed to spectroscopic surveys.Notable projects in this regard include eBOSS (Alam et al. 2017), DESI (Adame et al. 2023), Vera Rubin (Abell et al. 2009), Euclid (Laureĳs et al. 2011), andSPHEREx (Doré et al. 2018).CMB measurements from the ongoing and upcoming missions include ACTPol (Thornton et al. 2016), Simons Observatory (Abitbol et al. 2019), SPT-3G (Sobrin et al. 2022) and CMB-S4 (Abazajian et al. 2016).

Modelling of the astrophysical population of GW source
To accurately assess deviations from the GR model, it is crucial to rely on real GW events.The total number of GW events is contingent upon the merger rates, and mass population of the GW sources.We will explain below the astrophysical models used in the analysis.
Merger rate: In this study we use a delay time model of the binary mergers (Ochsner et al. 2012;Dominik et al. 2015).In the delay time model, the merger rate is described in terms of the delay time distribution, denoted as   .The delay time refers to the elapsed time between the formation of stars that will eventually become black holes and the actual merging of these black holes.It is important to note that the time delay is not uniform across all binary black holes but instead follows a specific distribution.This distribution function accounts for the variations in the delay time and is defined as follows: The delay time is given by   =   −   , where   and   are the lookback times of merger and formation respectively (Karathanasis et al. 2023).So the merger rate at redshift z can be defined as The parameter  0 represents the local merger rate, indicating the frequency of mergers at a redshift of z=0.According to the study in Abbott et al. (2023), the estimated values of  0 for BNS systems vary between 10 Gpc −3 yr −1 and 1700 Gpc −3 yr −1 .For NSBH systems, the  0 values range from 7.8 Gpc −3 yr −1 to 140 Gpc −3 yr −1 .
In this study, our analysis centers exclusively on the propagation of gravitational waves, employing the delay time merger rate model to sample binary systems.The considered modified theory of gravity does not influence binary system formation, ensuring the merger rates remain unaffected.
In our approach, we adopt a standard local merger rate of  0 = 20 Gpc −3 yr −1 for NSBH systems.For BNS systems, we consider four different scenarios with  0 values of 100, 200, 300, and 500 Gpc −3 yr −1 respectively.This methodology allows us to examine the effects of different local merger rates on our results.
The numerator of the expression involves the integration over redshift   from  to infinity, where   (  |   ,    , ) is the delay time distribution,    (  ) is the star formation rate, and    is the jacobian of the transformation.The star formation rate(   ()), is determined using the Madau & Dickinson (2014) star formation rate.
The total number of compact binary coalescing events per unit redshift is estimated as where   indicates the total observation time,    corresponds to the comoving volume element, and () denotes the merger rate (Karathanasis et al. 2022).We consider the delay time merger rate with a specific delay time   = 500Myrs and a power-law exponent of  = 1.However, it is regrettable to note that not all of these events can be detected due to the limitations imposed by the sensitivity of our current detectors.To determine which events are detectable, the calculation of the matched filtering signal-to-noise ratio (SNR) plays a crucial role.The SNR serves as a measure of the strength of the gravitational wave signal relative to the background noise.Only those events with a matched filtering SNR greater than or equal to a predetermined threshold SNR ( TH ) can be reliably detected.(Maggiore 2007).
For a GW emitted by an optimally oriented binary system, the optimized SNR, denoted as , is defined as follows (Sathyaprakash & Dhurandhar 1991;Cutler & Flanagan 1994;Balasubramanian et al. 1996;Nissanke et al. 2010) Here,   (  ) represents the power spectral density of the detector.The function ℎ(  ) corresponds to the Gravitational Wave (GW) strain in the restricted post-Newtonian approximation and is defined for plus (+) and cross (×) polarization as (Ajith et al. 2008): In this expression, the symbol  represents the symmetric mass ratio.The term   signifies the chirp mass of the system.The variable   denotes the luminosity distance.The constant c represents the speed of light in a vacuum.I + = (1 + cos 2 )/2 and I × = cos  depends on the inclination angle .Finally, Ψ(  ) stands for the phase of the waveform.
For the inspiral phase, this formula has proven to be accurate, particularly given the moderate mass ratio in our case.We employed this waveform model exclusively for SNR calculations.Given the moderate mass ratio, the impact of higher-order modes on SNR is minimal.Since SNR is tied to the number of events, the inclusion of higher-order Post-Newtonian corrections does not affect the individual events' error on F ().The parameter estimation for the selected sources, conducted with Bilby (Ashton et al. 2019), utilized the IMRPhenomHM waveform model (Kalaghatgi et al. 2020).
However, the signal detected by a GW detector h det is a complex interplay of several variables, including the detection antenna functions ( + ,  × ), and can be expressed as: here  + and  × are the antenna functions defined as follows (Finn & Chernoff 1993) where  and  define the location of the source in the sky, and  is related to the orientation of the binary system with respect to the detector.Consequently, the matched filtering SNR () takes the form (Finn 1996) where Θ 2 ≡ 4  2 + (1 + cos 2 ) 2 + 4 2 × cos 2  .Averaging over many binaries inclination angle and sky positions, Finn (1996) showed that Θ follows a distribution otherwise.As we have primarily focused on BNS and NSBHs, we did not include LISA in our analysis.Despite its ability to reach deeper into the redshift and provide more precise parameter measurements, the range of masses of compact objects considered in our analysis is comparatively lighter than the masses detectable from LISA.The LISA bright siren analysis will be presented in a future work.
For the LVK system, a threshold SNR of  Th = 12 is chosen.For CE & ET, the threshold SNR is adjusted to enable the detection of events up to a redshift of 3 (which corresponds to an SNR of 25 and 55 for BNS and NSBH respectively.).A redshift bin size of Δ = 0.025 is adopted, and within this bin, the total number of events is computed using Equation 9. Events are generated based on distances inferred from redshifts, the parameter Θ, and object masses from their respective distributions.
For NSBH systems, the LVK's 5-year operation is projected to yield approximately 30 detectable events (all of them need not have a detectable EM counterpart).However, with a one-year observation period, the CE & ET are expected to identify around 600 events.In the case of BNS systems, the LVK's 5-year operation is anticipated to result in approximately 25, 39, and 60 detectable events for local merger rates ( 0 ) of 200, 300, and 500 Gpc −3 yr −1 respectively.On the other hand, CE & ET are expected to detect approximately 4800 events in a one-year observation period with  0 =100.In Figure 3 and Figure 4, we present the projected total and detectable coalescing events involving BNS and NSBH systems within specific redshift intervals (Δ).To obtain these estimates, we first determine the number of merging events for each redshift bin using Equation 9.For each event, we employ the inverse transform method to derive the primary mass, secondary mass, and Θ from their respective probability distributions.The primary mass is generated within the range of 10 ⊙ to 20 ⊙ , the secondary mass within 1 ⊙ to 2 ⊙ , and Θ within the interval 0 to 4.
Incorporating redshift information for each event, we calculate the corresponding luminosity distance and SNR using Equation 14 for a single detector, which takes into account the orientation and inclination angle.This approach enhances the realism of our analysis by considering not only optimally oriented binary systems but also accounting for various orientations and inclination angles, providing a more comprehensive understanding of the system dynamics.However, the parameter estimation of the sources is performed using Bilby.
The total SNR denoted as  total , is then determined by combining individual SNRs from individual detectors using the formula This detection capability is evaluated using both the CE & ET and LVK network of GW detectors over a defined observation period,   .The depicted curve in the graph illustrates the distribution of these events up to a redshift of z = 4.This representation offers valuable insights into the temporal occurrence of these mergers throughout cosmic history and the likelihood of their successful detection using current and future planned instrumentation.

GW Source Parameter Estimation using Bilby
We initiate the parameter estimation process for these discernible sources by employing the Bilby (Ashton et al. 2019) package, which yields realistic posterior distributions on the GW luminosity distance marginalized over the other source parameters.The mass of the GW sources and the number of sources with redshift are drawn as described in the last section.For the remaining source parameters, such as the inclination angle (), polarization angle (), GW phase (), right ascension (RA), and declination (Dec), we sample uniformly, considering a non-spinning system.With these injected parameters, we proceed to generate a GW signal using the IMRPhenomHM (Kalaghatgi et al. 2020) waveform model.This model incorporates higher-order modes that can alleviate the degeneracy between the luminosity distance(  ) and inclination angle() for unequal mass sources.By constraining the priors of all other parameters to delta functions 2 , we generate posterior distributions for  1 ,  2 ,   and .Among these, the posterior of   is particularly crucial for our study, as it will be used for the inference of F ().Additionally, the detection probability of the electromagnetic counterpart is contingent on the value of the inclination  (Chase et al. 2022;Ronchini et al. 2022).

BARYON ACOUSTIC OSCILLATION SCALE FROM GALAXY POWER SPECTRUM
Inference of BAO from galaxy two-point angular correlation function: The BAO arises from the intricate interplay between radiation pressure and gravity during the early stages of the Universe.Previous works on this and corresponding measurements can be found here (Peebles 1973;Crocce et al. 2011). As where, the coefficients   represent the parameters of the model,   denotes the width of the BAO feature, and    corresponds to the best-fit value of the angular BAO scale.This signature serves as a robust standard ruler, enabling independent measurements of the angular diameter distance denoted as   () (Peebles & Hauser 1974;Davis & Peebles 1983;Hewett 1982;Hamilton 1993;Landy & Szalay 1993).
Modelling of the 2PACF: Theoretically, one can model the 2PACF () as where   is the Legendre polynomial of the   ℎ order and   the angular power spectrum, which can be expressed in terms of the three-dimensional matter power spectrum P () as follows (Peebles 1973;Crocce et al. 2011) 2 We have assumed here that the RA and Dec are inferred from the EM counterpart.where the quantity   () is defined as: where () is the galaxy selection function and   is the spherical Bessel function of the l th order, and r(z) is the comoving distance.
To improve the convergence of the integral for the angular power spectrum, a damping factor  (− 2 / 2    ) is introduced (Anderson et al. 2012).Throughout the calculations, we adopt a fixed value of 1/    = 3 Mpc/h, which has no significant impact on the angular power spectrum for the scales of interest.
In this study, we calculate the matter power spectrum from the module CAMB (Lewis & Challinor 2011), and the selection function on redshift is formulated as a normalized Gaussian function.The standard deviation of this Gaussian function is set to be 5% of the mean value.The anticipated photo-z error for the Vera Rubin Observatory is described by the expression 0.03(1 + ), as referenced in (Mandelbaum et al. 2018).Figure 5 visually demonstrates how the product of  2 () changes with angular separation at a specific redshift  = 0.125.The prominent peak at  ∼ 16.0 is a significant observation, highlighting the characteristic scale associated with BAO, which is crucial for our study.
The covariance of (), denoted as    ′ , is modelled as Here, 1/ is the shot noise which is related to the number density of galaxies  or, more precisely, the number of objects per steradian, and    represents the fraction of the sky that is covered by the survey or observation (Cabré et al. 2007;Dodelson & Schmidt 2020).Figure 6 illustrates how the covariance matrix, derived from the linear matter power spectrum, changes concerning angular separations at a specific redshift.The increase in the covariance matrix as both  and  ′ rise suggests a correlation between these angular scales.Figure 7 provides a comprehensive view of how the BAO scale and its error (calculated by using Equation 15) evolve with redshift, emphasizing the impact of photo-z inaccuracies.The consistent relative error observed across measurements serves as a key finding, shedding light on the reliability of BAO scale determinations in the given redshift range.In future  work on the application of this technique to data, one can estimate the covariance matrix from simulations.Several detectors are set to measure the BAO scale with impressive accuracy, delving into the depths of the universe's redshift.Among these detectors, DESI (Adame et al. 2023) stands out as a groundbreaking 5-year experiment on the ground, specifically designed to investigate BAO and the changes in cosmic structures as seen through redshift-space distortions.Covering a vast area of 14,000 square degrees in the sky (  sky = 0.3), DESI aims to achieve extremely precise measurements of the BAO feature.The goal is to surpass an impressive precision of 0.5%, focusing on redshift intervals 0.0 <  < 1.1, 1.1 <  < 1.9, and 1.9 <  < 3.7 (Adame et al. 2023).Additionally, Euclid (Laureĳs et al. 2011) and the Vera Rubin Observatory (Abell et al. 2009;Ivezić et al. 2019), covering a vast sky area of 18,000 square degrees will be able to make precise measurements of BAO up to high redshift.

Case with the fixed Hubble constant 𝐻 0
The intricate connection between the EM wave luminosity distance at a specific redshift (z) and key cosmological parameters, including the BAO scale (  ) and the sound horizon (  ), is succinctly captured by the equation This equation serves as a focal point for the frictional component under consideration.In the context of (F ()), it establishes a connection with three distinct length scales (i)GW luminosity distance: The luminosity distance for detectable GW events which is determined using Bilby.(ii) BAO scale(  ): We utilize the BAO scale, derived from the 2PACF using the Equation 15.This method is not bound to any particular survey, affording us the flexibility to apply our methodology across a range of surveys without being restricted to a specific one.

𝑙
|F ()  ,    ,   ,   ) corresponds to the likelihood function.(  ) and (   ) denote the prior probability distributions for the sound horizon (  ) and BAO scale respectively.The term Π(F ()) represents the prior distribution for F (), which encapsulates our knowledge or assumptions about the frictional term before incorporating the measured data.We adopt a flat prior on the frictional term, F () from 0.1 to 2. For this investigation, we make specific assumptions regarding the probability distribution of redshift(z) and prior to the frictional term F ().In this analysis, we consider bright sirens only.So, we assume that the host galaxy redshift of the GW sources is measured spectroscopically.In this case, it implies that we have precise knowledge or assume that the redshift of the sources is concentrated at a particular value.
In Figure 8, we presented the measurement of F (), with a fixed value of  0 for both the CE & ET and the LVK system, focusing on a single event.This plot illustrates a notable trend where the measurement error increases with higher redshifts.Additionally, it discerns that, at the same redshift, the error on F () is more pronounced for BNS systems compared to NSBH systems, attributed to the lower SNR observed for BNS events.

Case with varying the Hubble constant 𝐻 0
In the preceding sections, our analysis primarily focused on evaluating the frictional term denoted as F () in Equation 20, assuming a fixed Hubble constant ( 0 ).However, in this section, we jointly infer both F () and  0 .This extension enables the simultaneous estimation of both the frictional term F () and  0 , derived from sound horizon measurements that remain constant across redshift, serving as an overall normalization on F () measurements.
In this study, our primary focus is to constrain the frictional term F ().Simultaneously, while measuring this term, we also estimate the value of the  0 .Leveraging the capabilities of our future planned detectors, we can extend our observations to deeper redshifts, allowing us to measure  0 either in conjunction with F () or independently.This innovative approach holds the potential to address the Hubble tension (Abdalla et al. 2022), an observed discrepancy between the locally measured value of the  0 (Riess et al. 2022) and the value inferred from CMB observations (Aghanim et al. 2020b).
Employing a hierarchical Bayesian analysis, our unified framework efficiently estimates these two parameters by incorporating diverse observational datasets, including GW measurements, CMB measurements, and large-scale structure surveys, as discussed in previous sections.The hierarchical organization of these datasets ensures robust integration, effectively handling their unique uncertainties and systematics.The resulting correlated constraints significantly contribute to our understanding of cosmological phenomena, addressing potential discrepancies between datasets and revealing insightful connections among fundamental cosmological parameters.
The equation representing the joint posterior distribution on F (z) and  0 for these sources is given as (Mukherjee et al. 2021a) In this equation, P(F (z), H 0 |{d GW }, {z}) denotes the joint posterior probability density Function of F (z) and  0 .The term L (    |F ()  ,   0 ,     ,   ,   ) corresponds to the likelihood function.The term Π( 0 ) denotes the prior distribution for  0 , encapsulating our knowledge or assumptions regarding the Hubble constant before incorporating measured data.All other terms represent the same entities as defined in Section 5.1.In our hierarchical Bayesian framework, we adopt a flat prior assumption, indicating a lack of prior knowledge about the parameters for both the frictional term (F ()) and the cosmic expansion rate ( 0 ).The parameter ranges are defined from 0.1 to 2.0 for F () and from 40 km/s/Mpc to 100 km/s/Mpc for  0 .
In Figure 9, we present measurements of F () with varying values of  0 for both the CE & ET and the LVK system, focusing on the total number of events within each redshift bin (mentioned in the parenthesis in the upper x-axis in the plot).This plot reveals a discernible trend, indicating an increase in measurement errors with higher redshifts.Notably, due to the higher likelihood of BNS events compared to NSBH events, we can achieve a more accurate measurement of F () for BNS events than for NSBH events.
Detecting EM counterparts in BNS and NSBH mergers poses challenges related to prompt collapse, ejecta mass, and orientation.The visibility of collimated jets is influenced by their opening angle, with wider angles increasing detectability.Distance, instrument sensitivity, and rapid event localization further complicate detection.In astrophysics, terms like "structured jets" and "off-axis viewing" describe these phenomena (Rees & Mészáros 1998;Ramirez-Ruiz et al. 2001;Kumar & Piran 2000).Coordinated efforts between GW and EM observatories are vital for insights into these cosmic collisions.Despite progress, challenges persist, requiring proximity, orientation, and sensitivity considerations.Success relies on refining theoretical models, combining GW and EM observations, and technological advancements.Future endeavors aim to optimize multi-messenger observations, deepen theoretical comprehension, and amplify the detectability of NSBH merger EM counterparts.
In our study, we present cases for both the CE & ET and LVK systems.Specifically for CE & ET, we showcase outcomes for up to redshifts  = 2 for different numbers of detectable events.In the case of NSBH binaries, we consider a total of 30 events with a local merger rate ( 0 ) set at 20 Gpc −3 yr −1 in the LVK system.These events span redshifts from 0.125 to 0.425.Similarly, for BNS systems, we analyze approximately 25, 40, and 60 events with  0 set to 200 Gpc −3 yr −1 , 300 Gpc −3 yr −1 , and 500 Gpc −3 yr −1 , respectively.The redshifts of these events range from 0.025 to 0.225 over a 5-year observation period.In Figure 3 and Figure 4, the number of detectable events as a function of redshift with different  0 and   are presented for the LVK and CE & ET systems.
Figure 10 and Figure 11 provide a comprehensive comparison of F () measurements with fixed  0 and varying  0 across different scenarios for both LVK and CE & ET.These scenarios include the detection of a single event, a quarter of the events, half of the events, and the detection of all events.The aim is to demonstrate how the precision of measurements, ranging from percentages to sub-percentages, can be achieved using future detectors.This comparative analysis sheds light on the impact of event detection rates on the precision of F () measurements and underscores the potential of advanced detectors like CE & ET to provide increasingly accurate cosmological insights.In this section, we measure F () with a constant  0 that remains uniform across redshifts, serving as an overall normalization factor for F () measurements.Consequently, if an independent precise measurement of F () in our local universe is possible, it will allow us to normalize it to 1, eliminating the dependence on  0 .This normalization proves effective particularly for low redshifts, where the cumulative effects of modified gravity wave propagation have not had sufficient time to accumulate deviations from GR.
The error associated with the frictional term approximately scales as the inverse of the square root of the number of GW sources represented as 1/ √ N GW .This implies that as the number of GW sources increases, the error on the frictional term decreases, thereby improving the precision of the measurement.However, this improvement reaches a saturation point, beyond which an increase in the number of GW sources does not lead to a significant enhancement in the measurement's accuracy.This saturation is primarily due to the presence of error on the BAO scale.The BAO scale error acts as a limiting factor in the precision of the frictional term estimation.Therefore, to achieve further improvement in the estimation of the frictional term with redshift, a more precise measurement of the BAO scale is required.In essence, the precision of the frictional term estimation is a delicate balance between the number of GW sources and the ac-curacy of the BAO scale measurement and GW luminosity distance.Optimizing these factors can lead to significant advancements in our understanding of F () and discovering unexplored territories in fundamental physics.On the GW side, the measurements of F () are solely dependent on    ().Any improvement in the precision of    () will markedly enhance the accuracy of F (), as    () stands out as the predominant source of error for F ().Given that    () is degenerate with other GW parameters, most notably the inclination angle (), measuring the inclination angle from alternative sources (such as EM counterparts) can substantially reduce the error in    () and consequently enhance the precision of F () (Xie et al. 2023).For redshift inference from EM counterparts, the influence of peculiar velocity contamination is particularly significant at low redshifts ( < 0.03) but it is not important for high redshift (Mukherjee et al. 2021b;Nimonkar & Mukherjee 2024).As most of the detected sources are at high redshift, we have ignored this effect.Also, weak lensing of the GW sources can be a source of contamination to the luminosity distance up to a few percent at redshifts above  = 0.5 (Hirata et al. 2010).In comparison to the error bar on luminosity distance, the weak lensing error bar is small.As most of the sources with better luminosity distance measurement are mainly at low redshift ( < 0.5), we are getting maximum constraints around this redshift.So, the weak lensing uncertainty is not a major contamination for these sources.Moreover, weak lensing can provide additional information to constraint non-GR parameter The top and bottom plot relates to possible measurement using one BNS and NSBH source respectively jointly with the Hubble constant.The upper x-axis denotes the associated 1- error bars for each measurement, with the number in parentheses representing the total number of events in that specific redshift bin.The plot highlights the potential for more accurate measurements of F () at significantly deeper cosmic redshifts, particularly with future detectors like CE & ET, in contrast to the limitations posed by our current LVK detector, similar to the case with fixed Hubble constant.Additionally, it underscores that we can achieve more precise measurements of F () for BNS systems compared to NSBH events, given the higher likelihood of BNS events occurring.
space (Mukherjee et al. 2020a,b;Mpetha et al. 2023;Balaudo et al. 2023), which we have not considered in this analysis.

Reconstruction of F (𝑧) and 𝐻 0
By extension of the hierarchical Bayesian framework used in the previous section to measure F () with fixed  0 , we can jointly measure F () and  0 .In our framework, we measure the EM luminosity distance from different length scales, namely the BAO scale, sound horizon distance, and redshift.In Figure 10 and Figure 11, we offer a thorough comparison of F () measurements with both fixed and varying  0 across different scenarios for both the LVK and CE & ET systems.It is noteworthy that when allowing for variations in  0 , the errors in the measurements of F () increased by a factor ranging from 2 to 2.2 compared to measurements with a fixed  0 .
In Figure 1, we highlight the enhanced precision in measuring the redshift-dependent variation of the Planck mass (F ()) for NSBH and BNS systems.The figure emphasizes that the measurement of F () is notably improved in the redshift range of 0.2 to 1.0.The precision of estimating F () depends on accurate measurements of both the BAO scale and    ().Improved accuracy in the BAO scale measurements directly enhances precision in estimating F ().On the GW side, accurate measurements of F () rely solely on    ().Enhancing the precision of    () improves F (), as mentioned in the previous Subsection 6.1.There can also be peculiar velocity corrections and weak lensing contamination at low redshift as described in 6.1 which are not taken care of in this study.

CONCLUSION AND DISCUSSION
In this study, we present a model-independent approach to search for the frictional term, which can capture deviations at any redshift from GR.We propose a model-independent reconstruction of the frictional term as a function of redshift using a data-driven approach by amalgamating three length scales with redshift namely the luminosity distance from GW sources  GW  (), the angular scale of BAO scale from galaxy surveys  BAO (), and the sound horizon   from CMB observations.We show this provides a unique opportunity to scrutinize theories of gravity, specifically to test the presence of a frictional term.In the literature, various model-dependent searches for the frictional term have been conducted.These searches face challenges, especially near redshift z=0 and at large redshifts.Particularly noteworthy are the difficulties encountered when there is a potential loss of gravitons in extra dimensions.Furthermore, in the intermediate redshift range, the frictional term may display oscillatory be-  havior that cannot be adequately captured by these model-dependent searches.
We demonstrate that by integrating BAO measurements from galaxies, sound horizon distances from the CMB power spectrum, redshift information from the EM counterpart, and luminosity distance measurements from GW sources, one can perform a comprehensive measurement of cosmological parameters.Specifically, we highlight the potential to jointly determine the Hubble constant ( 0 ) and the frictional term through this integrated analysis.We have demonstrated the feasibility of measuring the frictional term using our ongoing detectors, such as LVK, across various redshifts.Additionally, we have explored the capabilities of upcoming detectors, specifically CE & ET, in quantifying the frictional term under various scenarios.Our analysis includes cases where a fraction of the event's EM counterpart is detectable, acknowledging the inherent difficulty in detecting EM counterparts.
We demonstrate the feasibility of achieving a precision measurement of the frictional term, ranging from a percent to sub-percent accuracy, as a function of cosmic redshift.This is achievable using data from LVK and, in the future, from CE & ET, despite the limitation of a finite number of bright sirens.The measurement of F () is particularly improved in the redshift range of 0.2 to 1.0.This improvement is attributed to two main factors: SNR and the number of events.The measurements from LVK can be further improved by including LIGO-India (Saleem et al. 2021)  3 .
In the (Ξ 0 , n) parametrizations of the frictional term, the fiducial value for Ξ 0 is set to 1.However, it is important to note that the parameter n is not defined at this fiducial value of Ξ 0 =1.Consequently, measuring the value of n becomes increasingly challenging as Ξ 0 approaches 1 (Belgacem et al. 2018b).By utilizing 1000 bright sirens up to redshift 8, Belgacem et al. (2018b) demonstrated that Ξ 0 can be measured with an accuracy of 0.8% using the Einstein Telescope for a fixed value of n.In our study, employing a model-independent parametrization, we found that with 1 year of observation using CE & ET, we can measure the frictional term with an accuracy of 1.62% and 2% for BNS and NSBH systems at redshift 0.5.Additionally, we can reconstruct the frictional term as a function of redshift without relying on any specific parametrization, as illustrated in Figure 1.This error scales inversely with the square root of the observation time, i.e., as 1/ √︁   / .Our proposed technique extends to dark sirens, where redshift inference becomes necessary.In such cases, we leverage the threedimensional clustering scale of GW sources with galaxies to infer redshift information (Mukherjee & Wandelt 2018;Mukherjee et al. 2020c;Abbott et al. 2021).Furthermore, this methodology holds promise for both bright and dark standard sirens in the mHz range, detectable by future space-based detectors like LISA.In this study, our exclusive focus is on bright sirens, and as a result, we exclude a prominent category of GW sources Binary Black Holes (BBH) due to their unlikely association with EM counterparts.Additionally, there are challenges associated with detecting EM counterparts for BNS and NSBH systems (Chase et al. 2022;Ronchini et al. 2022).One potential method to enhance the accuracy of frictional term measurements is to include dark sirens, which lack EM counterparts.In the future, the integration of LIGO-India and the LISA GW detector, in conjunction with projects such as eBOSS, DESI, Vera Rubin, and Euclid for BAO measurements, along with the application of this technique to dark sirens, holds the potential to substantially improve the precision of measuring the frictional term.

Figure 1 .
Figure1.This plot showcases the precision achievable in reconstructing the redshift variation of the Planck mass captured by F () for both NSBH and BNS systems, employing the current LVK detector with 5 years of observation time and the future CE & ET detectors with 1 year of observation time.For LVK, the BNS case considers the total events corresponding to a local merger rate  0 = 500 Gpc −3 yr −1 , while for NSBH,  0 = 20Gpc −3 yr −1 .For CE & ET, the BNS scenario assumes  0 = 100 Gpc −3 yr −1 , and for NSBH,  0 = 20 Gpc −3 yr −1 are used.The plot illustrates the enhanced accuracy in measuring the function F () using future detectors compared to our current instruments is possible in the redshift range 0.2 to 1 which is an interesting cosmic epoch when the acceleration of the Universe becomes important.

Figure 2 .
Figure 2.This is a schematic diagram demonstrates the measurement of three distinct length scales: the sound horizon (  ), BAO scale (   ), and the GW luminosity distance (

Figure 3 .
Figure 3.Total number of injected events (in blue) and the number of detectable GW events as a function of redshift for BNS sources are shown for CE & ET (in green) and LVK (in black).The total observation periods are mentioned in parentheses.

Figure 4 .
Figure 4. Total number of injected events (in blue) and the number of detectable GW events as a function of redshift for NSBH sources are shown for CE & ET (in green) and LVK (in black).The total observation periods are mentioned in parentheses.
photons and baryons decouple, a specific length known as the sound horizon at the drag epoch (  ), represented by   , leaves a distinct signature on the distribution of galaxies and the power spectrum of CMB anisotropies(Peebles & Yu 1970;Sunyaev & Zeldovich 1970;Bond & Efstathiou 1987;Hu & Dodelson 2002;Blake & Glazebrook 2003).By examining the two-point angular correlation function (2PACF) of galaxy catalogs denoted by (), one can identify the BAO signature, which manifests as a prominent feature resembling a bump.The BAO scale can be derived by fitting the angular correlation function () with a model defined as a combination of a power-law and a Gaussian component as(Xu et al. 2012;Carvalho et al. 2016;Alam et al. 2017)

Figure 5 .
Figure 5.This plot illustrates the relationship between the 2PACF( 2  (  )) is shown as a function of  for a specific redshift, z=0.125 calculated using the linear matter power spectrum.Notably, a distinct bump is observed around  ∼ 16.0, indicating the presence of the BAO characteristic scale   .

Figure 6 .
Figure6.This plot displays the covariance matrix(Cov   ′ ) calculated using linear matter power spectrum as a function of  and  ′ for a specific redshift, z=0.125.It is evident that as both  and  ′ increase, the value of the covariance matrix also increases indicating larger error at the large angular scales.

Figure 7 .
Figure 7.We depict the relationship between the BAO scale along with the associated error (in red), primarily arising from photo-z inaccuracies, as a function of redshift (z) up to z = 0.5.This presentation offers crucial insights into how the BAO scale and its error vary with z, highlighting a consistent relative error across measurements.For visual clarity, the error values on the plot have been magnified by a factor of 2.
(iii) Sound Horizons distance:The sound horizon is measured by analyzing acoustic oscillations in the CMB radiation using several missions such WMAP(Hinshaw et al. 2013), Planck(Aghanim et al. 2020a), ACTPol(Thornton et al. 2016), SPT-3G(Sobrin et al. 2022), and CMB-S4(Abazajian et al. 2016).In particular, the first peak observed in the angular power spectrum of the CMB corresponds to the scale of the sound horizon during the process of recombination.The measured sound horizon value at the drag epoch is found to be 147.09± 0.26 Mpc, indicating a high level of precision from Planck(Aghanim et al. 2020a).We utilize a Hierarchical Bayesian framework to calculate the posterior distribution of F () as a function of redshift for a total of   GW sources.The equation representing the posterior distribution on F (z) for bright sirens is given as(Mukherjee et al. 2021a)(F ()|{  }, {}) ∝ Π(F ())   =1 ∭   (  )    ×    (   )L (    |F (  ),    ,   ,   ).(21)In this equation, P(F (z)|{d GW }, {z}) denotes the posterior probability density function of F (z) given the GW source data   and spectroscopic redshift from EM counterpart of these sources {} for  ∈   events.The measurements of BAO angular scale   at the redshift (z) of the GW source where an EM counterpart of GW source can be measured using the galaxy survey, and it is marginalized over.The term L (

Figure 8 .
Figure 8.This violin plot illustrates the posterior on the reconstruction of the non-GR parameter F () as a function of cosmic redshift is feasible from LVK (in green) and CE & ET (in blue).The top and bottom plot relates to possible measurement using one BNS and NSBH source respectively, both assuming a fixed Hubble constant.The upper x-axis denotes the associated 1- error bars for each measurement.It demonstrates the potential to measure F () with greater precision and at significantly deeper cosmic redshifts using future detectors such as CE & ET, as opposed to the limitations of our current LVK detector.

Figure 9 .
Figure9.This violin plot illustrates the posterior on the reconstruction of the non-GR parameter F () as a function of cosmic redshift is feasible from LVK (in green) and CE & ET (in blue).The top and bottom plot relates to possible measurement using one BNS and NSBH source respectively jointly with the Hubble constant.The upper x-axis denotes the associated 1- error bars for each measurement, with the number in parentheses representing the total number of events in that specific redshift bin.The plot highlights the potential for more accurate measurements of F () at significantly deeper cosmic redshifts, particularly with future detectors like CE & ET, in contrast to the limitations posed by our current LVK detector, similar to the case with fixed Hubble constant.Additionally, it underscores that we can achieve more precise measurements of F () for BNS systems compared to NSBH events, given the higher likelihood of BNS events occurring.

Figure 10 .
Figure10.We show the 1- error-bar in the measurements of F () as a function of redshift and the number of GW events for NSBH events.The comparison is made between scenarios with a constant Hubble constant ( 0 ) and for the case with  0 inferred jointly with F ().All the events at each redshift bin are mentioned in parentheses in the figure legend.

Figure 11 .
Figure11.We show the 1- error-bar in the measurements of F () as a function of redshift and the number of GW events for BNS events.The comparison is made between scenarios with a constant Hubble constant ( 0 ) and those with varying  0 .All the events at each redshift bin are mentioned in parentheses in the figure legend.