# Cross-correlating dark sirens and galaxies: measurement of $H_{0}$ from
GWTC-3 of LIGO-Virgo-KAGRA

Suvodip Mukherjee<EMAIL_ADDRESS>Perimeter Institute for
Theoretical Physics, 31 Caroline St. North, Waterloo, ON NL2 2Y5, Canada  Alex
Krolewski<EMAIL_ADDRESS>AMTD Fellow, Waterloo Centre for
Astrophysics, University of Waterloo, Waterloo ON N2L 3G1, Canada  Perimeter
Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON NL2
2Y5, Canada  Benjamin D. Wandelt<EMAIL_ADDRESS>Institut d’Astrophysique de
Paris, UMR 7095, CNRS, Sorbonne Université, 98bis Boulevard Arago, 75014
Paris, France  Center for Computational Astrophysics, Flatiron Institute, 162
5th Avenue, 10010, New York, NY, USA  Joseph Silk
<EMAIL_ADDRESS><EMAIL_ADDRESS>Institut d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université, 98bis Boulevard Arago, 75014 Paris,
France  Beecroft Institute for Cosmology and Particle Astrophysics, University
of Oxford, Keble Road, Oxford OX1 3RH, UK  The Johns Hopkins University,
Department of Physics & Astronomy, Bloomberg Center for Physics and Astronomy,
Room 366, 3400 N. Charles Street, Baltimore, MD 21218, USA

###### Abstract

We measure the Hubble constant of the Universe using spatial cross-correlation
between gravitational wave (GW) sources without electromagnetic counterparts
from the third GW Transient Catalog (GWTC-3), and the photometric galaxy
surveys 2MPZ and WISE-SuperCOSMOS. Using the eight well-localised GW events,
we obtain Hubble constant $H_{0}=68.2_{-6.2}^{+26.0}$ km/s/Mpc (median and
68.3$\%$ equal-tailed interval (ETI)) after marginalizing over the matter
density and the GW bias parameters. Though the constraints are weak due to a
limited number of GW sources and poor sky localization, they are not subject
to assumptions regarding the GW mass distribution. By combining this
measurement with the Hubble constant measurement from binary neutron star
GW170817, we find a value of Hubble constant $H_{0}=67.0_{-3.8}^{+6.3}$
km/s/Mpc (median and 68.3$\%$ ETI).

Introduction – Discovery of gravitational waves (GW) Abbott et al. (2016a) by
the LIGO-Virgo-KAGRA (LVK) collaboration Aasi et al. (2015); Abbott et al.
(2016b); Acernese et al. (2014, 2019); Akutsu et al. (2019, 2020) has opened a
new observational window to study the cosmos using transient sources such as
binary neutron stars (BNSs), binary black holes (BBHs) and neutron star black
holes (NSBHs). GW sources are uniquely accurate tracers of the luminosity
distance. They can therefore be used to measure the expansion history of the
Universe (Schutz (1986)). This fact has earned GW sources the name standard
sirens. However, one of the key ingredients required to measure the expansion
history using GW sources is an independent measurement (or inference) of the
GW source redshifts. In the absence of electromagnetic counterparts, a
promising way to infer the GW source redshifts is through spatial cross-
correlation of the GW sources with galaxies. By using the clustering redshift
of the GW sources (Oguri, 2016; Bera et al., 2020; Cañas-Herrera et al., 2021)
or a 3D-cross correlation generalizing the Alcock-Paczynski effect Mukherjee
and Wandelt (2018); Mukherjee et al. (2020); Mukherjee et al. (2021a, b);
Cigarrán Díaz and Mukherjee (2022), we can measure the cosmic expansion
history after marginalizing over the GW bias parameters. Apart from cross-
correlation techniques, statistical host identification Schutz (1986); MacLeod
and Hogan (2008); Del Pozzo (2012); Arabsalmani et al. (2013); Chen et al.
(2018); Gray et al. (2020); Fishbach et al. (2019); Abbott et al. (2019);
Soares-Santos et al. (2019); Abbott et al. (2020); Borhanian et al. (2020) and
GW mass distribution Taylor et al. (2012); Farr et al. (2019); Mastrogiovanni
et al. (2021); You et al. (2021); Mancarella et al. (2021); Mukherjee (2021);
Leyde et al. (2022); Ezquiaga and Holz (2022) can also be used to infer
redshifts for BBHs. However, the mass distribution of the BBHs can have
intrinsic redshift dependence that influences parameter estimation, if the
full mass distribution gets affected Mukherjee (2021); Ezquiaga and Holz
(2022).

LVK dark standard sirens have been used to measure the expansion history using
O1+O2 data Abbott et al. (2019) and O1+O2+O3 data Abbott et al. (2021a), in
tandem with GLADE Dálya et al. (2018) and GLADE+ Dálya et al. (2021) for
statistical host identification. The latest LVK measurement yields
$H_{0}=68_{-6}^{+8}$ km/s/Mpc (68.3% highest density interval (HDI)) after
combining with the bright siren GW170817 Abbott et al. (2021a). Other
independent measurements of $H_{0}$ using statistical host identification have
also been performed on the GW data Finke et al. (2021); Palmese et al. (2021).

In this paper, we make the first measurement of the Hubble constant $H_{0}$
using the cross-correlation between the GWTC-3 catalog of the LVK
collaboration Abbott et al. (2021b) and the photometric galaxy surveys 2MPZ
Bilicki et al. (2014) and WISE-SuperCOSMOS (WSC) Bilicki et al. (2016).111We
have not used DeCALS Zhou et al. (2021) in this analysis, as discussed later
Though currently we cannot detect clustering between GW sources and galaxies
due to the limited number of GW sources and poor sky localization error, this
is the first proof of principle of this technique on data. The current
measurement is limited by the lack of high redshift galaxies and the number of
well localised GW sources. The statistical power in the cross-correlation
technique with limited number of sources having poor sky localization error
will be similar to (but not worse than) the statistical host identification
technique. However, the cross-correlation measurement does not depend on
assumptions about the GW source population, and provides an independent
technique to infer the value of the Hubble constant.

Figure 1: The sky map in equatorial coordinates (top) and the luminosity
distance (below) of the eight selected GW sources from GWTC-3.

