# Bayesian framework to infer the Hubble constant from the cross-correlation of individual gravitational wave events with galaxies

Tathagata Ghosh, Surhud More, Sayantani Bera, Sukanta Bose

► To cite this version:

Tathagata Ghosh, Surhud More, Sayantani Bera, Sukanta Bose. Bayesian framework to infer the Hubble constant from the cross-correlation of individual gravitational wave events with galaxies. 2026. [hal-05489029]

HAL Id: hal-05489029

https://hal.science/hal-05489029v1

Preprint submitted on 2 Feb 2026

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Bayesian framework to infer the Hubble constant from the cross-correlation of
individual gravitational wave events with galaxies

Tathagata Ghosh ①,<fnref n="1" /> Surhud More ①,<fnref n="1" />,² Sayantani Bera ①,³,⁴ and Sukanta Bose ①,⁵,<fnref n="1" />

¹Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India

²Kavli Institute for the Physics and Mathematics of the Universe
(WPI), University of Tokyo, 5-1-5, Kashiwanoha, 2778583, Japan

³Aix-Marseille Université, Université de Toulon, CNRS, CPT, Marseille, France

⁴Departament de Física, Universitat de les Illes Balears, IAC3 - IEEC, Carretera de Valldemossa km 7.5, E-07122 Palma, Spain

⁵Department of Physics & Astronomy, Washington State University, 1245 Webster, Pullman, WA 99164-2814, USA

(Dated: March 10, 2025)

Gravitational waves (GWs) from the inspiral of binary compact objects offer a one-step mea-
surement of the luminosity distance to the event, which is essential for the measurement of the
Hubble constant, $H_0$, which characterizes the expansion rate of the Universe. However, unlike
binary neutron stars, the inspiral of binary black holes is not expected to be accompanied by elec-
tromagnetic radiation and a subsequent determination of its redshift. Consequently, independent
redshift measurements of such GW events are necessary to measure $H_0$. In this study, we present
a novel Bayesian approach to infer $H_0$ by measuring the overdensity of galaxies around individual
binary black hole merger events in configuration space. We model the measured overdensity using
the 3D cross-correlation between galaxies and GW events, explicitly accounting for the GW event
localization uncertainty. We demonstrate the efficacy of our method with 250 simulated GW events
distributed within 1 Gpc in colored Gaussian noise of Advanced LIGO and Advanced Virgo detectors
operating at O4 sensitivity. We show that such measurements can constrain the Hubble constant
with a precision of ≲ 8% (90% highest density interval). We highlight the potential improvements
that need to be accounted for in further studies before the method can be applied to real data.

# I. INTRODUCTION

The discovery of gravitational wave (GW) events opened an independent way to probe the expansion history of the Universe. The very first multimessenger approach to infer the Hubble constant was possible due to the detection of GW from the binary neutron star (BNS) merger, GW170817 [1] along with the electromagnetic (EM) counterpart [2]. The independent measurement of luminosity distance of the BNS merger from the GW data and its redshift from the EM data allows an inference of the Hubble constant $H_0 = 70_{-8.0}^{+12.0} \text{ km s}^{-1} \text{ Mpc}^{-1}$ [3]. However, most of the GW observations by the detectors in the LIGO-Virgo-KAGRA collaboration [4-6] are primarily binary black hole (BBH) mergers [7], which are not accompanied by EM counterparts. Such events are aptly called *dark sirens*. Though the unique identification of the host galaxy for such BBHs is not currently possible due to large sky localization errors, there are alternative approaches to infer the Hubble constant from such dark sirens. The idea initially proposed by Schutz [8] involves the use of galaxy catalogs for the identification of galaxies clustered within the sky localization of the GW events, which would provide noisy but useful information on the redshift of the events. His approach has been implemented as a statistical host identification method in a Bayesian framework (also sometimes termed the galaxy catalog method) by Refs. [9-21]. In another approach called the spectral siren method, the source redshift can be inferred by invoking the universality of the source frame mass distribution motivated from astrophysics [22-28]. Both of these methods, i.e.,

the galaxy catalog method and spectral siren method,
have been applied to 47 BBHs from GWTC-3 [7] in ad-
dition to the bright siren GW170817 and its EM counter-
part [29], to constrain the Hubble constant to $H_0 = 68_{-6}^{+8}$
$\text{km s}^{-1} \text{ Mpc}^{-1}$ and $H_0 = 68_{-8}^{+12} \text{ km s}^{-1} \text{ Mpc}^{-1}$ [30], re-
spectively<fnref n="1" />. There are also ongoing efforts [31-34] to
unify the spectral and galaxy methods. Furthermore, if
at least one of the components of the merger is a neutron
star, one can use the measurement of its tidal deformabil-
ity and knowledge of the neutron star equation of state to
infer the source redshift [35-41], and hence $H_0$ inference
by combining the redshift information with the measured
luminosity distance of the merger event.

The current state-of-the-art methods [16, 26] adopted
to estimate $H_0$ in the recent LIGO-Virgo-KAGRA (LVK)
cosmology analysis [30] are dominated by the assump-
tions related to the nature of the GW source popula-
tion [30]. However, it is important to note that these
two methods mentioned above do not utilize the informa-
tion on the spatial clustering of galaxies explicitly. Both
the GW source population and the galaxy population
are tracers of the same underlying large-scale structure
in the Universe. So, GW sources share the same large-
scale structure as galaxies at their redshifts and, hence,
are correlated to each other. Therefore, one can expect
that the cross-correlation function between galaxies with
known redshift and GW sources would be nonzero at the
true redshift of the GW sources and, thus, in conjunc-
tion with the luminosity distances of the GW events, can

<footnote n="1">¹ The errors denote 68% credible intervals in this case.</footnote>
be used to constrain the cosmological parameters such as the Hubble constant [42–49]. In this work, we demonstrate a novel Bayesian framework to infer $H_0$ from individual GW events by utilizing the cross-correlation between these GW events and the galaxy catalog. This approach models the galaxy overdensity for each GW event while also accounting for its volume uncertainty region when estimating the Hubble constant through this cross-correlation. This method, which can be termed as cross-correlation, does not directly rely on assumptions related to the intrinsic population parameters that describe the GW sources.

The paper is structured as follows. In Sec. II, we describe the methodology of using galaxy clustering information to estimate $H_0$ from individual GW events. After that, we elaborate on the simulations performed to demonstrate the methodology in Sec. III. In Sec. IV, we show the efficacy of our method in recovering $H_0$ and examine how various parameters impact the precision of the resulting constraints. Finally, we conclude our paper and discuss future directions in Sec. V.

## II. METHODOLOGY

In the standard model of cosmology (also known as the ΛCDM model), the matter distribution in the Universe forms a cosmic web of large-scale structure with clustering properties determined by the value of cosmological parameters. Dark matter halos form at the peaks of the matter density distribution. The clustering properties of these halos on large scales reflect the clustering properties of the matter distribution, modified with a multiplicative factor $b_h^2$, where $b_h$ is the linear bias that corresponds to the ratio of the halo overdensity ($\delta_h$) to the matter overdensity ($\delta_m$),

$$ \delta_h(\mathbf{x}, M) = \frac{n_h(\mathbf{x}, M)}{\bar{n}_h} - 1 = b_h(M)\delta_m(\mathbf{x}). \quad (1) $$

Here, $n_h(\mathbf{x}, M)$ denotes the number of halos at a given position $\mathbf{x}$ with halo mass $M$ and $\bar{n}_h$ is the average number of halos [50]. The bias of the dark matter halos increases with increasing halo mass [51, 52]. More luminous galaxies are expected to form within more massive dark matter halos [53–57], and on large scales, they inherit the bias of the dark matter halos they inhabit, such that

$$ b_g = \int dM P_h(M) b_h(M) \quad (2) $$

where $P_h(M)$ specifies the distribution of halo masses to which the galaxies belong [52, 58, 59] (for a detailed review on the halo and galaxy bias, one may refer to [60]).

In this section, we describe the methodology to utilize the cross-correlation between individual GW events and galaxy catalogs to constrain the Hubble constant. The spatial clustering between galaxies and GW events, separated by the comoving distance $\mathbf{r}$, can be characterized

by the 3D cross-correlation function, $\xi_{gw,g}(\mathbf{r})$, given by

$$ \xi_{gw,g}(\mathbf{r}) \equiv \langle \delta_{gw}(\mathbf{x})\delta_g(\mathbf{x}+\mathbf{r}) \rangle, \quad (3) $$

where $\mathbf{x}$ represents the location of GW sources; $\delta_{gw}$ and $\delta_g$ denote the overdensity of the GW source distribution and the galaxy distribution, respectively. As the Universe is assumed to be isotropic, the cross-correlation function depends on the scalar distance $r = |\mathbf{r}|$. Since GW events and galaxies are tracers of the same underlying matter distribution, both GW overdensity and galaxy overdensity are related to the matter overdensity $\delta_m$ on large scales as

$$ \delta_g = b_g \delta_m \quad (4a) $$

$$ \delta_{gw} = b_{gw} \delta_m \quad (4b) $$

Here, the proportionality constants $b_g$ and $b_{gw}$ are the linear biases of the population of galaxies and GW sources with respect to the matter distribution, respectively. A bias of unity would imply that the observed population follows the underlying matter distribution. In general, the bias parameters depend on different properties of the sources. The dependence of galaxy bias on redshift can be estimated from the measurement of autocorrelation functions of galaxies in a given redshift slice with themselves. The redshift-dependent GW bias needs to be assumed a parametric form and be marginalized over those parameters (see Ref. [46]).

