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

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 ~$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.


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  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  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  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  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. 2021Mukherjee et al. , 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  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  0 nearly independently of merging binary population models1 .
Since the various techniques of measuring  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(Riess et al. , 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. 2018Scelfo et al. , 2020;;Adhikari et al. 2020;Gagnon et al. 2023).With the third generation of gravitational wave detectors likely to bring in ∼ 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. 2021Libanore et al. , 2022;;Balaudo et al. 2023;Vĳaykumar 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 -nearest-neighbour distributions (Banerjee & Abel 2021a) as summary statistics.The nearest-neighbour measurements are sensitive to all -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 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×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 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.

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.

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 from this parent set of ∼ 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.5 2 .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 twodetector 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. 2 The classification probabilities were calculated using the GW package of PESummary (Hoy & Raymond 2021).

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 population3 .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 (U), which translates to drawing their right ascension () from U (0, 2) and sine of declination (sin ) from 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 package4 .These are as follows: (i) 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) (ii) 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  w.r.t. the plane of orbital angular momentum, polarisation angle  and phase at coalescence Φ  over their allowed physical ranges, i.e., ,  ∈ U (0, ) and Φ  ∈ U (0, 2), 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, BAYESTAR5 (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 package6 .Finally, we compute the (phase-maximised) network matchedfiltered signal-to-noise (SNR henceforth) for each event and classify the events with an SNR ≥ 10 as 'detections'7 .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 ∼ 3 deviation from the observed distribution in the second last histogram bin of figure 4a.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.

Galaxy Catalogue
We use galaxies and quasars from the publicly available WISE×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 (∼ ) 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 hemisphere8 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,  1 = 3.4,  2 = 4.6,  3 = 12, and  4 = 23, while SuperCOSMOS has data in three optical bands, , , and .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    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.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  (−) and/or on  1 − 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 ∼ 3 steradians in the sky, corresponding to a sky coverage of ∼ 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 effects9 .The sky distribution of the WSC galaxies and quasars is shown in figure 5.
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 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.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).
The final catalogue contains ∼ 1.7×10 7 sources, with ∼ 15 million classified as galaxies and ∼ 2 million as quasars.This translates to an average number density of more than 600 sources per sq.deg. in the sky.

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

Auto-clustering
Consider a set  of   discrete, point-like tracers or data points10 .
where   ,   represents the declination and right ascension of the  th tracer in celestial equatorial (J2000) coordinates.In polar coordi-nates, the position of the  th tracer is given by To study clustering, we need a metric for the distance between two points ( 1 ,  1 ) , ( 2 ,  2 ) in the sky.A natural choice is the greatcircle distance given by  = , where  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 ( 1 ,  1 ) , ( 2 ,  2 ) is computed using the haversine formula (RIOS 1795) where hav () ≜ sin 2 (/2) is the haversine function.

