# Spatial clustering of gravitational wave sources with $k$-nearest neighbour
distributions

Kaustubh Rajesh Gupta1 and Arka Banerjee1
1Department of Physics, Indian Institute of Science Education and Research,
Homi Bhabha Road, Pashan, Pune 411008, India
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a framework to quantify the clustering of gravitational wave (GW)
transient sources and measure their spatial cross-correlation with the large-
scale structure (LSS) of the universe using the $k$-nearest neighbour ($k$NN)
formalism. As a first application, we measure the nearest-neighbour
distributions of $53$ suitably selected Binary Black Hole (BBH) mergers
detected in the first three observation runs of LIGO-Virgo-KAGRA and cross-
correlate these sources with $\mathchar 21016\relax\,1.7\times 10^{7}$
galaxies and quasars from the WISE$\times$SuperCOSMOS all-sky catalogue. To
determine the significance of the clustering signal while accounting for
observational systematics in the GW data, we create 135 realisations of mock
BBHs that are statistically similar to the observed BBHs but spatially
unclustered. We find no evidence for spatial clustering or cross-correlation
with LSS in the data and conclude that the present sky localisation and number
of detections are insufficient to get a statistically significant clustering
signal. Looking forward, the statistically large number of detections and the
significant improvements in sky localisations expected from future observing
runs of LIGO (including LIGO India) and the next generation of GW detectors
will enable measurement of the BBH-LSS cross-correlation and open a new window
into cosmology.

###### keywords:

methods: statistical – cosmology: large-scale structure of Universe –
gravitational waves

††pubyear: 2024††pagerange: Spatial clustering of gravitational wave sources
with $k$-nearest neighbour distributions–D

##  1 Introduction

The clustering of objects that trace structure formation in the universe
contains a wealth of information that can be used to test the standard model
of cosmology and constrain its parameters (see, e.g., DES Collaboration et
al., 2022a, b, c, 2023; Dvornik et al., 2023; Amon et al., 2023; Miyatake et
al., 2023; Fumagalli, A. et al., 2024). The detection of gravitational waves
by the LIGO-Virgo-KAGRA (LVK) collaboration (LIGO Scientific Collaboration and
Virgo Collaboration et al., 2016) has unveiled potential new tracers of
structure formation in the form of merging binaries of compact stellar
remnants such as black holes and neutron stars (LIGO Scientific Collaboration
et al., 2023b).

Gravitational waves allow a direct measurement of the luminosity distance to
their sources (see, e.g., Schutz, 1986; Holz & Hughes, 2005; Holz et al.,
2018) without the need for a hierarchical distance ladder or an empirical
calibration process, with the only fundamental assumption being that general
relativity is valid, making merging compact binaries ‘standard sirens’. Hence,
gravitational waves provide a mechanism to study the expansion history of our
universe if the source redshifts can be estimated. This led Schutz (1986) to
suggest that merging compact binaries can be used to constrain the Hubble-
Lemaître parameter $H_{0}$, which characterises the present-day rate of
expansion of the universe.

Since the redshift of a merging binary cannot be inferred directly from its
gravitational waves, many techniques have been developed in the literature to
measure $H_{0}$ using additional astrophysical observations (see
Mastrogiovanni et al. (2024) for a recent review). For example, the ‘bright
siren’ method uses direct electromagnetic counterparts of the merger events to
obtain redshifts (Abbott et al., 2017). In contrast, the statistical dark
siren method (see Gair et al. (2023) for a review) uses galaxy surveys to
identify potential hosts of the merger events inside the localisation volumes
provided by gravitational wave observations (Soares-Santos et al., 2019;
Palmese et al., 2020; Abbott et al., 2021; Alfradique et al., 2024). MacLeod &
Hogan (2008) proposed a method that uses galaxy clustering to extract redshift
information for a sample of merger events in a statistical sense to estimate
$H_{0}$ without needing to identify host galaxies for individual merger
events. Methods have also been proposed that try to estimate the redshift of
gravitational wave sources by breaking the mass-redshift degeneracy in
gravitational wave analyses; this can be achieved, for example, by
constraining the neutron star tidal deformability or by combining features in
the mass distribution and redshift evolution of merger rate (the so-called
‘spectral siren method’), to measure the source masses (Farr et al., 2019a, b;
Mastrogiovanni et al., 2021; Ezquiaga & Holz, 2022; Mancarella et al., 2022;
Abbott et al., 2023).

All of these methods, however, have various drawbacks (Mastrogiovanni et al.,
2024). The bright siren method relies on detecting rare events accompanied by
electromagnetic counterparts and possibly only applies to binary neutron star
(BNS) mergers. The dark siren method suffers from difficulties due to large
localisation volumes of gravitational wave events and is susceptible to
potential biases in the inference of $H_{0}$ due to the incompleteness of
galaxy catalogues (Trott & Huterer, 2022). The mass-redshift degeneracy method
is model-dependent; uncertainties in the modelling of the neutron star
equation of state, or a wrong model of the binary merger rate, can introduce
systematic biases in the cosmological inference (see section 2.3.3 of
Mastrogiovanni et al. (2024) and references therein).

If star-forming regions follow the fluctuations in the underlying cosmological
matter field, merging compact binaries are expected to be inherently clustered
and spatially correlated with other tracers such as galaxies and galaxy
clusters (Scelfo et al., 2018). The strength of the cross-correlation between
sources of gravitational waves and the large-scale structure of the universe
is sensitive to cosmological parameters. It can, therefore, be used as an
independent probe of the Hubble-Lemaître constant (Oguri, 2016; Bera et al.,
2020; Mukherjee et al., 2021, 2022), after marginalizing over the other
relevant parameters. As long as the gravitational wave sources and the galaxy
sample trace fluctuations in the same underlying density field, a measurement
of their spatial cross-correlation does not require uniquely identifying the
source of each merger event (Fang et al. (2020) have a similar discussion in
the context of cross-correlations between high-energy neutrinos and large-
scale structure). Therefore, the cross-correlation method of determining
$H_{0}$ does not suffer from biases due to the incompleteness of the galaxy
catalogues used (see Bera et al., 2020, for a systematic study). Moreover,
this method does not require any direct assumptions about the population
properties of the merging binaries. Therefore, it measures $H_{0}$ nearly
independently of merging binary population models111It is to be noted,
however, that the clustering of gravitational wave sources is also expected to
be a function of bias parameters that model the tracer-matter connection. One
needs to marginalise over these parameters to obtain constraints on
cosmological parameters (see, e.g., Banerjee et al., 2022; Peron et al., 2023,
and references therein). Since the bias parameters are controlled by the
source population properties (Raccanelli et al., 2016; Adhikari et al., 2020;
Peron et al., 2023; Vijaykumar et al., 2023a), an indirect dependence on
population models is introduced in the measurement of $H_{0}$..

Since the various techniques of measuring $H_{0}$ using gravitational wave
sources discussed above do not utilise information in the temperature
fluctuations of the cosmic microwave background (CMB) (Ade et al., 2016;
Aghanim et al., 2020) or cosmic distance ladder distance measurements from
supernovae and other standard candles (Riess et al., 2018, 2019; Riess, 2020;
Wong et al., 2020; Riess et al., 2022), measuring the cross-correlation
between gravitational wave sources and the large-scale structure of the
universe is important in the context of the so-called Hubble Tension (see Hu &
Wang (2023) and Valentino et al. (2021) for comprehensive reviews). In
addition to cosmology, these cross-correlation measurements can also be used
to study the astrophysical origins and formation channels of gravitational
wave sources (see, e.g., Raccanelli et al., 2016; Scelfo et al., 2018, 2020;
Adhikari et al., 2020; Gagnon et al., 2023). With the third generation of
gravitational wave detectors likely to bring in $\mathchar
21016\relax\,10^{5}$ more detections per year (Iacovelli et al., 2022;
Borhanian & Sathyaprakash, 2022), measuring the clustering of these objects
and modelling it as a function of cosmological parameters will, therefore,
play an essential role for both precision cosmology and compact binary
astrophysics in the coming decade.

There have been a few attempts to measure the angular two-point correlation
function (Cavaglià & Modi, 2020; Zheng et al., 2023) and the angular power
spectrum (Zheng et al., 2023) of the currently detected LVK events; to
determine the level of anisotropy in their sky distribution (Payne et al.,
2020; Essick et al., 2023); and to measure their spatial cross-correlation
with galaxy catalogues (Mukherjee et al., 2022). However, a statistically
significant detection of anisotropy or clustering has not yet been achieved.
Studies using forecasts for future detectors have also been performed
(Namikawa et al., 2016; Scelfo et al., 2018; Calore et al., 2020; Libanore et
al., 2021, 2022; Balaudo et al., 2023; Vijaykumar et al., 2023b; Gagnon et
al., 2023), primarily focusing on two-point summary statistics.

In this paper, we present a framework to quantify the clustering of
gravitational wave sources and their spatial cross-correlation with large-
scale structure catalogues using the $k$-nearest-neighbour distributions
(Banerjee & Abel, 2021a) as summary statistics. The nearest-neighbour
measurements are sensitive to all $N$-point correlation functions of the
tracers and hence are a much more powerful probe of clustering, compared to
the two-point function, on scales where the underlying matter field is not
well-approximated as a Gaussian random field, and the effect of gravitational
nonlinearities cannot be neglected (Banerjee & Abel, 2021a, b, 2023).
Application of these statistics could, in principle, lead to a detection of
the clustering signal from the same datasets used in previous two-point
analyses. To enable this new analysis, we extend the $k$NN formalism,
originally presented for 3D clustering in cartesian coordinates, to angular
clustering in the sky. As a first application to data, we compute the auto-
correlation of a suitable subset of the binary black holes (BBHs) detected in
the first three observing runs of LVK and their cross-correlation with the
WISE$\times$SuperCOSMOS all-sky survey. We also compare the results of the
two-point and nearest-neighbour analyses. Although we focus on BBHs in this
work, our framework can easily be extended to study the clustering of other
gravitational wave transients like binary neutron stars and neutron star-black
hole binaries.

The rest of the paper is structured as follows. We describe the data used in
this study in section 2. In section 3, we develop the mathematical formalism
for clustering statistics and discuss how to compute them numerically. At the
end of section 3, we present an illustrative example that demonstrates the
potential boost in the clustering signal of sparsely sampled tracers expected
from the nearest-neighbour measurements on small spatial scales over the two-
point summary statistics. Having motivated the $k$NN formalism, in section 4,
we describe our methods: the null hypothesis, our procedure to test the null
hypothesis, how we determine the statistical significance of the clustering
signal, and our strategy to deal with observational selection biases in the
data. We present our results in section 5 and discuss some interesting aspects
of this study. Finally, we summarise, draw conclusions, and discuss possible
future directions in section 6. Some additional material is presented in the
appendix.

##  2 Data

In this section, we discuss the data used in this study. We describe the
gravitational wave events selected for this work in section 2.1. Typically, a
mock catalogue of unclustered data points (known as ‘randoms’ in the
literature) is required to get a reliable measurement of the statistical
significance of the clustering signal in the presence of observational
selection biases (see., e.g., Wang et al., 2022). In section 2.2, we motivate
this requirement for the specific case of gravitational wave data, discuss the
procedure to create the unclustered catalogue, and present the resulting mock
data. Finally, we describe the large-scale structure catalogue used for cross-
correlating the BBHs in section 2.3.

###  2.1 Gravitational Wave Events

This work uses the compact binary merger events detected in the first three
observing runs of LIGO-Virgo-KAGRA, as reported in LIGO Scientific
Collaboration et al. (2023b). Following Zheng et al. (2023), from this parent
set of $\mathchar 21016\relax\,80$ events, we select the events detected with
a false alarm rate (FAR) less than 1 per year and crossed a detection
threshold of network matched-filtered signal-to-noise ratio (SNR) greater than
10. Since we are interested in binary black holes (BBHs), we further restrict
our sample to those events that have a probability of being a BBH merger
greater than 0.5222The classification probabilities were calculated using the
GW package of PESummary (Hoy & Raymond, 2021)..

Zheng et al. (2023) restrict their sample to events detected in all three
detectors that were in science mode during the LVK observing period, namely
LIGO Livingston, LIGO Hanford (LIGO Scientific Collaboration et al., 2015) and
Virgo (Accadia et al., 2012), to get better-localised events. However, this
step removes a significant fraction of the BBHs selected above. Since it is
not clear a priori whether the resulting gain in sky localisation accuracy
would compensate for the reduction in the BBH sample size, we keep two-
detector events in our final sample. To ensure better homogeneity in the sky
localisations of the BBHs, we remove all events from our sample detected
before the Virgo detector joined the observing run. Note that a non-detection
in one of the detectors does carry some information about the location of the
merger event in the sky. Hence, the two-detector events not detected in Virgo
when it is in science mode are still expected to be better localised than
those observed in the absence of Virgo.

Finally, we are left with 53 BBHs that constitute our observed catalogue.
Using the parameter estimation posterior samples on declination and right
ascension made publicly available by the LVK collaboration, we generate
skymaps for each event representing the uncertainty in their localisation in
the sky. Figure 1 shows a combined skymap of all events generated by stacking
the individual skymaps. We summarise the properties of the observed BBHs that
are most relevant for clustering, namely the distribution of sky localisation
areas and observed luminosity distances, in Figure 2. Distributions of other
BBH properties are summarised in appendix A.

Figure 1: Mollweide projection of the combined skymap of the 53 observed
events that constitute our BBH catalogue, in equatorial (J2000) coordinates.
Each banana-shaped cloud represents a single BBH. Skymaps were generated
through the Healpy package using the parameter estimation posterior samples
provided by the LVK collaboration for each event. The colour represents the
number of posterior samples per pixel in a logarithmic scale. The HEALPix
NSIDE for this map is 64.

(a) 68% credible area (b) Luminosity distance

Figure 2: Distribution of the $1\sigma$ sky localisation uncertainty areas
(left) and luminosity distances (right) of the observed BBH catalogue. The
credible areas are computed using the rapid Bayesian localisation code
BAYESTAR (Singer & Price, 2016).

###  2.2 Mock BBH Catalog

