Cross-correlating dark sirens and galaxies: measurement of H0 from GWTC-3 of

LIGO-Virgo-KAGRA

Suvodip Mukherjee,1, ∗ Alex Krolewski,2, 1, † Benjamin D. Wandelt,3, 4, ‡ and Joseph Silk3, 5, 6, §

1
Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON NL2 2Y5, Canada
2
AMTD Fellow, Waterloo Centre for Astrophysics,
University of Waterloo, Waterloo ON N2L 3G1, Canada
3
Institut d’Astrophysique de Paris, UMR 7095, CNRS,
Sorbonne Université, 98bis Boulevard Arago, 75014 Paris, France
4
Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA
5
Beecroft Institute for Cosmology and Particle Astrophysics,
University of Oxford, Keble Road, Oxford OX1 3RH, UK
6
The Johns Hopkins University, Department of Physics & Astronomy,
arXiv:2203.03643v1 [astro-ph.CO] 7 Mar 2022

Bloomberg Center for Physics and Astronomy, Room 366,

3400 N. Charles Street, Baltimore, MD 21218, USA
(Dated: March 9, 2022)
We measure the Hubble constant of the Universe using spatial cross-correlation between grav-
itational 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 H0 = 68.2+26.0
−6.2 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 H0 = 67.0+6.3
−3.8 km/s/Mpc (median and 68.3% ETI).