Angular Power Spectrum
The sky positions of the tracers  can be used to define a number density field   (, ) in the sky, such that 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 n =   4  .The overdensity field is given by The overdensity field contains all the information about the autoclustering of the tracer set .By expanding   (, ) in terms of spherical harmonics  ℓ (, ), one can derive the angular power spectrum C , ℓ , a widely used two-point summary statistic for clustering: where ℓ goes from 0 to ∞ and  takes values from −ℓ to ℓ.In practice, the summation is cut off at some ℓ max determined by the resolution of the numerical grid on which the field is defined.The power spectrum  , ℓ at a particular value of ℓ quantifies the clustering strength at an angular scale corresponding roughly to  ≈ /ℓ.In this study, we use the HEALPix11 scheme (Górski et al. 2005) as implemented in the python library healpy12 (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-ℓ cutoff of 3×NSIDE -1 in our analysis.
An equivalent measure of spatial clustering of a set of tracers is the two-point autocorrelation function,  , (), which captures the excess probability of finding two data points separated by an angular distance of  in the sky over finding two points drawn from a random (Poisson) distribution in the sky.Mathematically, where Ω1 and Ω2 are unit vectors separated by an angular distance  in the sky such that Ω1 • Ω2 = cos , and the angular brackets denote an average over all such configurations of Ω1 and Ω2 .It can be shown that the angular power spectrum is related to the two-point function of the tracers in the following way: where  ℓ (cos ) denotes the Legendre polynomial of order ℓ and argument cos .Therefore, the two-point function and the angular power spectrum encode the same physical information about the clustering of the tracer set ; the angular power spectrum is a twopoint clustering statistic.

Nearest-neighbour distributions: 𝑘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  = 2(1 − cos ) to study clustering at an angular scale 13 .
Suppose we are given the positions of a set  of discrete tracers.Given the discussion above, the fundamental quantity that quantifies their clustering at spatial scale  is the probability P  |  of finding  data points of  in a randomly placed spherical cap of area  in the sky.P  |  can be written in terms of a generating function (| ) as For the case of spherical caps, it can be shown that the generating function is given by (Banerjee & Abel 2021a)14 where  (  ) are the -point correlation functions of the underlying field of the tracers , with  (0) = 0 and  (1) = 1 by definition15 .
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 P > |  of having more than  tracers in a randomly chosen spherical cap of .P > |  can also be expressed in terms of a generating function  (| ) given by (Banerjee & Abel 2021a) such that 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 P > |  directly from the positions of the tracers.
The count distributions P  |  can then be calculated trivially using Consider a set of   area-filling, randomly distributed query points in the sky, such that   ≫   .Each query point will have a data point in  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 -nearest-neighbour data point, or NN-CDF, is precisely equal to P >−1|  .
To understand this connection, let us examine the case of  = 1.Consider   spherical caps of area  = 2 (1 − cos ), 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 .The nearest-neighbour CDF at angular scale  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 .This argument can be easily generalised if we consider the nearest-neighbour instead of the first nearest-neighbour.Therefore, we conclude In practice, the 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 -nearest-neighbour data point of each query point.The nearestneighbour 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  neighbouring data points for all query points simultaneously.Sorting the computed distances for each neighbour index  immediately gives the empirical CDF of the -nearest-neighbour distances over a range of spatial scales.The empirical CDF converges to the true NN-CDF in the limit of a large number of query points.In this study, we 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 ℓ max ∼ 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.
utilise the sklearn.neighbors.BallTree routine from the library scikit-learn16 (Pedregosa et al. 2012) with the haversine distance metric for our purposes.
It is evident from equation 11 that the count distributions P > |  , and hence the NN-CDFs, are formally sensitive to integrals of all 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 NN-CDFs over two-point clustering statistics.

Cross-clustering
There are two possible approaches to quantify the cross-clustering between the BBHs and the WSC catalogue: (i) directly cross-correlate the WSC source positions with the BBH positions (ii) 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 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 tracerfield cross-correlation formalism.
We now describe the summary statistics for measuring the spatial cross-correlation between a set of discrete tracers, , and a continuous field   .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.

Cross Angular Power Spectrum
As discussed in section 3.1.1,we can define an overdensity field   for the discrete tracers , and expand both   and   in spherical harmonics The cross angular power spectrum between  and  is defined as We compute the cross angular power spectrum in a similar manner to the auto angular power spectrum, using healpy's anafast routine.

Nearest-neighbour distributions
Banerjee & Abel (2023) generalised the 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 NN-CDFs.In other words, we must investigate the behaviour of the NN-CDFs of a set of discrete tracers of the field   as their average number density n tends to infinity.If the tracers represent a local Poisson process on the field   , Banerjee & Abel (2023) showed that the continuum version of the NN measurements at a spatial scale  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, where    represents the continuous field smoothed on the scale , (   ) represents the probability density function (PDF) of the smoothed field.The nearest-neighbour index  for discrete data maps to a threshold density  * on the smoothed continuous field, such that n  (1 +  * ) = 17 .Now that we have a continuum version of the NN-CDFs, we discuss the characterisation of spatial cross-correlations using the NN formalism.Following Banerjee & Abel (2023), we take the joint probability P ≥ ,>  * () of finding at least  tracers and the smoothed continuous field    to cross threshold  * in spherical caps of angular radius  as a measure of the spatial cross-correlations between tracers and a continuous field.Assuming that the tracers  represent a local Poisson process on the overdensity field   , we have, similar to Banerjee & Abel (2023), where (   ,    ) is the joint probability distribution of the two fields when smoothed on angular scale .The quantity of our interest, P ≥ ,>  * , can be written in terms of P ,>  * as Suppose the tracers  are completely uncorrelated and statistically independent of the continuous field   .In that case, the joint distribution function can be factored into a product of the individual PDFs of the smoothed fields, i.e., (   ,    ) ∝ (   ) (   ).In this case, we have (Banerjee & Abel 2023) P ≥ ,>  * = P ≥  × P >  * .This fact provides a convenient way to define a summary statistic that measures the excess cross-correlation between the tracers and the continuous field This is a useful quantity to measure, as a positive (negative) measurement for  ,  * − 1 would indicate that the tracer  is correlated (anti-correlated) with the field   , while  ,  * −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  and the field  is contained in the joint distribution (   ,    ).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 NN formalism provides a powerful way to characterise cosmological cross-correlations.
The joint probability distributions P ≥ ,>  * and excess crosscorrelations  ,  * are simple to compute numerically.We follow the procedure laid out in Banerjee & Abel (2023) (i) Create a set of area-filling query points by creating a finelyspaced HEALPix grid in the sky, such that the number of pixels  pix is far greater than the number of data points.
(ii) Build a Ball tree from the set of tracer positions and estimate the angular distances to the -nearest neighbour data points from each query point.For each , sort the distances to produce the empirical NN-CDF over a range of angular scales .In the limit of large  pix , the empirical CDF approaches P ≥  .
(iii) Smooth the continuous field   on an angular scale  using a top-hat filter.The smoothing is done in harmonic space using the {  ℓ } 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).
(iv) For a given  and threshold  * , compute the fraction of query points for which the  th nearest-neighbour lies at an angular distance less than  and the smoothed field, interpolated to that grid point, exceeds  * .In the limit of large  pix , this fraction approaches P ≥ ,>  * .
(v) Compute the fraction of query points for which the smoothed field, interpolated to the grid point, exceeds  * .In the limit of large  pix , this fraction approaches P >  * .
(vi) From the quantities calculated above, compute the excess cross-correlation using equation 19.
(vii) Repeat steps (iii) to (vi) for different values of the angular scale .
In this study, we choose the constant percentile threshold described in Banerjee & Abel (2023) to define the threshold value  * for the continuous field.Specifically, we set  * =  75 , the value of the 75 th percentile of    .This choice implies that P >  * = 0.25 irrespective of the smoothing scale .

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, nearestneighbour 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 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 footprint18 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 ∼ 1 • to ∼ 35 • , equivalent to ℓ = 6 to ℓ = 180.This choice ensures that (i) we have sufficient sampling in the nonlinear regime, where the nearest-neighbour distributions can capture information not accessible through two-point statistics (ii) the measurements are not affected by the lack of sampling towards the right tail of the auto NN-CDF (see Banerjee & Abel (2021a) for a discussion) These angular distances correspond to projected transverse distance scales of ∼ 15 to ∼ 400 Mpc for a median redshift of ∼ 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 highestdensity locations, the shaded band represents the cosmic variance, and the dash-dot line represents the expected value in the absence of cross-correlation.
The excess cross-correlation, as measured by the first nearestneighbour distribution, lies well outside the shaded region representing the 3 cosmic variance in the randoms at all angular scales smaller than ∼ 4 • .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.

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.

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.

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 (, ).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.

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 (i) for each BBH, draw an (RA, Dec) pair from the sky location posterior (ii) using the drawn samples as the 'true' locations of the events, compute the auto-clustering statistics as defined in section 3.1 (iii) 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.Vĳaykumar et al. (2023b) took a similar approach in a recent clustering analysis with forecast BBH data for the third generation of gravitational wave detectors.

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  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: (i) Draw 1000 samples of (RA, Dec) pairs from the sky location posteriors of each of the BBHs (ii) 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 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.
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 den- sity 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 nearestneighbour 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 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 crosscorrelation.Note that these steps are not taken when computing the BBH auto-clustering since we have BBH data on the entire sky.

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, (), which could be the angular power spectrum, the two-point function or the nearest-neighbour distribution, either as a measure of the BBH autocorrelation or BBH-Galaxy cross-correlation.Let the scales considered in the analysis be { 1 , ...,   }.We define the data vector   as the summary statistic  evaluated on the observed BBH catalogue at angular scale   .Similarly, we define a mock vector    as (  ) evaluated on the  th realisation of the mock BBH catalogue for each of the  realisations.Note that () 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 where the angular brackets denote an average over the  realisations of the mock catalogue.The covariance matrix is, by definition, a ×  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 Once we have the corrected inverse covariance matrix, we characterise the signal-to-noise for clustering using the  2 statistic.For the observed BBHs and each realisation of the mock BBHs, we define the  2 value as where  ≜ { 1 ,  2 , ...,   } and   ≜ {  1 ,   2 , ...,    }.The distribution of  2   (henceforth the null distribution) represents the signal-to-noise expected from data consistent with the null hypothesis, and  2  represents the signal-to-noise measured from the data.A larger value for  2  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 -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  2   >  2  .If the signal is strong enough that none of the mock realisations have a larger  2 than the data, then one must fit a  2 distribution to the null distribution to compute the -value.The null hypothesis is ruled out if the -value is smaller than a chosen detection threshold.