In this section, we describe the procedure for creating the mock BBH catalogue
that will serve as the set of ‘randoms’ used for the clustering analysis. At
first, it appears that simply distributing points uniformly in the sky should
be sufficient, as the resulting data set would be unclustered and uncorrelated
with large-scale structure. This would be a valid approach if we had perfect
observations of a sample of BBHs representative of the entire BBH population
in the universe. Unfortunately, due to the limited sensitivity of the current
gravitational wave detectors, the data are plagued with selection biases and
systematic effects; the BBHs selected for this study do not constitute a
representative sample of the population333It should be noted, however, that
the third generation of gravitational wave detectors is expected to detect
almost all BBH mergers in the universe up to a very high redshift (Hall &
Evans, 2019; Iacovelli et al., 2022; Borhanian & Sathyaprakash, 2022)..
Furthermore, the detectors are not equally sensitive to all regions in the
sky; each detector is most sensitive to the merger events that go off directly
on top of it, i.e., perpendicular to the plane of detector arms. This is a
consequence of the transverse nature of gravitational waves. As a result,
there is a selection function in the sky for the detector network as a whole
(see Chen et al. (2017) for example). These observational systematics have to
be carefully folded into the clustering analysis to avoid getting biased or
spurious signals.

One way of mitigating the selection biases outlined above is to create
realistic mock BBHs that reproduce the properties of the observed BBH sample
in a statistical sense but which are inherently unclustered and spatially
uncorrelated with the large-scale structure of the universe. Such a mock data
set allows us to naturally incorporate the effects of observational biases on
the clustering measurements.

We follow a procedure outlined in Zheng et al. (2023) to create our mock
catalogue. First, we distribute a population of BBH merger events
isotropically in the sky by sampling their locations from a uniform
distribution ($\mathcal{U}$), which translates to drawing their right
ascension ($\alpha$) from $\mathcal{U}(0,2\pi)$ and sine of declination
($\sin\delta$) from $\mathcal{U}(-1,1)$. We next draw their source parameters
from the population distributions inferred by the LVK collaboration (LIGO
Scientific Collaboration et al., 2023a) as implemented by the GWPopulation
package444https://colmtalbot.github.io/gwpopulation/. These are as follows:

  1. 1.

Power Law + Peak model for the mass of the primary (heavier) BBH and a power
law distribution for the ratio of component masses (Talbot & Thrane, 2018)

  2. 2.

power law distribution for redshift evolution of merger rate per unit comoving
volume per unit source-frame time (Fishbach et al., 2018)

The mathematical details of these models are discussed in appendix C. We
assume uniform distributions for inclination angle $\iota$ w.r.t. the plane of
orbital angular momentum, polarisation angle $\psi$ and phase at coalescence
$\Phi_{c}$ over their allowed physical ranges, i.e.,
$\iota,\psi\in\mathcal{U}(0,\pi)$ and $\Phi_{c}\in\mathcal{U}(0,2\pi)$, and
uniformly sample the BBH merger time during the LIGO observation period after
Virgo started taking data. For simplicity, we set the black hole spins
identically to zero since we do not expect them to affect the clustering
properties or the sky localisation uncertainties, which are most relevant to
us. We have further checked that including the spins does not affect our
analysis; the final mock BBH catalogues with and without spins turned on are
statistically similar in all aspects.

After creating the mock BBH population, we determine which events can be
‘detected’ by the current gravitational wave detectors to reproduce the
selection biases in the data. Here, we consider a detector network consisting
of LIGO Livingston, LIGO Hanford and Virgo, the same network that collected
the data for our observational sample. We conduct the following gravitational
wave data analysis using the rapid Bayesian localisation code for
gravitational wave events, BAYESTAR555We follow a similar procedure to the one
outlined in https://lscsoft.docs.ligo.org/ligo.skymap/quickstart/bayestar-
injections.html. (Singer & Price, 2016). First, we simulate the gravitational
wave signals for each BBH using the IMRPhenomXPHM model (Pratten et al.,
2021), which is the same waveform used in the LVK analysis of the data. Next,
we inject the simulated signals in stationary Gaussian noise created using
analytic estimates for the third observing run power spectral densities for
the LIGO Livingston, LIGO Hanford and Virgo detectors, as provided by the
PyCBC package666http://pycbc.org/pycbc/latest/html/index.html. Finally, we
compute the (phase-maximised) network matched-filtered signal-to-noise (SNR
henceforth) for each event and classify the events with an SNR $\geq 10$ as
‘detections’777Technically, the false alarm rate (FAR) is a better measure of
whether an event should be considered detectable, but computing FARs is
computationally expensive, as it requires doing parameter estimation for the
full set of injected events. An SNR cutoff of 10 is a reasonable proxy (see
(LIGO Scientific Collaboration et al., 2023b) or (Essick, 2023) for more
details).. Once we have the selected events, we use BAYESTAR to localise them.
BAYESTAR also returns the estimated luminosity distances and credible
intervals for the area of sky localisation uncertainty.

Using the procedure outline above, we generate a mock catalogue of $135$
realisations of $53$ BBHs each. The combined skymaps of $4$ sample
realisations are shown in figure 3, with the colour palate and resolution
identical to figure 1 for ease of comparison. The skymaps of the observed and
mock BBHs are visually similar. To investigate if the mock catalogue is
statistically similar to the observational data, we compare their sky
localisation uncertainty and luminosity distance distributions, averaged over
the $135$ realisations, with the corresponding distributions measured in the
data. The results are shown in figure 4. Similar plots for other source
properties are presented in appendix B.

Within the limit of sample variance across the realisations, the mock
catalogue is reasonably statistically similar to the observed BBH sample. It
does includes more events with higher uncertainty in the sky localisation, as
seen from the $\mathchar 21016\relax\,3\sigma$ deviation from the observed
distribution in the second last histogram bin of figure 4(a). However, this is
not expected to bias the clustering results since it would only lead to
slightly larger measurement errors for the clustering statistics of each mock
realisation. As we discuss in section 4.3, the relevant quantity for measuring
the significance of the clustering signal is the variance across the
realisations, which is independent of the measurement errors on the individual
realisations.

(a) Sample realisation 1 (b) Sample realisation 2

(c) Sample realisation 3 (d) Sample realisation 4

Figure 3: Mollweide projection of the combined skymap of the 4 sample
realisations of the mock BBH catalogue in equatorial (J2000) coordinates. Each
banana-shaped cloud represents a single BBH. These are visually similar to
figure 1. Skymaps were generated using BAYESTAR. As before, the colour
represents the number of posterior samples per pixel in a logarithmic scale.
The HEALPix NSIDE for these maps is 64.

(a) 68% credible area (b) Luminosity distance

Figure 4: Distribution of the $1\sigma$ sky localisation uncertainty areas
(left) and luminosity distances (right) of the observed (bold histograms) and
mock (light histograms) BBHs. The error bars on the mock histograms show the
variance over 135 realisations of the mock catalogue. Except for one or two
bins, the mock and the observed data agree reasonably within the error bars,
signifying that the mock dataset statistically reproduces the observations.

###  2.3 Galaxy Catalogue

We use galaxies and quasars from the publicly available
WISE$\times$SuperCOSMOS (hereafter WSC) catalogue (Bilicki et al., 2016),
which is a cross-match between two parent full-sky catalogues: the AllWISE
release (Cutri et al., 2013) from the Wide-field Infrared Survey Explorer
(WISE) (Wright et al., 2010), a mid-infrared ($\lambda\mathchar
21016\relax\,\mu m$) space survey; and the SuperCOSMOS Sky Survey (Hambly et
al., 2001), consisting of data from digitised optical photographic plates
taken by the United Kingdom Schmidt Telescope (UKST) in the southern
hemisphere888https://www.roe.ac.uk/ifa/wfau/ukstu/telescope.html and the
Palomar Observatory Sky Survey-II (POSS-II), in the northern hemisphere (Reid
et al., 1991). WISE, a NASA space-based mission, surveyed the entire sky in
four bands, $W_{1}=3.4\mu m$, $W_{2}=4.6\mu m$, $W_{3}=12\mu m$, and
$W_{4}=23\mu m$, while SuperCOSMOS has data in three optical bands, $B$, $R$,
and $I$. The interested reader is referred to Bilicki et al. (2016) for more
details on the WISE and SuperCOSMOS surveys and the cross-matching procedure.
We use photometric redshifts for the WSC catalogue provided by Bilicki et al.
(2016), which have been estimated using the artificial neural network code
ANNz (Collister & Lahav, 2004).

We work with the ‘SVM’ release of the WSC catalogue (Krakowski et al., 2016),
which classifies sources into galaxies, stars and quasars using a support
vector machines (SVM) learning algorithm. The reason for choosing the SVM
catalogue is as follows: in creating the original WSC catalogue, colour cuts
were placed that already removed the quasars. Since we expect quasars to trace
the large-scale fluctuations in the universe alongside galaxies, removing
quasars is unnecessary for a cross-correlation study like ours. We remove the
sources classified as ‘stars’ in the SVM catalogue and select the remaining
objects to form our raw catalogue. Finally, we remove sources lying in
problematic regions in the sky with unreliable data using the publicly
available WSC mask to create the final catalogue cross-correlated with the BBH
catalogue created in section 2.1. This process removes regions such as those
obscured by the plane of the Milky Way and by the Large and Small Magellanic
Clouds (SMC and LMC) and the areas with high stellar contamination (see
Bilicki et al. (2016) for a detailed description of the masking procedure).
Some authors (for example, Mukherjee et al., 2022) impose additional colour
cuts on $E(B-V)$ and/or on $W_{1}-W_{2}$ to mitigate dust extinction and
further reduce stellar contamination. However, this increases the purity of
the galaxy sample at the cost of completeness, which is undesirable for cross-
correlation studies. Since we do not expect nearby stars to correlate with the
extragalactic BBHs, we do not impose any additional colour cuts.

After masking out regions in the sky with unreliable data, we are left with a
catalogue that covers $\mathchar 21016\relax\,3\pi$ steradians in the sky,
corresponding to a sky coverage of $\mathchar 21016\relax\,68\%$. This makes
this catalogue suitable for a cross-correlation study with BBH catalogues,
which are inherently all sky since gravitational waves are not susceptible to
medium propagation effects999Note, however, that gravitation waves, like
light, are indeed susceptible to weak gravitational lensing (Oguri, 2016;
Meena & Bagla, 2019).. The sky distribution of the WSC galaxies and quasars is
shown in figure 5.

Figure 5: Mollweide projection of the skymap of $\mathchar
21016\relax\,1.7\times 10^{7}$ galaxies and quasars in the WSC catalogue, in
equatorial (J2000) coordinates. Skymaps were generated from the sky locations
of the sources using the Healpy package. The colour represents the number of
sources per pixel on a linear scale. The colour bar has been limited to a
maximum of 50 samples per pixel to enhance contrast, which makes it easier to
visualise the large-scale structure in the distribution of the galaxies and
quasars. The empty navy regions represent regions in the sky with unreliable
data and have been masked out. The HEALPix NSIDE for this map is 256.

Furthermore, the redshift distribution of the WSC sources significantly
overlaps with that of the BBHs selected for this study, as shown in figure 6.
This is important for cross-correlation studies since we do not expect
cosmological fluctuations at different redshifts to be correlated. Note that
the redshifts of the BBHs are not direct observables but are computed from the
luminosity distances, assuming a cosmological model for the expansion history
of the universe. For this dataset, the LVK collaboration (LIGO Scientific
Collaboration et al., 2023b) assumed a cosmological model consistent with the
Planck 2015 results (Ade et al., 2016).

Figure 6: A comparison of the redshift distributions of the observed BBHs (red
histogram) and the WSC catalogue sources (blue histogram). There is a
significant overlap between the redshifts of the two datasets, which is
crucial for conducting cross-correlation studies since we do not expect
cosmological fluctuations at different redshifts to be correlated. Also
plotted are filled histograms showing the distribution of galaxy (orange) and
quasar (green) redshifts separately. The quasars are at a higher redshift on
average.

The final catalogue contains $\mathchar 21016\relax\,1.7\times 10^{7}$
sources, with $\mathchar 21016\relax\,15$ million classified as galaxies and
$\mathchar 21016\relax\,2$ million as quasars. This translates to an average
number density of more than $600$ sources per sq. deg. in the sky.

##  3 Mathematical formalism

This section describes the summary statistics used to quantify clustering
strength. For each statistic, we give a mathematical definition followed by a
computational recipe to calculate the clustering strength for given data using
that particular statistic. We discuss the auto-clustering of a set of tracers
in section 3.1 and the cross-clustering between two sets of tracers in section
3.2. In section 3.3, we present an illustrative example that demonstrates the
potential of the nearest-neighbour formalism to measure the clustering of a
sparse sample of tracers and motivates their application to a BBH-galaxy
cross-correlation study.

###  3.1 Auto-clustering

Consider a set $X$ of $N_{X}$ discrete, point-like tracers or data
points101010We use tracer and data point interchangeably in this paper..

$X\equiv\left\\{\left(\delta_{1},\alpha_{1}\right),\left(\delta_{2},\alpha_{2}\right),...,\left(\delta_{N_{X}},\alpha_{N_{X}}\right)\right\\}$
(1)

where $\delta_{i}$, $\alpha_{i}$ represents the declination and right
ascension of the $i^{\rm th}$ tracer in celestial equatorial (J2000)
coordinates. In polar coordinates, the position of the $i^{\rm th}$ tracer is
given by

$\displaystyle\theta_{i}$ $\displaystyle=\delta_{i}-\pi/2$ (2b)
$\displaystyle\phi_{i}$ $\displaystyle=\alpha_{i}$ (2a)

To study clustering, we need a metric for the distance between two points
$\left(\delta_{1},\alpha_{1}\right),\left(\delta_{2},\alpha_{2}\right)$ in the
sky. A natural choice is the great-circle distance given by $d=\theta$, where
$\theta$ is the central angle between the two points on the sphere. Note that
we treat the sky as a unit sphere; hence, the radius term that usually
multiplies the angle to get the distance is absent. The central angle between
$\left(\delta_{1},\alpha_{1}\right),\left(\delta_{2},\alpha_{2}\right)$ is
computed using the haversine formula (RIOS, 1795)

