A Model-Independent Precision Test of General Relativity using LISA Bright Standard Sirens

The formation and evolution of MBBHs and their co-evolution with host galaxies pose fundamental questions in astrophysics, with implications for understanding the regulation of star formation and the growth of galaxies [65, 66]. While active galactic nuclei (AGN), powered by growing MBBHs, are believed to play a crucial role in these processes, many aspects of MBBHs growth and AGN feedback mechanisms remain poorly understood [67, 68]. One significant challenge in studying MBBHs is related to the difficulty of observing their early evolution at high redshifts. Traditional observational methods are limited to detecting only the most luminous and rapidly growing MBBHs. GWs provide a promising avenue for probing the history of MBBHs, as their propagation through the universe is essentially unobstructed [69]. The future space-based GWs detector LISA will be particularly instrumental in this area. These binaries lie within the LISA band, and LISA’s capabilities will allow for observations deep into high redshifts, shedding light on the formation and evolution of MBBHs [70]. As a consequence, our understanding of the formation and evolution of MBBHs currently heavily relies on theoretical models. Semi-analytic models of galaxy and MBBHs formation and evolution serve as powerful tools for this[71, 72, 73, 74, 75]. These models link the hierarchical formation of dark matter halos to the evolution of galaxies and their MBBHs, capturing a wide range of galaxy formation physics [76, 77, 72, 78]. In these models, the formation of MBBHs begins with the merger of two galaxies embedded in their dark matter halos. The subsequent evolution of MBBHs proceeds through various stages, from separations of 10s-100s kpc down to a few hundred pc, and eventually to separations of $\sim$ 0.001–0.01 pc, where GWs emission drives the MBBHs to coalesce. The dynamical evolution of MBH binaries is influenced by interactions with their stellar and gaseous environments, as well as with other MBBHs through three/four-body interactions [79, 80]. However, fully hydrodynamic simulations, which incorporate the complex physics involved in MBBHs evolution, often suffer from poor resolution and are constrained to relatively small volumes, resulting in limited statistical power [81, 82, 83, 84, 85, 86, 87, 88]. Additionally, these simulations are computationally very demanding. While cosmological simulations have been employed to study the evolution of unresolved MBBHs, they generally necessitate significant post-processing and involve considerable simplifications. The formation and evolution of MBBHs and their co-evolution with host galaxies rely heavily on several key ingredients. These include the initial mass function of the MBH seeds, the delays between galaxy mergers and MBBHs mergers. These factors introduce complexities in modeling the dynamics and growth patterns of MBBHs. The literature presents a variety of seed models, delay models, reflecting the diverse theoretical approaches to these phenomena [71, 72, 73, 74, 75, 89, 90].

Seed Model: There are mainly two classes of seed models in the semi-analytical modeling of MBBHs formation: the Light-Seed (LS) and Heavy-Seed (HS) models. The LS model suggests that MBH seeds are remnants of Population III stars from high-redshift, low-metallicity environments, with masses centered around 300 solar masses and a standard deviation of 0.2 dex. These seeds are located in the most massive halos formed from the collapse of 3.5-sigma peaks in the matter density field at redshifts greater than 15, allowing for mildly super-Eddington accretion rates. In contrast, the HS model posits that MBH seeds form through the collapse of protogalactic disks due to bar instabilities, typically reaching masses of approximately $\mathrm{10^{5}}$ solar masses, and are found in similar high-redshift halos but with accretion rates capped at the Eddington limit. The choice between these models affects the formation and evolution of MBBHs, impacting the observed event rates and characteristics of MBBHs mergers. This differentiation is crucial for understanding the astrophysical processes leading to MBBHs and their subsequent mergers [91, 92, 69, 93, 88]. Recent observations by the James Webb Space Telescope (JWST)[94] provide compelling evidence favoring the HS model [95]. JWST has detected extremely luminous quasars at redshifts greater than 7, suggesting the presence of very massive black holes already in the early universe[96]. The inferred masses of these early quasars are more consistent with the higher initial masses predicted by the HS model, as the rapid formation and growth required to reach such masses challenge the LS model’s predictions. These findings indicate that massive black hole seeds might have formed through more efficient, direct collapse mechanisms as proposed by the HS model[97].

Delay Models: In the study of MBBHs formation, time delays play a crucial role by affecting the evolution and eventual merger of these systems. These delays stem from various mechanisms, including dynamical friction, galaxy mergers, and the internal dynamical evolution of MBBHs within galaxies [72, 73, 98, 20]. The impact of these time delays is modeled through different scenarios, each with distinct implications for merger rates and characteristics[99, 69, 100, 101].

All of these factors and their combinations affect the event rate and the mass model of these MBBHs, which are crucial for understanding the range of scenarios detectable by LISA (for more details, see reference [69, 72]). These elements collectively influence our predictions of the frequency and characteristics of MBBHs merger events. The actual mechanisms of MBH mergers and the corresponding mass models remain largely unknown. The total number of compact binary coalescences per unit redshift is estimated by the equation