Method – The compact objects of both astrophysical and primordial origin will
exhibit spatial correlation with galaxies due to the underlying distribution
of dark matter. This spatial correlation can be used to infer the clustering
redshift of the dark standard sirens, by cross-correlating with photometric
(or spectroscopic) galaxy samples, as shown in Mukherjee and Wandelt (2018);
Mukherjee et al. (2021a, b). The clustering redshift technique is also used
widely in other cosmological studies Newman (2008); Menard et al. (2013);
Schmidt et al. (2013). From the observed spatial distribution of GW sources
(or galaxies) $n_{X}(\theta,\phi)$, we can construct a density map as
$\delta^{X}(\theta,\phi)=\frac{n_{X}(\theta,\phi)}{\bar{n}_{X}}-1$ (where
$\bar{n}_{X}$ is the mean density and $X\in\\{\text{GW, g}\\}$). The angular
correlation between a galaxy map and a GW map in the spherical harmonic
basis222$Y_{\ell m}(\theta,\phi)$ denote the spherical harmonic basis.
($\delta^{X}(\theta,\phi)=\sum_{\ell m}\delta^{X}_{\ell m}Y_{\ell
m}(\theta,\phi)$) can be written as
$\hat{C}^{XY}_{\ell}=\frac{\sum_{m}\tilde{\delta}^{X}_{\ell
m}\tilde{\delta}^{*Y}_{\ell m}}{2l+1}$. $\hat{C}^{XY}_{\ell}$ denotes the
pseudo auto (for $X=Y$) and cross (for $X\neq Y$) angular power spectrum
obtained from the masked density maps, denoted by $\tilde{\delta}$. One can
construct the binned average power spectrum as
$\hat{C}^{XY}_{\ell_{b}}=\sum_{\ell\in l_{b}}W_{\ell}C^{XY}_{\ell}$, where
$W_{\ell}$ denotes the normalised window function. The corresponding
covariance matrix for the angular power spectrum in the Gaussian limit can be
written as

$\Sigma_{C^{XY}_{\ell}C^{XY}_{\ell^{\prime}}}=\delta_{\ell\ell^{\prime}}\frac{(C^{XX}_{\ell}+n^{X}_{\ell})(C^{YY}_{\ell}++n^{Y}_{\ell})+C^{XY}_{\ell}}{(2\ell+1)f_{\rm{sky}}\Delta\ell},$
(1)

where $n^{X}_{\ell}$ denotes the shot noise for tracer $X$, equal to the
inverse of its number density. For the cross power spectrum, the shot noise is
zero. $f_{\rm sky}\equiv\Omega_{s}/4\pi$ denotes the overlapping sky fraction
between GW sources and galaxy catalog, and $\Delta\ell$ denotes the bin width
in $\ell$-space over which we estimate the band-averaged power spectrum.
Choosing a large bin width makes it possible to reduce the correlation between
different multipoles.

We estimate the cosmological parameters, Hubble constant $H_{0}$ and matter
density $\Omega_{m}$, along with the GW bias parameter
$b_{GW}(z)=b_{GW}(1+z)^{\alpha}$ using a Bayesian framework.333Mukherjee et
al. (2021a); Cigarrán Díaz and Mukherjee (2022) give a detailed derivation of
the Bayesian framework. THe framework is written for a spectroscopic redshift
survey, so an additional term to include to error on the redshift is included.
The posterior on the parameters given the GW data $\bm{\vartheta}_{GW}$ with
NGW sources and galaxy data $\bm{d}_{g}$ with photometric redshifts $z_{c}$
can be written, after marginalizing over the nuisance parameters
$\Theta_{n}\in\\{b_{GW},\alpha\\}$, as Mukherjee et al. (2021a); Cigarrán Díaz
and Mukherjee (2022)

$\displaystyle\mathcal{P}(\Theta_{c}|\bm{\vartheta}_{GW},\bm{d}_{g},z_{c})$
$\displaystyle\propto\iint
d\Theta_{n}\,dz\,\prod_{i=1}^{N_{GW}}\Pi(z)\Pi(\Theta_{n})\Pi(\Theta_{c})$
$\displaystyle\times\mathcal{P}(\bm{\vartheta}_{GW}|\\{C_{\ell}^{gg}(z_{c})\\},\Theta_{n},\bm{d}_{g}(z_{c}))$
$\displaystyle\times\mathcal{P}(\bm{d}_{g}(z_{c})|\\{C_{\ell}^{gg}(z_{c})\\})\mathcal{P}(z|z_{c})$
$\displaystyle\times\mathcal{P}({\\{d^{i}_{\ell}\\}}_{GW}|z,\Theta_{c},\\{\hat{\theta}^{i},\,\hat{\phi}^{i}\\}_{GW}),$
(2)

where the likelihood
$\mathcal{P}(\bm{\vartheta}_{GW}|\\{C_{\ell}^{gg}(z_{c})\\},\Theta_{n},\bm{d}_{g}(z_{c}))$
is written as

$\displaystyle\mathcal{P}(\bm{\vartheta}_{GW}|$
$\displaystyle\\{C_{\ell}^{gg}(z_{c})\\},\Theta_{n},\bm{d}_{g}(z_{c}))\propto$
(3) $\displaystyle\exp\bigg{(}-0.5\sum^{\ell_{\rm
max}}_{\ell_{b},\ell^{\prime}_{b}}{{D}(\ell_{b},z_{c}){\Sigma^{-1}_{C^{XY}_{\ell_{b}}C^{XY}_{\ell^{\prime}_{b}}}}{D}(\ell^{\prime}_{b},z_{c})}\bigg{)},$

Here, $\ell_{\rm max}$ denotes the maximum value of the multipoles that can be
explored (which depends on the sky localization error) and
${D}(\ell,z_{c})=\hat{C}_{\ell}^{GW\,g}(z_{c})-C_{\ell}^{GW\,g}(z_{c})$. The
angular cross-correlation power spectrum $\hat{C}_{\ell}^{GW\,g}(z_{c})$ is
obtained from cross-correlating GW sources detected above a network matched
filtering SNR with galaxy catalogs $\mathbf{d}_{g}(z_{c})$. The theoretical
angular cross-power spectrum is written in terms of the galaxy auto-power
spectrum $C_{\ell}^{gg}(z_{c})$, the galaxy bias $b_{g}(z_{c}$) and GW bias
$b_{GW}(z_{c})$ as
$C_{\ell}^{GW\,g}(z_{c})=b_{GW}(z_{c})/b_{g}(z_{c})C_{\ell}^{gg}(z_{c})$. The
term $\mathcal{P}(\bm{d}_{g}(z_{c})|\\{C_{\ell}^{gg}(z_{c})\\})$ denotes the
galaxy density field given the auto-power spectrum $C_{\ell}^{gg}(z_{c})$, and
$\mathcal{P}(z|z_{c})$ is the probability distribution of true redshifts
within $z_{c}$, capturing the photometric redshift error on the galaxies. The
likelihood on the luminosity distance given the cosmological parameters
$\Theta_{c}$ and redshift is denoted by
$\mathcal{P}({\\{d^{i}_{\ell}\\}}_{GW}|z,\Theta_{c},\\{\hat{\theta}^{i},\,\hat{\phi}^{i}\\}_{GW})$,
and the prior on the redshift, cosmological parameters, and nuisance
parameters are denoted by $\Pi(z)$, $\Pi(\Theta_{c})$, and $\Pi(\Theta_{n})$
respectively.