$\text{hav}\left(\theta\right)=\text{hav}\left(\delta_{2}-\delta_{1}\right)+\cos\delta_{1}\cos\delta_{2}\text{hav}\left(\alpha_{2}-\alpha_{1}\right)$
(3)

where $\text{hav}\left(\theta\right)\triangleq\sin^{2}\left(\theta/2\right)$
is the haversine function.

####  3.1.1 Angular Power Spectrum

The sky positions of the tracers $X$ can be used to define a number density
field $n_{X}(\theta,\phi)$ in the sky, such that

$\int_{\rm Allsky}n_{X}(\theta,\phi)\sin\theta d\theta d\phi=N_{X}$ (4)

This can be done numerically, for example, by dividing the sky into equal-area
pixels using some pixelation scheme, counting the number of tracers in each
pixel, and dividing the number counts by the area of each pixel. Let the
average tracer number density in the sky be $\bar{n}_{X}=\frac{N_{X}}{4\pi}$.
The overdensity field is given by

$\delta_{X}(\theta,\phi)=\frac{n_{X}(\theta,\phi)}{\bar{n}_{X}}-1$ (5)

The overdensity field contains all the information about the auto-clustering
of the tracer set $A$. By expanding $\delta_{X}(\theta,\phi)$ in terms of
spherical harmonics $Y_{\ell m}(\theta,\phi)$, one can derive the angular
power spectrum $\mathcal{C}^{X,X}_{\ell}$, a widely used two-point summary
statistic for clustering:

$\displaystyle\delta_{X}(\theta,\phi)=\sum_{\ell m}\alpha^{X}_{\ell m}Y_{\ell
m}(\theta,\phi)$ (6)
$\displaystyle\mathcal{C}^{X,X}_{\ell}=\frac{1}{2\ell+1}\sum_{m=-\ell}^{\ell}\lvert\alpha_{\ell
m}^{X}\rvert^{2}$ (7)

where $\ell$ goes from 0 to $\infty$ and $m$ takes values from $-\ell$ to
$\ell$. In practice, the summation is cut off at some $\ell_{\rm max}$
determined by the resolution of the numerical grid on which the field is
defined. The power spectrum $C^{X,X}_{\ell}$ at a particular value of $\ell$
quantifies the clustering strength at an angular scale corresponding roughly
to $\theta\approx\pi/\ell$. In this study, we use the
HEALPix111111https://healpix.sourceforge.io/ scheme (Górski et al., 2005) as
implemented in the python library
healpy121212https://healpy.readthedocs.io/en/latest/ (Zonca et al., 2019) to
compute the overdensity field on a grid of equal-area pixels in the sky. We
use the healpy’s anafast routine to compute the power spectra. We choose the
default high-$\ell$ cutoff of 3$\times$NSIDE \- 1 in our analysis.

An equivalent measure of spatial clustering of a set of tracers is the two-
point autocorrelation function, $w^{X,X}(\theta)$, which captures the excess
probability of finding two data points separated by an angular distance of
$\theta$ in the sky over finding two points drawn from a random (Poisson)
distribution in the sky. Mathematically,

$w^{X,X}(\theta)=\bigg{\langle}\delta_{X}(\boldsymbol{\hat{\Omega}}_{1})\delta_{X}(\boldsymbol{\hat{\Omega}}_{2})\bigg{\rangle}_{\boldsymbol{\hat{\Omega}}_{1}\cdot\boldsymbol{\hat{\Omega}}_{2}=\cos\theta}$
(8)

where $\boldsymbol{\hat{\Omega}}_{1}$ and $\boldsymbol{\hat{\Omega}}_{2}$ are
unit vectors separated by an angular distance $\theta$ in the sky such that
$\boldsymbol{\hat{\Omega}}_{1}\cdot\boldsymbol{\hat{\Omega}}_{2}=\cos\theta$,
and the angular brackets denote an average over all such configurations of
$\boldsymbol{\hat{\Omega}}_{1}$ and $\boldsymbol{\hat{\Omega}}_{2}$. It can be
shown that the angular power spectrum is related to the two-point function of
the tracers in the following way:

$w^{X,X}(\theta)=\frac{1}{4\pi}\sum_{\ell}\left(1+2\ell\right)\mathcal{C}^{X,X}_{\ell}P_{\ell}\left(\cos\theta\right)$
(9)

where $P_{\ell}\left(\cos\theta\right)$ denotes the Legendre polynomial of
order $\ell$ and argument $\cos{\theta}$. Therefore, the two-point function
and the angular power spectrum encode the same physical information about the
clustering of the tracer set $X$; the angular power spectrum is a two-point
clustering statistic.

####  3.1.2 Nearest-neighbour distributions: $k$NN-CDFs

The nearest-neighbour distributions as a measure of spatial clustering in 3D
were introduced in Banerjee & Abel (2021a). Here, we briefly summarise the
idea behind these statistics and extend the mathematical formalism to 2D
clustering in the sky using angular coordinates.

The key idea that motivates the nearest-neighbour clustering framework is as
follows: all the physical information about the clustering of a set of
discrete tracers is contained in the distribution of their number counts,
i.e., the number of tracers enclosed inside a randomly chosen spatial region
of a given spatial extent. The spatial regions can have an arbitrary
geometrical shape as long as there is a way of assigning a spatial extent to
them. Since we are concerned with tracers in the sky, which is represented by
the surface of a $3$-sphere, we choose to work with spherical caps of area
$A=2\pi(1-\cos\theta)$ to study clustering at an angular scale
$\theta$131313Henceforth, we use $\theta$ and $A$ interchangeably..

Suppose we are given the positions of a set $X$ of discrete tracers. Given the
discussion above, the fundamental quantity that quantifies their clustering at
spatial scale $\theta$ is the probability $\mathcal{P}_{k|A}$ of finding $k$
data points of $X$ in a randomly placed spherical cap of area $A$ in the sky.
$\mathcal{P}_{k|A}$ can be written in terms of a generating function $P(z|A)$
as

$\mathcal{P}_{k|A}=\frac{1}{k!}\left[\left(\frac{d}{dz}\right)^{k}P(z|A)\right]_{z=0}$
(10)

For the case of spherical caps, it can be shown that the generating function
is given by (Banerjee & Abel, 2021a)141414Banerjee & Abel (2021a) derived this
expression for 3D clustering in cartesian coordinates. The argument is similar
for 2D clustering in angular coordinates.

$\displaystyle P(z|A)=$
$\displaystyle\exp\bigg{[}\sum_{k=1}^{\infty}\frac{\bar{n}_{X}^{k}\left(z-1\right)^{k}}{k!}$
$\displaystyle\times\int_{A}...\int_{A}d\boldsymbol{\hat{\Omega}}_{1}...d\boldsymbol{\hat{\Omega}}_{k}\omega^{(k)}\left(\boldsymbol{\hat{\Omega}}_{1},...\boldsymbol{\hat{\Omega}}_{k}\right)\bigg{]}$
(11)

where $\omega^{(N)}$ are the $N$-point correlation functions of the underlying
field of the tracers $X$, with $\omega^{(0)}=0$ and $\omega^{(1)}=1$ by
definition151515Note that the $N$-point functions are defined analogously to
the two-point auto-correlation function defined earlier. In fact,
$\omega^{(2)}\equiv w$. However, for $N>2$, $\omega^{(N)}$ are defined as
functions of $N$ unit vectors
$\left\\{\boldsymbol{\hat{\Omega}}_{1},...,\boldsymbol{\hat{\Omega}}_{N}\right\\}$
instead of a single angular separation $\theta$, and the average is performed
over all possible configurations preserving the polyhedron formed by the $N$
vectors.. Equation 11 shows the connection between the number count
distribution and the correlation functions usually employed to measure
clustering.

An equivalent measure of clustering is the cumulative distribution of the
tracer number counts, which represents the probability $\mathcal{P}_{>k|A}$ of
having more than $k$ tracers in a randomly chosen spherical cap of $A$.
$\mathcal{P}_{>k|A}$ can also be expressed in terms of a generating function
$C(z|A)$ given by (Banerjee & Abel, 2021a)

$\displaystyle C(z|A)=\frac{1-P(z|A)}{1-z}$ (12)

such that

$\mathcal{P}_{>k|A}=\frac{1}{k!}\left[\left(\frac{d}{dz}\right)^{k}C(z|A)\right]_{z=0}$
(13)

From equations 11 to 13, it is not clear how to compute these distributions
without first computing all higher-order correlation functions. Following
Banerjee & Abel (2021a), we now discuss another interpretation of the
cumulative count distributions that allows us to compute $\mathcal{P}_{>k|A}$
directly from the positions of the tracers. The count distributions
$\mathcal{P}_{k|A}$ can then be calculated trivially using
$\mathcal{P}_{k|A}=\mathcal{P}_{>k-1|A}-\mathcal{P}_{>k|A}$.

Consider a set of $N_{r}$ area-filling, randomly distributed query points in
the sky, such that $N_{r}\gg N_{X}$. Each query point will have a data point
in $X$ that is nearest to it, a data point that is second-nearest to it and so
on. The distributions of the distances to these neighbouring data points, over
all query points in the sky, are directly connected to the count distributions
discussed above. We argue that the cumulative distribution function (CDF) of
the distances from the query points to the $k$-nearest-neighbour data point,
or $k$NN-CDF, is precisely equal to $\mathcal{P}_{>k-1|A}$.

To understand this connection, let us examine the case of $k=1$. Consider
$N_{r}$ spherical caps of area $A=2\pi\left(1-\cos\theta\right)$, the centres
of which are distributed randomly in the sky. The fraction of such caps
enclosing at least 1 data point is equal to that of cap centres with the
angular distance to their nearest neighbour less than $\theta$. The nearest-
neighbour CDF at angular scale $\theta$ is the precise measure of the fraction
of query points (equivalent to centres of the spherical caps) for which the
nearest data point is at a distance less than $\theta$. This argument can be
easily generalised if we consider the $k$-nearest-neighbour instead of the
first nearest-neighbour. Therefore, we conclude

$\mathcal{P}_{>k-1|A}=\text{CDF}_{k\rm NN}(\theta)$ (14)

In practice, the $k$NN-CDFs are simple to calculate in a computationally
efficient manner. We start by creating a HEALPix grid of query points with a
sufficiently high value of NSIDE, such that the resolution of the query grid
is much finer than the smallest angular scale at which we want to study
spatial clustering. As noted in Banerjee & Abel (2021a), placing query points
on a finely spaced grid gives the same results as randomly distributed query
points, as long as the grid separation is much smaller than the mean
interparticle separation of the data. Next, we compute the distances to the
$k$-nearest-neighbour data point of each query point. The nearest-neighbour
search is carried out very efficiently by constructing a Ball tree structure
(Omohundro, 2009) on the data points. Once a tree is built, it can be used to
calculate the distances to the first $k$ neighbouring data points for all
query points simultaneously. Sorting the computed distances for each neighbour
index $k$ immediately gives the empirical CDF of the $k$-nearest-neighbour
distances over a range of spatial scales. The empirical CDF converges to the
true $k$NN-CDF in the limit of a large number of query points. In this study,
we utilise the sklearn.neighbors.BallTree routine from the library scikit-
learn161616https://scikit-learn.org/ (Pedregosa et al., 2012) with the
haversine distance metric for our purposes.

It is evident from equation 11 that the count distributions
$\mathcal{P}_{>k|A}$, and hence the $k$NN-CDFs, are formally sensitive to
integrals of all $N$-point correlation functions of the underlying tracer
field. This makes these summary statistics extremely powerful probes of
clustering on small spatial scales where the higher-order correlation
functions contribute significantly. We refer the interested reader to Banerjee
& Abel (2021a) for a detailed study of the gain in clustering measurements as
well as cosmological constraints achieved using the $k$NN-CDFs over two-point
clustering statistics.

###  3.2 Cross-clustering

There are two possible approaches to quantify the cross-clustering between the
BBHs and the WSC catalogue:

  1. 1.

directly cross-correlate the WSC source positions with the BBH positions

  2. 2.

cross-correlate fluctuations in the WSC source number density field and the
BBH positions.

The number of sources in the WSC catalogue is a few orders of magnitude larger
than the number of BBH events under consideration, making a tracer-tracer
cross-correlation difficult to compute using the $k$NN formalism. Since the
source number density is large enough, on the scales that we will consider in
the analysis, the WSC catalogue can be effectively treated as a ‘continuous’
field. Furthermore, as shown in figure 7, the Poisson sampling noise, or shot
noise, is subdominant to the angular power spectrum at all scales of interest.
This means that the underlying galaxy density field is well-sampled by WSC
source positions. Therefore we adopt approach (ii) in this paper and compute
the BBH-Galaxy cross-clustering using a tracer-field cross-correlation
formalism.

Figure 7: The angular power spectrum for the WSC sources computed using
healpy’s anafast routine, with the shot (Poisson sampling) noise plotted for
reference. The power spectrum is well above the shot noise up to an $\ell_{\rm
max}\mathchar 21016\relax\,10^{3}$, corresponding to spatial scales much
smaller than those considered in this analysis. Therefore, the treatment of
the WSC sources as a continuous field is a reasonable approximation.

We now describe the summary statistics for measuring the spatial cross-
correlation between a set of discrete tracers, $X$, and a continuous field
$\delta_{Y}$. We only consider continuous fields in the form of dimensionless
fluctuations in a quantity, i.e., fields that are bounded below by -1 and
average to 0. The continuous field could be an overdensity field derived from
another sample of tracers, such as the galaxies relevant for this study, or an
actual continuum field like the CMB temperature fluctuations.

####  3.2.1 Cross Angular Power Spectrum

As discussed in section 3.1.1, we can define an overdensity field $\delta_{X}$
for the discrete tracers $X$, and expand both $\delta_{X}$ and $\delta_{Y}$ in
spherical harmonics

$\displaystyle\delta_{X}(\theta,\phi)$ $\displaystyle=\sum_{\ell
m}\alpha^{X}_{\ell m}Y_{\ell m}(\theta,\phi)$
$\displaystyle\delta_{Y}(\theta,\phi)$ $\displaystyle=\sum_{\ell
m}\alpha^{Y}_{\ell m}Y_{\ell m}(\theta,\phi)$