$$\mathrm{\frac{dN_{GW}(z)}{dzdt}=\frac{R(z)}{1+z}\frac{dV_{c}}{dz}},$$

(3.1)

where $\mathrm{\frac{dV_{c}}{dz}}$ represents the comoving volume element, and $\mathrm{R(z)}$ indicates the merger rate. To capture the wide range of possible merger rate scenarios, we parameterize the merger rate $\mathrm{R(z)}$ using a combination of two Gaussian distributions and a power law, as described by the following equation.

$$\mathrm{R(z)=A}\begin{cases}\mathrm{\exp\left(-\frac{(z-\mu_{L})^{2}}{2\cdot\sigma_{L}^{2}}\right)}&\mathrm{\text{if }z<\mu_{L}}\\[10.0pt] \mathrm{\exp\left(-\alpha\cdot(z-\mu_{L})\right)}&\mathrm{\text{if }\mu_{L}\leq z\leq\mu_{R}}\\[10.0pt] \mathrm{\exp\left(-\frac{(z-\mu_{R})^{2}}{2\cdot\sigma_{R}^{2}}\right)\cdot\exp\left(-\alpha\cdot(\mu_{R}-\mu_{L})\right)}&\mathrm{\text{if }z>\mu_{R}}\end{cases}$$

(3.2)

where $\mathrm{A}$ is the overall normalization, $\mathrm{\mu_{L}}$ is the mean of the left Gaussian distribution, $\mathrm{\sigma_{L}}$ is the standard deviation of the left Gaussian distribution, $\mathrm{\mu_{R}}$ is the mean of the right Gaussian distribution, $\mathrm{\sigma_{R}}$ is the standard deviation of the right Gaussian distribution, and $\mathrm{\alpha}$ is the decay parameter of the exponential function. This approach is motivated by the need to cover all possible variations in the merger rate. Different scenarios can lead to different distributions of mergers, and a flexible parametric form allows us to represent this diversity. The Gaussian distributions help capture localized features in the merger rate, while the power law accounts for the broader trends over redshift. By configuring these components with distinct parameters, we can match a variety of theoretical models and observational constraints. Similarly, the mass model is described using variable parameters that account for different theoretical combinations of the aforementioned factors. For each unique parameter set in the merger and mass models, the predicted number of observable events varies. By integrating Equation 3.1, we calculate the total number of events, allowing us to predict the observable merger events for LISA under different theoretical scenarios. Figure 2 displays the merger rates for three distinct MBBH models using LISA. The Light Seed model involves light MBBH seeds derived from Light Seed stars, considering the temporal delays between MBBHs and galaxy mergers, thereby providing a more conservative estimate. The Heavy Seed+delay model assumes heavy MBBH seeds originating from the collapse of protogalactic disks with similar delays, while Heavy Seed omits these delays, offering an optimistic scenario for LISA event rates.

Each model is fitted using two Gaussian distributions and a power law function, as defined in Equation 3.2, and is depicted with dashed lines in Figure 2. The dotted points in Figure 2 represent data from semi-analytical simulations, as referenced in[69, 72]¹¹1https://people.sissa.it/~barausse/catalogs/readme.pdf. The solid lines are obtained by fitting these data using Equation 3.2, and all parameter values are listed in Table 1. This figure also illustrates the corresponding total redshifted mass $\mathrm{M_{z}\equiv(1+z)(m_{1}+m_{2})}$ models for these scenarios in terms of the individual masses $\mathrm{m_{1}}$ and $\mathrm{m_{2}}$. In addition to the total mass, another crucial parameter for binary systems is the mass ratio ($\mathrm{q\equiv m_{1}/m_{2}}$, where $\mathrm{m_{1}\leq m_{2}}$). In this analysis we are interested in measuring any deviation from GR on the cosmological scale. So, we particularly focus on the high matched filtering signal to noise ratio (SNR) events in this analysis so that the uncertainty on any inferred deviation from GR is minimal. So, we consider events above SNR of fifty and only a sub-set of high SNR events are assumed to have EM counterparts. As the events with highest matched filtering SNR are those with mass ratio ($\mathrm{q=1}$) equal to one, we limit our analysis only for $\mathrm{q=1}$ in this analysis. Also, as it is currently unclear on how many LISA sources can have EM counterparts, we obtain forecast on the parameter $\mathcal{F}(z)$ for different scenarios of number of events with EM counterpart in this analysis, ranging from about $4\%$ to $75\%$ events in four years with matched filtering SNR above fifty.

