Two-dimensional correlation function of binary black hole coalescences

We compute the two-dimensional correlation functions of the binary black hole coalescence detections in LIGO-Virgo first and second observation runs. The observed correlation function is compared to two reference functions obtained from artificial maps to test for any possible excess correlation of the binary black hole coalescence events at different angular scales. No excess correlation at any angular scale is found. The power-law slope of the correlation function is estimated to be $\gamma= 2.24\pm 0.33$ at 3-$\sigma$ confidence level, a value consistent with a random distribution of sources.


Introduction
On September 14, 2015, researchers from the Laser Interferometer Gravitational-wave Observatory (LIGO) [1] Scientific Collaboration (LSC) and the European Virgo Collaboration [2] made the first direct detection of Gravitational Waves (GWs) from a pair of coalescing black holes [3]. Less than two years after that first announcement, LIGO and Virgo observed GWs from the merger of two neutron stars [4], an event which was rapidly followed by the Fermi Gamma-ray Space Telescope's detection of a gamma-ray flash, and eventually by optical, infrared, radio, and x-ray observations by hundreds of telescopes around the world in what became the most observed event in the history of modern astronomy [5].
Nowadays, GW astronomy is a well-established scientific discipline. In the first two observing runs of Advanced LIGO and Virgo, O1 and O2, LSC and Virgo collaboration researchers observed ten Binary Black Hole (BBH) coalescence detections and one binary neutron star coalescence detection. The third observation run, O3, brought us detections on a weekly basis [6], enabling a plethora of novel astrophysical and theoretical investigations.
The next decade will see GW astronomy further expand its reach in frontier scientific research. Japan's KAGRA detector [7] has joined the international network of GW ground-based observatories. India has established the LIGO-India Scientific Collaboration (LISC) and finalized plans for the construction of the LIGO-India detector [8]. The European space-based LISA mission [9], slated to launch in 2034, will greatly improve detection capabilities and localizations of astrophysical sources. The International Pulsar Timing Array project will detect ultra-low frequency GWs within ten years [10]. Optical, particle and GW astronomy will together explore the Universe through complementary physical carriers.
The publication of LIGO-Virgo's first catalog of compact binary merger signals [11] has shown that GW astrophysics is a powerful tool for population and source property studies of compact objects, tests of General Relativity, and large-scale cosmological measurements. However, many open questions still remain. For example, tests of GR have returned a null result [12], the formation channels of black hole binaries [13] and the physics of EM-bright mergers [14] are still unclear, as well as the determination of the Hubble constant from GW sirens [15]. Despite LIGO and Virgo running all-sky, unmodeled searches [16,17], no GW signal has been detected that cannot be modeled as a compact binary coalescence. Other sources of multi-messenger signals such as isolated compact objects [18], CCSN [19] and magnetars [20] have not been observed in the GW domain.
The rapid growth in the number of BBH coalescence detections and the dramatic improvement in their sky localizations are turning GW astrophysics into a precision observational science like large-scale structure astrophysics and early-universe cosmology. One important physical concept in large-scale structure investigations and observational cosmology is that of the Correlation Function (CF) [21]. The (two-point) (auto-)CF describes the excess probability of finding pairs of points at a given separation. In large-scale astrophysics and observational cosmology, the CF (or its Fourier transform, i.e., the power spectrum) is commonly used to describe the spatial distribution of galaxies or the density fluctuations observed in the cosmic microwave background. The CF from galaxy surveys, for example, allows astronomers to estimate the distance scales of galaxy clustering and gain information about the origin and evolution of the Universe's large-scale structures.
The purpose of this short article is to introduce the concept of the CF for BBH coalescence events. We use the public BBH coalescence detections in LIGO-Virgo's O1 and O2 runs to compute the two-dimensional CF for the population of these objects. Ten detections with sky localizations ranging from 39 square degrees to 1666 square degrees are clearly not sufficient to draw any meaningful conclusion on the spatial distribution of BBH coalescences. However, this calculation shows that the CF can be used to investigate the statistical properties of the population of these objects. We illustrate the method by comparing the two-dimensional CF obtained from the LIGO-Virgo O1-O2 BBH detections to a CF obtained by a random distribution of the same detections. The result shows that the two-dimensional spatial distribution of the detections is consistent with an isotropic distribution, as reported in Ref. [22] by implementing a pixelization-based method for the O1-O2 BBH detections. We also confirm this conclusion by comparing the CF to a synthetic CF obtained by simulating a number of BBH detections with sky localization error regions consistent with those of the LIGO-Virgo sample.