Angular scales
Similar to section 3.3, we choose 10 log-spaced angular bins from ∼ 1 • to ∼ 35 • for the clustering analysis, which leads to a Hartlap factor of 0.92.For the angular power spectrum, we choose 10 linearlyspaced bins between ℓ = 6 and ℓ = 180.
For computing the overdensity fields and the query points for the nearest-neighbour measurements, we use an NSIDE = 256 HEALPix grid with ∼ 7.8 × 10 5 pixels and an angular resolution of ∼ 0.22 • .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 • of the WSC mask boundaries for computing nearest-neighbour excess cross-correlation to avoid any biases due to the presence of the WSC mask.

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

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 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/( nBBH ).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.
The  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  2 for the observed BBHs.The data is consistent with the null hypothesis with a -value of 0.061, calculated by fitting a  2 function to the null distribution.Note that the relatively small -value is most likely caused by the large ( 3) fluctuation at ℓ∼ 80.As evident from figure 10 values over 135 realisations of the mock catalogue, and the curve enveloping it represents the best-fit  2 distribution.The vertical line represents the measured  2 value for the data.The data is consistent with the null hypothesis at a -value of 0.061, calculated using the CDF of the best  2 fit.There is no evidence for a statistically significant clustering signal in the present data.
we conclude that the angular power spectrum does not capture a statistically significant clustering signal in the presently available BBH 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  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 -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.