GW catalog and selection function : In this analysis, we use the publicly
available GW catalog GWTC-3 detected by the LVK collaboration Abbott et al.
(2021b). As the most constraining estimations of cosmological parameters can
be made from sources with high matched filtering SNR, we select samples from
GWTC-3 with SNR $\geq$ 11\. Also, as the cross-correlation technique is most
effective for sources with better sky localization error, we further select
sources with sky localization error $\Delta\Omega\leq 30$ sq. deg at $68.3\%$
CI. These two selections lead to a total of eight GW events, namely GW170818,
GW190412, GW190814, GW190701$\\_$203306, GW190720$\\_$000836,
GW200129$\\_$0065458, GW200224$\\_$222234, and GW200311$\\_$115853\. The
posteriors on the sky position and the luminosity distance are shown in Fig.
1.

Galaxy catalog and selection function : We use galaxies from the 2MASS
Photometric Redshift catalog (2MPZ, Fig. 2; Bilicki et al., 2014) and WISE
cross SuperCOSMOS Photometric Redshift catalog (WSC, Fig. 2; Bilicki et al.,
2016).444The details of both these catalogs can be found in Bilicki et al.
(2014) and Bilicki et al. (2016). 2MPZ is derived from the all-sky 2MASS near-
infrared extended source catalog (XSC) (Jarrett et al., 2000; Skrutskie et
al., 2006), cross-matched to the infrared AllWISE Wright et al. (2010) and
optical SuperCOSMOS catalogs Hambly et al. (2001a, b, c); Peacock et al.
(2016). WSC is constructed in a similar way, but cross-matching AllWISE and
SuperCOSMOS only. For WSC, we further apply a color cut of
$\rm{W}1-\rm{W}2>0.2$ to the publicly available sample to reduce stellar
contamination and increase uniformity (Xavier et al., 2019). The galaxy
samples are further described and validated in Krolewski et al. (2022–),
including the impact of changing the selection criteria.

Figure 2: The sky map along with mask (gray) in equatorial coordinates of 2MPZ
(top) and WSC (bottom).

The choice of mask555The details for the construction of the mask are
mentioned in the Appendix. for the 2MPZ Bilicki et al. (2014) and WSC Bilicki
et al. (2016) galaxy surveys left nearly $65\%$ of the sky. At higher redshift
(and with lower stellar contamination than WSC), DeCALS Zhou et al. (2021)
covers about $50\%$ of the sky area. Although the sky fraction reduction is
modest for DeCALS, it unfortunately misses four out of the eight GW events. As
a result, we use only 2MPZ and WSC in this paper.

The redshift distributions of 2MPZ and WSC are shown in Fig. 3 in orange and
blue respectively. At $z<0.1$ we use 2MPZ despite its lower number density, as
it has more precise photometric redshifts and far less stellar contamination.
For $z>0.1$, we use WSC exclusively. At $0.3<z<0.4$, WSC clustering is shot
noise dominated at $\ell>60$, and is shot noise dominated at all $\ell$ for
$z>0.4$. We still measure galaxy clustering at $0.4<z<0.5$, albeit with
increased errors (total SNR $\sim 4.7$, after subtracting shot noise, over the
relevant scales $10<\ell<40$), but exclude $z>0.5$ where there are very few
WSC galaxies. For both surveys, photometric redshifts are trained using the
ANNz algorithm Collister and Lahav (2004), yielding typical redshift errors
$\sigma_{z}=0.015$ for 2MPZ and $\sigma_{z}/(1+z)=0.033$ for WSC.

Figure 3: The redshift distribution of the 2MPZ (in orange) and WSC (in blue)
galaxies. The shaded region shows the redshift range that is shot noise
dominated.

Results – We have adopted the following uniform prior ranges:
$\Pi(H_{0})=\mathcal{U}[20,120]$ km/s/Mpc,
$\Pi(\Omega_{m})=\mathcal{U}[0.1,0.4]$, $\Pi(b_{GW})=\mathcal{U}[0.1,6]$, and
$\Pi(\alpha)=\mathcal{U}[-2,2]$. We use a redshift bin-width of $\Delta z=0.1$
which is nearly three times the WSC photo-$z$ error. For the cross-correlation
we include only the lower multipoles $\ell\leq 30$ and consider different
choices of bins $\Delta\ell=5,10,15$ in the analysis. We do not use the first
$\ell$-bin in the analysis to minimize low-$\ell$ contaminations.

From the auto-power spectrum, we infer the galaxy bias $b_{g}(z)$ by fitting a
simple linear bias times the nonlinear “Halofit” matter power spectrum model
(Takahashi et al., 2012). We use the NaMaster code (Hivon et al., 2002; Alonso
et al., 2019) to measure pseudo-$C_{\ell}$ for each redshift slice, applying a
$1^{\circ}$ apodization (“$C^{1}$ apodization”; Grain et al., 2009) to the
galaxy mask. For WSC, we additionally deproject the Schlegel-Finkbeiner-Davis
(Schlegel et al., 1998) dust extinction map and a stellar density map from
Gaia (Brown et al., 2016) to reduce the impact of contamination, following
Koukoufilippas et al. (2020). We fit the one-parameter bias model to the data
in the range $10<\ell<40$, with shot noise fixed. For WSC, we additionally
allow for systematic variations in the number density from variations in the
zeropoint between SuperCOSMOS plates. We add a template to the model
$C_{\ell}^{\rm plate}=A\exp{[-2(\ell\theta_{\rm plate})^{2}/12]}$, where
$\theta_{\rm plate}$ is the plate scale, 5∘ Koukoufilippas et al. (2020).
Finally, we fix the shot noise to the inverse of the angular number density
(in steradians) except for the $0.4<z<0.5$ bin, where we adjust it downwards
by 5% to match the high-$\ell$ power of $C_{\ell}^{gg}$. For the other bins,
we check that $1/\bar{n}_{g}$ matches the high-$\ell$ power in
$C_{\ell}^{gg}$, and the discrepancies are small compared to the clustering
amplitude at $10<\ell<40$.

To model the redshift distribution, we convolve the observed photometric
redshift distribution with a Gaussian for 2MPZ (Balaguera-Antolínez et al.,
2018) and a generalized Lorentzian for WSC, $P(\delta
z)\propto\left(1+\frac{\delta z^{2}}{2as^{2}}\right)^{-a}$ (Peacock and
Bilicki, 2018). The width evolves as a function of redshift, for the Gaussian
following $\sigma=0.027\tanh{(-20.78z_{p}^{2}+7.76z_{p}+0.05)}/(1+z_{p})$,
i.e. increasing from 0.0013 at $z_{p}=0$ to 0.013 at $z_{p}=0.1$; and for the
Lorentzian, $a(z_{c})=-4z_{c}+3$ and $s(z_{c})=0.04z_{c}+0.02$, where $z_{c}$
is the midpoint of each redshift bin. Other choices for the redshift error
(e.g. redshift-independent modified Lorentzian for 2MPZ in Bilicki et al.
(2014) and Peacock and Bilicki (2018)) yield very similar results.