The cross-correlation function $\xi_{gw,g}$, defined in Eq. (3) can be written as,

$$ \xi_{gw,g} = b_{gw} b_g \xi_m. \quad (5) $$

Here, $\xi_m$ is the nonlinear matter power spectrum [61], which we compute under the assumption of the standard ΛCDM model of cosmology. In the flat ΛCDM Universe, the luminosity distance $d_L$ is related to the redshift $z$ as follows:

$$ \begin{aligned} d_L &= \frac{c(1+z)}{H_0} \int_0^z \frac{dz'}{H(z')} \\ &= d_c(1+z) \end{aligned} \quad (6) $$

Here, $d_c$ denotes the comoving distance, and $H(z)$ is the Hubble parameter, expressed in terms of the Hubble constant $H_0$ and matter density parameter $\Omega_m$ as:

$$ H(z) = H_0 E(z) = H_0 \sqrt{\Omega_m(1+z)^3 + (1-\Omega_m)}. \quad (7) $$

Since the correlation function $\xi_{gw,g}(r)$ quantifies the excess probability over a random Poisson distribution of finding a galaxy in a volume element $dV$, which is separated by a comoving distance $\mathbf{r}$ from a GW event located at position $\mathbf{x}$, we have

$$ n_g(\mathbf{x}+\mathbf{r})dV = \bar{n}_V [1+\xi_{gw,g}(r)] dV, \quad (8) $$

where $\bar{n}_V$ is the average number density of galaxies in comoving volume. So, the galaxy overdensity can be calculated by integrating Eq. (8) over some finite volume if
<figure><img src="image_2.png" /><figcaption>FIG. 1. An illustration of the geometry in the three-dimensional Cartesian coordinate system (X Y Z) in order to compute the modeled galaxy overdensity δmod g following Eq. (9). The horizontal arrow denotes the line of sight from the observer towards the maximum posterior probability position of the GW event (αgw m, δgw m), the conical volume has an angular radius of θmax around this position, and the two vertical planes indicate the redshift range of width Δz, centered at z within which we compute this modeled galaxy overdensity.</figcaption></figure>

we have the theoretical prediction of the nonlinear matter power spectrum ξm.

Consider a GW event located at its true position r gw denoted by redshift z gw = z gw (d gw L , H$_{0}$ ), right ascension α gw , and declination δ gw . The modeled galaxy overden-

sity, δmod g , in an angle θmax around the maximum a pos- teriori probability (MAP) location (αgw m, δgw m) and within a redshift bin centered around z $^{′}$, can be computed by carrying out the following volume integral (see Fig. 1 for the geometry) assuming the MAP is along the line of sight:

$$
\delta_g^{\text{mod}}(z | \mathbf{r}^{\text{gw}}) = \frac{\int_{z-\Delta z/2}^{z+\Delta z/2} \int_0^{\theta_{\text{max}}} \int_0^{2\pi} \bar{n}_V(z') \xi_{\text{gw,g}}(|\mathbf{r}(\mathbf{r}^{\text{gw}}, \mathbf{r}')|) \sin \theta' d_c'^2 J\left(\frac{d_c'}{z'}\right) d\alpha' d\theta' dz'}{\int_{z-\Delta z/2}^{z+\Delta z/2} \int_0^{\theta_{\text{max}}} \int_0^{2\pi} \bar{n}_V(z') \sin \theta' d_c'^2 J\left(\frac{d_c'}{z'}\right) d\alpha' d\theta' dz'} . \quad (9)
$$

Here, r denotes the comoving distance (in h − 1 unit) from r gw to a position at r ′ = ( z $^{′}$, α$^{′}$, δ$^{′}$). Without loss of generality, we define our coordinate system such that δ gw m = 0. For any other orientation, we can apply an appropriate rotation of the coordinate axis. Based on this geometry, in Eq. (9), we use θ ′ = π/ 2 − δ $^{′}$.

The true GW event redshift can be inferred from the true luminosity distance by inverting Eq. (6). The quan- tity d ′ c in the above equation denotes the comoving dis- tance corresponding to a redshift z ′ within the redshift bin enclosed by the vertical planes and J ( d ′ c z $^{′}$) corre- sponds to the Jacobian of the coordinate transformation from comoving distance to redshift, which can be written

as [see Eq. (6)]

$$
J\left(\frac{d_c}{z}\right) = \frac{\partial d_c}{\partial z} = \frac{c}{H(z)}. \qquad (10)
$$

The modeled galaxy overdensity thus depends on the choice of H$_{0}$ .

The observed galaxy overdensity within the redshift bin, centered around z , can be measured within the an- gular distance θmax around the line of sight towards the MAP position (αgw m, δgw m), as

$$
\delta_g^{\text{obs}}(z) = \frac{n(z | \delta_m^{\text{gw}}, \alpha_m^{\text{gw}})}{\bar{n}(z)} - 1, \quad (11)
$$
where $n(z | \delta_m^{gw}, \alpha_m^{gw})$ denotes the number of galaxies within the redshift bin. For a galaxy survey with uniform depth, the average number of galaxies $\bar{n}(z)$ expected within the redshift bin can be calculated by averaging the observed number of galaxies within the same angular radius $\theta_{max}$ around random lines of sight.

The observed overdensity $\delta_g^{obs}$ can be calculated in different redshift bins of galaxies. We call this the galaxy data vector $\mathbf{d}_g^{obs}$. This galaxy data vector can be compared with the model vector of $\delta_g^{mod}$ for the same redshift bins of galaxies, denoted by $\mathbf{d}_g^{mod}$ to construct a galaxy

overdensity likelihood. Note that the model vector is conditioned on the true location (both on the sky and in redshift) of a particular GW event. The comparison of the two should provide a constraint on the Hubble constant, $H_0$.

The true position on which the modeled galaxy overdensity is conditioned has its own uncertainty as the sky localization of the GW event is not precisely known but known as a posterior distribution over $\mathbf{d}_{gw} = \{z^{gw}(d_L^{gw}, H_0), \alpha^{gw}, \delta^{gw}\}$ inferred from the gravitational wave strain data $\mathbf{d}_{strain}$. The posterior for the Hubble constant, in that case, is given by,

$$ 
\begin{aligned}
p(H_0 | \mathbf{d}_{\text{strain}}, \mathbf{d}_g^{\text{obs}}) &= \int p(H_0, \mathbf{d}_{gw} | \mathbf{d}_{\text{strain}}, \mathbf{d}_g^{\text{obs}}) d\mathbf{d}_{gw} \\
&\propto \int \mathcal{L}(\mathbf{d}_{\text{strain}}, \mathbf{d}_g^{\text{obs}} | H_0, \mathbf{d}_{gw}) P(H_0, \mathbf{d}_{gw}) d\mathbf{d}_{gw} \\
&\propto \int \mathcal{L}(\mathbf{d}_g^{\text{obs}} | H_0, \mathbf{d}_{gw}) \mathcal{L}(\mathbf{d}_{\text{strain}} | H_0, \mathbf{d}_{gw}) P(H_0, \mathbf{d}_{gw}) d\mathbf{d}_{gw},
\end{aligned} 
\quad (12) $$

where we have applied Bayes' theorem and assumed independence of the galaxy overdensity and the strain data from gravitational wave events. The integration in the last equation is performed over the localization uncertainty of the GW event. The first likelihood in Eq. (12) corresponds to the galaxy overdensity likelihood, and we assume it to be a Gaussian such that

$$ \mathcal{L}(\mathbf{d}_g^{\text{obs}} | H_0, \mathbf{d}_{gw}) \propto \exp \left[ -\frac{1}{2} (\mathbf{d}_g^{\text{obs}} - \mathbf{d}_g^{\text{mod}})^{\text{T}} \mathbf{C}^{-1} (\mathbf{d}_g^{\text{obs}} - \mathbf{d}_g^{\text{mod}}) \right] \quad (13) $$

Here, $C$ denotes the covariance of the overdensity field, and we assume it to be diagonal with individual elements $\sigma_{\delta}^2$ corresponding to the variance in the observed galaxy overdensity in a given redshift bin within the angle $\theta_{max}$. The procedure for calculating $\sigma_{\delta}$ using random lines of sight is described in Sec. III.

The second likelihood term in the last integral of Eq. (12) corresponds to the strain data. Often, the posterior distributions of $P(\mathbf{d}_{gw})$ are computed assuming certain priors. Therefore, the likelihood can be obtained by taking the ratio of the posteriors for source parameters in $\mathbf{d}_{gw}$ and the priors that are used to infer these parameters. This allows us to separately impose our own priors on the quantities $P(H_0, \mathbf{d}_{gw})$ in Eq. (12).

Since GW events are independent, the final joint posterior of the Hubble constant can be obtained by combining the independent likelihoods of all the detected sources,

$$ 
\begin{aligned}
P(H_0 | \{\mathbf{d}_{\text{strain}}\}, \{\mathbf{d}_g^{\text{obs}}\}) &\propto P(H_0) \prod_i \mathcal{L}(\mathbf{d}_{\text{strain}_i}, \mathbf{d}_{g_i}^{\text{obs}} | H_0) \\
&\propto P(H_0) \prod_i \int \mathcal{L}(\mathbf{d}_{\text{strain}_i}, \mathbf{d}_{g_i}^{\text{obs}} | H_0, \mathbf{d}_{gw}) P(\mathbf{d}_{gw}) d\mathbf{d}_{gw}
\end{aligned} 
\quad (14) $$