The cross angular power spectrum between $X$ and $Y$ is defined as

$\mathcal{C}^{X,Y}_{\ell}=\frac{1}{2\ell+1}\sum_{m=-\ell}^{\ell}\left\\{\alpha^{X}_{\ell
m}\right\\}^{*}\alpha^{Y}_{\ell m}$ (15)

We compute the cross angular power spectrum in a similar manner to the auto
angular power spectrum, using healpy’s anafast routine.

####  3.2.2 Nearest-neighbour distributions

Banerjee & Abel (2023) generalised the $k$NN formalism to study tracer-field
cross-correlations in 3D using the nearest-neighbour distributions. Here, we
summarise the main points in the context of 2D angular clustering.

To perform nearest-neighbour measurements on a continuous field, we must
define the continuum version of the $k$NN-CDFs. In other words, we must
investigate the behaviour of the $k$NN-CDFs of a set of discrete tracers of
the field $\delta_{Y}$ as their average number density $\bar{n}_{Y}$ tends to
infinity. If the tracers represent a local Poisson process on the field
$\delta_{Y}$, Banerjee & Abel (2023) showed that the continuum version of the
$k$NN measurements at a spatial scale $\theta$ are thresholded evaluations of
the CDF of the smoothed continuous field (see also Anbajagane et al. (2023)
for an application of the CDF formalism to the lensing convergence fields from
the first three years of the Dark Energy Survey, where a similar approach was
taken). Mathematically,

$\mathcal{P}_{\geq
k|A}\to\mathcal{P}_{>\delta^{*}}(\theta)=\int_{\delta^{*}}^{\infty}\phi(\delta^{\theta}_{Y})d\delta^{\theta}_{Y}=1-\text{CDF}(\delta^{*})$
(16)

where $\delta^{\theta}_{Y}$ represents the continuous field smoothed on the
scale $\theta$, $\phi(\delta^{\theta}_{Y})$ represents the probability density
function (PDF) of the smoothed field. The nearest-neighbour index $k$ for
discrete data maps to a threshold density $\delta^{*}$ on the smoothed
continuous field, such that
$\bar{n}_{Y}A\left(1+\delta^{*}\right)=k$171717Note that at fixed $k$,
$\delta^{*}$ is also a function of $\theta$, as evident from its definition.

Now that we have a continuum version of the $k$NN-CDFs, we discuss the
characterisation of spatial cross-correlations using the $k$NN formalism.
Following Banerjee & Abel (2023), we take the joint probability
$\mathcal{P}_{\geq k,>\delta^{*}}(\theta)$ of finding at least $k$ tracers and
the smoothed continuous field $\delta^{\theta}_{Y}$ to cross threshold
$\delta^{*}$ in spherical caps of angular radius $\theta$ as a measure of the
spatial cross-correlations between tracers and a continuous field. Assuming
that the tracers $X$ represent a local Poisson process on the overdensity
field $\delta_{X}$, we have, similar to Banerjee & Abel (2023),

$\mathcal{P}_{k,>\delta^{*}}(\theta)=\int_{\delta^{*}}^{\infty}\frac{\left[\lambda(\delta^{\theta}_{X})\right]^{k}}{k!}e^{-\lambda(\delta^{\theta}_{X})}\phi(\delta^{\theta}_{X},\delta^{\theta}_{Y})d\delta^{\theta}_{X}d\delta^{\theta}_{Y}$
(17)

where $\phi(\delta^{\theta}_{X},\delta^{\theta}_{Y})$ is the joint probability
distribution of the two fields when smoothed on angular scale $\theta$. The
quantity of our interest, $\mathcal{P}_{\geq k,>\delta^{*}}$, can be written
in terms of $\mathcal{P}_{k,>\delta^{*}}$ as

$\mathcal{P}_{\geq
k,>\delta^{*}}=\mathcal{P}_{>\delta^{*}}-\sum_{j<k}\mathcal{P}_{j,>\delta^{*}}$
(18)

Suppose the tracers $X$ are completely uncorrelated and statistically
independent of the continuous field $\delta_{Y}$. In that case, the joint
distribution function can be factored into a product of the individual PDFs of
the smoothed fields, i.e.,
$\phi(\delta^{\theta}_{X},\delta^{\theta}_{Y})\propto\phi(\delta^{\theta}_{X})$
$\phi(\delta^{\theta}_{Y})$. In this case, we have (Banerjee & Abel, 2023)
$\mathcal{P}_{\geq k,>\delta^{*}}=\mathcal{P}_{\geq
k}\times\mathcal{P}_{>\delta^{*}}$. This fact provides a convenient way to
define a summary statistic that measures the excess cross-correlation between
the tracers and the continuous field

$\psi_{k,\delta*}\triangleq\mathcal{P}_{\geq
k,>\delta^{*}}/\left(\mathcal{P}_{\geq
k}\times\mathcal{P}_{>\delta^{*}}\right)$ (19)

This is a useful quantity to measure, as a positive (negative) measurement for
$\psi_{k,\delta*}-1$ would indicate that the tracer $X$ is correlated (anti-
correlated) with the field $\delta_{Y}$, while $\psi_{k,\delta*}-1=0$ would
indicate that the there is no spatial cross-correlation between them.

All the physical information about the spatial cross-correlation between
fluctuations in the sky distribution of tracer $X$ and the field $Y$ is
contained in the joint distribution
$\phi(\delta^{\theta}_{X},\delta^{\theta}_{Y})$. Therefore, it is clear from
equations 17 to 19 that the excess cross-correlation, as defined using the
nearest-neighbour distributions, will be sensitive not just to the linear or
Gaussian correlations in the density fluctuations of the tracers and the
continuous field, but to correlations in fluctuations at all orders (See
Banerjee & Abel (2023) for a detailed demonstration). Therefore, the $k$NN
formalism provides a powerful way to characterise cosmological cross-
correlations.

The joint probability distributions $\mathcal{P}_{\geq k,>\delta^{*}}$ and
excess cross-correlations $\psi_{k,\delta*}$ are simple to compute
numerically. We follow the procedure laid out in Banerjee & Abel (2023)

  1. 1.

Create a set of area-filling query points by creating a finely-spaced HEALPix
grid in the sky, such that the number of pixels $N_{\rm{pix}}$ is far greater
than the number of data points.

  2. 2.

Build a Ball tree from the set of tracer positions and estimate the angular
distances to the $k$-nearest neighbour data points from each query point. For
each $k$, sort the distances to produce the empirical $k$NN-CDF over a range
of angular scales $\theta$. In the limit of large $N_{\rm{pix}}$, the
empirical CDF approaches $\mathcal{P}_{\geq k}$.

  3. 3.

Smooth the continuous field $\delta_{Y}$ on an angular scale $\theta$ using a
top-hat filter. The smoothing is done in harmonic space using the
$\\{\alpha^{Y}_{\ell m}\\}$ of the field, computed via spherical harmonic
transforms, to speed up the computation time (see appendix D for details).
Interpolate the smoothed field on the query grid defined in step (i).

  4. 4.

For a given $k$ and threshold $\delta^{*}$, compute the fraction of query
points for which the $k^{\rm th}$ nearest-neighbour lies at an angular
distance less than $\theta$ and the smoothed field, interpolated to that grid
point, exceeds $\delta^{*}$. In the limit of large $N_{\rm{pix}}$, this
fraction approaches $\mathcal{P}_{\geq k,>\delta^{*}}$.

  5. 5.

Compute the fraction of query points for which the smoothed field,
interpolated to the grid point, exceeds $\delta^{*}$. In the limit of large
$N_{\rm{pix}}$, this fraction approaches $\mathcal{P}_{>\delta^{*}}$.

  6. 6.

From the quantities calculated above, compute the excess cross-correlation
using equation 19.

  7. 7.

Repeat steps (iii) to (vi) for different values of the angular scale $\theta$.

In this study, we choose the constant percentile threshold described in
Banerjee & Abel (2023) to define the threshold value $\delta^{*}$ for the
continuous field. Specifically, we set $\delta^{*}=\delta_{75}$, the value of
the $75^{\rm th}$ percentile of $\delta^{\theta}_{Y}$. This choice implies
that $\mathcal{P}_{>\delta^{*}}=0.25$ irrespective of the smoothing scale
$\theta$.

###  3.3 Example application of the tracer-field formalism

As discussed in the previous section, the nearest-neighbour distributions are
sensitive to all higher-order cross-correlation functions of the discrete
tracers and the continuous field. Consequently, nearest-neighbour measurements
are expected to measure a stronger clustering signal than the two-point
summary statistics at small angular scales where the underlying fields are
non-Gaussian, and these higher-order correlation functions cannot be
neglected. In this section, we demonstrate the gains in clustering power
obtained using the $k$NN tracer-field formalism compared to the angular power
spectrum using an illustrative example.

We create a set of discrete tracers by picking the centres of $36$ pixels in
the WSC footprint181818This is chosen to match the average number of mock BBHs
in the WSC footprint, see figure 9 for details. that have the highest number
density of sources in the sky, and compute their spatial cross-correlations
with the full WSC overdensity field using both the nearest-neighbour
measurements and the angular power spectrum. We further assume that we know
the locations of these tracers perfectly. To estimate the cosmic variance
associated with the clustering measurements, we also compute these clustering
statistics for 100 realisations of $36$ randomly chosen points in the sky. We
compute the cross-correlations on angular distance scales from $\mathchar
21016\relax\,1^{\circ}$ to $\mathchar 21016\relax\,35^{\circ}$, equivalent to
$\ell=6$ to $\ell=180$. This choice ensures that

  1. 1.

we have sufficient sampling in the nonlinear regime, where the nearest-
neighbour distributions can capture information not accessible through two-
point statistics

  2. 2.

the measurements are not affected by the lack of sampling towards the right
tail of the auto $k$NN-CDF (see Banerjee & Abel (2021a) for a discussion)

These angular distances correspond to projected transverse distance scales of
$\mathchar 21016\relax\,15$ to $\mathchar 21016\relax\,400$ Mpc for a median
redshift of $\mathchar 21016\relax\,0.2$ for the WSC catalogue.

The results for the first nearest-neighbour distribution are shown in the left
panel of figure 8, and those for the power spectrum are shown in the right
panel. The inset in the right panel presents a zoomed-in version of the full
subplot focusing on the smaller scales. In each plot, the solid line shows the
excess cross-correlation between the highest-density locations, the shaded
band represents the cosmic variance, and the dash-dot line represents the
expected value in the absence of cross-correlation.

Figure 8: The excess cross-correlation between 36 locations with the highest
density of WSC sources and the WSC overdensity map as measured by the first
nearest-neighbour distribution (left) and the cross angular power spectrum
(right). The solid lines represent the measurement on the tracers, the shaded
band represents the 3$\sigma$ variance in the same measurement performed on
randomly chosen points in the sky, and the dash-dot line represents the
expected value in the absence of cross-correlation. The inset on the right
panel shows a zoomed-in view of the power spectrum at smaller scales.

The excess cross-correlation, as measured by the first nearest-neighbour
distribution, lies well outside the shaded region representing the $3\sigma$
cosmic variance in the randoms at all angular scales smaller than $\mathchar
21016\relax\,4^{\circ}$. In contrast, the power spectrum fluctuates within the
shaded band, barely crossing the detection threshold on one or two bins. Thus,
figure 8 demonstrates that the nearest-neighbour measurements can capture a
statistically significant clustering signal on small, nonlinear scales for a
well-localised sample of rare, highly biased tracers, whereas the power
spectrum can not. We would like to draw the attention of the reader to the
following important implication of this finding: with a set of $<50$ well-
localised events, the nearest-neighbour measurements on small spatial scales
are statistically robust enough to investigate whether BBHs reside in highly
biased environments in the universe. It should be noted, however, that on
large scales, the nearest-neighbour distributions do not perform any better
than the power spectrum, and there is no detectable signal in either
statistic. This is because, on large scales, the galaxy density field is well
approximated by a Gaussian random field, and the higher-order correlation
functions are negligible.

##  4 Methods

In this section, we describe the procedure to quantify the clustering of the
BBHs and their spatial cross-correlations with the WSC sources. We address
this problem using a hypothesis-testing approach: We consider a null
hypothesis, which proposes that there is no statistical significance for a
clustering signal in the data, and attempt to rule it out by investigating the
likelihood of reproducing the observed data, assuming the null hypothesis to
be true. We describe our null hypothesis in section 4.1. To test the null
hypothesis, we require a control dataset consistent with its premise. We
already described the procedure for creating a catalogue of unclustered mock
BBHs in section 2.2. This mock catalogue automatically serves as a control set
to test the null hypothesis. In section 4.2, we describe how we apply the
mathematical formalism motivated in section 3 to compute clustering statistics
for the specific data considered in this study. We discuss the method to
calculate the statistical significance of the clustering signal in section
4.3. Finally, we describe the angular distance scales used in our analysis in
section 4.4.

###  4.1 Null Hypothesis

Our null hypothesis is as follows

> The BBHs currently detected by the LVK collaboration are spatially
> unclustered, distributed uniformly (isotropically) in the sky and are not
> spatially correlated with other tracers of the large-scale structure of the
> universe, such as galaxies and quasars.

Any dataset consistent with this hypothesis would not contain a statistically
significant clustering signal. In the following sections, we discuss the
process of testing the null hypothesis using the observed and mock BBH
catalogues described in section 2.

###  4.2 Application of clustering formalism to data

The formalism to study clustering developed in section 3 is applicable to
discrete tracers that can be treated as point objects in the sky, i.e.,
objects that can be localised to a single sky position $(\delta,\alpha)$.
However, as shown in section 2, due to the limited resolving power of
gravitational wave detectors, we can only assign each event with a probability
distribution over an extended region in the sky (see figure 1 for example).
Moreover, the computational procedure outlined in section 3 gives unbiased
measurements for the tracer-field cross-correlations only when the continuous
field is defined on the entire sky. However, as discussed in section 2,
certain regions in the sky do not have reliable galaxy data. Hence, we can not
define the galaxy overdensity field there. In this section, we discuss our
strategy to deal with these challenges, namely the uncertain sky localisation
of the BBHs and the complicating effects due to the presence of the WSC mask.