We then use the inferred galaxy auto power spectrum to model the cross-power
spectrum between the GW sources and galaxy as
$C^{gwg}_{\ell}(z)=b_{gw}(z)/b_{g}(z)C^{gg}_{\ell}(z)$. On the GW side, we
construct three GW maps from the selected GW samples composed of Set-1
(GW190814) , Set-2 (GW170818, GW1901412, GW190720_000836, GW2001129_065458,
GW200311_115853), and Set-3 (GW190701_203306, GW200224_222234). These maps are
constructed on the basis of their luminosity distance distribution. Sources
with a similar maximum value of the posterior distribution are combined to
enhance the cross-correlation signal. However, due to fewer sources, the shot
noise for GW sources is very large compared to the galaxies (by nearly 4-6
orders of magnitude depending on the redshift bin). So, the GW auto-
correlation signal is completely dominated by shot noise.

The joint estimation of the Hubble constant along with the matter density and
GW bias parameters are shown in Fig. 4. Constraints on the Hubble constant are
bimodal, with the median value $H_{0}=68.2_{-6.2}^{+26.0}$ km/s/Mpc (the upper
and the lower limit indicates the $68.3\%$ equal-tailed interval (ETI)). At
the lower end of the $H_{0}$ constraints, the constraints arise from the
absence of the angular cross-correlation signal between the GW sources and the
galaxies in the first $z$-bin. Due to the limited number of GW sources, we are
not able to detect the cross-correlation signal with galaxies (see Fig. 6 in
the appendix for further details.). The parameter constraints are driven by
the structure in the cross-correlation signal and its covariance for different
redshifts. Also, there are few galaxies at high redshift and the galaxy
clustering signal is shot noise-dominated there. As a result, the galaxy-GW
source cross-correlation does not add any additional information at high
redshift. This leads to weak constraints at high $H_{0}$ and is one of the
major sources of uncertainty. Particularly for Set-3, the constraints at the
high $H_{0}>90$ km/s/Mpc are affected due unavailability of catalog above
$z=0.5$. Along with this, the uncertainties on the galaxies’ photometric
redshifts and GW luminosity distance, and the limited number of sources, play
a crucial role in the weak constraints. The value of the matter density
$\Omega_{m}$ and bias parameter $b_{GW}(z)$ are not constrained. The estimated
parameters obtained for different values of the bin width are consistent with
each other. For large choices of the $\ell$-bin ($\Delta\ell\geq 10$), we are
able to mitigate the effect from off-diagonal terms. The $\Delta\ell=5$ case
shows a stronger bi-modal distribution. This also implies any systematic error
associated with the covariance estimation is not causing any major systematic
error in the inferred value of the Hubble constant for $\Delta\ell\geq 10$.

By combining the bright standard siren measurement from GW170817 with a better
measurement of peculiar velocity Mukherjee et al. (2021c), we show the
corresponding posterior on $H_{0}$ in Fig. 5 with the median value of
$H_{0}=67.0_{-3.8}^{+6.3}$ km/s/Mpc (68.3$\%$ ETI). This provides tighter
constraints in the value of the Hubble constant than previous measurements. In
Fig. 5 we compare the GW measurements of $H_{0}$ with the measurements from
Planck, $H_{0}=67.4^{+0.5}_{-0.5}$ km/s/Mpc Aghanim et al. (2020) and with the
measurement of $H_{0}=73.04^{+1.05}_{-1.05}$ km/s/Mpc from SH0ES Riess et al.
(2021).666The ACT and WMAP measurement is $H_{0}=67.6^{+1.1}_{-1.1}$ km/s/Mpc
Aiola et al. (2020). The measurement using TRGB as a calibrator is
$H_{0}=69.8\pm 0.6$ (stat) $\pm 1.6$ (sys) Freedman (2021). The strong lensing
measurement from H0LiCOW is $H_{0}=73.3_{-1.8}^{+1.7}$ km/s/Mpc Wong et al.
(2020). The current measurements from the dark sirens are not sufficiently
constraining yet to resolve the tension in the value of $H_{0}$ Verde et al.
(2019); Di Valentino et al. (2021). Though the systematic uncertainties in our
measurement of $H_{0}$ are smaller than the statistical uncertainties, in
future with more GW sources and better galaxy catalog, we will be able to
better assess the influence of any systematic uncertainties.

This measurement of the expansion history does not depend directly on the
choice of GW mass distribution nor on whether the mass distribution follows a
power-law mass distribution or a power-law + Gaussian peak model Abbott et al.
(2021c). It only depends on the maximum allowed mass of individual BHs. This
is primarily because the maximum mass of the BHs determines the maximum
luminosity distance up to which a source can be detected from a given detector
sensitivity. This luminosity distance threshold, in combination with the
allowed priors on the cosmological parameters, sets the maximum redshift in
the prior on GW source redshift out to which one needs to explore the cross-
correlation signal. However, with the currently existing galaxy surveys, we
are not able to go beyond $z=0.5$ due to limited sky coverage and the limited
redshift reach of the galaxy catalogs.

The GW bias parameters encode information about the GW source population and
their connection to galaxies. Hence we marginalize over the GW bias parameters
to infer cosmological parameters correctly and mitigate population
assumptions. At large angular scales ($\ell<30$), the GW bias parameter is
expected to follow the dark matter distribution in the linear regime, and will
not exhibit scale dependence. So, in this analysis, we have not considered any
scale dependence in the GW bias parameter.

Even though this is the first application of the cross-correlation technique
on GW data, the measurements of $H_{0}$ presented in this work are from a low
number of GW sources. As a result, the clustering of the GW sources is not
measured with any statistical significance. So, the estimates of $H_{0}$
presented in this work are weak and subject to systematic uncertainties
associated with small number statistics.777The measurement from other
techniques such as statistical host identification and from the GW mass
distribution are also affected by poor number statistics. To test the
robustness of our results, we have checked the following aspects, (i) randomly
varied the galaxy bias parameter within their error bars, (ii) enhanced the
covariance matrix by a factor of four, (iii) changed the cosmological
parameters such as $S_{8}\equiv\sigma_{8}\sqrt{\Omega_{m}}$ that is used to
fit the galaxy power spectrum $C_{\ell}^{gg}$ from $S_{8}=0.832$ (Planck-2018
Aghanim et al. (2018)) to a lower value $S_{8}=0.75$ as indicated by the KiDS
Collaboration Heymans et al. (2021), (iv) changed the value of $H_{0}=67$
km/s/Mpc to $H_{0}=74$ km/s/Mpc Riess et al. (2019) to estimate the galaxy
bias parameters, (v) changed the galaxy sample selection by additionally
removing WSC sources with a low probability of being galaxies using the SVM
catalog of (Krakowski et al., 2016), requiring $p_{\rm gal}>0.67$, and (vi)
changed the redshift bin width $\Delta z$ to 0.05 which is comparable to the
photo-z errors. The posterior on $H_{0}$ did not show any significant
variation for (i)–(v) cases. For the scenario (vi), the $H_{0}$ posterior show
some variation and $68.3\%$ ETI gets bigger than the estimates presented with
$\Delta z=0.1$. This is because the galaxy redshift kernels begin to overlap
due to photo-$z$ errors, violating our assumption that the GW cross-
correlations in neighboring bins are uncorrelated.