However, it is important to note that LISA binaries can detectable for mass ratios less than one as well and a fraction of them can also have EM counterpart. In this analysis we consider only high SNR events and make a pessimistic choice by not considering EM counterpart to $q<1$ sources. For q within the range of $\mathrm{10^{-3}}$ to $\mathrm{1}$, the probability on $\mathrm{q}$ defined as $\mathrm{P(q)}$ in the range of q values $\mathrm{0.5}$ to $\mathrm{1}$ is always at least 50% for flat ($\mathrm{P(q)=}$ constant) and for positive power law distributions with $\mathrm{P(q)\propto q^{\alpha},\,\alpha>0}$ [102]. This percentage decreases if we assume a power law with a negative exponent in q. For other $\mathrm{q}$ distributions, we have shown the results in Appendix B. To select a fraction of sources with the value of mass ratio close to one having EM counterpart, in this analysis we consider scenarios, namely $\mathrm{75\%}$ (fiducial case) and $\mathrm{50\%}$ of the total events after applying the SNR cutoff. We also consider a pessimistic case with only six events detected up to redshift $z=6$ with EM counterpart in four years, which corresponds to a scenario of about $4\%$ events for the Light Seed model and Heavy Seed model, and about $30\%$ events for the Heavy Seed+ delay model.

Galaxy simulations to forecast the EM counterparts for massive binary mergers typically assume a well-understood evolution of galaxies and star formation rates, calibrated from large-scale surveys. However, simulating the formation and distribution of galaxies, especially at high redshifts, comes with significant uncertainties[103, 104, 105]. One of the major challenges arises from the assumptions about gas accretion, feedback processes, and the precise merger history, which can vary significantly across different simulation models. Recent observations from the JWST have further complicated our understanding of the high-redshift universe. JWST has revealed a higher abundance of massive galaxies and earlier star formation than previously predicted, challenging existing models of galaxy formation and evolution. JWST has also detected a variety of galaxies with different dust and gas content, which directly affects the probability of producing an EM counterpart. Galaxies rich in gas are more likely to produce bright EM counterparts during binary mergers. Given these uncertainties, we assume that approximately 75% and 50% of massive binary mergers may have observable EM counterparts, depending on their redshift and environment. We have also shown the results for scenarios with $25\%$ and $10\%$ bright sirens in the appendix A. This assumption reflects the current understanding of galaxy formation processes and the likelihood that such mergers would occur in gas-rich environments, which are more conducive to generating EM radiation. While these fractions are estimates, they align with predictions from theoretical models and observational data, with the emerging JWST results suggesting a broader range of environments that could host potential EM counterparts. The models considered here are related to different theoretical approaches to MBBHs formation and evolution, significantly affecting the predicted total mass distributions and influencing the detectability of merger events by LISA. The plots effectively demonstrate how variations in initial assumptions about MBBHs seed formation and merger delays impact the mass spectrum observable by future missions. It is important to that all of these events cannot be detected due to the limitations imposed by the sensitivity of detectors. The calculation of the matched filtering SNR is essential to determine which events are detectable.

The matched filtering SNR, denoted as $\mathrm{\rho}$, is crucial for detecting GW events as it quantifies the GW signal’s strength relative to the background noise. An event is deemed detectable if its SNR surpasses a predefined threshold, $\mathrm{\rho_{TH}}$. The optimal SNR for a GW emitted by a binary system is defined by the equation [106, 107, 108]

$$\mathrm{\rho^{2}=4\int_{f_{\text{min}}}^{f_{\text{max}}}\frac{|\tilde{h}(f)|^{2}}{S_{n}(f)}\,df},$$

(3.3)

where $\mathrm{\tilde{h}(f)}$ represents the Fourier transform of the GWs strain, while $\mathrm{S_{n}(f)}$ denotes the power spectral density of the detector noise.

For $\mathrm{S_{n}(f)}$, we adopt the semi-analytical forms from [109], as shown in Figure3, and $\mathrm{f_{min}}$ and $\mathrm{f_{max}}$ are the boundaries for the frequency integral. In our study, for the computation of the SNR, we set $\mathrm{f_{max}}$ to the frequency of the innermost stable circular orbit ($\mathrm{f_{ISCO}}$) and $\mathrm{f_{min}}$ to the frequency one year prior to $\mathrm{f_{ISCO}}$, considering our observation period is one year. We do not include the Galactic background noise in the power spectral density. The Galactic background can have an impact on the overall noise, which in turn affects the number of detectable events.

Specifically, $\mathrm{\tilde{h}(f)}$ is expressed as

$$\mathrm{\tilde{h}(f)=h(f)Q},$$

(3.4)

with $\mathrm{h(f)}$ denoting the strain amplitude

$$\mathrm{h(f)=\sqrt{\frac{5\eta}{24}}\frac{(GM_{\text{chirp}})^{5/6}}{D_{L}\pi^{2/3}c^{3/2}}f^{-7/6}e^{-i\psi(f)}},$$

(3.5)

where $\mathrm{M_{chirp}}$ is the chirp mass, $\mathrm{\eta}$ is the mass ratio, $\mathrm{D_{L}}$ is the luminosity distance to the source, $\mathrm{c}$ is the speed of light, and $\mathrm{\psi(f)}$ is the phase of the waveform [110]. The factor $\mathrm{Q}$ depends on the inclination angle and the detection antenna functions $\mathrm{F_{+}}$ and $\mathrm{F_{\times}}$

$$\mathrm{Q=\left[F_{+}\left(\frac{1+\cos^{2}(\iota)}{2}\right)+\iota F_{\times}\cos(\iota)\right]}.$$