####  4.2.1 Strategy to deal with the uncertainty in BBH sky localisations

For each BBH, we have a probability distribution in the sky. Since we need a
single location for the BBHs, one possible approach is to assign each BBH the
most probable position in its sky localisation area, i.e., the position where
the probability distribution is maximised. However, the sky posteriors of many
events show signs of bimodality, and assigning a single representative
location to them is problematic. Furthermore, in reducing a probability
distribution to a single point, we lose a lot of information contained in the
shapes of the posteriors. For example, utilising the full sky distribution can
also allow us to characterise the measurement errors on the clustering
strength. Therefore, we adopt a different strategy in this work, which is as
follows

  1. 1.

for each BBH, draw an (RA, Dec) pair from the sky location posterior

  2. 2.

using the drawn samples as the ‘true’ locations of the events, compute the
auto-clustering statistics as defined in section 3.1

  3. 3.

repeat steps (i) and (ii) for 1000 draws from the posteriors

The average over the 1000 draws gives the estimated value of the clustering
strength, while the variance over the 1000 draws gives the estimated
measurement error due to the uncertainty in the sky localisation of the BBHs.
Therefore, our strategy naturally preserves the information present in the
full sky distribution of each BBH while giving us an estimate of the
measurement errors on the clustering strength. Vijaykumar et al. (2023b) took
a similar approach in a recent clustering analysis with forecast BBH data for
the third generation of gravitational wave detectors.

####  4.2.2 Strategy to deal with complications due to the WSC Mask

Since we do not have reliable data in regions outside the WSC mask, the galaxy
overdensity field is ill-defined there. How do we compute cross-correlations
with the BBHs in this scenario? One possible approach is to assign a
$\delta_{\rm Gal}=0$ to all pixels in the masked region since that is the
expected average value of an overdensity field. Although this would not bias
the results since the BBHs outside the mask would not contribute to the cross-
correlation signal, it would unnecessarily add to the noise budget. Instead,
we take the approach of removing the BBH events that lie outside the mask.
However, since the BBHs are not perfectly localised, we must be careful when
handling events whose sky localisation areas are partly inside the mask. We
follow the following strategy:

  1. 1.

Draw 1000 samples of (RA, Dec) pairs from the sky location posteriors of each
of the BBHs

  2. 2.

for each sample, remove the sky locations that lie outside the WSC mask

By keeping the posterior samples for events which are partly outside the mask,
our method not only preserves the information in the sky distribution of the
BBHs but also leads to more number of BBHs contributing to the analysis than
simply removing all BBHs whose localisation area intersects the mask would.
However, this process leads to different tracers in each sample. Since the
$k$NN-CDFs are highly sensitive to the number density of the tracers (see
equation 11), care needs to be taken to ensure that averaging the CDFs over
samples with different number densities does not lead to any issues. An
important check is whether the distribution of the number of events inside the
mask over the 1000 samples for the mock catalogue generated in section 2.2 is
statistically similar to that of the data. Figure 9 shows that, indeed, that
is the case, and hence any systematics that arise due to this would affect the
clustering of the observed and mock BBHs equally.

Figure 9: Distribution of the number of unmasked events over 1000 samples
drawn from the sky distributions of the observed (bold histogram) and mock
(light histogram) BBHs. The error bars on the mock histogram show the
$1\sigma$ variance over 135 realisations of the mock catalogue. There is an
excellent match between the two distributions within error bars. Therefore,
any systematics that arise due to differing number densities between different
samples would affect the clustering of the observed and mock BBHs equally.

We now have positions for the BBHs, from which we can compute the overdensity
fields needed to compute the cross angular power spectrum. Since there is no
data outside the mask, the value of the overdensity field computed there would
be artificially low (negative). To account for this, we set the BBH
overdensity fields outside the mask to zero before calculating their spherical
transforms.

To compute the nearest-neighbour measurements of the BBH-Galaxy excess cross-
correlation, we need to smooth the galaxy density field on various spatial
scales. To minimise the effects due to the mask, we set the field outside the
mask to zero before smoothing. We restrict the query points to inside the
mask. This is needed because the query points outside the mask would have
artificially large nearest-neighbour distances and would skew the measured
distributions. We also remove all query points within a certain angular
distance from the mask boundaries to ensure that the smoothed density field
interpolated at the query points is not affected by spurious contributions
from the regions outside the mask. Wang et al. (2022) followed a similar
procedure to analyse the spatial clustering of SDSS clusters using $k$NN-CDFS
in a recent study. In practice, we observe that a threshold distance of
roughly half the maximum angular scale used in the analysis leads to an
unbiased measurement of the excess cross-correlation. Note that these steps
are not taken when computing the BBH auto-clustering since we have BBH data on
the entire sky.

###  4.3 Statistical significance

In this section, we describe the way we compute the statistical significance
of the clustering measurements once the summary statistics have been computed
for the observed and mock BBHs. Consider a summary statistic as a function of
angular scale, $S(\theta)$, which could be the angular power spectrum, the
two-point function or the nearest-neighbour distribution, either as a measure
of the BBH auto-correlation or BBH-Galaxy cross-correlation. Let the scales
considered in the analysis be $\\{\theta_{1},...,\theta_{p}\\}$. We define the
data vector $D_{a}$ as the summary statistic $S$ evaluated on the observed BBH
catalogue at angular scale $\theta_{a}$. Similarly, we define a mock vector
$M^{i}_{b}$ as $S(\theta_{b})$ evaluated on the $i^{\rm th}$ realisation of
the mock BBH catalogue for each of the $n$ realisations. Note that $S(\theta)$
represents the summary statistic already averaged over 1000 samples drawn from
the sky distribution of the BBHs, as prescribed in section 4.2.1. To
characterise the noise properties of the measurement, we compute the
covariance matrix

$\boldsymbol{\Sigma}^{{}^{\prime}}_{ab}=\bigg{\langle}\left(M^{i}_{a}-\langle
M_{a}\rangle\right)\left(M^{i}_{b}-\langle M_{b}\rangle\right)\bigg{\rangle}$
(20)

where the angular brackets denote an average over the $n$ realisations of the
mock catalogue. The covariance matrix is, by definition, a $p\times p$ matrix.
The object that is relevant for the statistical calculations is the inverse of
the covariance matrix, which is multiplied by the Hartlap correction factor
(Hartlap, J. et al., 2007) to get an unbiased estimate

$\boldsymbol{\Sigma}^{-1}=\frac{n-p-2}{n-1}\left(\boldsymbol{\Sigma}^{{}^{\prime}}\right)^{-1}$
(21)

Once we have the corrected inverse covariance matrix, we characterise the
signal-to-noise for clustering using the $\chi^{2}$ statistic. For the
observed BBHs and each realisation of the mock BBHs, we define the $\chi^{2}$
value as

$\displaystyle\chi^{2}_{D}$
$\displaystyle=\big{(}\boldsymbol{D}-\langle\boldsymbol{M}\rangle\big{)}^{T}\boldsymbol{\Sigma}^{-1}\big{(}\boldsymbol{D}-\langle\boldsymbol{M}\rangle\big{)}$
(22) $\displaystyle\chi^{2}_{M^{i}}$
$\displaystyle=\big{(}\boldsymbol{M}^{i}-\langle\boldsymbol{M}\rangle\big{)}^{T}\boldsymbol{\Sigma}^{-1}\big{(}\boldsymbol{M}^{i}-\langle\boldsymbol{M}\rangle\big{)}$
(23)

where $\boldsymbol{D}\triangleq\\{D_{1},D_{2},...,D_{p}\\}$ and
$\boldsymbol{M}^{i}\triangleq\\{M^{i}_{1},M^{i}_{2},...,M^{i}_{p}\\}$. The
distribution of $\chi^{2}_{M^{i}}$ (henceforth the null distribution)
represents the signal-to-noise expected from data consistent with the null
hypothesis, and $\chi^{2}_{D}$ represents the signal-to-noise measured from
the data. A larger value for $\chi^{2}_{D}$ relative to the null distribution
implies a stronger statistical significance of the clustering signal.

From the null distribution and the measured signal-to-noise, we compute the
$p$-value, or probability of reproducing the observations assuming the null
hypothesis is true, by estimating the area enclosed under the (normalised)
null distribution curve after it crosses the measured signal-to-noise. In
practice, this can be estimated by counting the fraction of mock realisations
with $\chi^{2}_{M^{i}}>\chi^{2}_{D}$. If the signal is strong enough that none
of the mock realisations have a larger $\chi^{2}$ than the data, then one must
fit a $\chi^{2}$ distribution to the null distribution to compute the
$p$-value. The null hypothesis is ruled out if the $p$-value is smaller than a
chosen detection threshold.

###  4.4 Angular scales

Similar to section 3.3, we choose $10$ log-spaced angular bins from
$\theta\mathchar 21016\relax\,1^{\circ}$ to $\theta\mathchar
21016\relax\,35^{\circ}$ for the clustering analysis, which leads to a Hartlap
factor of $0.92$. For the angular power spectrum, we choose $10$ linearly-
spaced bins between $\ell=6$ and $\ell=180$.

For computing the overdensity fields and the query points for the nearest-
neighbour measurements, we use an NSIDE = 256 HEALPix grid with $\mathchar
21016\relax\,7.8\times 10^{5}$ pixels and an angular resolution of $\mathchar
21016\relax\,0.22^{\circ}$. As required for the nearest-neighbour analysis,
the number of query pixels is much larger than the number of data points, and
the query grid has sufficient resolution to sample the smallest spatial scales
analysed.

As discussed in section 4.2.2, we remove all query points within $20^{\circ}$
of the WSC mask boundaries for computing nearest-neighbour excess cross-
correlation to avoid any biases due to the presence of the WSC mask.

##  5 Results

In this section, we present the results of our clustering analysis, with
section 5.1 devoted to the auto and cross angular power spectrum and section
5.2 to the nearest-neighbour measurements. We discuss the implications of our
findings in section 5.3.

###  5.1 Angular power spectrum

Figure 10 shows the angular power spectrum of the BBHs. Filled circles
represent the power spectrum of the observed BBHs, and the error bars are
measurement errors, defined as 3 times the standard deviation across 1000
samples drawn from the BBH skymaps. The bold line and shaded band show the
mean angular power spectrum and $3\sigma$ variation around the mean for 135
realisations of the mock BBH catalogue. The dash-dot line represents the
Poisson sampling noise (shot noise) corresponding to the number density of the
BBH sample, which is equal to $1/(\bar{n}_{\rm BBH})$. Since the error bars
make it difficult to visualise the shaded band, we plot a zoomed-in version of
the full figure in the inset.

The mean power spectrum for the mock catalogue is very close to the shot
noise, as should be the case for a dataset consistent with the null
hypothesis. The variance over realisations of the mock catalogue gives an
estimate of the cosmic variance expected in the power spectrum if the null
hypothesis holds. It is evident from the figure that with the present number
of BBH detections, the angular power spectrum is shot noise-dominated and
cannot capture a clustering signal. As can be seen from the figure, the
measurement errors on the data are extremely large and exceed the cosmic
variance. This is a consequence of the significant uncertainties in the sky-
localisation of the BBHs.

Figure 10: Results of auto-clustering analysis conducted using the angular
power spectrum as the summary statistic, with a zoomed-in view provided in the
inset for better visibility. Filled circles represent the measured angular
power spectrum of the observed BBHs, with error bars representing variance
across 1000 samples drawn from the BBH skymaps. The bold line represents the
angular power spectrum averaged over 135 realisations of the mock BBH
catalogue. The shaded band, which shows the variance of the power spectrum
across realisations, represents the cosmic variance. The dash-dot line shows
the shot noise. All errors are displayed at the $3\sigma$ level. It is evident
from the figure that with the present number of BBH detections, the angular
power spectrum is shot noise-dominated; visually, there are no signs of a
clustering signal.

The $\chi^{2}$ significance test results for the angular power spectrum are
presented in figure 11, with the histogram representing the null distribution
and the solid vertical line representing the measured $\chi^{2}$ for the
observed BBHs. The data is consistent with the null hypothesis with a
$p$-value of 0.061, calculated by fitting a $\chi^{2}$ function to the null
distribution. Note that the relatively small $p$-value is most likely caused
by the large ($~{}3\sigma$) fluctuation at $\ell\mathchar 21016\relax\,80$. As
evident from figure 10, this fluctuation does not indicate a clustering
signal. Hence, we conclude that the angular power spectrum does not capture a
statistically significant clustering signal in the presently available BBH
data.

Figure 11: Results of the statistical significance test for the auto angular
power spectrum of the BBHs. The histogram represents the distribution of
$\chi^{2}$ values over 135 realisations of the mock catalogue, and the curve
enveloping it represents the best-fit $\chi^{2}$ distribution. The vertical
line represents the measured $\chi^{2}$ value for the data. The data is
consistent with the null hypothesis at a $p$-value of 0.061, calculated using
the CDF of the best $\chi^{2}$ fit. There is no evidence for a statistically
significant clustering signal in the present data.

Figure 12 shows the cross angular power spectrum measurements between the BBHs
and the WSC sources, with the plotting scheme identical to figure 10. Again,
we plot a zoomed-in version in the inset for better visibility of the shaded
band. The mean power spectrum for the mock catalogue is very close to 0 at all
scales, as expected from the null hypothesis, which stipulates that the BBHs
and WSC sources are spatially uncorrelated (the cross power spectrum is
unaffected by shot noise). Even for cross-correlation, the measurement errors
on the data are extremely large and exceed the cosmic variance. The $\chi^{2}$
significance test for the cross angular power spectrum, summarised in figure
13, confirms the visual impression given by figure 12: the data is consistent
with the null hypothesis at a $p$-value of 0.572, implying that the cross
angular power spectrum does not capture statistically significant evidence for
spatial cross-correlation between the presently observed BBHs and the WSC
sources.