The measurement in this work agrees with the dark siren measurement of
$H_{0}=67_{-12}^{+13}$ km/s/Mpc (68.3% HDI) by the LVK collaboration Abbott et
al. (2021a). Though the current constraints on $H_{0}$ from LVK are driven by
population assumptions, the MAP value of the distribution agrees with this
measurement. This is an independent validation of the LVK population
assumptions. The constraints on the higher value of $H_{0}$ from the LVK
analysis arise from the empty catalog component (which is driven by the
population assumption). So, the upper limit on $H_{0}$ is looser for the
cross-correlation technique than for the LVK analysis. However, at lower
values of $H_{0}$, a tighter constraint is obtained.

Figure 4: The joint constraints on $H_{0}=100h_{0}$ km/s/Mpc, $\Omega_{m}$,
and $b_{GW}(z)=b_{GW}(1+z)^{\alpha}$ for different choices of bin-width
$\Delta l$. Figure 5: Hubble constant $H_{0}$ measurement from GWTC-3 dark
sirens, bright siren GW170817, and combining the both are shown along with the
mean and the standard deviation on the measurements from Planck-2018 Aghanim
et al. (2018) and SH0ES Riess et al. (2021).

Conclusion and future outlook – We present the first measurement of the Hubble
constant $H_{0}$ from dark standard sirens using the cross-correlation
technique. The cross-correlation technique uses the spatial clustering of GW
sources with galaxies and includes information beyond statistical host
identification. With the best eight sources available from GWTC-3, we obtain a
median value of Hubble constant $H_{0}=68.2_{-6.2}^{+26.0}$ km/s/Mpc ($68.3\%$
ETI). Due to the limited number of GW sources and absence of galaxy samples at
high redshift, the cross-correlation signal is not detected, leading to only a
mild improvement from the previous constraints using galaxy catalog Abbott et
al. (2021a). In the future, with the availability of $z<0.8$ spectroscopic
galaxy catalogs such as DESI Aghamousa et al. (2016) and SPHEREx Dore et al.
(2018) (supplemented by $z>0.8$ spectroscopy from Euclid Blanchard et al.
(2020) and photometric redshifts from Vera Rubin Observatory LSST Science
Collaboration et al. (2009)), cross-correlation of the GW sources with
galaxies will be a powerful technique to measure the expansion history
Mukherjee and Wandelt (2018); Mukherjee et al. (2021a); Cigarrán Díaz and
Mukherjee (2022) and testing the general theory of relativity Mukherjee et al.
(2021b); Mukherjee et al. (2022–).

Acknowledgements – Authors are thankful to Gergely Dalya for carefully
reviewing the manuscript and providing useful comments during the LVK internal
review. Authors are also thankful to Maciej Bilicki for providing insightful
suggestions on the paper. Research at Perimeter Institute is supported in part
by the Government of Canada through the Department of Innovation, Science and
Economic Development and by the Province of Ontario through the Ministry of
Colleges and Universities. SM is supported by the Simons Foundation. AK thanks
the AMTD Foundation for support. This analysis is carried out at the Symmetry
computing facility of the Perimeter Institute and the Infinity cluster hosted
by Institut d’Astrophysique de Paris. We thank Stephane Rouberol for smoothly
running the Infinity cluster. We acknowledge the use of following packages in
this work: Astropy Astropy Collaboration et al. (2013, 2018), emcee: The MCMC
Hammer Foreman-Mackey et al. (2013), Giant-Triangle-Confusogram Bocquet and
Carter (2016), healpyGórski et al. (2005); Zonca et al. (2019), IPython Pérez
and Granger (2007), Matplotlib Hunter (2007), NaMaster Alonso et al. (2019),
NumPy van der Walt et al. (2011), and SciPy Jones et al. (2001–). This
research has made use of data obtained from the SuperCOSMOS Science Archive,
prepared and hosted by the Wide Field Astronomy Unit, Institute for Astronomy,
the University of Edinburgh, which is funded by the UK Science and Technology
Facilities Council. The authors would like to thank the LIGO/Virgo/KAGRA
scientific collaboration for providing the data. LIGO is funded by the U.S.
National Science Foundation. Virgo is funded by the French Centre National de
Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica
Nucleare (INFN), and the Dutch Nikhef, with contributions by Polish and
Hungarian institutes. This material is based upon work supported by NSF’s LIGO
Laboratory which is a major facility fully funded by the National Science
Foundation.