(3.6)

The antenna functions are defined by[109]

	$\displaystyle\mathrm{F_{+}}$	$\displaystyle=\mathrm{\frac{1}{2}(1+\cos^{2}\theta)\cos 2\phi\cos 2\psi-\cos\theta\sin 2\phi\sin 2\psi},$		(3.7)
	$\displaystyle\mathrm{F_{\times}}$	$\displaystyle=\mathrm{\frac{1}{2}(1+\cos^{2}\theta)\cos 2\phi\sin 2\psi+\cos\theta\sin 2\phi\cos 2\psi},$		(3.8)

where $\mathrm{\theta}$ and $\mathrm{\phi}$ determine the source’s location in the sky, and $\psi$ relates to the orientation of the binary system with respect to the detector [111]. Consequently, the SNR formula becomes

$$\mathrm{\rho=\frac{\Theta}{4}\left[4\int_{f_{\text{min}}}^{f_{\text{max}}}\frac{h(f)^{2}}{S_{n}(f)}df\right]^{1/2}},$$

(3.9)

where $\mathrm{\Theta^{2}\equiv 4\left(F_{+}^{2}(1+\cos^{2}i)^{2}+4F_{\times}^{2}\cos^{2}i\right)}$. According to [112], averaging over many binaries, the inclination angle, and sky positions, $\mathrm{\Theta}$ follows the distribution

$$\mathrm{P(\Theta)}=\begin{cases}\mathrm{5\Theta(4-\Theta)^{3}/256}&\text{if }\mathrm{0<\Theta<4},\\ \mathrm{0}&\text{otherwise}.\end{cases}$$

(3.10)

In this study, we investigate various merger and mass models, recognizing that our understanding of these models is still incomplete. By examining different scenarios, we demonstrate how different assumptions about the formation and evolution of MBBHs impact the measurements of $\mathrm{\mathcal{F}(z)}$ using LISA sources. To estimate the total number of events, we integrate using Equation 3.1 and determine the number of merging events within each redshift bin, selecting a bin size of $\mathrm{\Delta z=0.4}$. Figure 4 displays the number of events binned by redshift for the Light Seed, Heavy Seed, and Heavy Seed + Delay models. The total number of events for each model, along with the observation years, is summarized in Table 2. For each event, we use the inverse transform method to derive the total mass of the binary system and the parameter $\mathrm{\Theta}$ from their probability distributions. We then generate an equal-mass binary system and calculate the SNR using Equation 3.9.

The events with an SNR exceeding the threshold are selected. Since we only consider a mass ratio $\mathrm{q}$ of 1 when generating the binary sources, we retain $\mathrm{75\%}$ of the events after the SNR cutoff. The number of events selected if only $\mathrm{50\%}$ are retained is also noted and listed in Table 2. In Sec. A, we have also shown the results for different fractions of bright standard sirens detectable from LISA. For comparison, in [102], it was found that, for a mission duration of 5 years with an 80% duty cycle up to a redshift of 10 and an SNR cutoff of 10, there were a few hundred events for both light seed and heavy seed scenarios, and about 30 events for the heavy seed with delay scenario. Here, we consider a more pessimistic choice for SNR and redshift when selecting EM bright sources, as we focus only on high SNR golden events for testing GR (see Sec. 7 and appendix A). The detection probability for the Light Seed model is lower compared to the other two models because its mass distribution peaks at much lower values, resulting in reduced signal strength and, consequently, lower SNR. Figure 5 illustrates the total mass of the two components of the binary system as a function of the redshift of the selected events, with an SNR threshold of fifty.

LISA will consist of a constellation of three spacecraft forming a giant equilateral triangle in space, each containing free-falling test masses. Lasers measure the minute changes in distance between these masses, changes induced by passing gravitational waves. However, a significant challenge in detecting these waves is the overwhelming noise from laser frequency fluctuations [24]. To combat this, LISA uses Time-Delayed Interferometry (TDI), a sophisticated data processing technique that combines the laser measurements in a way that reduces noise by several orders of magnitude, making the subtle signals of gravitational waves detectable. The parameter estimation process in LISA involves several complex steps. First, the raw data from the TDI observables specifically the second-generation observables X, Y, and Z, which accommodate the fact that the spacecraft formation allows for non rigid arm lengths are transformed into the A, E, and T channels [113, 114, 115, 116, 117, 118, 119]. These channels are designed to be approximately uncorrelated in terms of their noise properties, improving the clarity and reliability of the data for further analysis.