where $P(H_0)$ represents the prior over $H_0$, the index $i$ denotes individual GW events, and $\mathbf{d}_{g_i}^{\text{obs}}$ is the galaxy overdensity likelihood. Note that the model vector is conditioned on the true location (both on the sky and in redshift) of a particular GW event. The comparison of the two should provide a constraint on the Hubble constant, $H_0$.
density data vector corresponding to the *i*th GW event.

It is important to highlight that although the Bayesian formalism proposed in this work may look structurally similar to population inference techniques using GWs [62, 63] and other cosmological inference methods such as spectral siren method and galaxy catalog method, it does not share the same selection biases as mentioned in Ref. [64]. Unlike those methods, the present work does not explicitly rely on mass-model assumptions. Rather, it measures and models observed galaxy overdensity associated with the detected GW events, determined by an SNR threshold. The observable here is the galaxy overdensity around the location of GW events, which we have assumed to be independent of the population properties of the GW events. Eq. (14) signifies that we are only taking into account the detected events and not using the events that are below our SNR threshold. While this introduces a bias for population inference, our method is expected to be free from such biases, since the overdensity maps (or the 3-D power spectrum) of GWs, for both detected or undetected events, are the same by our construction of the mock events. Hence we are not introducing any bias by leaving out subthreshold events. Even in the case that there are correlations between GW event properties and their host galaxy properties, those will be absorbed by the GW overdensity bias, $b_{GW}$, and its redshift evolution. For our simulated catalog, we have currently set it to unity. Similarly, the effect of flux limited galaxy surveys will have to be accounted for through the galaxy bias term $b_g$. We plan to examine such effects and their impact on cosmological inference in future work.

When constructing the posterior in Eq. (12), we take the flat $\Lambda$CDM cosmology to be the true cosmological model, with $H_0$ as a free parameter and the matter density $\Omega_m$ kept fixed at 0.307. In the flat $\Lambda$CDM Universe, the relation between luminosity distance $d_L$ and the redshift $z$ in Eq. (6) is used to convert the measured luminosity distance of each GW event to its redshift for a given $H_0$ whenever required while performing the integral over $d_{gw}$ in Eq. (12).

We would like to emphasize here that, although this work stems from the same core idea (i.e., leveraging the clustering of GW events with other large-scale structure tracers) as outlined in Refs. [8, 44], the methodology presented here is not only significantly different but also incorporates several improvements over the previous work. As an example, consider Ref. [44], which uses the angular cross-correlation between GW events and galaxies in different redshift bins to infer the Hubble constant in two steps. In the first step, the GW events are binned in multiple luminosity distance bins based on the median luminosity distance estimate from the posterior distribution of parameters given the GW data. The GW events in each bin are cross-correlated with galaxies in narrow redshift bins in order to estimate the mean redshift of the GW events in that bin. The peak of the angular cross-correlation technique was used to infer the luminosity distance-redshift diagram. In the second step, the

luminosity distance-redshift relation thus obtained was fit in the context of an assumed model ($\Lambda$CDM in this case) to constrain the Hubble constant using the standard Bayesian inference technique. In contrast, the method we present in this work involves modeling the galaxy overdensity for each GW event while accounting for its volume uncertainty region and estimating the Hubble constant by modeling the 3D cross-correlation. We provide a detailed comparison of these two methodologies in the Appendix.

### III. SIMULATION

We test the framework discussed in Sec. II with a mock galaxy and GW catalog in order to assess how well the expansion rate can be recovered given realistic localization errors in the positions of GW events. We consider dark matter halos from the Big MultiDark Planck (BigMDPL) cosmological N-body simulation [65] from the CosmoSim database<fnref n="2" /> to construct the galaxy catalog. The BigMDPL simulation follows the evolution of the matter distribution in a flat $\Lambda$CDM cosmology with the following parameters: Hubble parameter $h = H_0/(100 \text{ km s}^{-1} \text{ Mpc}^{-1}) = 0.6777$, the matter density parameter $\Omega_m = 0.307$, the amplitude of density fluctuations characterized by $\sigma_8 = 0.823$, and the power spectrum slope of initial density fluctuations $n_s = 0.96$. The simulation is set up as 3840³ collisionless particles in a cubic box with comoving side length $2.5h^{-1}\text{Gpc}$ with periodic boundary conditions and corresponds to a mass resolution of $2.359 \times 10^{10}h^{-1}M_{\odot}$.

We only consider well-resolved dark matter halos at the $z=0$ snapshot of BigMDPL, with masses $\ge M_{\text{th}} = 10^{12}h^{-1}M_{\odot}$. For simplicity, we assume that all galaxies are central galaxies and thus place each galaxy at the center of these dark matter halos. We place an observer at the center of the box and compute the sky positions of each of these mock galaxies with respect to the center, which is within the sphere of radius $1.25h^{-1}\text{Gpc}$. We utilize the cosmological parameters of the simulation to determine the redshift of each galaxy from its comoving distance from the central observer. We use this mock galaxy catalog for further analysis. For simplicity, we have ignored any redshift space distortion effects. In practice, galaxy catalogs obtained from real galaxy surveys suffer from incompleteness issues due to the sky-position dependent flux limits. Flux-limited galaxy catalogs must be considered to understand the effect of the limited number of observed fainter galaxies at higher redshift. In such cases, one has to account for how the flux limit introduces an evolution in the galaxy bias factor. This is expected to broaden the $H_0$ posterior compared to that obtained using a complete galaxy catalog. Similar effects

<footnote n="2">² This is publicly available at https://www.cosmosim.org/.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Probability distribution of observed overdensity divided by the corresponding fluctuation in different redshift bins, mentioned at the top of each panel, considering $\theta_{max} = 0.03$ radian.</figcaption></figure>

are also expected for galaxy catalogs that have limited sky coverage. For such cases, the effect of appropriate sky masks needs to be considered while modeling the measurements. For simplicity, we consider a complete galaxy catalog with full sky coverage in this study.

Once we have the mock galaxy catalog, we randomly subselect galaxies to be the hosts of GW events without any regard to the properties of these halos or, equivalently, the properties of the galaxies in the catalog. The true redshifts and sky positions of the GW events are determined from their respective host galaxies in the cosmological box. We use the same cosmological parameters of the simulation box in order to compute the luminosity distance of GW events. The galaxy catalog we use extends up to a comoving distance of $1.25h^{-1}$Gpc given the size of the box before we run into incompleteness issues. This comoving distance corresponds to redshift $z \sim 0.46$, assuming the true cosmology. In this work, however, we consider uniform prior over $H_0$ between 20 and 120 km s$^{-1}$ Mpc$^{-1}$. Consequently, a higher value of $H_0$ can lead to an increase in the redshift of a GW event at a given luminosity distance. Thus, the GW event can have redshift up to $z \sim 0.36$ when $H_0 = 120$ km s$^{-1}$ Mpc$^{-1}$ and yet be detectable in the box. Therefore, we restrict the injection of GW events within a luminosity distance of 1 Gpc to avoid any biases creeping in due to the above selection effects.

The distribution of the source-frame masses of the BBHs follows the Power Law+Peak model [66] between minimum mass $m_{min} = 5M_{\odot}$ and maximum mass $m_{max} = 60M_{\odot}$. The injected values of other mass hyperparameters (see Appendix B of Ref. [62] for the descriptions of the hyperparameters) are $\lambda_{peak} = 0.06$, $\alpha = 3.5$, $\delta_m = 4$, $\mu_m = 34$, $\sigma_m = 5$, and $\beta = 1.1$, which are

consistent with the previous observations [62, 63]. However, note that the detected masses are redshifted and are expected to differ from the source-frame mass, $m$, as $m_z = m(1+z)$, where $z$ is the redshift of the event. The other BBH parameters (see Table 1 of Ref. [67] for the list of parameters) are randomly chosen from the default priors of *bilby* (see Table 1 of Ref. [67] for the default priors). The simulated BBH signal is added to stationary Gaussian noise, colored with the O4 design sensitivity [68] of LIGO and Virgo detectors. We analyze the strain signal between 20 and 2048 Hz to infer the source parameters. We consider the same non-precessing waveform model *IMRPhenomD* [69] for both the signal injection and recovery $^3$. We have used the nested sampler *dynesty* [72] implemented in *bilby* [67] to infer the posterior distribution of all the parameters given the waveform. We consider uniform prior over observed chirp mass within the range $M_{c,inj}^{z} \pm 5 M_{\odot}$, with the constraint that the component mass remains above 2.5 $M_{\odot}$ and over mass-ratio $q = m_2/m_1 \in [0.001, 1]$. These choices of mass parameters and their specified prior ranges optimize sampling efficiency while lowering computational costs. We employ cosmology independent $d_L^2$ prior to the luminos-