## References

  * Abbott et al. (2016a) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016a), https://link.aps.org/doi/10.1103/PhysRevLett.116.061102.
  * Aasi et al. (2015) J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), eprint 1411.4547.
  * Abbott et al. (2016b) B. P. Abbott et al., Phys. Rev. D 93, 112004 (2016b), [Addendum: Phys.Rev.D 97, 059901 (2018)], eprint 1604.00439.
  * Acernese et al. (2014) F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, et al., Classical and Quantum Gravity 32, 024001 (2014), ISSN 1361-6382, http://dx.doi.org/10.1088/0264-9381/32/2/024001.
  * Acernese et al. (2019) F. Acernese, M. Agathos, L. Aiello, A. Allocca, A. Amato, S. Ansoldi, S. Antier, M. Arène, N. Arnaud, S. Ascenzi, et al. (Virgo Collaboration), Phys. Rev. Lett. 123, 231108 (2019), https://link.aps.org/doi/10.1103/PhysRevLett.123.231108.
  * Akutsu et al. (2019) T. Akutsu et al. (KAGRA), Nat. Astron. 3, 35 (2019), eprint 1811.08079.
  * Akutsu et al. (2020) T. Akutsu et al. (KAGRA), arXiv:2005.05574 (2020), eprint 2005.05574.
  * Schutz (1986) B. F. Schutz, Nature 323, 310 (1986).
  * Oguri (2016) M. Oguri, Phys. Rev. D 93, 083511 (2016), https://link.aps.org/doi/10.1103/PhysRevD.93.083511.
  * Bera et al. (2020) S. Bera, D. Rana, S. More, and S. Bose, Astrophys. J. 902, 79 (2020), eprint 2007.04271.
  * Cañas-Herrera et al. (2021) G. Cañas-Herrera, O. Contigiani, and V. Vardanyan, ApJ 918, 20 (2021), eprint 2105.04262.
  * Mukherjee and Wandelt (2018) S. Mukherjee and B. D. Wandelt, arXiv:1808.06615 (2018), eprint 1808.06615.
  * Mukherjee et al. (2020) S. Mukherjee, B. D. Wandelt, and J. Silk, Mon. Not. Roy. Astron. Soc. 494, 1956 (2020), eprint 1908.08951.
  * Mukherjee et al. (2021a) S. Mukherjee, B. D. Wandelt, S. M. Nissanke, and A. Silvestri, Phys. Rev. D 103, 043520 (2021a), eprint 2007.02943.
  * Mukherjee et al. (2021b) S. Mukherjee, B. D. Wandelt, and J. Silk, Mon. Not. Roy. Astron. Soc. 502, 1136 (2021b), eprint 2012.15316.
  * Cigarrán Díaz and Mukherjee (2022) C. Cigarrán Díaz and S. Mukherjee, MNRAS 511, 2782 (2022), eprint 2107.12787.
  * MacLeod and Hogan (2008) C. L. MacLeod and C. J. Hogan, Phys. Rev. D 77, 043512 (2008), eprint 0712.0618.
  * Del Pozzo (2012) W. Del Pozzo, Phys. Rev. D 86, 043011 (2012), eprint 1108.1317.
  * Arabsalmani et al. (2013) M. Arabsalmani, V. Sahni, and T. D. Saini, Phys. Rev. D 87, 083001 (2013), eprint 1301.5779.
  * Chen et al. (2018) H.-Y. Chen, M. Fishbach, and D. E. Holz, Nature 562, 545 (2018), eprint 1712.06531.
  * Gray et al. (2020) R. Gray, I. M. n. Hernandez, H. Qi, A. Sur, P. R. Brady, H.-Y. Chen, W. M. Farr, M. Fishbach, J. R. Gair, A. Ghosh, et al., Phys. Rev. D 101, 122001 (2020), https://link.aps.org/doi/10.1103/PhysRevD.101.122001.
  * Fishbach et al. (2019) M. Fishbach et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 871, L13 (2019), eprint 1807.05667.
  * Abbott et al. (2019) B. Abbott et al. (LIGO Scientific, Virgo), arXiv:1908.06060 (2019), eprint 1908.06060.
  * Soares-Santos et al. (2019) M. Soares-Santos et al. (DES, LIGO Scientific, Virgo), Astrophys. J. Lett. 876, L7 (2019), eprint 1901.01540.
  * Abbott et al. (2020) R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. 896, L44 (2020), eprint 2006.12611.
  * Borhanian et al. (2020) S. Borhanian, A. Dhani, A. Gupta, K. Arun, and B. Sathyaprakash (2020), eprint 2007.02883.
  * Taylor et al. (2012) S. R. Taylor, J. R. Gair, and I. Mandel, Phys. Rev. D 85, 023535 (2012), eprint 1108.5161.
  * Farr et al. (2019) W. M. Farr, M. Fishbach, J. Ye, and D. Holz, Astrophys. J. Lett. 883, L42 (2019), eprint 1908.09084.
  * Mastrogiovanni et al. (2021) S. Mastrogiovanni, K. Leyde, C. Karathanasis, E. Chassande-Mottin, D. A. Steer, J. Gair, A. Ghosh, R. Gray, S. Mukherjee, and S. Rinaldi, arXiv:2103.14663 (2021), eprint 2103.14663.
  * You et al. (2021) Z.-Q. You, X.-J. Zhu, G. Ashton, E. Thrane, and Z.-H. Zhu, Astrophys. J. 908, 215 (2021), eprint 2004.00036.
  * Mancarella et al. (2021) M. Mancarella, E. Genoud-Prachex, and M. Maggiore (2021), eprint 2112.05728.
  * Mukherjee (2021) S. Mukherjee (2021), eprint 2112.10256.
  * Leyde et al. (2022) K. Leyde, S. Mastrogiovanni, D. A. Steer, E. Chassande-Mottin, and C. Karathanasis (2022), eprint 2202.00025.
  * Ezquiaga and Holz (2022) J. M. Ezquiaga and D. E. Holz (2022), eprint 2202.08240.
  * Abbott et al. (2021a) R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA) (2021a), eprint 2111.03604.
  * Dálya et al. (2018) G. Dálya, G. Galgóczi, L. Dobos, Z. Frei, I. S. Heng, R. Macas, C. Messenger, P. Raffai, and R. S. de Souza, Mon. Not. Roy. Astron. Soc. 479, 2374 (2018), eprint 1804.05709.
  * Dálya et al. (2021) G. Dálya et al. (2021), eprint 2110.06184.
  * Finke et al. (2021) A. Finke, S. Foffa, F. Iacovelli, M. Maggiore, and M. Mancarella, arXiv:2101.12660 (2021), eprint 2101.12660.
  * Palmese et al. (2021) A. Palmese, C. R. Bom, S. Mucesh, and W. G. Hartley (2021), eprint 2111.06445.
  * Abbott et al. (2021b) R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA) (2021b), eprint 2111.03606.
  * Bilicki et al. (2014) M. Bilicki, T. H. Jarrett, J. A. Peacock, M. E. Cluver, and L. Steward, Astrophys. J. Suppl. 210, 9 (2014), eprint 1311.5246.
  * Bilicki et al. (2016) M. Bilicki et al., Astrophys. J. Suppl. 225, 5 (2016), eprint 1607.01182.
  * Zhou et al. (2021) R. Zhou et al., Mon. Not. Roy. Astron. Soc. 501, 3309 (2021), eprint 2001.06018.
  * Newman (2008) J. A. Newman, Astrophys. J. 684, 88 (2008), eprint 0805.1409.
  * Menard et al. (2013) B. Menard, R. Scranton, S. Schmidt, C. Morrison, D. Jeong, T. Budavari, and M. Rahman, arXiv:1303.4722 (2013), eprint 1303.4722.
  * Schmidt et al. (2013) S. Schmidt, B. Menard, R. Scranton, C. Morrison, and C. McBride, Mon. Not. Roy. Astron. Soc. 431, 3307 (2013), eprint 1303.0292.
  * Jarrett et al. (2000) T. H. Jarrett, T. Chester, R. Cutri, S. Schneider, M. Skrutskie, and J. P. Huchra, AJ 119, 2498 (2000), eprint astro-ph/0004318.
  * Skrutskie et al. (2006) M. F. Skrutskie, R. M. Cutri, R. Stiening, M. D. Weinberg, S. Schneider, J. M. Carpenter, C. Beichman, R. Capps, T. Chester, J. Elias, et al., AJ 131, 1163 (2006).
  * Wright et al. (2010) E. L. Wright, P. R. M. Eisenhardt, A. K. Mainzer, M. E. Ressler, R. M. Cutri, T. Jarrett, J. D. Kirkpatrick, D. Padgett, R. S. McMillan, M. Skrutskie, et al., AJ 140, 1868 (2010), eprint 1008.0031.
  * Hambly et al. (2001a) N. C. Hambly, H. T. MacGillivray, M. A. Read, S. B. Tritton, E. B. Thomson, B. D. Kelly, D. H. Morgan, R. E. Smith, S. P. Driver, J. Williamson, et al., MNRAS 326, 1279 (2001a), eprint astro-ph/0108286.
  * Hambly et al. (2001b) N. C. Hambly, M. J. Irwin, and H. T. MacGillivray, MNRAS 326, 1295 (2001b), eprint astro-ph/0108290.
  * Hambly et al. (2001c) N. C. Hambly, A. C. Davenhall, M. J. Irwin, and H. T. MacGillivray, MNRAS 326, 1315 (2001c), eprint astro-ph/0108291.
  * Peacock et al. (2016) J. A. Peacock, N. C. Hambly, M. Bilicki, H. T. MacGillivray, L. Miller, M. A. Read, and S. B. Tritton, MNRAS 462, 2085 (2016), eprint 1607.01189.
  * Xavier et al. (2019) H. S. Xavier, M. V. Costa-Duarte, A. Balaguera-Antolínez, and M. Bilicki, JCAP 2019, 037 (2019), eprint 1812.08182.
  * Krolewski et al. (2022–) A. Krolewski et al., _Under preparation_ (2022–).
  * Collister and Lahav (2004) A. A. Collister and O. Lahav, PASP 116, 345 (2004), eprint astro-ph/0311058.
  * Takahashi et al. (2012) R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, ApJ 761, 152 (2012), eprint 1208.2701.
  * Hivon et al. (2002) E. Hivon, K. M. Górski, C. B. Netterfield, B. P. Crill, S. Prunet, and F. Hansen, ApJ 567, 2 (2002), eprint astro-ph/0105302.
  * Alonso et al. (2019) D. Alonso, J. Sanchez, and A. Slosar (LSST Dark Energy Science), Mon. Not. Roy. Astron. Soc. 484, 4127 (2019), eprint 1809.09603.
  * Grain et al. (2009) J. Grain, M. Tristram, and R. Stompor, Phys.Rev.D 79, 123515 (2009), eprint 0903.2350.
  * Schlegel et al. (1998) D. J. Schlegel, D. P. Finkbeiner, and M. Davis, ApJ 500, 525 (1998), eprint astro-ph/9710327.
  * Brown et al. (2016) A. G. A. Brown, A. Vallenari, T. Prusti, J. H. de Bruijne, F. Mignard, R. Drimmel, C. Babusiaux, C. A. Bailer-Jones, U. Bastian, M. Biermann, et al., Astronomy & Astrophysics 595, A2 (2016), ISSN 1432-0746, http://dx.doi.org/10.1051/0004-6361/201629512.
  * Koukoufilippas et al. (2020) N. Koukoufilippas, D. Alonso, M. Bilicki, and J. A. Peacock, MNRAS 491, 5464 (2020), eprint 1909.09102.
  * Balaguera-Antolínez et al. (2018) A. Balaguera-Antolínez, M. Bilicki, E. Branchini, and A. Postiglione, MNRAS 476, 1050 (2018), eprint 1711.04583.
  * Peacock and Bilicki (2018) J. A. Peacock and M. Bilicki, MNRAS 481, 1133 (2018), eprint 1805.11525.
  * Mukherjee et al. (2021c) S. Mukherjee, G. Lavaux, F. R. Bouchet, J. Jasche, B. D. Wandelt, S. M. Nissanke, F. Leclercq, and K. Hotokezaka, Astron. Astrophys. 646, A65 (2021c), eprint 1909.08627.
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
  * Riess et al. (2021) A. G. Riess et al. (2021), eprint 2112.04510.
  * Aiola et al. (2020) S. Aiola et al. (ACT), JCAP 12, 047 (2020), eprint 2007.07288.
  * Freedman (2021) W. L. Freedman, Astrophys. J. 919, 16 (2021), eprint 2106.15656.
  * Wong et al. (2020) K. C. Wong et al., Mon. Not. Roy. Astron. Soc. 498, 1420 (2020), eprint 1907.04869.
  * Verde et al. (2019) L. Verde, T. Treu, and A. Riess, in _Tensions between the Early and the Late Universe_ (2019), eprint 1907.10625.
  * Di Valentino et al. (2021) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021), eprint 2103.01183.
  * Abbott et al. (2021c) R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA) (2021c), eprint 2111.03634.
  * Aghanim et al. (2018) N. Aghanim et al. (Planck), arXiv: 1807.06209 (2018), eprint 1807.06209.
  * Heymans et al. (2021) C. Heymans, T. Tröster, M. Asgari, C. Blake, H. Hildebrandt, B. Joachimi, K. Kuijken, C.-A. Lin, A. G. Sánchez, J. L. van den Busch, et al., A&A 646, A140 (2021), eprint 2007.15632.
  * Riess et al. (2019) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Astrophys. J. 876, 85 (2019), eprint 1903.07603.
  * Krakowski et al. (2016) T. Krakowski, K. Małek, M. Bilicki, A. Pollo, A. Kurcz, and M. Krupa, A&A 596, A39 (2016), eprint 1607.01188.
  * Aghamousa et al. (2016) A. Aghamousa et al. (DESI) (2016), eprint 1611.00036.
  * Dore et al. (2018) O. Dore et al., arXiv (2018), eprint 1805.05489.
  * Blanchard et al. (2020) A. Blanchard et al. (Euclid), Astron. Astrophys. 642, A191 (2020), eprint 1910.09273.
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, P. A. Abell, J. Allison, S. F. Anderson, J. R. Andrew, J. R. P. Angel, L. Armus, D. Arnett, S. J. Asztalos, T. S. Axelrod, et al., ArXiv e-prints (2009), eprint 0912.0201.
  * Mukherjee et al. (2022–) S. Mukherjee, A. Krolewski, et al., _Under preparation_ (2022–).
  * Astropy Collaboration et al. (2013) Astropy Collaboration, T. P. Robitaille, E. J. Tollerud, P. Greenfield, M. Droettboom, E. Bray, T. Aldcroft, M. Davis, A. Ginsburg, A. M. Price-Whelan, et al., A&A 558, A33 (2013), eprint 1307.6212.
  * Astropy Collaboration et al. (2018) Astropy Collaboration, A. M. Price-Whelan, B. M. Sipőcz, H. M. Günther, P. L. Lim, S. M. Crawford, S. Conseil, D. L. Shupe, M. W. Craig, N. Dencheva, et al., AJ 156, 123 (2018), eprint 1801.02634.
  * Foreman-Mackey et al. (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, PASP 125, 306 (2013), eprint 1202.3665.
  * Bocquet and Carter (2016) S. Bocquet and F. W. Carter, The Journal of Open Source Software 1 (2016), http://dx.doi.org/10.21105/joss.00046.
  * Górski et al. (2005) K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelman, Astrophys. J. 622, 759 (2005), eprint astro-ph/0409513.
  * Zonca et al. (2019) A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon, and K. Gorski, Journal of Open Source Software 4, 1298 (2019), https://doi.org/10.21105/joss.01298.
  * Pérez and Granger (2007) F. Pérez and B. E. Granger, Computing in Science and Engineering 9, 21 (2007), ISSN 1521-9615, http://ipython.org.
  * Hunter (2007) J. D. Hunter, Computing In Science & Engineering 9, 90 (2007).
  * van der Walt et al. (2011) S. van der Walt, S. C. Colbert, and G. Varoquaux, Computing in Science and Engineering 13, 22 (2011), eprint 1102.1523.
  * Jones et al. (2001–) E. Jones, T. Oliphant, P. Peterson, et al., _SciPy: Open source scientific tools for Python_ (2001–), [Online; accessed ¡today¿], http://www.scipy.org/.
  * Alonso et al. (2015) D. Alonso, A. I. Salvador, F. J. Sánchez, M. Bilicki, J. García-Bellido, and E. Sánchez, MNRAS 449, 670 (2015), eprint 1412.5151.
  * Rafiei-Ravandi et al. (2021) M. Rafiei-Ravandi, K. M. Smith, D. Li, K. W. Masui, A. Josephy, M. Dobbs, D. Lang, M. Bhardwaj, C. Patel, K. Bandura, et al., ApJ 922, 42 (2021), eprint 2106.04354.
  * Krolewski et al. (2020) A. Krolewski, S. Ferraro, E. F. Schlafly, and M. White, JCAP 2020, 047 (2020), eprint 1909.07412.
  * Novaes et al. (2018) C. P. Novaes, A. Bernui, H. S. Xavier, and G. A. Marques, MNRAS 478, 3253 (2018), eprint 1805.04078.

##  Appendix A Cross-correlation signal between GW sources and galaxies

The cross-correlation between the GW sources and galaxies are dominated by
poor number statistics and we do not measure any cross-correlation clustering
signal with statistical significance. We show the measured band average cross-
correlation signal $\hat{C}_{\ell}^{GWg}$ as a function of redshift in Fig. 6
for $\Delta\ell=15$ with the three GW maps constructed from the eight events.
The value of maximum multipole $\ell_{\rm max}=30$ is chosen to take into
account the poor sky localization error of the GW sources for sources with sky
localization error $\Delta\Omega\leq 30$ sq. deg. The diagonal error bars are
shown on the measured values. The band average signal shows a signal that is
consistent with zero at all the redshift bins. The constraints on the
clustering signal at low redshift $z$ for all the three cases drives the
constraints on the lower value of $H_{0}\leq 50$ km/s/Mpc. For the higher
redshift bins, a non-zero mean value causes a non-zero posterior distribution
on $H_{0}$ which is different from the prior (taken as $\mathcal{U}[20,120]$
km/s/Mpc. This indicates that the constraints are weak and are affected by the
systematic associated with the non-detection of the cross-correlation signal.
At high redshift ($z>0.5$) the absence of sufficient galaxies leads to an
abrupt cut in the galaxy distribution and provides no information on higher
values of $H_{0}$. Due to the absence of detection of the clustering signal,
the statistical power of measuring $H_{0}$ is similar to the statistical host
identification technique.

Figure 6: The band average cross-correlation power spectrum between GW sources
and galaxies with a $\Delta\ell=15$ bin-width is shown for three different
maps of the GW sources Set-1 (GW190814) , Set-2 (GW170818, GW1901412,
GW190720_000836, GW2001129_065458, GW200311_115853), and Set-3
(GW190701_203306, GW200224_222234) composed from the selected eight events, as
a function of the median value of the redshift bin.

##  Appendix B Construction of galaxy mask for 2MPZ and WSC

The galaxy masks were carefully constructed to remove areas with large numbers
of stars or other systematics that could affect galaxy clustering, either by
direct stellar contamination or by correlations, e.g. suppressed galaxy
density in regions of high stellar density or extinction. We follow Balaguera-
Antolínez et al. (2018) to construct the 2MPZ mask, starting by masking low
Galactic latitudes ($|b|<10^{\circ}$, areas of high galactic extinction
($E(B-V)>0.3$ from Schlegel et al. (1998), and areas of high stellar density
as estimated from the 2MASS Point Source Catalog ($\log{n_{\rm
star}}>3.5$)888https://www.ipac.caltech.edu/2mass/releases/allsky/
doc/sec4_5c.html. We further include manual cutouts around the LMC and SMC,
excluding $275.47<{\rm RA}<285.47$ and $-37.89<{\rm DEC}<-27.89$, and
$300.81<{\rm RA}<304.81$ and $-46.33<{\rm DEC}<-42.33$. Finally, we mask
additional areas with low completeness in 2MPZ, determined by comparing the
number counts of 2MPZ sources and 2MASS XSC sources (with $K_{s}<13.9$) in
NSIDE=64 HEALPixels. We remove pixels with $<85\%$ completeness, mostly
corresponding to areas of lower depth around bright stars. We test variations
in the masking procedure (i.e. additionally multiplying by the WSC mask,
following Alonso et al. (2015), or changing the completeness threshold to 80%
or 90%) and find minimal changes in results.

For WSC, we follow the masking procedure of Xavier et al. (2019). We start
with the mask distributed with the WSC data release Bilicki et al.
(2016)999http://ssa.roe.ac.uk/WISExSCOSmask.fits.gz. We additionally mask
regions with high extinction ($E(B-V)>0.10$) and high stellar density (density
of stars from GAIA greater than 7 times the mean). We additionally test
several variations in the masking procedure, adding an additional mask at low
Galactic latitudes following Rafiei-Ravandi et al. (2021); adding a WISE
bright stars mask Krolewski et al. (2020); and adding a mask of regions in
WISE with high moon contamination, as determined by HEALPix pixels in which
GLADE+ Dálya et al. (2021) is incomplete compared to WSC. We also test
variations in the sample-selection procedure, i.e. additionally using the SVM
catalog of Krakowski et al. (2016) to restrict to likely galaxies Rafiei-
Ravandi et al. (2021); Novaes et al. (2018). We find that these variations
generally lead to a scale-independent shift in the amplitude of
$C_{\ell}^{gg}$, either corresponding to a change in galaxy bias due to
differing populations, or a change in the stellar contamination fraction,
which is entirely degenerate with bias at $\ell>10$ where the stellar power
spectrum is small. In this regime, the effect of changing stellar
contamination is degenerate with bias in both the galaxy auto spectrum and the
galaxy cross-spectrum with GW sources, so it will cause systematic errors in
our modeling.

We assume that the galaxy bias is redshift-independent in each bin, and obtain
best-fit values of $b_{g}(z_{c}=0.05)=1.18$, $b_{g}(z_{c}=0.15)=0.66$,
$b_{g}(z_{c}=0.25)=1.35$, $b_{g}(z_{c}=0.35)=1.76$, and
$b_{g}(z_{c}=0.45)=2.33$. The very low value of $b_{g}$ in the second bin is
driven by the SuperCOSMOS plate template, which is degenerate with the
cosmological contribution due to the limited multipole range considered
($10<\ell<40$).