To study the GW events detected by the LISA, it is essential to estimate the parameters of these sources. The process of parameter estimation for LISA sources is quite different from that used for ground-based observatories like LIGO, Virgo, and KAGRA [120]. For our analysis of LISA sources, we employ the BBHx Python package[121] , which is specifically designed for this purpose. Our methodology involves generating waveforms using BBHx [121], which incorporates the IMRPhenomD [20, 122] waveform approximation model. We obtain the waveform for a one-year observation period of the LISA sources. Subsequently, we use BBHx [121] to calculate the likelihood of the generated waveforms. For the estimation of the posterior distributions of key parameters for this study specifically, the masses of the black holes ($\mathrm{m_{1}}$ and $\mathrm{m_{2}}$), the luminosity distance ($\mathrm{D_{L}}$), and the inclination angle ($\mathrm{i}$) we constrain the priors of all other parameters to delta functions. Finally, for the parameter estimation, we employ the Dynesty [123] sampler. Among these, the posterior of $\mathrm{D_{L}}$ is particularly crucial for our study, as it will be used for the inference of $\mathrm{\mathcal{F}(z)}$.

The accuracy of luminosity distance measurements in GW astronomy is influenced not only by the sensitivity and configuration of the detector but also significantly by weak lensing, especially for sources detected at high redshifts by LISA [124]. Gravitational lensing affects GWs similarly to EM radiation. For GW events expected from large redshifts ($\mathrm{z>1}$), weak lensing is a common occurrence, along with the occasional strongly-lensed sources. A lens with magnification $\mu$ modifies the observed luminosity distance to $\mathrm{\frac{D_{L}}{\sqrt{\mu}}}$, introducing a systematic error $\mathrm{\Delta D_{L}/D_{L}=1-\frac{1}{\sqrt{\mu}}}$ [125]. This lensing error, when convolved with the expected magnification distribution $\mathrm{p(\mu)}$ from a standard $\mathrm{\Lambda}$CDM model, significantly influences the overall measurement accuracy. In this study, we utilize the following model to estimate the error due to weak lensing [126]

$$\mathrm{\frac{\sigma_{WL}}{D_{L}}=\frac{0.096}{2}\left(\frac{1-(1+z)^{-0.62}}{0.62}\right)^{2.36}}.$$

(3.11)

The total uncertainty in the measured luminosity distance, therefore, is expressed as

$$\mathrm{\sigma_{D_{L}}^{2}=\sigma_{GW}^{2}+\sigma_{WL}^{2}},$$

(3.12)

where $\mathrm{\sigma_{GW}}$ is the uncertainty derived from the detector’s setup during parameter estimation, and $\mathrm{\sigma_{WL}}$ accounts for the lensing effects. This comprehensive approach ensures a more accurate understanding of the intrinsic properties of GW sources.

The figure 6 shows the estimated values and uncertainties for component masses, luminosity distance, and inclination angle of LISA source at a redshift of z=2.4. This highlights the potential of LISA to offer precise insights into the cosmological properties of binary systems from the early universe. The figure 7 presents the luminosity distance errors as a function of redshift for LISA sources, showcasing the variability and precision of measurements across different redshifts. Each point in the plot represents the measured luminosity distance at various redshifts, with total uncertainty, derived from two specific sources of errors: GW source parameter estimation error and WL error. This figure underscores the advanced capabilities of LISA in probing the universe’s most distant reaches, illustrating the impact of different sources of uncertainty on measurement accuracy, and paving the way for profound cosmological discoveries.

The imprint manifests as a distinct peak within the two-point angular correlation function (2PACF) of galaxies, $\mathrm{w(\theta)}$, which offers a method to measure the universe’s expansion rate across its history. To detect the BAO feature within $\mathrm{w(\theta)}$, it is modeled using a combination of a power-law and a Gaussian component, represented by the following equation[153, 154, 155]

$$\mathrm{w(\theta,z)=a_{1}+a_{2}\theta^{k}+a_{3}\exp\left(-\frac{(\theta-\theta_{\text{FIT}}(z))^{2}}{\sigma_{\theta}^{2}}\right)},$$

(5.1)

where $a_{i}$ are parameters representing the fit, $\mathrm{\sigma_{\theta}}$ is the width of the BAO peak, and $\mathrm{\theta_{FIT}}$ indicates the angular position of the BAO feature. This method offers a robust standard ruler that enables independent estimates of the angular diameter distance, $\mathrm{D_{a}(z)}$, further elucidating the dynamics of the universe’s expansion[156, 157, 158, 159, 160].

The theoretical representation of the 2PACF, $\mathrm{w(\theta)}$, is formulated as

$$\mathrm{w(\theta)=\sum_{l\geq 0}\left(\frac{2l+1}{4\pi}\right)P_{l}(\cos(\theta))C_{l}},$$

(5.2)

where $\mathrm{P_{l}}$ denotes the Legendre polynomial of the $\mathrm{l^{th}}$ order, and $\mathrm{C_{l}}$ represents the angular power spectrum. The latter is derived from the three-dimensional matter power spectrum $\mathrm{\mathcal{P}(k)}$ through the integral[60, 147]

$$\mathrm{C_{l}=\frac{1}{2\pi^{2}}\int 4\pi k^{2}\mathcal{P}(k)\psi_{l}^{2}(k)e^{-k^{2}/k_{eff}^{2}}},$$

(5.3)

where $\mathrm{\psi_{l}(k)}$ is defined as:

$$\mathrm{\psi_{l}(k)=\int dz\phi(z)j_{l}(kr(z))},$$

(5.4)

incorporating $\mathrm{\phi(z)}$, the galaxy selection function, and $\mathrm{j_{l}}$, the spherical Bessel function of the $\mathrm{l^{th}}$ order, with $\mathrm{r(z)}$ being the comoving distance. An exponential damping factor $\mathrm{e^{-k^{2}/k_{eff}^{2}}}$ is included to enhance the integration’s convergence in calculating the angular power spectrum, where $\mathrm{1/k_{eff}=3\text{ Mpc/h}}$ is consistently used across all calculations, proven to have negligible effects on the scales under consideration[161]. For deriving the matter power spectrum, we employ the CAMB[162] module, configuring the galaxy selection function as a normalized Gaussian distribution. We set the standard deviation of this distribution to $\mathrm{0.03(1+z)\%}$ of its mean value, aligning it with the anticipated photo-z error from the Vera Rubin Observatory [139, 163, 164] .

Figure 8 showcases the variation of $\mathrm{\theta^{2}w(\theta)}$ with angular separation at a fixed redshift $\mathrm{z=1.0}$. A marked peak at approximately $\mathrm{\theta\sim 2.5}$ reveals the BAO scale, central to our analysis. Furthermore, the covariance of $\mathrm{w(\theta)}$, expressed as $\mathrm{Cov_{\theta\theta^{\prime}}}$, is computed as:

$$\mathrm{Cov_{\theta\theta^{\prime}}=\frac{2}{f_{sky}}\sum_{l\geq 0}\frac{2l+1}{(4\pi)^{2}}P_{l}(\cos(\theta))P_{l}(\cos(\theta^{\prime}))\left(C_{l}+\frac{1}{n}\right)^{2}},$$

(5.5)

where $\mathrm{1/n}$ signifies the shot noise linked to the galaxy’s number density per steradian, and $\mathrm{f_{sky}}$ is the fraction of the sky surveyed[165, 166]. A depiction (Figure 9) outlines how the covariance matrix evolves with angular separations at a given redshift(z=1), indicating a growing correlation as $\mathrm{\theta}$ and $\mathrm{\theta^{\prime}}$ increase. We have included the complete covariance matrix in the analysis, as there are sufficient correlations between different scales.

Figure 10 elaborates on the BAO scale’s variation and associated error across different redshifts, emphasizing the influence of photo-z errors. For LSST’s planned photo-z error, we illustrate the error on BAO measurements from $\mathrm{z=1}$ to $\mathrm{z=3}$ in cyan, showing how the expected errors will affect the precision of BAO measurements. Vera Rubin Observatory[139, 163], along with projects like Euclid [167], will significantly enhance our understanding of BAO by surveying a total area of 18,000 square degrees, extending precise measurements to higher redshifts. The figure also depicts the BAO measurement error anticipated from Dark Energy Spectroscopic Instrument(DESI) [168] up to redshift 3.7, shown in black, demonstrating the expected performance of DESI in this redshift range. The DESI is a notable five-year terrestrial survey focusing on exploring BAO and the evolution of cosmic structures through redshift-space distortions. Covering 14,000 square degrees of the sky, corresponding to a sky fraction ($\mathrm{f_{\text{sky}}}$) of 0.3, DESI aims for an accuracy better than 0.5% within the redshift intervals $\mathrm{0.0<z<3.7}$. Despite the lack of currently planned instruments to measure BAO beyond redshift 3.7, we extend the BAO measurements up to redshift 6 using the same photo-z error anticipated for LSST, shown in red. This extended range provides a hypothetical look at the potential capabilities of future surveys and their ability to probe deeper into the universe. This visualization provides essential insights into the variability of the BAO scale and its error with increasing redshift, highlighting a rise in relative error as the standard deviation of the selection function increases with redshift. With the capability of LISA in probing the high redshift Universe, this plot shows that future galaxy surveys with capabilities of measuring BAO beyond $\mathrm{z=3.5}$ can bring interesting tests of fundamental physics explored in this analysis.

The correlation between the EM luminosity distance at a specific redshift (z) and critical cosmological parameters, such as the BAO scale ($\mathrm{\theta_{BAO}}$) and the sound horizon ($\mathrm{r_{s}}$), is succinctly expressed by the Equation 2.7. This equation is central to the analysis of the frictional term, linking it to three significant cosmological lengths: (i)GW Luminosity Distance: The distance for detectable GW events is estimated through the parameter estimation of GW events as described in Subsection 3.3. (ii)BAO Scale ($\mathrm{\theta_{BAO}}$): Derived from the 2PACF and expressed in Equation 5.1, this scale is utilized in a versatile manner, allowing for application across various surveys. (iii)Sound Horizon Distance: Estimated by analyzing CMB acoustic oscillations through missions like WMAP [63], Planck [64], ACTPol [169], SPT-3G [170], and CMB-S4 [171]. The first peak in the angular power spectrum at recombination reflects the sound horizon scale, measured precisely by Planck as $\mathrm{147.09\pm 0.26}$ Mpc at the drag epoch and redshift: Inferred from the EM counterparts.