<footnote>$^3$ Note that we consider these simplified waveform model, as we would like to demonstrate a proof-of-concept for the cross-correlation based framework that we present in this paper. Analysis of real events will involve upgraded waveform models, potentially including higher-order modes [70] and precession [71], which can better constrain the luminosity distance by breaking the degeneracy between luminosity distance and inclination angle. In this context, the $H_0$ measurement can be regarded as a conservative estimate for the simulated events considered in this work.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Observed overdensity fluctuation as a function of redshift over $10^6$ random line of sights.</figcaption></figure>

ity distance of the source between 10 Mpc and 5 Gpc. The priors over the rest of the parameters (except for the mass parameters and luminosity distance) are the default priors implemented in bilby (refer to Table 1 in Ref. [67] for the default priors).

The parameters thus inferred are utilized to constrain $H_0$ by cross-correlating the localization of these individual events with galaxies in the mock catalog. We are primarily interested in the extrinsic source parameters, such as $d_L^{gw}$, $\alpha^{gw}$, and $\delta^{gw}$ in this work. So, we compute the GW likelihood $\mathcal{L}_{gw}$ by marginalizing over the rest of the intrinsic and extrinsic parameters. The galaxy overdensity likelihood $\mathcal{L}_g$ is also calculated for different sky coordinates within the localization of the GW events as a function of redshift in order to perform the integration of Eq. (12). The galaxy overdensity likelihood includes the computation of the observed and modeled galaxy overdensity fields as a function of redshift. The modeled galaxy overdensity is calculated from the cross-correlation function $\xi$ (see, Eq. (9)), which depends on the separation (comoving distance in $h^{-1}$ scale) between the GW event and the redshift where the observed galaxy overdensity is also measured. We have used the publicly available code AUM$^4$ to compute $\xi$ [73] from the nonlinear matter power spectrum [61] assuming the flat $\Lambda$CDM cosmology. To study the effect of the cross-correlation length scale on the $H_0$ estimation, we assume three different values of the angular radius $\theta_{max}$ (elaborated in Sec. IV) along the line of sight to calculate the observed and modeled galaxy overdensity fields.

In order to compute the likelihood of observing the galaxy overdensity given the modeled overdensity, we quantify the variance of the galaxy overdensity field [denoted as $\sigma_\delta^2$ in Eq. (12)] around random lines of sight

at each redshift bin. The variance of the overdensity depends upon the specific value of $\theta_{max}$ over which we average. We use $10^6$ random lines of sight within the mock galaxy catalog used in this work. The probability distribution of $\delta/\sigma^5$ for eight representative redshift bins of width $\Delta z = 0.01$ up to $z = 0.41$ are shown in Fig. 2, for $\theta_{max} = 0.03$ radian. The figure clearly shows that $\delta/\sigma$ becomes non-Gaussian at lower redshifts. Despite this, we assume a Gaussian likelihood for galaxy overdensity for simplicity and find that this assumption does not affect $H_0$ estimation in this study, given the sensitivity of the GW detectors considered. Developing a more rigorous framework for modeling the galaxy overdensity likelihood can be explored in future work. The corresponding standard deviation of the overdensity field is shown as a function of redshift as solid lines in Fig. 3, where the colors signify different values of $\theta_{max}$. As expected, $\sigma$ decreases with redshift for any $\theta_{max}$ due to the increase in cosmological volume probed within the angular region per redshift interval. Similarly, a higher value of $\theta_{max}$ also results in a smaller standard deviation for the overdensity at any given redshift. Notably, the redshift width of $\Delta z = 0.01$ for calculating galaxy overdensity has been carefully selected to ensure a sufficient number of galaxies in each redshift bin. We empirically choose the optimal value of $\Delta z = 0.01$ by experimenting with various choices, optimizing the impact of Poisson noise on the overdensity measurement for smaller $\Delta z$ and the diminished cross-correlation signal for larger $\Delta z$. In a realistic scenario, this choice needs to be adjusted depending on the galaxy number densities used for the cross-correlation. We use this variance in order to evaluate the likelihood of the galaxy overdensity given the modeled overdensity and a true position of the GW event, which, among other parameters, depends upon the Hubble constant. In order to infer the posterior distribution of the Hubble constant given the data, we assume a uniform prior in $H_0 \in U[20, 120]$ km $s^{-1}$ Mpc$^{-1}$. We finally infer the joint posterior for $H_0$ for multiple events as given by Eq. (14).

## IV. RESULTS

We apply our method as mentioned in Sec. II to three different mock GW catalogs for three different choices of $\theta_{max} = \{0.01, 0.02, 0.03\}$ radian to study the statistical robustness. The GW catalogs are random realizations of the same population, as discussed in Sec. III. Each GW catalog consists of 250 simulated GW events that are detected with a network signal-to-noise ratio (SNR) $\rho_{net} \ge 12$ in the three LIGO-Virgo detectors, with O4 design sensitivity [68]. We also study $H_0$ posteriors for two subsets of 150 and 200 GW events from those three sets

<footnote>$^4$ http://surhudm.github.io/aum/index.html</footnote>

<footnote>$^5$ We are omitting superscripts and subscripts.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Estimation of Hubble constant from BBH events whose network SNR $\rho_{\text{net}} \ge 12$, as detected by LIGO-Virgo detector operating at O4 sensitivity, for different choices of $\theta_{\text{max}}$ (in radian). The black vertical dashed line in each plot corresponds to the injected value of $H_0 = 67.77 \text{ km s}^{-1}\text{Mpc}^{-1}$.</figcaption></figure>

of 250 detected GW events. The threshold on network SNR in this work is chosen based on the earlier LVK cosmology papers [14, 30]. Sources with low network SNR would not contribute significantly to the improvement in the inference of the Hubble constant due to the poor sky localization and larger uncertainty in distance measurement. Figure 4 shows the $H_0$ posteriors for three different sets of GW events with varying numbers of detected events. Each row of Fig. 4 corresponds to a particular value of $\theta_{\text{max}}$, mentioned in the right of the figure for different numbers of GW events. The columns of Fig. 4 correspond to $H_0$ posteriors for a fixed number of GW events, mentioned at the top of that panel, with varying

$\theta_{\text{max}}$. The black vertical dashed lines correspond to the injected value of $H_0 = 67.77 \text{ km s}^{-1}\text{Mpc}^{-1}$ in Fig. 4.

We also summarize the uncertainty corresponding to 90% highest density interval in constraining $H_0$ in Fig. 6 for a different number of detected GW events. The different panels correspond to varying numbers of detected events $N_{\text{GW}}$. In each panel, the error bar of $H_0$ uncertainty for three different sets is depicted with three different colors. Set 1 and 3 have offsets of -0.001 and 0.001 radian, respectively, corresponding to the true $\theta_{\text{max}}$ for clarity in the representation. In Fig. 6, the black horizontal dashed lines represent the injected value of $H_0 = 67.77 \text{ km s}^{-1}\text{Mpc}^{-1}$.
<figure><img src="image_6.png" /><figcaption>FIG. 5. Posterior distribution of the Hubble constant inferred using $\theta_{\text{max}} = 0.02$ radian from one realization of 50 randomly chosen BBH events, each with a network SNR $\rho_{\text{net}} \ge 12$, as detected by the LIGO-Virgo detectors operating at O4 sensitivity. The accuracy with which $H_0$ is inferred in this realization is representative amongst the numerous random realizations we constructed. The median and 90% highest density interval are indicated at the top of the plot. The black vertical dashed line in the plot represents the injected value of $H_0 = 67.77 \text{ km s}^{-1} \text{ Mpc}^{-1}$.</figcaption></figure>

As expected, the posterior of $H_0$ starts to converge to the true value of $H_0 = 67.77 \text{ km s}^{-1} \text{ Mpc}^{-1}$ with an increasing number of GW events. However, for the smaller number of GW events $\sim 150$, the $H_0$ posteriors can show a wide variation. We speculate that the non-Gaussian nature of galaxy overdensity fluctuation $\sigma_\delta$ (see Fig 2) and not accounting for the covariance between different redshift bins is likely the reason for this wide variation. Though this is not an issue for unbiased estimation of $H_0$ with larger GW events, it is an aspect that needs to be carefully explored in the future. We have studied $H_0$ posteriors for larger $\theta_{\text{max}} > 0.03$ radian and found the posteriors degrade even when the number of GW events is significant, i.e., $\gtrsim 150$. This is likely related to the fact that large values of $\theta_{\text{max}}$ would decrease the amplitude of the correlation signal, reducing the total signal-to-noise ratio.

On the other hand, we have observed it becomes difficult to capture the correlation signal for smaller values of $\theta_{\text{max}} \lesssim 0.005$ radian due to the significant fluctuation in measuring the clustering for a smaller number of galaxies, especially at low redshift. So, it is crucial to explore the optimal value of $\theta_{\text{max}}$ to perform this method for unbiased estimation of the Hubble constant. From Fig. 6, we observe that $\theta_{\text{max}} = 0.02$ radian may be the optimal choice in our study. However, the optimal value of $\theta_{\text{max}}$ may vary for different GW events. In this work, we keep the values of $\theta_{\text{max}}$ fixed for all GW events for simplicity. We leave the exercise of finding the optimal $\theta_{\text{max}}$ given the localization area on the sky to investigations in the future.