Introduction – Discovery of gravitational waves (GW) highest density interval (HDI)) after combining with the
[1] by the LIGO-Virgo-KAGRA (LVK) collaboration [2– bright siren GW170817 [35]. Other independent mea-
7] has opened a new observational window to study the surements of H0 using statistical host identification have
cosmos using transient sources such as binary neutron also been performed on the GW data [38, 39].
stars (BNSs), binary black holes (BBHs) and neutron In this paper, we make the first measurement of the
star black holes (NSBHs). GW sources are uniquely Hubble constant H0 using the cross-correlation between
accurate tracers of the luminosity distance. They can the GWTC-3 catalog of the LVK collaboration [40] and
therefore be used to measure the expansion history of the photometric galaxy surveys 2MPZ [41] and WISE-
the Universe (Schutz [8]). This fact has earned GW SuperCOSMOS (WSC) [42].1 Though currently we can-
sources the name standard sirens. However, one of the not detect clustering between GW sources and galaxies
key ingredients required to measure the expansion his- due to the limited number of GW sources and poor sky
tory using GW sources is an independent measurement localization error, this is the first proof of principle of
(or inference) of the GW source redshifts. In the ab- this technique on data. The current measurement is lim-
sence of electromagnetic counterparts, a promising way ited by the lack of high redshift galaxies and the number
to infer the GW source redshifts is through spatial cross- of well localised GW sources. The statistical power in
correlation of the GW sources with galaxies. By using the cross-correlation technique with limited number of
the clustering redshift of the GW sources [9–11] or a 3D- sources having poor sky localization error will be similar
cross correlation generalizing the Alcock-Paczynski effect to (but not worse than) the statistical host identification
[12–16], we can measure the cosmic expansion history af- technique. However, the cross-correlation measurement
ter marginalizing over the GW bias parameters. Apart does not depend on assumptions about the GW source
from cross-correlation techniques, statistical host identi- population, and provides an independent technique to
fication [8, 17–26] and GW mass distribution [27–34] can infer the value of the Hubble constant.
also be used to infer redshifts for BBHs. However, the Method – The compact objects of both astrophysical
mass distribution of the BBHs can have intrinsic redshift and primordial origin will exhibit spatial correlation with
dependence that influences parameter estimation, if the galaxies due to the underlying distribution of dark mat-
full mass distribution gets affected [32, 34]. ter. This spatial correlation can be used to infer the
LVK dark standard sirens have been used to mea- clustering redshift of the dark standard sirens, by cross-
sure the expansion history using O1+O2 data [23] and correlating with photometric (or spectroscopic) galaxy
O1+O2+O3 data [35], in tandem with GLADE [36] and
GLADE+ [37] for statistical host identification. The latest
1
LVK measurement yields H0 = 68+8 −6 km/s/Mpc (68.3%
We have not used DeCALS [43] in this analysis, as discussed later
 2

We estimate the cosmological parameters, Hubble con-

stant H0 and matter density Ωm , along with the GW
bias parameter bGW (z) = bGW (1 + z)α using a Bayesian
framework.3 The posterior on the parameters given
the GW data ϑGW with NGW sources and galaxy data
dg with photometric redshifts zc can be written, af-
ter marginalizing over the nuisance parameters Θn ∈
{bGW , α}, as [14, 16]
ZZ NY
GW

P(Θc |ϑGW , dg , zc ) ∝ dΘn dz Π(z)Π(Θn )Π(Θc )

i=1
× P(ϑGW |{C`gg (zc )}, Θn , dg (zc ))
× P(dg (zc )|{C`gg (zc )})P(z|zc )
× P({di` }GW |z, Θc , {θ̂i , φ̂i }GW ),
(2)

where the likelihood P(ϑGW |{C`gg (zc )}, Θn , dg (zc )) is

written as

FIG. 1: The sky map in equatorial coordinates (top) and P(ϑGW |{C`gg (zc )}, Θn , dg (zc )) ∝ (3)
the luminosity distance (below) of the eight selected GW  `X
max 
sources from GWTC-3. exp − 0.5 D(`b , zc )Σ−1 D(`0b , zc ) ,
C XY C XY `b `0
`b ,`0b b

samples, as shown in [12, 14, 15]. The clustering redshift

Here, `max denotes the maximum value of the multipoles
technique is also used widely in other cosmological stud-
that can be explored (which depends on the sky local-
ies [44–46]. From the observed spatial distribution of GW
ization error) and D(`, zc ) = Ĉ`GW g (zc ) − C`GW g (zc ).
sources (or galaxies) nX (θ, φ), we can construct a density
The angular cross-correlation power spectrum Ĉ`GW g (zc )
map as δ X (θ, φ) = nXn̄(θ,φ) − 1 (where n̄X is the mean
X is obtained from cross-correlating GW sources detected
density and X ∈ {GW, g}). The angular correlation
above a network matched filtering SNR with galaxy
between a galaxy map and a P GW map in the spherical
catalogs dg (zc ). The theoretical angular cross-power
harmonic basis2 (δ X (θ, φ) = `m δ`m X
Y`m (θ, φ)) can be
P
δ̃ X δ̃ ∗Y
spectrum is written in terms of the galaxy auto-power
written as Ĉ`XY = m2l+1`m `m
. Ĉ`XY denotes the pseudo spectrum C`gg (zc ), the galaxy bias bg (zc ) and GW bias
auto (for X = Y ) and cross (for X 6= Y ) angular power bGW (zc ) as C`GW g (zc ) = bGW (zc )/bg (zc )C`gg (zc ). The
spectrum obtained from the masked density maps, de- term P(dg (zc )|{C`gg (zc )}) denotes the galaxy density
noted by δ̃. One can construct the binned average power field given the auto-power spectrum C`gg (zc ), and P(z|zc )
spectrum as Ĉ`XY XY
P
b
= `∈lb W ` C` , where W` denotes is the probability distribution of true redshifts within zc ,
the normalised window function. The corresponding co- capturing the photometric redshift error on the galax-
variance matrix for the angular power spectrum in the ies. The likelihood on the luminosity distance given the
Gaussian limit can be written as cosmological parameters Θc and redshift is denoted by
(C`XX + nX` )(C`
YY
+ +nY` ) + C`XY P({di` }GW |z, Θc , {θ̂i , φ̂i }GW ), and the prior on the red-
ΣC`XY C XY = δ``0 , shift, cosmological parameters, and nuisance parameters
` 0
(2` + 1)fsky ∆`
(1) are denoted by Π(z), Π(Θc ), and Π(Θn ) respectively.
where nX` denotes the shot noise for tracer X, equal to GW catalog and selection function : In this analysis,
the inverse of its number density. For the cross power we use the publicly available GW catalog GWTC-3
spectrum, the shot noise is zero. fsky ≡ Ωs /4π denotes detected by the LVK collaboration [40]. As the most
the overlapping sky fraction between GW sources and constraining estimations of cosmological parameters can
galaxy catalog, and ∆` denotes the bin width in `-space be made from sources with high matched filtering SNR,
over which we estimate the band-averaged power spec- we select samples from GWTC-3 with SNR ≥ 11. Also,
trum. Choosing a large bin width makes it possible to
reduce the correlation between different multipoles.
3 [14, 16] give a detailed derivation of the Bayesian framework.
THe framework is written for a spectroscopic redshift survey, so
2 Y`m (θ, φ) denote the spherical harmonic basis. an additional term to include to error on the redshift is included.
 3

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

the sky fraction reduction is modest for DeCALS, it un-

fortunately 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
FIG. 2: The sky map along with mask (gray) in equatorial shown in Fig. 3 in orange and blue respectively. At
coordinates of 2MPZ (top) and WSC (bottom). z < 0.1 we use 2MPZ despite its lower number density, as
it has more precise photometric redshifts and far less stel-
lar contamination. For z > 0.1, we use WSC exclusively.
as the cross-correlation technique is most effective for At 0.3 < z < 0.4, WSC clustering is shot noise dominated
sources with better sky localization error, we further at ` > 60, and is shot noise dominated at all ` for z > 0.4.
select sources with sky localization error ∆Ω ≤ 30 sq. We still measure galaxy clustering at 0.4 < z < 0.5, albeit
deg at 68.3% CI. These two selections lead to a total with increased errors (total SNR ∼ 4.7, after subtracting
of eight GW events, namely GW170818, GW190412, shot noise, over the relevant scales 10 < ` < 40), but
GW190814, GW190701 203306, GW190720 000836, exclude z > 0.5 where there are very few WSC galaxies.
GW200129 0065458, GW200224 222234, and For both surveys, photometric redshifts are trained using
GW200311 115853. The posteriors on the sky po- the ANNz algorithm [56], yielding typical redshift errors
sition and the luminosity distance are shown in Fig. 1. σz = 0.015 for 2MPZ and σz /(1 + z) = 0.033 for WSC.
Galaxy catalog and selection function : We use galaxies Results – We have adopted the following uniform
from the 2MASS Photometric Redshift catalog [2MPZ, prior ranges: Π(H0 ) = U[20, 120] km/s/Mpc, Π(Ωm ) =
Fig. 2; 41] and WISE cross SuperCOSMOS Photometric U[0.1, 0.4], Π(bGW ) = U[0.1, 6], and Π(α) = U[−2, 2]. We
Redshift catalog [WSC, Fig. 2; 42].4 2MPZ is derived use a redshift bin-width of ∆z = 0.1 which is nearly three
from the all-sky 2MASS near-infrared extended source times the WSC photo-z error. For the cross-correlation
catalog (XSC) [47, 48], cross-matched to the infrared All- we include only the lower multipoles ` ≤ 30 and consider
WISE [49] and optical SuperCOSMOS catalogs [50–53]. different choices of bins ∆` = 5, 10, 15 in the analysis.
WSC is constructed in a similar way, but cross-matching We do not use the first `-bin in the analysis to minimize
AllWISE and SuperCOSMOS only. For WSC, we fur- low-` contaminations.
ther apply a color cut of W1 − W2 > 0.2 to the publicly
From the auto-power spectrum, we infer the galaxy
available sample to reduce stellar contamination and in-
bias bg (z) by fitting a simple linear bias times the non-
crease uniformity [54]. The galaxy samples are further
linear “Halofit” matter power spectrum model [57]. We
described and validated in [55], including the impact of
use the NaMaster code [58, 59] to measure pseudo-C`
changing the selection criteria.
for each redshift slice, applying a 1◦ apodization [“C 1
The choice of mask5 for the 2MPZ [41] and WSC [42]
apodization”; 60] to the galaxy mask. For WSC, we
galaxy surveys left nearly 65% of the sky. At higher red-
additionally deproject the Schlegel-Finkbeiner-Davis [61]
shift (and with lower stellar contamination than WSC),
dust extinction map and a stellar density map from Gaia
DeCALS [43] covers about 50% of the sky area. Although
[62] to reduce the impact of contamination, following [63].
We fit the one-parameter bias model to the data in the
range 10 < ` < 40, with shot noise fixed. For WSC, we
4 The details of both these catalogs can be found in [41] and [42]. additionally allow for systematic variations in the num-
5 The details for the construction of the mask are mentioned in ber density from variations in the zeropoint between Su-
the Appendix. perCOSMOS plates. We add a template to the model
 4

C`plate = A exp [−2(`θplate )2 /12], where θplate is the plate uncertainty. Particularly for Set-3, the constraints at the
scale, 5◦ [63]. Finally, we fix the shot noise to the inverse high H0 > 90 km/s/Mpc are affected due unavailability
of the angular number density (in steradians) except for of catalog above z = 0.5. Along with this, the uncer-
the 0.4 < z < 0.5 bin, where we adjust it downwards tainties on the galaxies’ photometric redshifts and GW
by 5% to match the high-` power of C`gg . For the other luminosity distance, and the limited number of sources,
bins, we check that 1/n̄g matches the high-` power in play a crucial role in the weak constraints. The value of
C`gg , and the discrepancies are small compared to the the matter density Ωm and bias parameter bGW (z) are
clustering amplitude at 10 < ` < 40. not constrained. The estimated parameters obtained for
To model the redshift distribution, we convolve the different values of the bin width are consistent with each
observed photometric redshift distribution with a Gaus- other. For large choices of the `-bin (∆` ≥ 10), we are
sian for 2MPZ [64] and a generalized Lorentzian for able to mitigate the effect from off-diagonal terms. The
 −a ∆` = 5 case shows a stronger bi-modal distribution. This
δz 2
WSC, P (δz) ∝ 1 + 2as 2 [65]. The width evolves also implies any systematic error associated with the co-
as a function of redshift, for the Gaussian following variance estimation is not causing any major systematic
σ = 0.027 tanh (−20.78zp2 + 7.76zp + 0.05)/(1 + zp ), i.e. error in the inferred value of the Hubble constant for
increasing from 0.0013 at zp = 0 to 0.013 at zp = 0.1; ∆` ≥ 10.
and for the Lorentzian, a(zc ) = −4zc + 3 and s(zc ) = By combining the bright standard siren measurement
0.04zc + 0.02, where zc is the midpoint of each redshift from GW170817 with a better measurement of pecu-
bin. Other choices for the redshift error (e.g. redshift- liar velocity [66], we show the corresponding poste-
independent modified Lorentzian for 2MPZ in [41] and rior on H0 in Fig. 5 with the median value of H0 =
[65]) yield very similar results. 67.0+6.3
−3.8 km/s/Mpc (68.3% ETI). This provides tighter
We then use the inferred galaxy auto power spectrum constraints in the value of the Hubble constant than pre-
to model the cross-power spectrum between the GW vious measurements. In Fig. 5 we compare the GW mea-
sources and galaxy as C`gwg (z) = bgw (z)/bg (z)C`gg (z). surements of H0 with the measurements from Planck,
On the GW side, we construct three GW maps from H0 = 67.4+0.5−0.5 km/s/Mpc [67] and with the measure-
the selected GW samples composed of Set-1 (GW190814) ment of H0 = 73.04+1.05 6
−1.05 km/s/Mpc from SH0ES [68].
, Set-2 (GW170818, GW1901412, GW190720 000836, The current measurements from the dark sirens are not
GW2001129 065458, GW200311 115853), and Set-3 sufficiently constraining yet to resolve the tension in the
(GW190701 203306, GW200224 222234). These maps value of H0 [72, 73]. Though the systematic uncertainties
are constructed on the basis of their luminosity distance in our measurement of H0 are smaller than the statisti-
distribution. Sources with a similar maximum value of cal uncertainties, in future with more GW sources and
the posterior distribution are combined to enhance the better galaxy catalog, we will be able to better assess the
cross-correlation signal. However, due to fewer sources, influence of any systematic uncertainties.
the shot noise for GW sources is very large compared to This measurement of the expansion history does not
the galaxies (by nearly 4-6 orders of magnitude depend- depend directly on the choice of GW mass distribution
ing on the redshift bin). So, the GW auto-correlation nor on whether the mass distribution follows a power-
signal is completely dominated by shot noise. law mass distribution or a power-law + Gaussian peak
The joint estimation of the Hubble constant along with model [74]. It only depends on the maximum allowed
the matter density and GW bias parameters are shown mass of individual BHs. This is primarily because the
in Fig. 4. Constraints on the Hubble constant are bi- maximum mass of the BHs determines the maximum lu-
modal, with the median value H0 = 68.2+26.0 −6.2 km/s/Mpc minosity distance up to which a source can be detected
(the upper and the lower limit indicates the 68.3% equal- from a given detector sensitivity. This luminosity dis-
tailed interval (ETI)). At the lower end of the H0 con- tance threshold, in combination with the allowed priors
straints, the constraints arise from the absence of the on the cosmological parameters, sets the maximum red-
angular cross-correlation signal between the GW sources shift in the prior on GW source redshift out to which one
and the galaxies in the first z-bin. Due to the limited needs to explore the cross-correlation signal. However,
number of GW sources, we are not able to detect the with the currently existing galaxy surveys, we are not
cross-correlation signal with galaxies (see Fig. 6 in the able to go beyond z = 0.5 due to limited sky coverage
appendix for further details.). The parameter constraints and the limited redshift reach of the galaxy catalogs.
are driven by the structure in the cross-correlation sig-
nal and its covariance for different redshifts. Also, there
are few galaxies at high redshift and the galaxy cluster-
ing signal is shot noise-dominated there. As a result, the 6 The ACT and WMAP measurement is H0 = 67.6+1.1 −1.1 km/s/Mpc
galaxy-GW source cross-correlation does not add any ad- [69]. The measurement using TRGB as a calibrator is H0 =
ditional information at high redshift. This leads to weak 69.8±0.6 (stat) ±1.6 (sys) [70]. The strong lensing measurement
constraints at high H0 and is one of the major sources of from H0LiCOW is H0 = 73.3+1.7 −1.8 km/s/Mpc [71].
 5

The GW bias parameters encode information about

the GW source population and their connection to galax- =
ies. Hence we marginalize over the GW bias parame- =
ters to infer cosmological parameters correctly and mit- =
igate population assumptions. At large angular scales

1.6 0.8 0.0 0.8 1.6 1.5 3.0 4.5 0.120.180.240.300.36

(` < 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 H0 presented in this work are from a low number
of GW sources. As a result, the clustering of the GW

0.4
0.6
0.8
0.10
0.12
0.28
0.34
0.30
6
1.5
3.0
4.5
1.6
0.8
0.0
0.8
1.6
sources is not measured with any statistical significance.

1.
So, the estimates of H0 presented in this work are weak h
and subject to systematic uncertainties associated with
small number statistics.7 To test the robustness of our FIG. 4: The joint constraints on H0 = 100h0 km/s/Mpc,
results, we have checked the following aspects, (i) ran- Ωm , and bGW (z) = bGW (1 + z)α for different choices of
bin-width ∆l.
domly varied the galaxy bias parameter within their er-
ror bars, (ii) enhanced the covariance matrix by a factor
of four, (iii)√changed the cosmological parameters such
as S8 ≡ σ8 Ωm that is used to fit the galaxy power
spectrum C`gg from S8 = 0.832 (Planck-2018 [75]) to a
lower value S8 = 0.75 as indicated by the KiDS Collabo-
ration [76], (iv) changed the value of H0 = 67 km/s/Mpc
to H0 = 74 km/s/Mpc [77] to estimate the galaxy bias
parameters, (v) changed the galaxy sample selection by
additionally removing WSC sources with a low proba-
bility of being galaxies using the SVM catalog of [78],
requiring pgal > 0.67, and (vi) changed the redshift bin
width ∆z to 0.05 which is comparable to the photo-z er-
rors. The posterior on H0 did not show any significant FIG. 5: Hubble constant H0 measurement from GWTC-3
variation for (i)–(v) cases. For the scenario (vi), the H0 dark sirens, bright siren GW170817, and combining the both
posterior show some variation and 68.3% ETI gets big- are shown along with the mean and the standard deviation
ger than the estimates presented with ∆z = 0.1. This is on the measurements from Planck-2018 [75] and SH0ES [68].
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 LVK analysis. However, at lower values of H0 , a
The measurement in this work agrees with the dark tighter constraint is obtained.
siren measurement of H0 = 67+13 −12 km/s/Mpc (68.3% Conclusion and future outlook – We present the first
HDI) by the LVK collaboration [35]. Though the cur- measurement of the Hubble constant H0 from dark stan-
rent constraints on H0 from LVK are driven by popu- dard sirens using the cross-correlation technique. The
lation assumptions, the MAP value of the distribution cross-correlation technique uses the spatial clustering of
agrees with this measurement. This is an independent GW sources with galaxies and includes information be-
validation of the LVK population assumptions. The con- yond statistical host identification. With the best eight
straints on the higher value of H0 from the LVK analysis sources available from GWTC-3, we obtain a median
arise from the empty catalog component (which is driven value of Hubble constant H0 = 68.2+26.0 −6.2 km/s/Mpc
by the population assumption). So, the upper limit on (68.3% ETI). Due to the limited number of GW sources
H0 is looser for the cross-correlation technique than for 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 [35]. In the future, with the availability
7 The measurement from other techniques such as statistical host of z < 0.8 spectroscopic galaxy catalogs such as DESI
identification and from the GW mass distribution are also af- [79] and SPHEREx [80] (supplemented by z > 0.8 spec-
fected by poor number statistics. troscopy from Euclid [81] and photometric redshifts from
 6

Vera Rubin Observatory [82]), cross-correlation of the 024001 (2014), ISSN 1361-6382, .
GW sources with galaxies will be a powerful technique [5] F. Acernese, M. Agathos, L. Aiello, A. Allocca, A. Am-
to measure the expansion history [12, 14, 16] and testing ato, S. Ansoldi, S. Antier, M. Arène, N. Arnaud, S. As-
the general theory of relativity [15, 83]. cenzi, et al. (Virgo Collaboration), Phys. Rev. Lett. 123,
231108 (2019), .
Acknowledgements – Authors are thankful to Gergely [6] T. Akutsu et al. (KAGRA), Nat. Astron. 3, 35 (2019),
Dalya for carefully reviewing the manuscript and pro- 1811.08079.
viding useful comments during the LVK internal review. [7] T. Akutsu et al. (KAGRA), arXiv:2005.05574 (2020),
2005.05574.
Authors are also thankful to Maciej Bilicki for providing
[8] B. F. Schutz, Nature 323, 310 (1986).
insightful suggestions on the paper. Research at Perime- [9] M. Oguri, Phys. Rev. D 93, 083511 (2016), .
ter Institute is supported in part by the Government [10] S. Bera, D. Rana, S. More, and S. Bose, Astrophys. J.
of Canada through the Department of Innovation, Sci- 902, 79 (2020), 2007.04271.
ence and Economic Development and by the Province [11] G. Cañas-Herrera, O. Contigiani, and V. Vardanyan, ApJ
of Ontario through the Ministry of Colleges and Uni- 918, 20 (2021), 2105.04262.
versities. SM is supported by the Simons Foundation. [12] S. Mukherjee and B. D. Wandelt, arXiv:1808.06615
(2018), 1808.06615.
AK thanks the AMTD Foundation for support. This
[13] S. Mukherjee, B. D. Wandelt, and J. Silk, Mon. Not. Roy.
analysis is carried out at the Symmetry computing fa- Astron. Soc. 494, 1956 (2020), 1908.08951.
cility of the Perimeter Institute and the Infinity clus- [14] S. Mukherjee, B. D. Wandelt, S. M. Nissanke, and A. Sil-
ter hosted by Institut d’Astrophysique de Paris. We vestri, Phys. Rev. D 103, 043520 (2021), 2007.02943.
thank Stephane Rouberol for smoothly running the In- [15] S. Mukherjee, B. D. Wandelt, and J. Silk, Mon. Not. Roy.
finity cluster. We acknowledge the use of following Astron. Soc. 502, 1136 (2021), 2012.15316.
packages in this work: Astropy [84, 85], emcee: The [16] C. Cigarrán Dı́az and S. Mukherjee, MNRAS 511, 2782
(2022), 2107.12787.
MCMC Hammer [86], Giant-Triangle-Confusogram [87],
[17] C. L. MacLeod and C. J. Hogan, Phys. Rev. D 77, 043512
healpy[88, 89], IPython [90], Matplotlib [91], NaMaster (2008), 0712.0618.
[59], NumPy [92], and SciPy [93]. This research has [18] W. Del Pozzo, Phys. Rev. D 86, 043011 (2012),
made use of data obtained from the SuperCOSMOS Sci- 1108.1317.
ence Archive, prepared and hosted by the Wide Field [19] M. Arabsalmani, V. Sahni, and T. D. Saini, Phys. Rev.
Astronomy Unit, Institute for Astronomy, the University D 87, 083001 (2013), 1301.5779.
of Edinburgh, which is funded by the UK Science and [20] H.-Y. Chen, M. Fishbach, and D. E. Holz, Nature 562,
545 (2018), 1712.06531.
Technology Facilities Council. The authors would like to
[21] R. Gray, I. M. n. Hernandez, H. Qi, A. Sur, P. R.
thank the LIGO/Virgo/KAGRA scientific collaboration Brady, H.-Y. Chen, W. M. Farr, M. Fishbach, J. R. Gair,
for providing the data. LIGO is funded by the U.S. Na- A. Ghosh, et al., Phys. Rev. D 101, 122001 (2020), .
tional Science Foundation. Virgo is funded by the French [22] M. Fishbach et al. (LIGO Scientific, Virgo), Astrophys.
Centre National de Recherche Scientifique (CNRS), the J. Lett. 871, L13 (2019), 1807.05667.
Italian Istituto Nazionale della Fisica Nucleare (INFN), [23] B. Abbott et al. (LIGO Scientific, Virgo),
and the Dutch Nikhef, with contributions by Polish and arXiv:1908.06060 (2019), 1908.06060.
[24] M. Soares-Santos et al. (DES, LIGO Scientific, Virgo),
Hungarian institutes. This material is based upon work
Astrophys. J. Lett. 876, L7 (2019), 1901.01540.
supported by NSF’s LIGO Laboratory which is a major [25] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J.
facility fully funded by the National Science Foundation. 896, L44 (2020), 2006.12611.
[26] S. Borhanian, A. Dhani, A. Gupta, K. Arun, and
B. Sathyaprakash (2020), 2007.02883.
[27] S. R. Taylor, J. R. Gair, and I. Mandel, Phys. Rev. D
85, 023535 (2012), 1108.5161.
∗
smukherjee1@perimeterinstitute.ca [28] W. M. Farr, M. Fishbach, J. Ye, and D. Holz, Astrophys.
†
akrolews@uwaterloo.ca J. Lett. 883, L42 (2019), 1908.09084.
‡
wandelt@iap.fr [29] S. Mastrogiovanni, K. Leyde, C. Karathanasis,
§
joseph.silk@physics.ox.ac.uk, silk@iap.fr E. Chassande-Mottin, D. A. Steer, J. Gair, A. Ghosh,
[1] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Aber- R. Gray, S. Mukherjee, and S. Rinaldi, arXiv:2103.14663
nathy, F. Acernese, K. Ackley, C. Adams, T. Adams, (2021), 2103.14663.
P. Addesso, R. X. Adhikari, et al. (LIGO Scientific Col- [30] Z.-Q. You, X.-J. Zhu, G. Ashton, E. Thrane, and Z.-H.
laboration and Virgo Collaboration), Phys. Rev. Lett. Zhu, Astrophys. J. 908, 215 (2021), 2004.00036.
116, 061102 (2016), . [31] M. Mancarella, E. Genoud-Prachex, and M. Maggiore
[2] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32, (2021), 2112.05728.
074001 (2015), 1411.4547. [32] S. Mukherjee (2021), 2112.10256.
[3] B. P. Abbott et al., Phys. Rev. D 93, 112004 (2016), [33] K. Leyde, S. Mastrogiovanni, D. A. Steer, E. Chassande-
[Addendum: Phys.Rev.D 97, 059901 (2018)], 1604.00439. Mottin, and C. Karathanasis (2022), 2202.00025.
[4] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Alle- [34] J. M. Ezquiaga and D. E. Holz (2022), 2202.08240.
mandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, [35] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA)
G. Ballardin, et al., Classical and Quantum Gravity 32, (2021), 2111.03604.
 7

[36] G. Dálya, G. Galgóczi, L. Dobos, Z. Frei, I. S. [61] D. J. Schlegel, D. P. Finkbeiner, and M. Davis, ApJ 500,
Heng, R. Macas, C. Messenger, P. Raffai, and R. S. 525 (1998), astro-ph/9710327.
de Souza, Mon. Not. Roy. Astron. Soc. 479, 2374 (2018), [62] A. G. A. Brown, A. Vallenari, T. Prusti, J. H. de Bruijne,
1804.05709. F. Mignard, R. Drimmel, C. Babusiaux, C. A. Bailer-
[37] G. Dálya et al. (2021), 2110.06184. Jones, U. Bastian, M. Biermann, et al., Astronomy &
[38] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore, and Astrophysics 595, A2 (2016), ISSN 1432-0746, .
M. Mancarella, arXiv:2101.12660 (2021), 2101.12660. [63] N. Koukoufilippas, D. Alonso, M. Bilicki, and J. A. Pea-
[39] A. Palmese, C. R. Bom, S. Mucesh, and W. G. Hartley cock, MNRAS 491, 5464 (2020), 1909.09102.
(2021), 2111.06445. [64] A. Balaguera-Antolı́nez, M. Bilicki, E. Branchini, and
[40] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA) A. Postiglione, MNRAS 476, 1050 (2018), 1711.04583.
(2021), 2111.03606. [65] J. A. Peacock and M. Bilicki, MNRAS 481, 1133 (2018),
[41] M. Bilicki, T. H. Jarrett, J. A. Peacock, M. E. Clu- 1805.11525.
ver, and L. Steward, Astrophys. J. Suppl. 210, 9 (2014), [66] S. Mukherjee, G. Lavaux, F. R. Bouchet, J. Jasche,
1311.5246. B. D. Wandelt, S. M. Nissanke, F. Leclercq, and
[42] M. Bilicki et al., Astrophys. J. Suppl. 225, 5 (2016), K. Hotokezaka, Astron. Astrophys. 646, A65 (2021),
1607.01182. 1909.08627.
[43] R. Zhou et al., Mon. Not. Roy. Astron. Soc. 501, 3309 [67] N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6
(2021), 2001.06018. (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)],
[44] J. A. Newman, Astrophys. J. 684, 88 (2008), 0805.1409. 1807.06209.
[45] B. Menard, R. Scranton, S. Schmidt, C. Morrison, [68] A. G. Riess et al. (2021), 2112.04510.
D. Jeong, T. Budavari, and M. Rahman, arXiv:1303.4722 [69] S. Aiola et al. (ACT), JCAP 12, 047 (2020), 2007.07288.
(2013), 1303.4722. [70] W. L. Freedman, Astrophys. J. 919, 16 (2021),
[46] S. Schmidt, B. Menard, R. Scranton, C. Morrison, and 2106.15656.
C. McBride, Mon. Not. Roy. Astron. Soc. 431, 3307 [71] K. C. Wong et al., Mon. Not. Roy. Astron. Soc. 498, 1420
(2013), 1303.0292. (2020), 1907.04869.
[47] T. H. Jarrett, T. Chester, R. Cutri, S. Schneider, [72] L. Verde, T. Treu, and A. Riess, in Tensions between the
M. Skrutskie, and J. P. Huchra, AJ 119, 2498 (2000), Early and the Late Universe (2019), 1907.10625.
astro-ph/0004318. [73] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang,
[48] M. F. Skrutskie, R. M. Cutri, R. Stiening, M. D. Wein- A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class.
berg, S. Schneider, J. M. Carpenter, C. Beichman, Quant. Grav. 38, 153001 (2021), 2103.01183.
R. Capps, T. Chester, J. Elias, et al., AJ 131, 1163 [74] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA)
(2006). (2021), 2111.03634.
[49] E. L. Wright, P. R. M. Eisenhardt, A. K. Mainzer, [75] N. Aghanim et al. (Planck), arXiv: 1807.06209 (2018),
M. E. Ressler, R. M. Cutri, T. Jarrett, J. D. Kirkpatrick, 1807.06209.
D. Padgett, R. S. McMillan, M. Skrutskie, et al., AJ 140, [76] C. Heymans, T. Tröster, M. Asgari, C. Blake, H. Hilde-
1868 (2010), 1008.0031. brandt, B. Joachimi, K. Kuijken, C.-A. Lin, A. G.
[50] N. C. Hambly, H. T. MacGillivray, M. A. Read, S. B. Sánchez, J. L. van den Busch, et al., A&A 646, A140
Tritton, E. B. Thomson, B. D. Kelly, D. H. Morgan, (2021), 2007.15632.
R. E. Smith, S. P. Driver, J. Williamson, et al., MNRAS [77] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and
326, 1279 (2001), astro-ph/0108286. D. Scolnic, Astrophys. J. 876, 85 (2019), 1903.07603.
[51] N. C. Hambly, M. J. Irwin, and H. T. MacGillivray, MN- [78] T. Krakowski, K. Malek, M. Bilicki, A. Pollo, A. Kurcz,
RAS 326, 1295 (2001), astro-ph/0108290. and M. Krupa, A&A 596, A39 (2016), 1607.01188.
[52] N. C. Hambly, A. C. Davenhall, M. J. Irwin, and [79] A. Aghamousa et al. (DESI) (2016), 1611.00036.
H. T. MacGillivray, MNRAS 326, 1315 (2001), astro- [80] O. Dore et al., arXiv (2018), 1805.05489.
ph/0108291. [81] A. Blanchard et al. (Euclid), Astron. Astrophys. 642,
[53] J. A. Peacock, N. C. Hambly, M. Bilicki, H. T. A191 (2020), 1910.09273.
MacGillivray, L. Miller, M. A. Read, and S. B. Tritton, [82] LSST Science Collaboration, P. A. Abell, J. Allison, S. F.
MNRAS 462, 2085 (2016), 1607.01189. Anderson, J. R. Andrew, J. R. P. Angel, L. Armus, D. Ar-
[54] H. S. Xavier, M. V. Costa-Duarte, A. Balaguera- nett, S. J. Asztalos, T. S. Axelrod, et al., ArXiv e-prints
Antolı́nez, and M. Bilicki, JCAP 2019, 037 (2019), (2009), 0912.0201.
1812.08182. [83] S. Mukherjee, A. Krolewski, et al., Under preparation
[55] A. Krolewski et al., Under preparation (2022–). (2022–).
[56] A. A. Collister and O. Lahav, PASP 116, 345 (2004), [84] Astropy Collaboration, T. P. Robitaille, E. J. Tollerud,
astro-ph/0311058. P. Greenfield, M. Droettboom, E. Bray, T. Aldcroft,
[57] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Davis, A. Ginsburg, A. M. Price-Whelan, et al., A&A
M. Oguri, ApJ 761, 152 (2012), 1208.2701. 558, A33 (2013), 1307.6212.
[58] E. Hivon, K. M. Górski, C. B. Netterfield, B. P. Crill, [85] Astropy Collaboration, A. M. Price-Whelan, B. M.
S. Prunet, and F. Hansen, ApJ 567, 2 (2002), astro- Sipőcz, H. M. Günther, P. L. Lim, S. M. Crawford,
ph/0105302. S. Conseil, D. L. Shupe, M. W. Craig, N. Dencheva, et al.,
[59] D. Alonso, J. Sanchez, and A. Slosar (LSST Dark Energy AJ 156, 123 (2018), 1801.02634.
Science), Mon. Not. Roy. Astron. Soc. 484, 4127 (2019), [86] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Good-
1809.09603. man, PASP 125, 306 (2013), 1202.3665.
[60] J. Grain, M. Tristram, and R. Stompor, Phys.Rev.D 79, [87] S. Bocquet and F. W. Carter, The Journal of Open
123515 (2009), 0903.2350. Source Software 1 (2016), .
 8