Figure 12: Results of BBH-Galaxy cross-clustering analysis conducted using the
cross angular power spectrum as the summary statistic, with a zoomed-in view
provided in the inset for better visibility. The plotting scheme is identical
to figure 10, except the dash-dot line, $y=0$, represents the cross power
spectrum of two uncorrelated data sets. With the present number of BBH
detections, no visual evidence exists for a clustering signal in the cross
angular power spectrum. Figure 13: Results of the statistical significance
test for the BBH-Galaxy cross angular power spectrum. The plotting scheme is
the same as figure 11. The data is consistent with the null hypothesis at a
$p$-value of 0.572, calculated using the CDF of the best $\chi^{2}$ fit. No
statistical evidence exists for spatial cross-correlation between the
presently detected BBHs and the galaxies and quasars from the WSC catalogue.

###  5.2 Nearest-neighbour measurements

Figure 14 shows the first two $k$NN-CDFs of the BBHs in blue and orange
standing for the first and second neighbours, respectively. The plotting
scheme is similar to figure 10, except the dash-dot line represents the
expectation for the $k$NN-CDFs of an unclustered, Poisson-distributed dataset
in the sky. The analytic expression for the $k$NN-CDFs of Poisson distributed
points is only a function of $\bar{n}_{\rm BBH}A$, as can be seen from
equation 11 (see also Banerjee & Abel (2021a)).

The left panel shows that $\rm CDF_{\rm 1NN}$ is smaller than $\rm CDF_{\rm
2NN}$ at all scales. This is intuitive since, for each query point, the
distance to the $k^{\rm th}$ nearest-neighbour is always smaller than the
distance to the $(k+1)^{\rm th}$ nearest-neighbour. As a result, for a given
scale $\theta$, the fraction of query points with the first nearest neighbour
at a distance less than $\theta$ would be larger than those with the second
nearest neighbour. This generally holds, i.e., the $k$NN-CDFs always shift
towards larger angular scales with increasing $k$.

Since the CDFs span a large range relative to the error bars, it is difficult
to visualise the results from the left panel. We plot the CDFs normalised by
the mean of the mock catalogue in the right panel of figure 14. It is clear
from these plots that the mean of the mock catalogue is consistent with the
expected value for Poisson distributed data. As a consequence of the
uncertainties in the sky localisation of the BBHs, the measurement errors on
the CDFs of the observed BBHs are large, even exceeding the variance across
the mock catalogue, similar to what was observed for the angular power
spectrum in figure 10.

Note how the variance across the mock realisations for the $\rm{CDF}_{\rm
1NN}$ becomes vanishingly small as we approach the smallest angular scales.
This happens because at spatial scales much smaller than the mean tracer-
tracer separation, the $\rm{CDF}_{\rm 1NN}$ approaches the expected value for
a Poisson distribution regardless of the positions of the tracers.
Mathematically, as the area $A$ goes to 0, the leading order term in the
summation inside the exponential in equation 11 dominates. Hence, any
deviations from the Poisson expression, either due to a clustering signal or
due to sampling noise, are exponentially suppressed (see also Banerjee & Abel,
2021a).

Figure 14: Left panel: Results of auto-clustering analysis conducted using the
first and second nearest-neighbour cumulative distribution functions as the
summary statistics. The plotting scheme is similar to figure 10, except the
dash-dot line represents the analytic expectation for the $k$NN-CDFs of an
unclustered, Poisson-distributed dataset in the sky. Different colours
represent different values of the nearest-neighbour index $k$, with blue and
orange standing for $k=1$ and $2$, respectively. Right panel: CDFs divided by
the mean over the mock catalogue to reduce the dynamic range of the plot. The
error bars for the bottom plot have not been shown to the full extent to make
the rest of the plot clearer. With the present number of BBH detections, there
is no visual evidence for a clustering signal in the $k$NN-CDFs.

The data seems to indicate no clustering signal in any of the $k$NN-CDFs. To
confirm the visual intuition, we perform a $\chi^{2}$ significance test using
measurements at 5 spatial bins each from the first and second nearest-
neighbour distributions and find that the data is consistent with the null
hypothesis at a $p$-value of 0.081. Note that in this case, since the null
distribution is heavy-tailed and poorly characterised by a $\chi^{2}$
function, we compute the $p$-value by counting the number of mock realisations
with $\chi^{2}_{M^{i}}>\chi^{2}_{D}$. These results are summarised in figure
15. We conclude that the nearest-neighbour distributions do not capture a
statistically significant clustering signal in the presently available BBH
data.

Figure 15: Results of the statistical significance test for the combined first
and second nearest-neighbour cumulative distribution functions. The plotting
scheme is the same as in figure 11. The data is consistent with the null
hypothesis at a $p$-value of 0.081, calculated by counting the number of mock
realisations with $\chi^{2}_{M^{i}}>\chi^{2}_{D}$ since the best chi-square
fit does not characterise the right tail of the distribution well. There is no
evidence for a statistically significant clustering signal in the present
data.

Figure 16 shows the excess cross-correlation between the BBHs and the quasars
and galaxies from the WSC catalogue, as measured by the first and second
nearest-neighbour measurements. The plotting and colour schemes are identical
to figure 14, except the dash-dot line represents the expected excess cross-
correlation between two spatially uncorrelated datasets and is identically
equal to 1 at all scales. The inset on the right shows a zoomed-in version of
the full subplot for better visibility. The mean of the mock catalogue is
consistent with 1, as expected from the null hypothesis. As was the case for
the other summary statistics, the measurement errors on the data are extremely
large even for the excess cross-correlation191919Note that the excess cross-
correlation, by definition, cannot take negative values; the error bars
extending to negative values is an effect of the choice to plot the mean $\pm$
3 $\times$ standard deviation, but the actual measurements are always
positive..

Interestingly, the nearest-neighbour distributions indicate a mild anti-
correlation between the BBHs and the WSC catalogue at all angular scales
considered, which is not picked up by the cross angular power spectrum. An
anti-correlation between BBHs and large-scale structure would imply that these
binary black holes primarily merge in highly isolated, pristine environments
such as cosmic voids or in small field galaxies that are not part of larger
gravitationally bound groups or clusters. However, even though the anti-
correlation seems to manifest in both the nearest-neighbours systematically,
extreme care needs to be taken while analysing such plots of the nearest-
neighbour distributions; since nearest-neighbour measurements are cumulative,
a noise-driven fluctuation on one spatial scale can affect the measurement at
nearby scales, and our visual intuition can not be trusted. We need to take
this into account by calculating the full covariance matrix before reaching
any conclusions. Moreover, the deviation from $\psi=1$ in each case is well
within the limits of cosmic variance as characterised by the mock catalogue.
This is likely an effect of sample variance due to the small number of
observed BBHs, and we will return to this point in section 5.3.

Figure 16: The excess BBH-Galaxy cross-correlation measured by the first
(left) and second (right) nearest-neighbour measurements. The plotting scheme
is identical to the right panel of figure 14, except the dash-dot line here
represents the expected excess cross-correlation between two spatially
uncorrelated datasets and is identically equal to 1. Since the error bars for
the right panel make visualising the rest of the figure difficult, a zoomed-in
view is provided in the inset. The plots visually indicate a mild anti-
correlation between the BBHs and the WSC catalogue at all angular scales.

We perform a $\chi^{2}$ significance test using 5 spatial bins each from the
first and second nearest-neighbour distributions and find that the data is
consistent with the null hypothesis at a $p$-value of 0.444. Similar to the
auto-CDFs, the null distribution here is also heavy-tailed. Hence, we compute
the $p$-value by counting the number of mock realisations with
$\chi^{2}_{M^{i}}>\chi^{2}_{D}$. Figure 17 shows the summary plot of this
analysis. We conclude that the nearest-neighbour measurements do not capture
statistically significant evidence for spatial cross-correlation between the
presently observed BBHs and the WSC sources.

Figure 17: Results of the statistical significance test for the excess BBH-
Galaxy cross-correlation as measured using a combination of the first and
second nearest-neighbour measurements. The plotting scheme is the same as
figure 11. The data is consistent with the null hypothesis at a $p$-value of
0.444, calculated by counting the number of mock realisations with
$\chi^{2}_{M^{i}}>\chi^{2}_{D}$ since the best chi-square fit does not
characterise the right tail of the distribution well. This implies that
despite a visual indication for an anti-correlation, there is no statistical
evidence for any spatial cross-correlation between the presently detected BBHs
and the galaxies and quasars from the WSC catalogue.

###  5.3 Discussion

As we saw in section 5, none of the summary statistics considered in this
study were able to capture a statistically significant signal, either for the
auto-clustering of BBHs or for spatial cross-correlations of BBHs with the
large-scale structure of the universe, in the presently available data. Our
results are consistent with previous attempts in the literature (see, for
example, Cavaglià & Modi, 2020; Mukherjee et al., 2022; Zheng et al., 2023).
What explains these results? If BBHs reside primarily in galaxies, where most
stars in the universe live and die, their locations are expected to be
clustered and spatially cross-correlated with the observed fluctuations in
large-scale structure surveys.

We believe two aspects of the data conspire to obscure the clustering signal:
first, the observed BBHs constitute a statistically small sample that is
susceptible to sample variance; with only $\mathchar 21016\relax\,50$ data
points, it is challenging to tell apart a clustered sample from a Poisson-
distributed one. Second, with the current sensitivities of the gravitational
wave detectors, there is considerable uncertainty in the sky localisation of
the BBHs, which tends to smear out any clustering signal at small scales that
the nearest-neighbour distributions are the most sensitive to. Of course,
there is always the possibility that the null hypothesis is true, in which
case, we would not see a clustering signal even if we had perfect observations
and a statistically large sample of BBHs.

The detection of a clustering signal in a small sample such as the one
selected for this study would imply that binary black holes reside in
extremely biased environments, such as highly dense nodes of the cosmic web or
huge cosmic voids. As discussed in section 3.3, the nearest-neighbour
distributions are statistically powerful enough to detect the clustering of
rare and highly biased tracers at nonlinear scales. Unfortunately, the non-
detection of a clustering signal in the current BBH data does not rule out the
possibility of BBHs being highly biased tracers because the uncertainty in the
sky localisations would completely wash it out even if a signal existed.
Repeating this analysis with better-localised events from future gravitational
wave observations would be worthwhile.

We briefly discussed in section 5 that the cross-correlation measurements
indicate a mild anti-correlation between the observed BBHs and the WSC
sources, albeit statistically insignificant. To investigate this further, we
plot, in figure 18, the most probable sky locations of the observed BBHs on
top of the galaxy overdensity field smoothed on $10^{\circ}$ scale using a
top-hat filter. Many of the observed BBHs appear to lie near large-scale
underdense regions. Due to the small sample size, this is picked up as a mild
anti-correlation in the nearest-neighbour distributions. However, this is not
evidence for anti-correlation, as established earlier using the $\chi^{2}$
test. We have further checked that a significant number of the mock
realisations show similar behaviour, which is most likely a result of sample
variance.

Figure 18: Mollweide projection of the overdensity field for the WSC catalogue
smoothed on a $10^{\circ}$ scale using a top-hat filter, with the superimposed
white dots representing the most probable positions of the observed BBHs.
Warmer colours represent a higher density of galaxies and quasars, while
cooler colours represent underdensities. Many of the observed points lie near
large-scale underdensities, leading to a slight anti-correlation in the cross-
clustering measurements. We believe that this may be due to sample variance
since the anti-correlation is not statistically significant, and a significant
number of the mock realisations show similar behaviour.

##  6 Conclusions and Future Directions

In this paper, we developed a framework for quantifying the spatial clustering
of sources of gravitational waves and their cross-correlation with the large-
scale structure of the universe, using the angular power spectrum and nearest-
neighbour distributions as summary statistics. We extended the $k$-nearest-
neighbour formalism, originally developed in Banerjee & Abel (2021a) and
Banerjee & Abel (2023) for 3D clustering, to angular clustering in the sky.
Our framework implements robust strategies to deal with the extended sky
localisation of sources and selection biases associated with gravitational
wave detections. It can handle observational systematics due to the presence
of masked regions in the sky with unreliable electromagnetic observations.

We illustrated the statistical power of the nearest-neighbour distributions as
measures of spatial clustering of sparsely sampled and highly biased tracers
by cross-correlating the overdensity field of the WISE$\times$SuperCOSMOS
(WSC) all-sky catalogue with $36$ tracers residing in the highest density
regions in the sky. Even with such a small sample size, the first nearest-
neighbour distribution captured a statically significant signal at small
scales where the angular power spectrum did not. This example demonstrated
that the nearest-neighbour distributions can access information in the higher-
order correlation functions at small scales where cosmological fluctuations
are non-Gaussian.

As a first application to data, we measured the angular power spectrum and
nearest-neighbour distributions of the Binary Black Hole (BBH) mergers
detected in the first three observation runs of LIGO-Virgo-KAGRA and cross-
correlated these sources with galaxies and quasars from the WSC catalogue. We
adopted a hypothesis-testing approach to determine the significance of the
clustering signal, with the null hypothesis stipulating that BBHs are
distributed uniformly in the sky. To mitigate observational biases in the BBH
data, we created a catalogue of mock BBHs that statistically reproduce the
observed properties of the detected BBHs but are spatially unclustered and
uncorrelated with the large-scale structure of the universe. This sample
served as a natural control set to compare with the data while testing the
null hypothesis.

Using chi-squared distributions to measure statistical deviations from the
null hypothesis, we found no evidence for spatial clustering of BBHs or their
cross-correlation with large-scale structure in the presently available data.
These results are consistent with similar studies in the literature (Cavaglià
& Modi, 2020; Mukherjee et al., 2022; Zheng et al., 2023). We discussed that
an absence of a clustering signal is not unexpected, given the small sample
size and considerable uncertainty in the sky localisation of the BBHs.

A detection of clustering with so few events would indicate that BBHs reside
in extremely biased environments in the universe, such as cosmic web nodes and
massive voids. However, a non-detection of this cross-correlation in currently
available data does not rule out this scenario since the sky localisation
uncertainty smears out the clustering signal at small scales where the
measurements are most sensitive. We demonstrated that with well-localised
BBHs, the $k$NN tracer-field formalism has the exciting potential to test the
possibility of BBHs being highly biased tracers of large-scale structures.