<table><caption>TABLE I. Summary of uncertainty (maximum a posteriori probability with 90% highest density interval) of inferred values of the Hubble constant from different sets of GW events with varying $\theta_{\text{max}}$.</caption><thead><tr><th rowspan="2">N<sub>GW</sub></th><th rowspan="2">θ<sub>max</sub> (radian)</th><th colspan="3">H<sub>0</sub>(km s<sup>-1</sup> Mpc<sup>-1</sup>)</th></tr><tr><th>Set 1</th><th>Set 2</th><th>Set 3</th></tr></thead><tbody><tr><td rowspan="3">150</td><td>0.01</td><td>67.30<sup>+15.69</sup><sub>-7.84</sub></td><td>65.97<sup>+18.19</sup><sub>-15.97</sub></td><td>71.41<sup>+10.59</sup><sub>-5.77</sub></td></tr><tr><td>0.02</td><td>68.10<sup>+5.74</sup><sub>-5.10</sub></td><td>60.90<sup>+12.26</sup><sub>-9.90</sub></td><td>69.57<sup>+13.50</sup><sub>-5.57</sub></td></tr><tr><td>0.03</td><td>68.65<sup>+5.91</sup><sub>-4.64</sub></td><td>58.73<sup>+9.27</sup><sub>-6.29</sub></td><td>90.28<sup>+12.72</sup><sub>-23.01</sub></td></tr><tr><td rowspan="3">200</td><td>0.01</td><td>68.76<sup>+6.71</sup><sub>-4.76</sub></td><td>55.24<sup>+13.59</sup><sub>-4.24</sub></td><td>68.02<sup>+4.97</sup><sub>-3.19</sub></td></tr><tr><td>0.02</td><td>68.37<sup>+4.09</sup><sub>-3.37</sub></td><td>63.34<sup>+6.58</sup><sub>-7.34</sub></td><td>67.64<sup>+4.35</sup><sub>-3.51</sub></td></tr><tr><td>0.03</td><td>68.85<sup>+4.14</sup><sub>-3.23</sub></td><td>61.42<sup>+11.07</sup><sub>-7.42</sub></td><td>68.22<sup>+4.78</sup><sub>-3.90</sub></td></tr><tr><td rowspan="3">250</td><td>0.01</td><td>69.66<sup>+7.41</sup><sub>-4.67</sub></td><td>68.53<sup>+4.92</sup><sub>-4.53</sub></td><td>67.82<sup>+4.84</sup><sub>-2.82</sub></td></tr><tr><td>0.02</td><td>69.68<sup>+3.83</sup><sub>-3.68</sub></td><td>66.07<sup>+4.93</sup><sub>-3.48</sub></td><td>66.07<sup>+4.92</sup><sub>-3.48</sub></td></tr><tr><td>0.03</td><td>69.84<sup>+4.15</sup><sub>-3.54</sub></td><td>66.30<sup>+6.18</sup><sub>-4.30</sub></td><td>68.00<sup>+4.00</sup><sub>-3.34</sub></td></tr></tbody></table>

We report the 90% highest density interval of the inferred values of $H_0$ corresponding to the three different sets in Table I for three different choices of $\theta_{\text{max}} = \{0.01, 0.02, 0.03\}$ radian. The corresponding $H_0$ posteriors are also shown in Fig. 6. We observe that we can constrain $H_0$ with a precision of $\lesssim 8\%$ (90% highest density interval) by combining the individual $H_0$ posteriors from 250 GW events. In the best-case scenario, $H_0$ can be measured as precisely as $\sim 5\%$. For a smaller number of 150 GW events, the uncertainty in measuring $H_0$ varies between $8\% - 14\%$ across different sets of GW events when considering the optimal $\theta_{\text{max}} = 0.02$ radian. However, the uncertainty in $H_0$ increases to $\sim 25\%$ with other choices of $\theta_{\text{max}}$. Following Ref. [44], we also investigate the impact of our method on a smaller sample of nearby events, specifically 50 events within 700 Mpc, and present a representative example of the $H_0$ posterior for one of the realizations with a median accuracy of about 11%, in Fig. 5. This precision is similar to that obtained in Ref. [44], even after our careful accounting for the localization uncertainty. We speculate that the uncertainties degrade only slightly due to our improved treatment of the cross-correlation of galaxies with individual GW events, which compensates for the possible degradation due to localization uncertainty. However, we note that, with 50 events, a minority of posteriors exhibit multimodal nature with a significantly increased uncertainty. We have cross-checked that by adding more events, the posterior converges to the correct mode. This suggests that measuring the Hubble constant from fewer events could lack statistical reliability, which could be remedied by increasing the sample of GW events to be as large as 250. It is important to note that we consider GW events within 1 Gpc as the mock galaxy catalog used in this work is volume limited, extending up to redshift $z \sim 0.46$, already detailed in Sec. III. The sky localizations of these simulated GW events are better than those from larger distances, which are observed by the current generation detectors. So, the broad sky localizations of GW events from significantly large distances may affect
<figure><img src="image_7.png" /><figcaption>FIG. 6. Uncertainty corresponding to 90% highest density interval in estimating $H_0$ from different subsets of GW events for various values of $\theta_{\max}$. The dashed horizontal line in each of the plots corresponds to the true value of $H_0 = 67.77 \text{ km s}^{-1} \text{ Mpc}^{-1}$. Each plot corresponds to three different sets of a fixed number of GW events, mentioned at the top of the corresponding plot. $H_0$ error bar plot of set 2 GW events corresponds to the correct $\theta_{\max}$. However, $H_0$ error bar plot of set 1 and set 3 GW events are plotted with the offsets of -0.001 and 0.001 radian in $\theta_{\max}$. These offsets are used only for illustration purposes.</figcaption></figure>

$H_0$ estimation and result in broadening the $H_0$ posterior.

V. DISCUSSION

Observations of GWs from coalescing binaries provide the direct measurement of the luminosity distance unaffected by systematics in various external astrophysical calibrations that rely on the cosmic distance ladder. However, GW estimates can be biased due to systematics inherent to its measurements, such as GW detector calibration [74, 75]. Assuming that the instrument can be calibrated sufficiently accurately [74, 75], GW sources can be used to estimate the Hubble constant and other cosmological parameters. However, redshift information is not available from GW observations. Additionally, most of the detected GW events are BBH mergers and are not accompanied by EM counterparts. In this paper, we successfully demonstrated the inference of the Hubble constant from individual GW events by cross-correlating with the galaxy catalog. This method is useful even if the mergers do not have any EM observations and can supplement the inference of the Hubble constant from BNS mergers with detected EM counterparts.

In this work, we have considered varying numbers of GW events in the range of [150, 250] to assess the performance of our method. The total number of GW events that we expect when we consider events from GWTC-3 as well as the ongoing fourth observation (O4) run of the LVK detectors[7, 63] are similar. However, they will be distributed in a much larger volume than the 1 Gpc in which we populate our simulated sources. We have restricted ourselves to the 1 Gpc due to the limitations of the simulation volume available to us. The method

that we show in our paper is, in principle, expected to be applicable to all GW events. The precision of the $H_0$ estimation depends on the population of GW events as the strength of the cross-correlation signal between GW events and galaxies is expected to diminish for GW events at larger distances due to the larger volume uncertainty. Given that the estimation of the Hubble constant (in particular, the credible intervals) is not robust (see Sec. IV), we refrain from analyzing the current GW events from GWTC-3 due to the small sample size. Once the new catalog of GW events, including those from O4, is available, we will explore the use of this method for the inference of the Hubble constant.

Our method consists of measuring the galaxy overdensity as a function of redshift at the maximum *a posteriori* probability in the sky localization of individual gravitational wave events. We compare the data vector of observed galaxy overdensity $d_g^{\text{obs}}$ with that of the modeled galaxy overdensity $d_g^{\text{mod}}$ which depends on the value of $H_0$, to define a likelihood and infer the posterior distribution of $H_0$ in a Bayesian analysis. As examples, we considered three different mock GW catalogs to infer $H_0$ using the method described in Sec. II and were able to constrain the Hubble constant with an accuracy of $\lesssim 8\%$ (90% highest density interval). The efficacy of this method is tested for the second-generation detector network comprised of the Advanced LIGO [4] and Advanced Virgo [5] detectors. With the improved sensitivity in the third-generation detectors like Cosmic Explorer [76] and Einstein Telescope [77], GW sources will be detected with much higher detection rate, and precise estimation of luminosity distance and sky localization [78–82]. We expect this method to constrain $H_0$ with higher precision with such detectors.
In our current analysis, we have only considered $H_0$ as a free parameter to be constrained from the GW and galaxy overdensity observations. However, it is straightforward to extend our analysis to constrain other cosmological parameters, such as the matter density at the current epoch and dark energy equation of state. We do not consider them here to avoid the extra computational costs involved in estimating cosmological parameter(s). With the second-generation GW detectors, cosmological parameters except for $H_0$ are expected to be weakly constrained [30].