[88] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt,

F. K. Hansen, M. Reinecke, and M. Bartelman, Astro-
phys. J. 622, 759 (2005), astro-ph/0409513.
[89] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset,
E. Hivon, and K. Gorski, Journal of Open Source Soft-
ware 4, 1298 (2019), .
[90] F. Pérez and B. E. Granger, Computing in Science and
Engineering 9, 21 (2007), ISSN 1521-9615, .
[91] J. D. Hunter, Computing In Science & Engineering 9, 90
(2007).
[92] S. van der Walt, S. C. Colbert, and G. Varoquaux,
Computing in Science and Engineering 13, 22 (2011),
1102.1523.
[93] E. Jones, T. Oliphant, P. Peterson, et al., SciPy: Open
source scientific tools for Python (2001–), [Online; ac-
cessed ¡today¿], .
[94] D. Alonso, A. I. Salvador, F. J. Sánchez, M. Bilicki,
J. Garcı́a-Bellido, and E. Sánchez, MNRAS 449, 670
(2015), 1412.5151.
[95] 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), 2106.04354.
[96] A. Krolewski, S. Ferraro, E. F. Schlafly, and M. White,
JCAP 2020, 047 (2020), 1909.07412.
[97] C. P. Novaes, A. Bernui, H. S. Xavier, and G. A. Mar-
ques, MNRAS 478, 3253 (2018), 1805.04078.

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 sig-
nal with statistical significance. We show the measured
band average cross-correlation signal Ĉ`GW g as a func-
tion of redshift in Fig. 6 for ∆` = 15 with the three GW FIG. 6: The band average cross-correlation power
maps constructed from the eight events. The value of spectrum between GW sources and galaxies with a ∆` = 15
maximum multipole `max = 30 is chosen to take into ac- bin-width is shown for three different maps of the GW
count the poor sky localization error of the GW sources sources Set-1 (GW190814) , Set-2 (GW170818, GW1901412,
for sources with sky localization error ∆Ω ≤ 30 sq. deg. GW190720 000836, GW2001129 065458,
GW200311 115853), and Set-3 (GW190701 203306,
The diagonal error bars are shown on the measured val-
GW200224 222234) composed from the selected eight
ues. The band average signal shows a signal that is events, as a function of the median value of the redshift bin.
consistent with zero at all the redshift bins. The con-
straints on the clustering signal at low redshift z for all
the three cases drives the constraints on the lower value
of H0 ≤ 50 km/s/Mpc. For the higher redshift bins, a Construction of galaxy mask for 2MPZ and WSC
non-zero mean value causes a non-zero posterior distri-
bution on H0 which is different from the prior (taken as
U[20, 120] km/s/Mpc. This indicates that the constraints The galaxy masks were carefully constructed to remove
are weak and are affected by the systematic associated areas with large numbers of stars or other systematics
with the non-detection of the cross-correlation signal. At that could affect galaxy clustering, either by direct stellar
high redshift (z > 0.5) the absence of sufficient galaxies contamination or by correlations, e.g. suppressed galaxy
leads to an abrupt cut in the galaxy distribution and pro- density in regions of high stellar density or extinction.
vides no information on higher values of H0 . Due to the We follow [64] to construct the 2MPZ mask, starting by
absence of detection of the clustering signal, the statis- masking low Galactic latitudes (|b| < 10◦ , areas of high
tical power of measuring H0 is similar to the statistical galactic extinction (E(B − V ) > 0.3 from [61], and ar-
host identification technique. eas of high stellar density as estimated from the 2MASS
 9

Point Source Catalog (log nstar > 3.5)8 . We further in- [96]; and adding a mask of regions in WISE with high
clude manual cutouts around the LMC and SMC, exclud- moon contamination, as determined by HEALPix pixels
ing 275.47 < RA < 285.47 and −37.89 < DEC < −27.89, in which GLADE+ [37] is incomplete compared to WSC.
and 300.81 < RA < 304.81 and −46.33 < DEC < We also test variations in the sample-selection procedure,
−42.33. Finally, we mask additional areas with low com- i.e. additionally using the SVM catalog of [78] to restrict
pleteness in 2MPZ, determined by comparing the number to likely galaxies [95, 97]. We find that these variations
counts of 2MPZ sources and 2MASS XSC sources (with generally lead to a scale-independent shift in the ampli-
Ks < 13.9) in NSIDE=64 HEALPixels. We remove pix- tude of C`gg , either corresponding to a change in galaxy
els with < 85% completeness, mostly corresponding to bias due to differing populations, or a change in the stel-
areas of lower depth around bright stars. We test vari- lar contamination fraction, which is entirely degenerate
ations in the masking procedure (i.e. additionally multi- with bias at ` > 10 where the stellar power spectrum is
plying by the WSC mask, following [94], or changing the small. In this regime, the effect of changing stellar con-
completeness threshold to 80% or 90%) and find minimal tamination is degenerate with bias in both the galaxy
changes in results. auto spectrum and the galaxy cross-spectrum with GW
sources, so it will cause systematic errors in our modeling.
For WSC, we follow the masking procedure of [54]. We assume that the galaxy bias is redshift-independent
We start with the mask distributed with the WSC data in each bin, and obtain best-fit values of bg (zc = 0.05) =
release [42]9 . We additionally mask regions with high ex- 1.18, bg (zc = 0.15) = 0.66, bg (zc = 0.25) = 1.35,
tinction (E(B −V ) > 0.10) and high stellar density (den- bg (zc = 0.35) = 1.76, and bg (zc = 0.45) = 2.33. The
sity of stars from GAIA greater than 7 times the mean). very low value of bg in the second bin is driven by the
We additionally test several variations in the masking SuperCOSMOS plate template, which is degenerate with
procedure, adding an additional mask at low Galactic lat- the cosmological contribution due to the limited multi-
itudes following [95]; adding a WISE bright stars mask pole range considered (10 < ` < 40).

8 https://www.ipac.caltech.edu/2mass/releases/allsky/ 9
doc/sec4 5c.html