Our framework provides a powerful means to study spatial cross-correlations
between continuous fields and rare transient events with uncertain sky
localisation in the presence of selection effects and observational
systematics. We focused on binary black holes in this work. However, our
methods can also be applied to study binary neutron star mergers or neutron
star black hole mergers with minor modifications. In addition to gravitational
wave sources, other astrophysical transients, such as gamma-ray bursts, are
often poorly localised (see Michael Burgess, J. et al. (2021) and references
therein). The methods presented in this paper will be useful for conducting
multi-messenger studies with these objects. Similarly, the tracer-field
correlation formalism discussed here can be applied to conduct cross-
correlation studies between large-scale structure and cosmological fields,
such as the cosmic microwave background and the cosmological 21 cm neutral
hydrogen signal.

While we could not measure the clustering signal in the presently available
data on binary black hole mergers, with a statistically significant population
of better-localised merger events expected to be detected in future observing
runs of LIGO-Virgo-KAGRA and the third generation of gravitational wave
detectors (Hall & Evans, 2019; Iacovelli et al., 2022; Borhanian &
Sathyaprakash, 2022), we would have access to even smaller scales where the
nearest-neighbour distributions are expected to offer significant gains over a
two-point analysis (Banerjee & Abel, 2021a, b, 2023). Moreover, the larger
sample size will also allow us to take a tomographic approach to compute
cross-correlations by dividing the data in redshift bins202020Such an approach
has been shown to increase the diagnostic power of clustering statistics (see
Calore et al., 2020, for a detailed discussion).. Hence, the techniques
developed in this paper would be crucial for measuring the clustering of
gravitational wave sources that will be detected in the coming decades.

In work under preparation, we are conducting forecast studies analysing the
angular clustering of BBHs expected to be detected in future observing runs of
LIGO and the third-generation detectors, as well as their cross-correlations
with forecast galaxy data for stage-IV large-scale surveys. Our preliminary
findings suggest that a future detector network including LIGO India (Saleem
et al., 2022) will detect $\mathchar 21016\relax\,1.6\times 10^{4}$ BBHs with
$1\sigma$ sky localisation area less than $50$ sq. deg. in 10 years of
observation, while the third-generation detectors are expected to detect an
order of magnitude more number of BBHs with even better sky localisation. We
will present the results of a clustering analysis on this forecast data in an
upcoming paper (Gupta & Banerjee 2024 in prep.).

## Acknowledgements

The authors would like to thank Aditya Vijaykumar and Susmita Adhikari for
helpful comments on the manuscript. We thank Aditya Vijakumar and Prayush
Kumar for useful discussions and for their help in creating the mock BBH
catalogue. The support and the resources provided by PARAM Brahma Facility
under the National Supercomputing Mission, Government of India at the Indian
Institute of Science Education and Research; Pune are gratefully acknowledged.
The iPython notebook environment JupyterLab
212121https://jupyterlab.readthedocs.io/en/latest/index.html and the python
libraries healpy1212footnotemark: 12, NumPy222222https://numpy.org/,
pandas232323https://pandas.pydata.org/, scikit-learn1616footnotemark: 16 and
SciPy242424https://scipy.org/ were used extensively in this work. All plots in
this paper were made using Matplotlib252525https://matplotlib.org.

## Data Availability

All data used in this paper are publicly available; the gravitational wave
parameter estimation data and skymaps for the parent BBH catalogue can be
found at https://zenodo.org/records/5546663 and the WSC SVM catalogue and mask
can be found at http://ssa.roe.ac.uk/WISExSCOS.html. The data generated in
this study, including the mock BBH catalogues, are available upon reasonable
request.

## References

  * Abbott et al. (2017) Abbott B. P., et al., 2017, Nature, 551, 85
  * Abbott et al. (2021) Abbott B. P., et al., 2021, The Astrophysical Journal, 909, 218
  * Abbott et al. (2023) Abbott R., et al., 2023, The Astrophysical Journal, 949, 76
  * Accadia et al. (2012) Accadia T., et al., 2012, Journal of Instrumentation, 7, P03012
  * Ade et al. (2016) Ade P. a. R., et al., 2016, Astronomy & Astrophysics, 594, A13
  * Adhikari et al. (2020) Adhikari S., Fishbach M., Holz D. E., Wechsler R. H., Fang Z., 2020, The Astrophysical Journal, 905, 21
  * Aghanim et al. (2020) Aghanim N., et al., 2020, Astronomy & Astrophysics, 641, A6
  * Alfradique et al. (2024) Alfradique V., et al., 2024, Monthly Notices of the Royal Astronomical Society, 528, 3249
  * Amon et al. (2023) Amon A., et al., 2023, Monthly Notices of the Royal Astronomical Society, 518, 477
  * Anbajagane et al. (2023) Anbajagane D., et al., 2023, Monthly Notices of the Royal Astronomical Society, 526, 5530
  * Balaudo et al. (2023) Balaudo A., Garoffolo A., Martinelli M., Mukherjee S., Silvestri A., 2023, J. Cosmology Astropart. Phys., 2023, 050
  * Banerjee & Abel (2021a) Banerjee A., Abel T., 2021a, Monthly Notices of the Royal Astronomical Society, 500, 5479
  * Banerjee & Abel (2021b) Banerjee A., Abel T., 2021b, Monthly Notices of the Royal Astronomical Society, 504, 2911
  * Banerjee & Abel (2023) Banerjee A., Abel T., 2023, Monthly Notices of the Royal Astronomical Society, 519, 4856
  * Banerjee et al. (2022) Banerjee A., Kokron N., Abel T., 2022, Monthly Notices of the Royal Astronomical Society, 511, 2765
  * Bera et al. (2020) Bera S., Rana D., More S., Bose S., 2020, The Astrophysical Journal, 902, 79
  * Bilicki et al. (2016) Bilicki M., et al., 2016, The Astrophysical Journal Supplement Series, 225, 5
  * Borhanian & Sathyaprakash (2022) Borhanian S., Sathyaprakash B. S., 2022, Listening to the Universe with Next Generation Ground-Based Gravitational-Wave Detectors, doi:10.48550/arXiv.2202.11048, http://arxiv.org/abs/2202.11048
  * Calore et al. (2020) Calore F., Cuoco A., Regimbau T., Sachdev S., Serpico P. D., 2020, Physical Review Research, 2, 023314
  * Cavaglià & Modi (2020) Cavaglià M., Modi A., 2020, Universe, 6, 93
  * Chen et al. (2017) Chen H.-Y., Essick R., Vitale S., Holz D. E., Katsavounidis E., 2017, The Astrophysical Journal, 835, 31
  * Collister & Lahav (2004) Collister A. A., Lahav O., 2004, Publications of the Astronomical Society of the Pacific, 116, 345
  * Cutri et al. (2013) Cutri R. M., et al., 2013, Technical report, Explanatory Supplement to the AllWISE Data Release Products, https://ui.adsabs.harvard.edu/abs/2013wise.rept....1C. https://ui.adsabs.harvard.edu/abs/2013wise.rept....1C
  * DES Collaboration et al. (2022a) DES Collaboration et al., 2022a, Physical Review D, 105, 023520
  * DES Collaboration et al. (2022b) DES Collaboration et al., 2022b, Physical Review D, 106, 043520
  * DES Collaboration et al. (2022c) DES Collaboration et al., 2022c, Physical Review D, 106, 103530
  * DES Collaboration et al. (2023) DES Collaboration et al., 2023, Physical Review D, 107, 083504
  * Devaraju (2015) Devaraju B., 2015, doctoralThesis, doi:10.18419/opus-3985, http://elib.uni-stuttgart.de/handle/11682/4002
  * Dvornik et al. (2023) Dvornik A., et al., 2023, Astronomy & Astrophysics, 675, A189
  * Essick (2023) Essick R., 2023, Physical Review D, 108, 043011
  * Essick et al. (2023) Essick R., Farr W. M., Fishbach M., Holz D. E., Katsavounidis E., 2023, Phys. Rev. D, 107, 043016
  * Ezquiaga & Holz (2022) Ezquiaga J. M., Holz D. E., 2022, Phys. Rev. Lett., 129, 061102
  * Fang et al. (2020) Fang K., Banerjee A., Charles E., Omori Y., 2020, The Astrophysical Journal, 894, 112
  * Farr et al. (2019a) Farr W. M., Fishbach M., Ye J., Holz D. E., 2019a, The Astrophysical Journal Letters, 883, L42
  * Farr et al. (2019b) Farr W. M., Fishbach M., Ye J., Holz D. E., 2019b, The Astrophysical Journal Letters, 883, L42
  * Fishbach et al. (2018) Fishbach M., Holz D. E., Farr W. M., 2018, The Astrophysical Journal Letters, 863, L41
  * Fumagalli, A. et al. (2024) Fumagalli, A. Costanzi, M. Saro, A. Castro, T. Borgani, S. 2024, A&A, 682, A148
  * Gagnon et al. (2023) Gagnon E. L., Anbajagane D., Prat J., Chang C., Frieman J., 2023, arXiv e-prints, p. arXiv:2312.16289
  * Gair et al. (2023) Gair J. R., et al., 2023, The Astronomical Journal, 166, 22
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, The Astrophysical Journal, 622, 759
  * Hall & Evans (2019) Hall E. D., Evans M., 2019, Classical and Quantum Gravity, 36, 225002
  * Hambly et al. (2001) Hambly N., et al., 2001, Monthly Notices of the Royal Astronomical Society, 326, 1279
  * Hartlap, J. et al. (2007) Hartlap, J. Simon, P. Schneider, P. 2007, A&A, 464, 399
  * Holz & Hughes (2005) Holz D. E., Hughes S. A., 2005, The Astrophysical Journal, 629, 15
  * Holz et al. (2018) Holz D. E., Hughes S. A., Schutz B. F., 2018, Physics Today, 71, 34
  * Hoy & Raymond (2021) Hoy C., Raymond V., 2021, SoftwareX, 15
  * Hu & Wang (2023) Hu J.-P., Wang F.-Y., 2023, Universe, 9, 94
  * Iacovelli et al. (2022) Iacovelli F., Mancarella M., Foffa S., Maggiore M., 2022, The Astrophysical Journal, 941, 208
  * Krakowski et al. (2016) Krakowski T., Małek K., Bilicki M., Pollo A., Kurcz A., Krupa M., 2016, Astronomy & Astrophysics, 596, A39
  * LIGO Scientific Collaboration and Virgo Collaboration et al. (2016) LIGO Scientific Collaboration and Virgo Collaboration et al., 2016, Physical Review Letters, 116, 061102
  * LIGO Scientific Collaboration et al. (2015) LIGO Scientific Collaboration a. J. A., et al., 2015, Classical and Quantum Gravity, 32, 074001
  * LIGO Scientific Collaboration et al. (2023a) LIGO Scientific Collaboration and KAGRA Collaboration V. C., et al., 2023a, Physical Review X, 13, 011048
  * LIGO Scientific Collaboration et al. (2023b) LIGO Scientific Collaboration and KAGRA Collaboration V. C., et al., 2023b, Physical Review X, 13, 041039
  * Libanore et al. (2021) Libanore S., et al., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 035
  * Libanore et al. (2022) Libanore S., Artale M., Karagiannis D., Liguori M., Bartolo N., Bouffanais Y., Mapelli M., Matarrese S., 2022, Journal of Cosmology and Astroparticle Physics, 2022, 003
  * MacLeod & Hogan (2008) MacLeod C. L., Hogan C. J., 2008, Phys. Rev. D, 77, 043512
  * Mancarella et al. (2022) Mancarella M., Genoud-Prachex E., Maggiore M., 2022, Physical Review D, 105, 064030
  * Mastrogiovanni et al. (2021) Mastrogiovanni S., et al., 2021, Phys. Rev. D, 104, 062009
  * Mastrogiovanni et al. (2024) Mastrogiovanni S., Karathanasis C., Gair J., Ashton G., Rinaldi S., Huang H.-Y., Dálya G., 2024, Annalen der Physik, 536, 2200180
  * Meena & Bagla (2019) Meena A. K., Bagla J. S., 2019, Monthly Notices of the Royal Astronomical Society, 492, 1127
  * Michael Burgess, J. et al. (2021) Michael Burgess, J. Cameron, Ewan Svinkin, Dmitry Greiner, Jochen 2021, A&A, 654, A26
  * Miyatake et al. (2023) Miyatake H., et al., 2023, Physical Review D, 108, 123517
  * Mukherjee et al. (2021) Mukherjee S., Wandelt B. D., Nissanke S. M., Silvestri A., 2021, Phys. Rev. D, 103, 043520
  * Mukherjee et al. (2022) Mukherjee S., Krolewski A., Wandelt B. D., Silk J., 2022, Cross-correlating dark sirens and galaxies: measurement of $H_0$ from GWTC-3 of LIGO-Virgo-KAGRA, doi:10.48550/arXiv.2203.03643, http://arxiv.org/abs/2203.03643
  * Namikawa et al. (2016) Namikawa T., Nishizawa A., Taruya A., 2016, Phys. Rev. Lett., 116, 121302
  * Oguri (2016) Oguri M., 2016, Phys. Rev. D, 93, 083511
  * Omohundro (2009) Omohundro S. M., 2009. https://api.semanticscholar.org/CorpusID:61067117
  * Palmese et al. (2020) Palmese A., et al., 2020, The Astrophysical Journal Letters, 900, L33
  * Payne et al. (2020) Payne E., Banagiri S., Lasky P. D., Thrane E., 2020, Phys. Rev. D, 102, 102004
  * Pedregosa et al. (2012) Pedregosa F., et al., 2012, Scikit-learn: Machine Learning in Python, https://arxiv.org/abs/1201.0490v4
  * Peron et al. (2023) Peron M., Libanore S., Ravenni A., Liguori M., Artale M. C., 2023, arXiv e-prints, p. arXiv:2305.18003
  * Pratten et al. (2021) Pratten G., et al., 2021, Physical Review D, 103, 104056
  * RIOS (1795) RIOS J. d. M. Y., 1795, Memoria sobre algunos Métodos nuevos de calcular la Longitud por las distancias lunares: y aplicacion de su teórica á la solucion de otras problemas de navegacion. (Tabla, etc.).. Imp. Real
  * Raccanelli et al. (2016) Raccanelli A., Kovetz E. D., Bird S., Cholis I., Muñoz J. B., 2016, Phys. Rev. D, 94, 023516
  * Reid et al. (1991) Reid I. N., et al., 1991, Publications of the Astronomical Society of the Pacific, 103, 661
  * Riess (2020) Riess A. G., 2020, Nature Reviews Physics, 2, 10
  * Riess et al. (2018) Riess A. G., et al., 2018, The Astrophysical Journal, 855, 136
  * Riess et al. (2019) Riess A. G., Casertano S., Yuan W., Macri L. M., Scolnic D., 2019, The Astrophysical Journal, 876, 85
  * Riess et al. (2022) Riess A. G., et al., 2022, The Astrophysical journal letters, 934, L7
  * Saleem et al. (2022) Saleem M., et al., 2022, Classical and Quantum Gravity, 39, 025004
  * Scelfo et al. (2018) Scelfo G., Bellomo N., Raccanelli A., Matarrese S., Verde L., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 039
  * Scelfo et al. (2020) Scelfo G., Boco L., Lapi A., Viel M., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 045
  * Schutz (1986) Schutz B. F., 1986, Nature, 323, 310
  * Singer & Price (2016) Singer L. P., Price L. R., 2016, Physical Review D, 93, 024013
  * Soares-Santos et al. (2019) Soares-Santos M., et al., 2019, The Astrophysical Journal Letters, 876, L7
  * Talbot & Thrane (2018) Talbot C., Thrane E., 2018, The Astrophysical Journal, 856, 173
  * Trott & Huterer (2022) Trott E., Huterer D., 2022, Challenges for the statistical gravitational-wave method to measure the Hubble constant (arXiv:2112.00241)
  * Valentino et al. (2021) Valentino E. D., et al., 2021, Classical and Quantum Gravity, 38, 153001
  * Vijaykumar et al. (2023a) Vijaykumar A., Fishbach M., Adhikari S., Holz D. E., 2023a, arXiv e-prints, p. arXiv:2312.03316
  * Vijaykumar et al. (2023b) Vijaykumar A., Saketh M., Kumar S., Ajith P., Choudhury T. R., 2023b, Physical Review D, 108, 103017
  * Wang et al. (2022) Wang Y., Banerjee A., Abel T., 2022, Monthly Notices of the Royal Astronomical Society, 514, 3828
  * Wong et al. (2020) Wong K. C., et al., 2020, Monthly Notices of the Royal Astronomical Society, 498, 1420
  * Wright et al. (2010) Wright E. L., et al., 2010, The Astronomical Journal, 140, 1868
  * Zheng et al. (2023) Zheng Y., Kouvatsos N., Golomb J., Cavaglià M., Renzini A. I., Sakellariadou M., 2023, Physical Review Letters, 131, 171403
  * Zonca et al. (2019) Zonca A., Singer L. P., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K. M., 2019, Journal of Open Source Software, 4, 1298