Our results show that the cross-correlation method can be utilized to infer $H_0$ from individuals from GW events as a proof of concept. We have made several simplifying assumptions in this work. The mock galaxy catalog used here is assumed to be volume limited with uniform depth over the entire sky and no limit on the galaxy fluxes. In contrast, real galaxy catalogs are flux-limited and also have varying depths in different directions. We defer the application of our method to galaxy catalogs, which include these imperfections in future work. We have also considered the idealistic scenario where the GW events are the unbiased tracers of underlying matter distribution. In reality, the GW events and the galaxies trace the large-scale structure with a redshift-dependent bias [83, 84]. This redshift-dependent bias could arise from the evolving way in which the gravitational wave events populate dark matter halos or due to selection effects. These effects need to be parametrized and marginalized to constrain the cosmological parameters from real data. We also have not accounted for any impact that weak gravitational lensing will have on the luminosity distance estimates due to the intervening large-scale structure. Weak lensing induces nontrivial spatial correlations on the sky with galaxies at redshifts lower than that of the GW events, albeit with a smaller amplitude than the physical cross-correlation with galaxies at the same redshift as the GW events. This effect could be additionally incorporated into our methodology while working with real data, especially for high redshift events, by adding in a cross-correlation term related to weak lensing [43]. Furthermore, weak lensing can also result in an increased covariance for the measurement of the overdensity. Our method of utilizing random lines of sight would naturally account for this increased covariance. The exact quantitative impact of such effects, however, will require full sky galaxy catalogs with weak lensing shears computed at multiple different lens planes and we defer such a study to future work.

## ACKNOWLEDGMENTS

The authors acknowledge Divya Rana for helpful discussions throughout this work. The authors would like to thank Simone Mastrogiovanni, Aditya Vijaykumar, and Anuj Mishra for carefully reading the manuscript and useful suggestions. We also thank the anonymous references for the careful review of the manuscript and the valuable comments.

T.G. acknowledges using the IUCAA LDG cluster Sarathi, accessed through the LVK collaboration, for the computation involved in this work. The work of S. Bera was supported by the Universitat de les Illes Balears (UIB); the Spanish Agencia Estatal de Investigación Grants No. PID2022-138626NB-I00, PID2019-106416GB-I00, RED2022-134204-E, RED2022-134411-T, funded by MCIN/AEI/10.13039/501100011033; the MCIN with funding from the European Union NextGenerationEU/PRTR (PRTR-C17.I1); Comunitat Autonòma de les Illes Balears through the Direcció General de Recerca, Innovació I Transformació Digital with funds from the Tourist Stay Tax Law (PDR2020/11-ITS2017-006), the Conselleria d'Economia, Hisenda i Innovació grant numbers SINCO2022/18146 and SINCO2022/6719, cofinanced by the European Union and FEDER Operational Program 2021-2027 of the Balearic Islands. S. Bera received partial support from the French government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix Marseille Université - A*MIDEX AMX-22-CEI-02. S. Bose acknowledges support from the NSF under Grant No. PHY-2309352.

The MultiDark Database used in this paper and the web application providing online access to it were constructed as part of the activities of the German Astrophysical Virtual Observatory as a result of a collaboration between the Leibniz-Institute for Astrophysics Potsdam (AIP) and the Spanish MultiDark Consolider Project No. CSD2009-00064. The Big MultiDark Planck simulation has been performed in the Supermuc supercomputer at LRZ using time granted by PRACE.

## Appendix A: DIFFERENCES AND IMPROVEMENTS OVER THE PREVIOUS WORK

A general overview of the differences between this study and Ref. [44], both of which utilize the clustering of GW events with galaxies, pertaining to their methodology has been highlighted at the end of Sec.II. Below, we provide more details comparing the two studies with regard to the simulated GW source population as well as the analysis technique.

To test the clustering method, in Ref. [44], the authors used simulated GW events up to a distance ~ 1300 Mpc (or, $900h^{-1}$ Mpc). They divided the GW events into bins of width 200 Mpc, based on the median luminosity distance of the posterior distribution of the luminosity distance. For a sample of GW events in a given luminosity distance bin, they then proceeded to compute the angular cross-correlation with galaxies [see Eq. (13) in Ref. [44] for details on the estimator of angular cross-correlation] using the maximum probable posterior sky position of each GW event. The uncertainties in the angular position of the GW event are expected to reduce the ampli-
tude of the cross-correlation signal, but this was never modeled in Ref. [44]. Instead, they utilized the fact that the angular cross-correlation for each distance bin would still peak at the correct redshift corresponding to the GW events in the luminosity distance bin and modeled it as a Gaussian distribution [Eq. (14) in Ref. [44]], to estimate the mean redshift and its error. For each bin, they assigned an error-weighted average luminosity distance based on the errors in the luminosity distance from the posterior distribution of parameters given the GW data. The resultant Hubble diagram (luminosity distance versus redshift) was then used to estimate $H_0$. Their study thus used a variety of simplifications to demonstrate a proof of concept of how the clustering of gravitational waves with galaxies can be used to infer the Hubble constant.

In this work, we incorporate several further improvements. In contrast with the earlier work, we adopt a single-step Bayesian inference method, eliminating the need to bin the GW events into luminosity distance bins. Instead of measuring the angular cross-correlation of a sample of GW events with galaxies, we measure the overdensity of galaxies around each GW event separately.

The biggest difference is in our modeling scheme. We model the 3D cross-correlation of GW events with galaxies, thus predicting both the amplitude and shape of how the overdensities around GW events depend upon the cosmological parameters. In our work, we also account for the posterior distribution of the sky localization as well as the luminosity distance while modeling the observables, thus naturally accounting for the smearing of the angular cross-correlation function. In this manner, any non-Gaussianities in the posterior distributions, as well as the correlated errors in the localization volume, are automatically incorporated in our study.

Finally, we also comment on the precision of the inference of the Hubble constant obtained in the two studies. In the earlier work, a large number of GW events were distributed between 200 and 1300 Mpc, out of which 50 and 500 (distributed up to ~700 Mpc) sources were selected to infer the Hubble constant with an accuracy of ~10% and 2%, respectively. In the present work, in addition to analyzing 250 events out to 1 Gpc, we have also repeated the exercise for 50 events and shown in the main text (see Fig. 5) that we get similar accuracy for the inference on $H_0$ as in the earlier study even after marginalizing over the localization uncertainty.

[1] B. P. Abbott et al. (LIGO Scientific, Virgo), “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc].

[2] B. P. Abbott et al. (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), “Multi-messenger Observations of a Binary Neutron Star Merger,” Astrophys. J. Lett. 848, L12 (2017), arXiv:1710.05833 [astro-ph.HE].

[3] B. P. Abbott et al. (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), “A gravitational-wave standard siren measurement of the Hubble constant,” Nature 551, 85–88 (2017), arXiv:1710.05835 [astro-ph.CO].

[4] J. Aasi et al. (LIGO Scientific), “Advanced LIGO,” Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc].

[5] F. Acernese et al. (VIRGO), “Advanced Virgo: a second-

generation interferometric gravitational wave detector,” Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc].

[6] T. Akutsu et al. (KAGRA), “KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector,” Nature Astron. 3, 35–40 (2019), arXiv:1811.08079 [gr-qc].

[7] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run,” Phys. Rev. X 13, 041039 (2023), arXiv:2111.03606 [gr-qc].

[8] B. F. Schutz, “Determining the Hubble constant from gravitational wave observations,” Nature (London) 323, 310–311 (1986).

[9] Walter Del Pozzo, “Inference of the cosmological parameters from gravitational waves: application to second generation interferometers,” Phys. Rev. D 86, 043011 (2012), arXiv:1108.1317 [astro-ph.CO].

[10] Hsin-Yu Chen, Maya Fishbach, and Daniel E. Holz, “A two per cent Hubble constant measurement from standard sirens within five years,” Nature 562, 545–547 (2018), arXiv:1712.06531 [astro-ph.CO].

[11] M. Fishbach et al. (LIGO Scientific, Virgo), “A Standard Siren Measurement of the Hubble Constant from GW170817 without the Electromagnetic Counterpart,” Astrophys. J. Lett. 871, L13 (2019), arXiv:1807.05667 [astro-ph.CO].

[12] M. Soares-Santos et al. (DES, LIGO Scientific, Virgo), “First Measurement of the Hubble Constant from a Dark Standard Siren using the Dark Energy Survey Galaxies and the LIGO/Virgo Binary-Black-hole Merger GW170814,” Astrophys. J. Lett. 876, L7 (2019), arXiv:1901.01540 [astro-ph.CO].
[13] Rachel Gray et al., "Cosmological inference using gravitational wave standard sirens: A mock data analysis," Phys. Rev. D 101, 122001 (2020), arXiv:1908.06050 [gr-qc].

[14] B. P. Abbott et al. (LIGO Scientific, Virgo, VIRGO), "A Gravitational-wave Measurement of the Hubble Constant Following the Second Observing Run of Advanced LIGO and Virgo," Astrophys. J. 909, 218 (2021), arXiv:1908.06060 [astro-ph.CO].

[15] R. Abbott et al. (LIGO Scientific, Virgo), "GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object," Astrophys. J. Lett. 896, L44 (2020), arXiv:2006.12611 [astro-ph.HE].

[16] Rachel Gray, Chris Messenger, and John Veitch, "A pixelated approach to galaxy catalogue incompleteness: improving the dark siren measurement of the Hubble constant," Mon. Not. Roy. Astron. Soc. 512, 1127-1140 (2022), arXiv:2111.04629 [astro-ph.CO].

[17] Andreas Finke, Stefano Foffa, Francesco Iacovelli, Michele Maggiore, and Michele Mancarella, "Cosmology with LIGO/Virgo dark sirens: Hubble parameter and modified gravitational wave propagation," JCAP 08, 026 (2021), arXiv:2101.12660 [astro-ph.CO].