Nearest-neighbour measurements
Figure 14 shows the first two 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 NN-CDFs of an unclustered, Poisson-distributed dataset in the sky.The analytic expression for the NN-CDFs of Poisson distributed points is only a function of nBBH , as can be seen from equation 11 (see also Banerjee & Abel (2021a)).
The left panel shows that CDF 1NN is smaller than CDF 2NN at all scales.This is intuitive since, for each query point, the distance to the  th nearest-neighbour is always smaller than the distance to the (+1) th nearest-neighbour.As a result, for a given scale , the fraction of query points with the first nearest neighbour at a distance less than  would be larger than those with the second nearest neighbour. .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 -value of 0.572, calculated using the CDF of the best  2 fit.No statistical evidence exists for spatial cross-correlation between the presently detected BBHs and the galaxies and quasars from the WSC catalogue.This generally holds, i.e., the NN-CDFs always shift towards larger angular scales with increasing .
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 CDF 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 CDF 1NN approaches the expected value for a Poisson distribution regardless of the positions of the tracers.Mathematically, as the area  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).
The data seems to indicate no clustering signal in any of the NN-CDFs.To confirm the visual intuition, we perform a  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 -value of 0.081.Note that in this case, since the null distribution is heavy-tailed and poorly characterised by a  2 function, we compute the -value by counting the number of mock realisations with  2   >  2  .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 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-correlation 19 .
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 nearestneighbour 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  = 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.
We perform a  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 -value of 0.444.Similar to the auto-CDFs, the null distribution here is also heavytailed.Hence, we compute the -value by counting the number of   >  2  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.mock realisations with  2   >  2  .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.

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 ∼ 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  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 anticorrelation, there is no statistical evidence for any spatial cross-correlation between the presently detected BBHs and the galaxies and quasars from the WSC catalogue.
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 in-Figure 18. Mollweide projection of the overdensity field for the WSC catalogue smoothed on a 10 • 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 crossclustering 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.vestigate 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 • 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  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.