To estimate the posterior distribution of $\mathrm{\mathcal{F}(z)}$ as a function of redshift for $\mathrm{n_{GW}}$ GW sources, we employ a Hierarchical Bayesian framework. The posterior distribution is represented as follows [57]

$$\mathrm{P(\mathcal{F}(z))\propto\Pi(\mathcal{F}(z))\prod_{i=1}^{n_{GW}}\iiint dr_{s}P(r_{s})dD_{L}^{GW^{i}}d\theta_{BAO}^{i}P(\theta_{BAO}^{i})\mathcal{L}(D_{L}^{GW^{i}}|\mathcal{F}(z^{i}),\theta_{BAO}^{i},r_{s},z^{i})}.$$

(6.1)

In this framework, $\mathrm{P(\mathcal{F}(z))}$ denotes the posterior probability density of $\mathrm{\mathcal{F}(z)}$. The term $\mathrm{\mathcal{L}(D_{L}^{\text{GW}^{i}}|\mathcal{F}(z^{i}),\theta_{\text{BAO}}^{i},r_{s},z^{i})}$ corresponds to the likelihood function, which represents the posterior distribution of the GW luminosity distance. In this expression, $\mathrm{P(r_{s})}$ and $\mathrm{P(\theta_{BAO}^{i})}$ representing the prior probabilities for the sound horizon and BAO scale, respectively. $\mathrm{\Pi(\mathcal{F}(z))}$ indicates the prior distribution for $\mathrm{\mathcal{F}(z)}$, reflecting pre-existing knowledge or assumptions about the frictional term before the data is considered. A flat prior from 0.1 to 2 is adopted for $\mathrm{\mathcal{F}(z)}$. This investigation focuses solely on bright sirens, assuming that the redshift of the host galaxies of the GW sources is measured spectroscopically. This assumption implies precise measurements of redshift for these sources, as the error in redshift measurements is very small compared to the errors associated with the BAO scale and distance measurements.

In previous discussions, we concentrated on analyzing the frictional term, denoted as $\mathrm{\mathcal{F}(z)}$ in Equation 2.7, under the assumption of a constant Hubble constant ($\mathrm{H_{0}}$). Moving forward, we expand our analysis to concurrently estimate both $\mathrm{\mathcal{F}(z)}$ and $\mathrm{H_{0}}$. Our primary objective remains the precise determination of $\mathrm{\mathcal{F}(z)}$, while concurrently assessing $\mathrm{H_{0}}$. This dual estimation strategy is especially pertinent with the advent of future detectors, which will enhance our capability to probe deeper into redshifts, thus enabling either combined or independent measurements of these parameters. This methodology is particularly promising for addressing the so-called Hubble tension (for a comprehensive review [172]), the notable discrepancy between the locally measured values of $\mathrm{H_{0}}$ and those inferred from CMB observations [173].

By employing a hierarchical Bayesian framework, we efficiently estimate these parameters using a comprehensive array of observational data, encompassing GW detections, CMB data, and large-scale structure surveys, as elaborated in earlier sections. The hierarchical integration of these datasets is meticulously designed to manage their distinct uncertainties and biases effectively. This integration produces correlated constraints that significantly enhance our understanding of cosmological dynamics, help resolve discrepancies among observational datasets, and uncover critical relationships between key cosmological parameters. The joint posterior distribution for $\mathrm{\mathcal{F}(z)}$ and $\mathrm{H_{0}}$, based on these considerations, is formulated as follows

$\displaystyle\mathrm{P(\mathcal{F}(z),H_{0})\propto\Pi(\mathcal{F}(z))\Pi(H_{0})\prod_{i=1}^{n_{GW}}\iiint}\mathrm{dr_{s}P(r_{s}^{i})dD_{L}^{GW^{i}}d\theta_{BAO}^{i}P(\theta_{BAO}^{i})\mathcal{L}(D_{L}^{GW^{i}}|\mathcal{F}(z^{i}),H_{0},\theta_{BAO}^{i},r_{s},z^{i})}.$

(6.2)

In this model, $\mathrm{P(\mathcal{F}(z),H_{0}})$ denotes the joint posterior probability density function for $\mathrm{\mathcal{F}(z)}$ and $\mathrm{H_{0}}$. The likelihood function is represented by $\mathrm{\mathcal{L}(D_{L}^{\text{GW}^{i}}|\mathcal{F}(z)^{i},H_{0}^{i},\theta_{\text{BAO}^{i}},r_{s},z^{i})}$. The terms $\mathrm{\Pi(\mathcal{F}(z))}$ and $\mathrm{\Pi(H_{0})}$ define the prior distributions for $\mathrm{\mathcal{F}(z)}$ and $\mathrm{H_{0}}$, respectively, encapsulating pre-existing beliefs or assumptions prior to analyzing the data. Within our hierarchical Bayesian framework, we utilize flat priors for both parameters, indicating minimal initial bias. The range for $\mathrm{\mathcal{F}(z)}$ is set from 0.1 to 2.0, and for $\mathrm{H_{0}}$, from 40 km/s/Mpc to 100 km/s/Mpc.