[18] Antonella Palmese, Clecio R. Bom, Sunil Mucesh, and William G. Hartley, "A Standard Siren Measurement of the Hubble Constant Using Gravitational-wave Events from the First Three LIGO/Virgo Observing Runs and the DESI Legacy Survey," Astrophys. J. 943, 56 (2023), arXiv:2111.06445 [astro-ph.CO].

[19] Jonathan R. Gair et al., "The Hitchhiker's Guide to the Galaxy Catalog Approach for Dark Siren Gravitational-wave Cosmology," Astron. J. 166, 22 (2023), arXiv:2212.08694 [gr-qc].

[20] V. Alfradique et al., "A dark siren measurement of the Hubble constant using gravitational wave events from the first three LIGO/Virgo observing runs and DELVE," Mon. Not. Roy. Astron. Soc. 528, 3249-3259 (2024), arXiv:2310.13695 [astro-ph.CO].

[21] W. Ballard et al. (DESI), "A Dark Siren Measurement of the Hubble Constant with the LIGO/Virgo Gravitational Wave Event GW190412 and DESI Galaxies," Res. Notes AAS 7, 250 (2023), arXiv:2311.13062 [astro-ph.CO].

[22] David F. Chernoff and Lee Samuel Finn, "Gravitational radiation, inspiraling binaries, and cosmology," Astrophys. J. Lett. 411, L5-L8 (1993), arXiv:gr-qc/9304020.

[23] Stephen R. Taylor, Jonathan R. Gair, and Ilya Mandel, "Hubble without the Hubble: Cosmology using advanced gravitational-wave detectors alone," Phys. Rev. D 85, 023535 (2012), arXiv:1108.5161 [gr-qc].

[24] Will M. Farr, Maya Fishbach, Jiani Ye, and Daniel Holz, "A Future Percent-Level Measurement of the Hubble Expansion at Redshift 0.8 With Advanced LIGO," Astrophys. J. Lett. 883, L42 (2019), arXiv:1908.09084 [astro-ph.CO].

[25] Zhi-Qiang You, Xing-Jiang Zhu, Gregory Ashton, Eric Thrane, and Zong-Hong Zhu, "Standard-siren cosmology using gravitational waves from binary black holes," Astrophys. J. 908, 215 (2021), arXiv:2004.00036 [astro-ph.CO].

[26] S. Mastrogiovanni, K. Leyde, C. Karathanasis, E. Chassande-Mottin, D. A. Steer, J. Gair, A. Ghosh, R. Gray, S. Mukherjee, and S. Rinaldi, "On the importance of source population models for gravitational-

wave cosmology," Phys. Rev. D 104, 062009 (2021), arXiv:2103.14663 [gr-qc].

[27] Michele Mancarella, Edwin Genoud-Prachex, and Michele Maggiore, "Cosmology and modified gravitational wave propagation from binary black hole population models," Phys. Rev. D 105, 064030 (2022), arXiv:2112.05728 [gr-qc].

[28] Jose María Ezquiaga and Daniel E. Holz, "Spectral Sirens: Cosmology from the Full Mass Distribution of Compact Binaries," Phys. Rev. Lett. 129, 061102 (2022), arXiv:2202.08240 [astro-ph.CO].

[29] B. P. Abbott et al. (LIGO Scientific, Virgo), "Properties of the binary neutron star merger GW170817," Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc].

[30] R. Abbott et al. (LIGO Scientific, Virgo,, KAGRA, VIRGO), "Constraints on the Cosmic Expansion History from GWTC-3," Astrophys. J. 949, 76 (2023), arXiv:2111.03604 [astro-ph.CO].

[31] Simone Mastrogiovanni, Danny Laghi, Rachel Gray, Giada Caneva Santoro, Archisman Ghosh, Christos Karathanasis, Konstantin Leyde, Daniele A. Steer, Stephane Perries, and Gregoire Pierra, "Joint population and cosmological properties inference with gravitational waves standard sirens and galaxy surveys," Phys. Rev. D 108, 042002 (2023), arXiv:2305.10488 [astro-ph.CO].

[32] Rachel Gray et al., "Joint cosmological and gravitational-wave population inference using dark sirens and galaxy catalogues," JCAP 12, 023 (2023), arXiv:2308.02281 [astro-ph.CO].

[33] Federico Stachurski, Christopher Messenger, and Martin Hendry, "Cosmological inference using gravitational waves and normalizing flows," Phys. Rev. D 109, 123547 (2024), arXiv:2310.13405 [gr-qc].

[34] Nicola Borghi, Michele Mancarella, Michele Moresco, Matteo Tagliazucchi, Francesco Iacovelli, Andrea Cimatti, and Michele Maggiore, "Cosmology and Astrophysics with Standard Sirens and Galaxy Catalogs in View of Future Gravitational Wave Observations," Astrophys. J. 964, 191 (2024), arXiv:2312.05302 [astro-ph.CO].

[35] C. Messenger and J. Read, "Measuring a cosmological distance-redshift relationship using only gravitational wave observations of binary neutron star coalescences," Phys. Rev. Lett. 108, 091101 (2012), arXiv:1107.5725 [gr-qc].

[36] Deep Chatterjee, Abhishek Hegade KR, Gilbert Holder, Daniel E. Holz, Scott Perkins, Kent Yagi, and Nicolás Yunes, "Cosmology with Love: Measuring the Hubble constant using neutron star universal relations," Phys. Rev. D 104, 083528 (2021), arXiv:2106.06589 [gr-qc].

[37] Shang-Jie Jin, Tian-Nuo Li, Jing-Fei Zhang, and Xin Zhang, "Prospects for measuring the Hubble constant and dark energy using gravitational-wave dark sirens with neutron star tidal deformation," JCAP 08, 070 (2023), arXiv:2202.11882 [gr-qc].

[38] Tathagata Ghosh, Bhaskar Biswas, and Sukanta Bose, "Simultaneous inference of neutron star equation of state and the Hubble constant with a population of merging neutron stars," Phys. Rev. D 106, 123529 (2022), arXiv:2203.11756 [astro-ph.CO].

[39] Banafsheh Shiralilou, Geert Raaijmakers, Bastien Duboeuf, Samaya Nissanke, Francois Foucart, Tanja Hinderer, and Andrew R. Williamson, "Measuring the Hubble Constant with Dark Neutron Star-Black Hole Merg-
ers," Astrophys. J. **955**, 149 (2023), arXiv:2207.11792 [astro-ph.CO].

[40] Arnab Dhani, Ssohrab Borhanian, Anuradha Gupta, and Bangalore Sathyaprakash, "Cosmography with bright and Love sirens," (2022), arXiv:2212.13183 [gr-qc].

[41] Tathagata Ghosh, Bhaskar Biswas, Sukanta Bose, and Shasvath J. Kapadia, "Joint Inference of Population, Cosmology, and Neutron Star Equation of State from Gravitational Waves of Dark Binary Neutron Stars," (2024), arXiv:2407.16669 [gr-qc].

[42] Chelsea L. MacLeod and Craig J. Hogan, "Precision of Hubble constant derived using black hole binary absolute distances and statistical redshift information," Phys. Rev. D **77**, 043512 (2008), arXiv:0712.0618 [astro-ph].

[43] Masamune Oguri, "Measuring the distance-redshift relation with the cross-correlation of gravitational wave standard sirens and galaxies," Phys. Rev. D **93**, 083511 (2016), arXiv:1603.02356 [astro-ph.CO].

[44] Sayantani Bera, Divya Rana, Surhud More, and Sukanta Bose, "Incompleteness Matters Not: Inference of $H_0$ from Binary Black Hole–Galaxy Cross-correlations," Astrophys. J. **902**, 79 (2020), arXiv:2007.04271 [astro-ph.CO].

[45] Suvodip Mukherjee, Benjamin D. Wandelt, and Joseph Silk, "Probing the theory of gravity with gravitational lensing of gravitational waves and galaxy surveys," Mon. Not. Roy. Astron. Soc. **494**, 1956–1970 (2020), arXiv:1908.08951 [astro-ph.CO].

[46] Suvodip Mukherjee, Benjamin D. Wandelt, Samaya M. Nissanke, and Alessandra Silvestri, "Accurate precision Cosmology with redshift unknown gravitational wave sources," Phys. Rev. D **103**, 043520 (2021), arXiv:2007.02943 [astro-ph.CO].

[47] Francesca Calore, Alessandro Cuoco, Tania Regimbau, Surabhi Sachdev, and Pasquale Dario Serpico, "Cross-correlating galaxy catalogs and gravitational waves: a tomographic approach," Phys. Rev. Res. **2**, 023314 (2020), arXiv:2002.02466 [astro-ph.CO].

[48] Cristina Cigarran Diaz and Suvodip Mukherjee, "Mapping the cosmic expansion history from LIGO-Virgo-KAGRA in synergy with DESI and SPHEREx," Mon. Not. Roy. Astron. Soc. **511**, 2782–2795 (2022), arXiv:2107.12787 [astro-ph.CO].

[49] Suvodip Mukherjee, Alex Krolewski, Benjamin D. Wandelt, and Joseph Silk, "Cross-correlating dark sirens and galaxies: measurement of $H_0$ from GWTC-3 of LIGO-Virgo-KAGRA," (2022), arXiv:2203.03643 [astro-ph.CO].

[50] Asantha Cooray and Ravi K. Sheth, "Halo Models of Large Scale Structure," Phys. Rept. **372**, 1–129 (2002), arXiv:astro-ph/0206508.