CONCLUSIONS AND FUTURE DIRECTIONS
In this paper, we developed a framework for quantifying the spatial clustering of sources of gravitational waves and their crosscorrelation with the large-scale structure of the universe, using the angular power spectrum and nearest-neighbour distributions as summary statistics.We extended the -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×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 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 crosscorrelations 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 astro-physical 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 tracerfield 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 bins20 .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 ∼ 1.6 × 10 4 BBHs with 1 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.).where  obs is the total observation time and  c  is the differential comoving volume.The probability distribution function for redshifts is given by normalising equation C6 (|,  max ) = where  max is the maximum redshift out to which the population has been created.Note that the constants R 0 and  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  max in LIGO analyses is typically taken to be 2.3, in this work, we assume a  max of 0.7 since assuming a lower value of  max is computationally less demanding.We have checked that inputting a higher  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  = 3 in our analysis.

APPENDIX D: SMOOTHING IN HARMONIC SPACE
Consider a field () defined in the sky.The expression for the field smoothed on an angular scale  is given by The smoothed field is equivalent to the field averaged over spherical caps of angular radius .Equation D1 can be re-written as This expression makes sense since computing the integral in equation D3 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 D7 would be different.See chapter 2 of Devaraju (2015) for a nice discussion on the topic.
In practice, we compute the { ℓ } for the unsmoothed field using the healpy package.The This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
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.

Figure 2 .
Figure 2. Distribution of the 1 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).

Figure 3 .
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.

Figure 4 .
Figure 4. Distribution of the 1 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.

Figure 5 .
Figure 5. Mollweide projection of the skymap of ∼ 1.7 × 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.

Figure 6 .
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.
Figure7.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 ℓ max ∼ 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.

Figure 8 .
Figure8.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 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.

Figure 9 .
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 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.
Figure10.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 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.

Figure 11 .
Figure10shows 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 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/( nBBH ).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.The  2 significance test results for the angular power spectrum are presented in figure11, with the histogram representing the null distribution and the solid vertical line representing the measured  2 for the observed BBHs.The data is consistent with the null hypothesis with a -value of 0.061, calculated by fitting a  2 function to the null distribution.Note that the relatively small -value is most likely caused by the large ( 3) fluctuation at ℓ∼ 80.As evident from figure10, this fluctuation does not indicate a clustering signal.Hence,

Figure 12 .
Figure12.Results of BBH-Galaxy cross-clustering analysis conducted using the cross angular power spectrum as the summary statistic, with a zoomedin view provided in the inset for better visibility.The plotting scheme is identical to figure 10, except the dash-dot line,  = 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 14 .Figure 15 .
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 NN-CDFs of an unclustered, Poisson-distributed dataset in the sky.Different colours represent different values of the nearest-neighbour index , with blue and orange standing for  = 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 NN-CDFs.

Figure 16 .Figure 17 .
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.

Figure A1 .
Figure A1.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 ∼ 30 ⊙ is clearly visible.

Figure B1 .
Figure B1.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.
( ′ , ) (D2)where   ( ′ , ) is the top-hat filter in configuration space, given by  ( ′ , ) = 1 2  (1−cos  ) arccos( Ω • Ω′ ) ≤  0 otherwise (D3)The integral in equation D3 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 ( Ω) = ∑︁ ℓ  ℓ  ℓ ( -hat filter represents a homogeneous and isotropic smoothing kernel that is only a function of the central angle  between Ω and Ω′ , it can be expanded in terms of the Legendre polynomials  ℓ (cos ) as:  = ∑︁ ℓ  ℓ  ℓ (cos ) (D6)It can be shown, by substituting equations D4 to D6 in equation D2, that the spherical harmonic expansion coefficients   ℓ of the smoothed field are given by the product of the  ℓ of the unsmoothed field and the Legedre expansion coefficients  ℓ of the top-hat filter(Devaraju 2015) 26 :