In cases where $\mathrm{H_{0}}$ is held constant, it acts as an overall normalization factor for $\mathrm{\mathcal{F}(z)}$ across different redshifts. This setup provides a unique opportunity to normalize $\mathrm{\mathcal{F}(z)}$ measurements to 1 if an independent and precise measurement of $\mathcal{F}(z)$ in our local universe is possible. Such normalization, particularly effective at low redshifts, helps eliminate the dependence on $\mathrm{H_{0}}$ and demonstrates the potential deviations from GR due to the cumulative effects of modified gravity wave propagation over cosmological distances. The uncertainty in the measurement of the frictional term, $\mathrm{\mathcal{F}(z)}$, is inversely proportional to the square root of the number of GW sources, denoted as $\mathrm{1/\sqrt{N_{GW}}}$. This relationship suggests that increasing the count of GW sources detected by LISA enhances the precision of $\mathrm{\mathcal{F}(z)}$ measurements. However, the accuracy improvement reaches a plateau beyond a certain number of sources, primarily due to errors associated with the BAO scale. The precision in estimating the BAO scale emerges as a critical factor limiting the refinement of $\mathrm{\mathcal{F}(z)}$ estimation accuracy. Thus, advancing the accuracy of BAO scale measurements is essential for further enhancements in determining the frictional term with redshift. Essentially, optimizing the number of GW observations, alongside improvements in BAO scale accuracy and GWs luminosity distance measurements, is crucial for advancing our understanding of $\mathrm{\mathcal{F}(z)}$ and probing new aspects of fundamental physics.

On the GWs front, the accuracy of $\mathrm{\mathcal{F}(z)}$ relies predominantly on the precision of the luminosity distance to GWs sources, $\mathrm{D_{L}^{GW}(z)}$. Enhancements in measuring $\mathrm{D_{L}^{GW}(z)}$ would significantly boost the fidelity of $\mathrm{\mathcal{F}(z)}$ estimations. Given the degeneracy of $\mathrm{D_{L}^{GW}(z)}$ with other GWs parameters, notably the inclination angle ($\mathrm{i}$), acquiring more accurate measurements of the inclination angle, potentially through EM counterparts, could substantially reduce uncertainties in $\mathrm{D_{L}^{GW}(z)}$ and, by extension, refine $\mathrm{\mathcal{F}(z)}$ estimates [174]. The effect of peculiar velocities on redshift inference from EM counterparts is notably critical at lower redshifts ($\mathrm{z<0.03}$) but becomes negligible at higher redshifts [175]. Given that the majority of GW sources detected by LISA are expected to be at higher redshifts, this effect is disregarded in our analysis.

Building upon the hierarchical Bayesian framework introduced earlier to estimate $\mathrm{\mathcal{F}(z)}$ with a fixed $\mathrm{H_{0}}$, we now extend this approach to concurrently measure $\mathrm{\mathcal{F}(z)}$ and the Hubble constant, $\mathrm{H_{0}}$. In this extended framework, we derive the EM luminosity distance by integrating various cosmological scales, including the BAO scale, sound horizon distance, and redshift. Figures 11 provide a detailed analysis of $\mathrm{\mathcal{F}(z)}$ measurements under scenarios of both fixed and varying $\mathrm{H_{0}}$ with only six events up to $z=6$. The uncertainty on the parameter $\mathrm{\mathcal{F}(z)}$ increased by approximately a factor of 2 to 2.5 when $\mathrm{H_{0}}$ is allowed to vary relative to scenarios with a fixed $\mathrm{H_{0}}$. This scenario with only six events corresponds to a pessimistic case with only $4\%$ EM bright events for the Light Seed model and the Heavy Seed model and about $30\%$ for the Heavy Seed+Delay model in four years of observation period. This corresponds to a single events at every redshift up to $z=6$. We have shown the constraints for a few scenarios of detectable bright sirens from LISA in appendix A. If the efficiency of LISA bright sirens are less than this, then one needs to consider a longer observation period of the LISA mission for at least one detection of bright standard siren.

In Figure 12, we specifically illustrate the improved precision in tracking the redshift-dependent variations of the Planck mass ($\mathrm{\mathcal{F}(z)}$) for MBBHs systems, emphasizing enhanced measurements in the redshift range from $z=1.2$ to $z=6.0$. Additionally, Figure 13 illustrates the optimal joint measurements of $\mathrm{\mathcal{F}(z)}$ and the Hubble constant ($\mathrm{H_{0}}$). In Figure 14, we show the measurement of $H_{0}$ marginalized over $\mathrm{\mathcal{F}(z)}$ for the three different population models, for all events up to a redshift of $\mathrm{z=6.0}$. The results on $\mathrm{\mathcal{F}(z)}$ and $H_{0}$ for scenario with MBBHs distribution having a mass ratio beyond $q=1$ is shown in the appendix B.