[51] Nick Kaiser, "On the Spatial correlations of Abell clusters," Astrophys. J. Lett. **284**, L9–L12 (1984).

[52] Jeremy L. Tinker, Brant E. Robertson, Andrey V. Kravtsov, Anatoly Klypin, Michael S. Warren, Gustavo Yepes, and Stefan Gottlöber, "The large-scale bias of dark matter halos: Numerical calibration and model tests," The Astrophysical Journal **724**, 878 (2010).

[53] Marcello Cacciato, Frank C. van den Bosch, Surhud More, Ran Li, H. J. Mo, and Xiaohu Yang, "Galaxy Clustering & Galaxy-Galaxy Lensing: A Promising Union to Constrain Cosmological Parameters," Mon. Not. Roy. Astron. Soc. **394**, 929–946 (2009), arXiv:0807.4932 [astro-ph].

[54] Idit Zehavi, Zheng Zheng, David H. Weinberg,

Michael R. Blanton, Neta A. Bahcall, Andreas A. Berlind, Jon Brinkmann, Joshua A. Frieman, James E. Gunn, Robert H. Lupton, Robert C. Nichol, Will J. Percival, Donald P. Schneider, Ramin A. Skibba, Michael A. Strauss, Max Tegmark, and Donald G. York, "Galaxy Clustering in the Completed SDSS Redshift Survey: The Dependence on Color and Luminosity," Astrophys. J. **736**, 59 (2011), arXiv:1005.2413 [astro-ph.CO].

[55] Surhud More, Frank C. van den Bosch, Marcello Cacciato, Ramin Skibba, H. J. Mo, and Xiaohu Yang, "Satellite kinematics - III. Halo masses of central galaxies in SDSS," Mon. Not. Roy. Astron. Soc. **410**, 210–226 (2011), arXiv:1003.3203 [astro-ph.CO].

[56] Surhud More, "How accurate is our knowledge of galaxy bias?" Astrophys. J. **741**, 19 (2011), arXiv:1107.1498 [astro-ph.CO].

[57] Peter Behroozi, Risa H. Wechsler, Andrew P. Hearin, and Charlie Conroy, "UNIVERSEMACHINE: The correlation between galaxy growth and dark matter halo assembly from z = 0–10," Mon. Not. Roy. Astron. Soc. **488**, 3143–3194 (2019), arXiv:1806.07893 [astro-ph.GA].

[58] Ravi K. Sheth and Giuseppe Tormen, "Large scale bias and the peak background split," Mon. Not. Roy. Astron. Soc. **308**, 119 (1999), arXiv:astro-ph/9901122.

[59] Titouan Lazeyras, Marcello Musso, and Fabian Schmidt, "Large-scale assembly bias of dark matter halos," JCAP **03**, 059 (2017), arXiv:1612.04360 [astro-ph.CO].

[60] Vincent Desjacques, Donghui Jeong, and Fabian Schmidt, "Large-scale galaxy bias," Physics Reports **733**, 1–193 (2018).

[61] Ryuichi Takahashi, Masanori Sato, Takahiro Nishimichi, Atsushi Taruya, and Masamune Oguri, "Revising the Halofit Model for the Nonlinear Matter Power Spectrum," Astrophys. J. **761**, 152 (2012), arXiv:1208.2701 [astro-ph.CO].

[62] R. Abbott et al. (LIGO Scientific, Virgo), "Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog," Astrophys. J. Lett. **913**, L7 (2021), arXiv:2010.14533 [astro-ph.HE].

[63] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), "Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3," Phys. Rev. X **13**, 011048 (2023), arXiv:2111.03634 [astro-ph.HE].

[64] Ilya Mandel, Will M. Farr, and Jonathan R. Gair, "Extracting distribution parameters from multiple uncertain observations with selection biases," Mon. Not. Roy. Astron. Soc. **486**, 1086–1093 (2019), arXiv:1809.02063 [physics.data-an].

[65] Anatoly Klypin, Gustavo Yepes, Stefan Gottlober, Francisco Prada, and Steffen Hess, "MultiDark simulations: the story of dark matter halo concentrations and density profiles," Mon. Not. Roy. Astron. Soc. **457**, 4340–4359 (2016), arXiv:1411.4001 [astro-ph.CO].

[66] Colm Talbot and Eric Thrane, "Measuring the binary black hole mass spectrum with an astrophysically motivated parameterization," Astrophys. J. **856**, 173 (2018), arXiv:1801.02699 [astro-ph.HE].

[67] Gregory Ashton, Moritz Hübner, Paul D. Lasky, Colm Talbot, Kendall Ackley, Sylvia Biscoveanu, Qi Chu, Atul Divakarla, Paul J. Easter, Boris Goncharov, and et al., "Bilby: A user-friendly bayesian inference library for gravitational-wave astronomy," The Astrophysical Journal Supplement Series **241**, 27 (2019).
[68] B. P. Abbott et al. (KAGRA, LIGO Scientific, Virgo, VIRGO), “Prospects for observing and localiz- ing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Rev. Rel. 21, 3 (2018), arXiv:1304.0670 [gr-qc].

[69] Sebastian Khan, Sascha Husa, Mark Hannam, Frank Ohme, Michael Pürrer, Xisco Jiménez Forteza, and Alejandro Bohé, “Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomeno- logical model for the advanced detector era,” Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc].

[70] Lionel London, Sebastian Khan, Edward Fauchon-Jones, Cecilio García, Mark Hannam, Sascha Husa, Xisco Jiménez-Forteza, Chinmay Kalaghatgi, Frank Ohme, and Francesco Pannarale, “First higher-multipole model of gravitational waves from spinning and coalescing black-hole binaries,” Phys. Rev. Lett. 120, 161102 (2018), arXiv:1708.00404 [gr-qc].

[71] Salvatore Vitale and Hsin-Yu Chen, “Measuring the Hub- ble constant with neutron star black hole mergers,” Phys. Rev. Lett. 121, 021303 (2018), arXiv:1804.07337 [astro- ph.CO].

[72] Joshua S. Speagle, “Dynesty: a dynamic nested sam- pling package for estimating bayesian posteriors and ev- idences,” Mon. Not. Roy. Astron. Soc. 493, 3132–3158 (2020), arXiv:1904.02180 [astro-ph.IM].

[73] Frank C. van den Bosch, Surhud More, Marcello Cac- ciato, Houjun Mo, and Xiaohu Yang, “Cosmological constraints from a combination of galaxy clustering and lensing - I. Theoretical framework,” Mon. Not. Roy. As- tron. Soc. 430, 725–746 (2013), arXiv:1206.6890 [astro- ph.CO].

[74] Ling Sun et al., “Characterization of systematic error in Advanced LIGO calibration in the second half of O3,” (2021), arXiv:2107.00129 [astro-ph.IM].

[75] D. Bhattacharjee, Y. Lecoeuche, S. Karki, J. Betzwieser, V. Bossilkov, S. Kandhasamy, E. Payne, and R. L. Savage, “Fiducial displacements with improved accuracy for the global network of gravitational wave detectors,”

Class. Quant. Grav. 38, 015009 (2021), arXiv:2006.00130 [astro-ph.IM].

[76] Matthew Evans et al., “A Horizon Study for Cosmic Ex- plorer: Science, Observatories, and Community,” (2021), arXiv:2109.09882 [astro-ph.IM].

[77] Michele Maggiore et al., “Science Case for the Ein- stein Telescope,” JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO].

[78] Ssohrab Borhanian and B. S. Sathyaprakash, “Listen- ing to the Universe with next generation ground-based gravitational-wave detectors,” Phys. Rev. D 110, 083040 (2024), arXiv:2202.11048 [gr-qc].

[79] Mauro Pieroni, Angelo Ricciardone, and Enrico Ba- rausse, “Detectability and parameter estimation of stel- lar origin black hole binaries with next generation grav- itational wave detectors,” Sci. Rep. 12, 17940 (2022), arXiv:2203.12586 [astro-ph.CO].

[80] Francesco Iacovelli, Michele Mancarella, Stefano Foffa, and Michele Maggiore, “Forecasting the Detection Ca- pabilities of Third-generation Gravitational-wave Detec- tors Using GWFAST,” Astrophys. J. 941, 208 (2022), arXiv:2207.02771 [gr-qc].

[81] Marica Branchesi et al., “Science with the Einstein Tele- scope: a comparison of different designs,” JCAP 07, 068 (2023), arXiv:2303.15923 [gr-qc].

[82] Ish Gupta et al., “Characterizing gravitational wave de- tector networks: from A to cosmic explorer,” Class. Quant. Grav. 41, 245001 (2024), arXiv:2307.10421 [gr- qc].

[83] Aditya Vijaykumar, M. V. S. Saketh, Sumit Kumar, Parameswaran Ajith, and Tirthankar Roy Choudhury, “Probing the large scale structure using gravitational- wave observations of binary black holes,” Phys. Rev. D 108, 103017 (2023), arXiv:2005.01111 [astro-ph.CO].

[84] S. Libanore, M. C. Artale, D. Karagiannis, M. Liguori, N. Bartolo, Y. Bouffanais, N. Giacobbo, M. Mapelli, and S. Matarrese, “Gravitational Wave mergers as trac- ers of Large Scale Structures,” JCAP 02, 035 (2021), arXiv:2007.06905 [astro-ph.CO].