Two-dimensional Correlation Function
In our analysis, we follow the customary definition for the two-dimensional (angular) CF of largescale astrophysics [21]. The two-dimensional CF of a population of objects describes the excess probability of finding two objects separated by the angular distance θ w.r.t. a uniform distribution. To compute the CF of the BBH population, we treat the sky localization error regions of the BBH detections as heat maps. Given the (normalized) sky localization error region map of the i-th BBH detection in the sample, M i (θ, ϕ), where θ and ϕ are angles on the celestial sphere, we define the sky localization probability density map of the sample as: where N is the number of BBH detections and F i are probability weights that depend on the GW detector network sensitivity. By expanding the sky localization map in spherical harmonics, the sky correlation function of the BBH sample can be defined as where the average is taken over the observed sky with angular separation θ 21 held fixed. Using the addition theorem of spherical harmonics, the CF can be written as where P l (cos θ) denotes the Legendre polynomial of order l and argument cos θ, and we have defined a 2 l = m |a lm | 2 . Note that the CF in Eq. (4) differs with the usual definition of angular power spectrum that is used in CMB cosmology, where a 2 l = (2l + 1)C l . As the map M (θ, ϕ) describes a probability density field, rather than the perturbation field of a physical quantity, in the following we will focus on the CF instead of the power spectrum which is the standard measure for fluctuation fields.
The quantities a 2 l are measured from the sky localization map M (θ, ϕ) and determine the 2dimensional angular distribution of the BBH sample. Comparison of the CF to theoretical models involves modification of Eq. (4) by multiplying the a 2 l coefficients by a window function W l to take into account experimental constraints in the observations. For example, the finite beam resolution of the detector introduces a high-l cutoff that can be modeled with a window function W l ∝ exp[−l(l + 1)σ 2 ], where σ is the detector resolution. If the object population cannot be observed across the full sky, a mask is required. In contrast with CMB observations, where the region of the sky along the galactic plane must be masked in CMB observations due to the impossibility of measuring temperature fluctuations along the galactic plane, the full sky is transparent to GWs and no mask is necessary. As the sky map in Eq. (1) is obtained by summing the sky localization error regions of the BBH detections, the angular resolution is determined by the diffraction-limited spot size of the LIGO detectors where d is the typical separation of the detectors in the network, c is the speed of light, and f is the frequency of the measurement. Assuming a typical frequency of 200 Hz for the detector sensitivity and a LIGO-Virgo detector distance d ∼ 7000 km, a crude estimate of the minimum map angular resolution is θ res ∼ 3 • , or σ ∼ 1/30, implying a high-l cut-off of l max ∼ 30. In the following analysis, for the sake of simplicity, we assume the probability weights in Eq. (1) to be constant, i.e., we assume that the sensitivity of the LIGO-Virgo detector network does not depend on the sky position. (See Ref. [22] for a more refined analysis and a discussion on the effects of detector sensitivity on isotropic test of GW detections.) An additional, possible modification of Eq. (4) is due to the different sensitivities of the GW detector network across the O1-O2 epochs and the varying number of detectors observing each BBH event in the sample. These systematics can be eliminated, at least partially, by comparing the observed CF C obs (θ) to a CF which is computed from reference maps M ref,i (θ, ϕ) obtained by uniformly distributing the observed BBH sky localization error regions in the sky. A more refined analysis could be performed by injecting a population of simulated BBH signals with a uniform angular distribution and then creating the reference map by recovering the sky localization error regions of these injections with the GW network in the same configuration as in the real case. While this procedure would produce a more rigorous CF estimate than the one considered here, we consider it beyond the scope of this paper due to the small sample of BBH detections and the illustrative purpose of our analysis. We plan to revisit this procedure in a future work.