##  Appendix A BBH properties

Figure 19 shows the distribution of component masses, chirp masses and network
matched-filtered SNRs of the observed BBH catalogue.

(a) Component mass (b) Chirp mass (c) SNR

Figure 19: Distribution of the component masses (left), chirp masses (middle)
and SNRs (right) of the observed BBH catalogue. In the left panel, the primary
(heavier) BBH mass distribution is shown in blue while the secondary (lighter)
BBH mass distribution is shown in orange. The characteristic peak at
$\mathchar 21016\relax\,30M_{\odot}$ is clearly visible.

##  Appendix B Observed vs. Mock BBHs

Figure 20 compares the distributions of the primary masses, chirp masses and
network matched-filtered SNRs of the observed and mock BBH catalogue. Within
error bars, the mock catalogue statistically reproduces the observed data.

(a) Primary mass (b) Chirp mass (c) SNR

Figure 20: A comparison of the distribution of the primary masses (left),
chirp masses (middle) and SNRs (right) of the observed and mock BBH catalogue.
The plotting scheme is the same as figure 4. Within error bars, the mock
catalogue statistically reproduces the observed data.

##  Appendix C Population Models

For creating the mock BBH catalogue, we assume the Power Law + Peak model for
the mass of the primary (heavier) BBH, a power law distribution for the ratio
of component masses, and a power law distribution for redshift evolution of
merger rate per unit comoving volume per unit source-frame time, as specified
in LIGO Scientific Collaboration et al. (2023a).

Let $m_{1}$ denote the mass of the primary black hole, and $q$ denote the mass
ratio, such that the mass of the secondary black hole is $qm_{1}$. The
population distribution model for $m_{1}$ is given by

$\displaystyle\centering\pi(m_{1}|\@add@centering$ $\displaystyle\lambda_{\rm
peak},\alpha,m_{\rm min},\delta_{m},m_{\rm max},\mu_{m},\sigma_{m})=$
$\displaystyle\bigg{[}(1-\lambda_{\rm peak})\mathcal{P}(m_{1}|-\alpha,m_{\rm
max})+\lambda_{\rm peak}G(m_{1}|\mu_{m},\sigma_{m})\bigg{]}$
$\displaystyle\times S(m_{1}|m_{\rm min},\delta_{m})$ (24)

where $\mathcal{P}(m_{1}|-\alpha,m_{\rm max})$ is a power law distribution
with spectral index $-\alpha$ and high-mass cutoff $m_{\rm max}$,
$G(m_{1}|\mu_{m},\sigma_{m})$ is a Gaussian distribution with mean $\mu_{m}$
and standard deviation $\sigma_{m}$, and $\lambda_{\rm peak}$ is the mixing
fraction that determines the relative importance of the power law and Gaussian
components. The lower mass end of the distribution is tapered using a
smoothing function $S(m_{1}|m_{\rm min},\delta_{m})$ which rises from 0 to 1
over the interval $(m_{\rm min},m_{\rm min}+\delta_{m})$, given by

$\displaystyle S(m_{1}|m_{\rm min},\delta_{m})$
$\displaystyle=\begin{cases}0&m<m_{\rm min}\\\ \left[f(m-m_{\rm
min},\delta_{m})+1\right]^{-1}&m_{\rm min}\leq m<m_{\rm min}+\delta_{m}\\\
1&m\geq m_{\rm min}+\delta_{m}\end{cases}$ (25)

with

$f(m,\delta_{m})=\exp\left(\frac{\delta_{m}}{m}+\frac{\delta_{m}}{m-\delta_{m}}\right)$
(26)

The conditional mass ratio distribution is given by a power law, also smoothed
at the lower mass end

$\pi(q|m_{1},\beta_{q},m_{\rm min},\delta_{m})\propto
q^{\beta_{q}}S(qm_{1}|m_{\rm min},\delta_{m})$ (27)

The values for the model parameters assumed in this work are taken from the
publicly available LVK population analysis results and are summarised in table
1.

Parameter | Value
---|---
$\alpha$ | 3.4
$\beta_{q}$ | 1.08
$m_{\rm min}$ | 5.08
$m_{\rm max}$ | 86.85
$\lambda_{\rm peak}$ | 0.04
$\mu_{m}$ | 33.73
$\sigma_{m}$ | 3.56
$\delta_{m}$ | 4.83
Table 1: Power Law + Peak model parameters

The power law redshift evolution model parameterises the merger rate density
per comoving volume and source-frame time as

$\mathcal{R}(z)=\frac{dN}{dV_{\rm c}dt_{\rm s}}=\mathcal{R}_{0}(1+z)^{\kappa}$
(28)

where $\mathcal{R}_{0}$ is the merger rate density at $z=0$, $t_{\rm s}$ is
the source-frame time, related to observer-frame time as $t_{\rm s}=t_{\rm
o}/(1+z)$ due to cosmological redshift. This implies that the observed
redshift distribution is

$\displaystyle\frac{dN}{dz}$ $\displaystyle=\int dt_{\rm o}\frac{dV_{\rm
c}}{dz}\mathcal{R}_{0}(1+z)^{\kappa-1}$ $\displaystyle=t_{\rm
obs}\mathcal{R}_{0}\frac{dV_{c}}{dz}(1+z)^{\kappa-1}$ (29)

where $t_{\rm obs}$ is the total observation time and $\frac{dV_{\rm c}}{dz}$
is the differential comoving volume. The probability distribution function for
redshifts is given by normalising equation 29

$\pi(z|\kappa,z_{\rm
max})=\frac{\frac{dV_{c}}{dz}(1+z)^{\kappa-1}}{\int_{0}^{z_{\rm
max}}dz\frac{dV_{c}}{dz}(1+z)^{\kappa-1}}$ (30)

where $z_{\rm max}$ is the maximum redshift out to which the population has
been created. Note that the constants $\mathcal{R}_{0}$ and $t_{\rm obs}$ drop
out of the expression since they are simply normalisation constants and do not
affect the shape of the distribution. However, they do indeed control the
total number of mock events to be drawn from the normalised distribution.
Although the $z_{\rm max}$ in LIGO analyses is typically taken to be 2.3, in
this work, we assume a $z_{\rm max}$ of 0.7 since assuming a lower value of
$z_{\rm max}$ is computationally less demanding. We have checked that
inputting a higher $z_{\rm max}$ value does not affect our mock catalogue.
This is due to the fact that a negligible fraction of injections outside this
redshift are detected by our assumed network. We assume $\kappa=3$ in our
analysis.

##  Appendix D Smoothing in harmonic space

Consider a field $\delta(\boldsymbol{\Omega})$ defined in the sky. The
expression for the field smoothed on an angular scale $\theta$ is given by

$\delta^{\theta}(\boldsymbol{\Omega})=\frac{1}{2\pi(1-\cos\theta)}\int_{\arccos(\boldsymbol{\hat{\Omega}}\cdot\boldsymbol{\hat{\Omega}}^{\prime})\leq\theta}d\boldsymbol{\hat{\Omega}}^{\prime}\delta(\boldsymbol{\hat{\Omega}}^{\prime})$
(31)

The smoothed field is equivalent to the field averaged over spherical caps of
angular radius $\theta$. Equation 31 can be re-written as

$\delta^{\theta}(\boldsymbol{\Omega})=\int_{\text{All
sky}}d\boldsymbol{\hat{\Omega}}^{\prime}\delta(\boldsymbol{\hat{\Omega}}^{\prime})W^{\theta}(\boldsymbol{\Omega}^{\prime},\boldsymbol{\Omega})$
(32)

where $W^{\theta}(\boldsymbol{\Omega}^{\prime},\boldsymbol{\Omega})$ is the
top-hat filter in configuration space, given by

$W^{\theta}(\boldsymbol{\Omega}^{\prime},\boldsymbol{\Omega})=\begin{cases}\frac{1}{2\pi(1-\cos\theta)}&\arccos(\boldsymbol{\hat{\Omega}}\cdot\boldsymbol{\hat{\Omega}}^{\prime})\leq\theta\\\
0&\text{otherwise}\end{cases}$ (33)

The integral in equation 33 is computationally expensive to perform in
configuration space, but one can use properties of the spherical harmonics and
massively speed up the calculation by going to harmonic space. Let the
expansions of the smoothed and unsmoothed fields in spherical harmonics be
given by

$\displaystyle\delta(\boldsymbol{\hat{\Omega}})$ $\displaystyle=\sum_{\ell
m}\alpha_{\ell m}Y_{\ell m}(\boldsymbol{\hat{\Omega}})$ (34)
$\displaystyle\delta^{\theta}(\boldsymbol{\hat{\Omega}})$
$\displaystyle=\sum_{\ell m}\alpha^{\theta}_{\ell m}Y_{\ell
m}(\boldsymbol{\hat{\Omega}})$ (35)

Since the top-hat filter represents a homogeneous and isotropic smoothing
kernel that is only a function of the central angle $\theta$ between
$\boldsymbol{\hat{\Omega}}$ and $\boldsymbol{\hat{\Omega}}^{\prime}$, it can
be expanded in terms of the Legendre polynomials $P_{\ell}(\cos\theta)$ as:

$W^{\theta}=\sum_{\ell}b_{\ell}P_{\ell}(\cos\theta)$ (36)

It can be shown, by substituting equations 34 to 36 in equation 32, that the
spherical harmonic expansion coefficients $\alpha^{\theta}_{\ell m}$ of the
smoothed field are given by the product of the $\alpha_{\ell m}$ of the
unsmoothed field and the Legedre expansion coefficients $b_{\ell}$ of the top-
hat filter (Devaraju, 2015)262626Note that we have an additional factor of
$4\pi$ not present in the expression derived by Devaraju (2015), since we use
a normalised definition for the top-hat filter that already incorporates an
extra factor of $1/4\pi$ that appears in their equivalent of equation 32.:

$\alpha^{\theta}_{\ell m}=4\pi\frac{b_{\ell}}{2\ell+1}\alpha_{\ell m}$ (37)

This expression makes sense since computing the integral in equation 33 is
equivalent to performing a spatial convolution on the surface of a sphere. A
convolution in configuration space corresponds to a product in harmonic space.
It is to be noted that for a general inhomogeneous smoothing kernel, however,
spatial smoothing is not equivalent to a convolution since each point on the
sphere has a different kernel. For such cases, equation 37 would be different.
See chapter 2 of Devaraju (2015) for a nice discussion on the topic.

In practice, we compute the $\\{\alpha_{\ell m}\\}$ for the unsmoothed field
using the healpy package. The $\\{b_{\ell}\\}$ for the top-hat function are
computed as follows. By definition,

$\displaystyle b_{\ell}$
$\displaystyle=\frac{2\ell+1}{2}\int_{-1}^{1}W^{\theta}P_{\ell}(\cos\theta)d(\cos\theta)$
$\displaystyle=\frac{2\ell+1}{2}\frac{1}{2\pi(1-\cos\theta)}\int_{\cos\theta}^{1}P_{\ell}(\cos\theta)d(\cos\theta)$
(38)

Using the recursion relation

$P_{\ell}(x)=\frac{1}{2\ell+1}\frac{d}{dx}\left[P_{\ell+1}(x)-P_{\ell-1}(x)\right]$
(39)

we get

$b_{\ell}=\frac{1}{4\pi(1-\cos\theta)}\left[P_{\ell-1}(\cos\theta)-P_{\ell+1}(\cos\theta)\right]$
(40)