Results
We use the public sky localizations of the O1-O2 LIGO-Virgo BBH detections from the GW Open Science Center [23] and the open source Healpy package [24] to compute the CF. The sky localization error regions of the BBH detections come with different resolutions. We first rescale each map to an NSIDE resolution of 256, corresponding to a pixel angular resolution of θ pix ∼ 0.23 • θ res , and then create the map M obs (θ, ϕ) in Eq. (1) by summing the sky localization error regions of each BBH event and normalizing to the number of detections, such as j M obs (p j ) = 1, where p j denotes the j-th pixel. A Mollweide representation of M obs (θ, ϕ) in Eq. (1) is shown in Fig. 1.
The M obs (θ, ϕ) map is treated as a heat map and the Healpy function map2alm is used to compute the a lm . The coefficients of the Legendre expansion in Eq. (4) are then obtained by summing the a lm in m. We follow the same procedure to compute the CF from reference maps M ref,i (θ, ϕ) used to test possible angular correlation signatures in the CF.
In our analysis, we compare the observed CF C obs (θ) to two reference CFs. The first CF (model A), C ref,A (θ), is obtained by averaging the CF of 1000 artificial maps, each obtained by randomly rotating the maps of each single BBH detection in the sky by arbitrary θ and ϕ angles. The second CF (model B), C ref,B (θ), is obtained by averaging 500 synthetic maps, each consisting of 10 elliptically-shaped sky localization error regions with random orientation and uniformly distributed in the sky. The sky localization areas of these artificial events are chosen such that their semi-axes are R · (x, 1/x), where x is uniformly distributed in (0, 10), and their area πR 2 is drawn from a lognormal distribution with mean (standard deviation) equal to the mean (standard deviation) of the sky localization areas of the observed events. Probability distribution contours of each of these artificial sky localization areas are simulated by superimposing 100 regions built as described above and radius decreasing as f n (R) = R ln(2)/ ln(2 + n), where n = 0 . . . 99. Both reference maps are normalized to the number of detections in the sample, following the same procedure used for the observed map. One example of a synthetic map is shown in Fig. 2. Figure 3 shows the observed CF compared to the reference CF for model A. be seen from Fig. 1, the LIGO-Virgo sky localization error regions of observed BBH events are not perfect ellipses. Even if they were, their ellipticity would not follow a uniform distribution in their semi-axis ratio. Finally, drawing samples from a lognormal distribution of sky localization areas does not accurately represent the observed distribution of sky localizations in O1 and O2. A much more accurate estimate of BBH events angular correlations could be obtained by simulating realistic sky maps by injecting, recovering and localizing events according to the actual sensitivity of the GW detector network. The CF in Eq. (4) can be interpreted as a weighted projection of the spatial two-point CF ξ(r). At small scales, the power-law behavior of the CF is expected to be where θ 0 is an angular correlation scale and γ is the power-law slope of the spatial two-point CF where r 0 is the spatial correlation length. The power-law slope of the BBH distribution can be obtained by fitting C obs (θ) at small angular scales. A weighted best fit of Eq. (6) from θ ∼ θ res to θ ∼ 18 • , where departures from the power-law behavior become evident, gives for the power-law slope γ = 2.24±0.33 at 3-σ c.l., a value consistent with a random distribution of objects, ξ(r) ∼ r −2 .
As a comparison, the power-law slope from Sloan Digital Sky Survey (SDSS) data is γ ∼ 1.8 over the range 0.005 • -10 • [25]. The VIMOS Public Extragalactic Redshift Survey (VIPERS) reports γ ∼ 1.7 − 1.8 for a broad range of galaxy luminosities and stellar masses in the redshift range 0.5 < z < 1.1 [26]. The VIMOS-VLT Deep Survey observes a significant redshift evolution of the luminosity dependence of power-law slope parameter with γ steepening from γ ∼ 1.7 at low redshift to γ ∼ 2.4 for z ∼ 0.9 and galaxies with high intrinsic luminosity [27]. In contrast to SDSS, VIPERS and VIMOS-VLT results, which point to galaxy clustering in the redshift range of BBH detections, z 0.5, our result shows no evidence of clustering at these distance scales. It will be interesting to test whether any evidence of clustering will appear in the data with more BBH detections and better sample statistic.

Conclusions
In this short article, we have computed the two-dimensional CF of BBH observations in the first and second observation runs of advanced LIGO and Virgo. The CF is commonly used in large-scale structure astrophysics and precision cosmology to quantify the spatial distribution of an object class population. Similarly, we have used the two-dimensional CF to measure the statistical properties of the BBH coalescence spatial distribution. By comparing the CF of the LIGO-Virgo detections to a simulated CF from a synthetic sample of sky localizations and a CF obtained by randomly re-orienting the BBH detections, we have shown that the distribution of O1-O2 BBH events in the sky is in agreement with a random distribution of sources, as previously reported in Ref. [22]. The power-law slope of the CF is found to be γ = 2.24 ± 0.33, a value consistent with the upper bound of the power-law slope from galaxy surveys at low redshift z.
While the limited number of O1-O2 detections with large sky localization error regions does not allow us to draw any significant physical conclusions, our work lays the formalism for computing the CF of a class of GW detections. Our analysis is clearly rudimental and can be improved in many ways. The extension to the tens of LIGO-Virgo detections in O3 is straightforward. A better estimate for the two-dimensional spatial distribution of BBH coalescence events could be obtained by comparing the detected CF to a synthetic CF from a realistic population of events as done in Ref. [22]. This could be done by testing the detected CF against a CF from injection sets consistent with the observed BBH coalescence population and detector network sensitivity. The existence of angular correlations in the spatial distribution of BBH coalescences could be tested by building CFs for events distributed isotropically in the sky, or at given angular scales. Comparisons to CFs of anisotropic models for the astrophysical GW background [28] and other astrophysical objects could be used to test correlations with the spatial distribution of these objects and test BBH population paradigms. Our method could also be extended to include information about the distances of the BBH sources by computing the three-dimensional CF [29]. The latter could be compared to CFs obtained from given models of population synthesis, as well as three-dimensional CFs of other astrophysical objects. With the anticipated higher rate of detections and more accurate sky localizations in future LIGO-Virgo observing runs, the CF of BBH and other GW-bright sources may prove itself as another useful tool for GW astronomy investigations.