The angular power spectrum of gravitational-wave transient sources as a probe of the large-scale structure

We present a new, simulation-based inference method to compute the angular power spectrum of the distribution of foreground gravitational-wave transient events. As a ﬁrst application of this method, we use the binary black hole mergers observed during the LIGO, Virgo, and KAGRA third observation run to test the spatial distribution of these sources. We ﬁnd no evidence for anisotropy in their angular distribution. We discuss further applications of this method to investigate other gravitational-wave source populations and their correlations to the cosmological large-scale structure

Introduction-Since the first detection of a gravitational wave (GW) signal from a binary black hole (BBH) coalescence in 2015 [1], LIGO, Virgo, and KAGRA (LVK) have detected dozens more such signals during the first three observation runs [2].At the end of the next two observation runs, the number of detections is expected to reach the thousands [3].This abundance of detected events will allow us to continuously refine our knowledge of the GW emitters.
In this context, an area of growing interest is the measurement of the spatial distribution of GW (SDGW) transient sources and its relation to the large-scale structure (LSS) of the universe [4][5][6][7].The SDGW provides a means to test the LSS that is complementary to electromagnetic measurements as well as dark siren analyses [8,9], which rely on cross-referencing GW detections with galaxy catalogs and are prone to complications such as catalog incompleteness and selection bias.Developing a scheme to accurately measure the SDGW constitutes one of the critical milestones towards precision cosmology with GWs [10].
In this paper, we present a novel, simulation-based inference method to test the SDGW that borrows from techniques used in electromagnetic precision cosmology, in particular the study of the cosmic microwave background radiation (CMB).Specifically, we show how to calculate the observed angular power spectrum of foreground GW events and use it to probe the SDGW.This technique provides complementary information to analogous studies based on the astrophysical GW background, where the angular power spectrum is derived from the clustering statistics of the BBH host galaxies [11][12][13][14][15].
As a first application of our method, we test the isotropic source distribution hypothesis for the confident BBH mergers observed during the third LVK observing run (O3).However, it should be stressed that our approach is not limited to this specific instance.The technique that we present here can be easily generalized to various GW sources, future GW searches with additional detections, and different test hypotheses on the SDGW and its correlation with the LSS.
In the next two sections, we discuss the basics of our method, the selection of GW events, the generation of synthetic signals to test the isotropic hypothesis, and the production of sky localization maps via parameter estimation.In the last two sections, we present the main results and discuss future extensions of this work.
Methodology-Our method probes the spatial distribution of BBH merger events by computing their observed angular power spectrum [16] and comparing it to a fiducial distribution.In this work, we select the isotropic distribution, which corresponds to testing whether BBHs are isotropically distributed in the local universe.First, we compute the power spectrum of observed BBH events from the LVK GW catalogs.We choose a suitable subset of these events by imposing the selection cuts detailed in the next section.Then, we compute the power spectra of a number of mock sets obtained by injecting synthetic signals into real detector data.We sample their parameters from the latest LVK population analysis posterior distributions [17] and inject the signals isotropically on the sky.We then select a subset of events by imposing the same selection cuts used for the observed BBH mergers.The synthetic power spectra are combined to produce a fiducial distribution of an isotropically distributed angular power spectrum as would be measured by the LVK detectors.Finally, we perform a statistical consistency test of the observed BBH angular distribution with the fiducial isotropic distribution; for each multipole component of the power spectrum, we compute the p-value that the observed multipole belongs to the fiducial distribution.
We consider the subset of BBH events detected during the LVK O3 observing run with a false alarm rate (FAR) smaller than 1 yr −1 as reported in Ref. [17].We further restrict our sample to three-detector events.This is required for the generation of a consistent fiducial angular distribution, as the accuracy of sky localizations depends on the number of detectors [18].These conditions restrict the sample of O3 events to 34.These events constitute our catalog of observed signals.To generate the synthetic signals, we draw their source parameters from their inferred median population distributions [17], assuming the Power Law + Peak model for the primary mass [19] with a power law on mass ratio, the Default spin model [20,21], and a power law model for redshift evolution [22].The phase and orientation parameters are sampled from distributions with isotropic orientations.We inject the signals into real detector data with an isotropic distribution in the sky.The times of the injections are uniformly sampled during O3.We then downselect these times to periods that do not overlap with known non-astrophysical transient noise [23] and GWTC-3 confident detections [2].The signals are simulated with the IMRPhenomPv2 [24,25] waveform model.Selecting the synthetic events based on their FAR is computationally expensive, as it requires doing PE for the full set of events.To avoid this computational cost, we substitute the FAR selection cut with a threshold on the optimal network signal-to-noise ratio (SNR) ρ N .We choose ρ N > 10, following the approximate threshold used for the semianalytic sensitivity estimates in Ref. [17].
We generate a catalog of 3,400 synthetic events.This allows us to produce meaningful statistical results while limiting the computational cost required to perform PE and generate the sky localization maps.We use the synthetic signals to create 100 random mock sets of 34 events each.These sets provide independent realizations of what the detectors would observe under the hypothesis that the events are isotropically distributed in the sky.We use these sets to generate the fiducial distribution.We perform PE of all observed and synthetic events with bilby pipe [26].We use the IMRPhenomPv2 waveform for the signal model and draw the samples from the posterior distribution with the nested sampler dynesty [27].
We adopt the standard LVK uniform priors on the mass ratio and chirp mass from Ref. [2].We restrict the chirp mass to a ±12M ⊙ range around the injected values of the synthetic events and the median values of the O3 observed events.Additionally, we constrain the priors on the primary and secondary masses to be within the interval [1,120] M ⊙ .The prior on all other parameters is chosen according to the uninformative priors adopted in standard LVK analyses [2].We then use the posterior samples for the declination and the right ascension to produce sky maps.Angular power spectrum-Following Ref. [16], we treat the event sky localization error regions as probability density heat maps.We generate the combined sky localization map of the observed GW events, M (χ, ϕ), by stacking the sky localization density maps of all events in the observed catalog.Here, χ and ϕ are the polar and azimuthal angles on the celestial sphere, respectively.Figure 1 shows the Mollweide representation of M (χ, ϕ).We repeat this procedure to obtain a cumulative sky localization map for each set of synthetic events.Figure 2 shows the combined sky localization map obtained by stacking the 100 synthetic maps, each made from 34 events.The map shows that the synthetic events are isotropically distributed in the sky.It also depicts what the GW sky would look like with 3400 foreground BBH events, a not too unrealistic scenario in a few years.
We then compute the angular power spectra of the combined sky localization maps by expanding each of them into spherical harmonics: The multipole components of the angular power spectrum are obtained by summing the absolute square of the α lm coefficients of the expansion over m: ( The physical information contained in the power spectrum can also be expressed in terms of the two-dimensional (angular) correlation function (CF).The CF describes the excess probability of finding two objects in the directions n1 and n2 and angular separation θ with respect to a uniform distribution.Given the cumulative sky localization map M (χ, ϕ), the CF is defined as where the average is taken over the observed sky with angular separation held fixed [16].The CF can be written in terms of the power spectrum as where P l (cos θ) denotes the Legendre polynomial of order l and argument cos θ.Typically, the finite beam resolution of the detectors leads to a high-l cutoff l max in Eq.( 3).This effect can be modeled by introducing a window function W l ∝ exp −l(l + 1)σ 2 res , where σ res is the detector resolution [30].
The diffraction-limited angular resolution of the LIGO-Virgo network determines the high-l cutoff as l max ∼ π/θ res , where θ res is the angular resolution.We estimate l max directly from the distributions of the skymaps.We fit the distribution of the observed skymap 90% contour regions as a proxy for the square angular resolution ∆Ω res = 2π[1 − cos(θ res /2)] with a gamma distribution.We then perform a one-tailed test and choose ∆Ω res such that 90% of the observed events have a larger localization area than that value.This provides an estimate for the angular resolution of θ res,o ∼ 6.95 • , corresponding to l max,o ∼ 26.We then repeat the procedure for the whole set of synthetic events.This yields θ res,s ∼ 4.83 • , corresponding to l max,s ∼ 37.The resolution of the simulated set is better than the resolution of the observed set.We expect this is due to the larger number of events in the simulated set compared to the observations.As a consistency check, we also estimate ∆Ω res using the theoretical estimate of Ref. [31].For a monochromatic GW at frequency f , the square angular resolution of a three-detector network is where A N is the triangular area formed by the three detector sites, i N is the angle between the wave direction and the three-detector plane, ρ N is the network optimal SNR of the GW signal, and ρ i (i =1,2,3) are the single-detector SNRs.We consider a triangular area A N = 10 17 cm 2 for the LIGO-Virgo network and a mean incidence angle of 45 • with the detector plane.We use the posterior sample median values to estimate the SNRs and approximate f with the ISCO frequency obtained from the posterior median chirp mass and mass ratio.Using the means of the SNRs and f in Eq. ( 4), we obtain the angular resolution θ res,o ∼ 4.04 • for the observed events and θ res,s ∼ 4.44 • for the synthetic events, corresponding to l max,o ∼ 45 and l max,s ∼ 41, respectively.The theoretical estimate gives higher bounds than the data sets.This is expected, as Eq. ( 4) is derived under optimal assumptions and a Fisher approximation.In the following, out of an excess of caution, we will use l max = 26 as a conservative upper bound.
Results-Figure 3 shows the power spectrum of the observed events (red curve) and the mean spectrum of the 100 synthetic sets (black curve) up to l max = 26.For each l, we fit the C l distribution from the synthetic sets with a gamma distribution.The three gray-filled areas in Fig. 3 (darker to lighter gray) denote the 1 -3σ confidence level regions from the mean.All observed C l values lie within the 2σ band.Therefore, we conclude that the observed angular distribution of observed BBH events shows no significant inconsistencies relative to an isotropic distribution.To quantify this statement, we performed two statistical tests.In the first test, we compute the cumulative distributions of p-values for the observed C l under the hypothesis that the BBH are distributed isotropically in the sky.The observed power spectrum of the O3 BBH events considered in the analysis (red curve) and the fiducial power spectrum obtained from the 100 synthetic sets under the isotropic hypothesis (black curve).The gray-filled regions denote 1 -3σ deviations from the mean.
Figure 4 shows the cumulative distributions of p-values (red dots).The expected distribution is represented by the black dashed line, with the gray-filled regions denoting the 1 -3σ confidence levels.All p-values lie within the 2σ region, in agreement with the results of Fig. 3.In the second test, we assess the goodness of fit of the observed power spectrum with the fiducial spectrum by performing a χ 2 test, which yields a p-value of 0.82, in agreement with the null isotropic hypothesis.Finally, we test the isotropy hypothesis with the CF. Figure 5 shows the CF for the observed set and the fiducial correlation function obtained from the 100 synthetic sets under the isotropic hypothesis, where we have set the window function resolution to σ res = l max .Consistent with the power spectrum result, the observed CF is in agreement with the fiducial isotropic distribution within 2σ.The behavior of the CF at small scales, C(θ) = (θ/θ 0 ) 1−γ , provides a test of isotropy [16].We first compute the power-law slope γ of each synthetic CF at the minimum angular resolution θ res,s with a log-log fit.Averaging the values, we obtain a fiducial value of γ s = 2.05 ± 0.35, which is consistent with an isotropic distribution (γ = 2).We then compute the power-law slope for the observed set at the same angular scale.The observed power-slope is γ o = 1.96.This is in agreement with the null isotropic hypothesis with a p-value of 0.45.
Conclusions-In this paper, we have developed a new, simulation-based inference framework to probe the spatial distribution of observed, foreground GW events.Our approach compares the power spectrum of observed GW signals to a fiducial power spectrum from a theoretical distribution.As an application of this method, we tested the isotropy hypothesis of the BBH mergers observed during the O3 LVK observing run.As foreseen [4][5][6][7], we found no evidence of anisotropy at the 2σ confidence level.
Our method provides a powerful framework for testing the universe's LSS that complements current GW back- ground searches [32,33].Due to the phase-coherence of matched-filter searches employed in GWTC-3 [2], we are able to access higher multipole moments than background searches [34].Relying on resolved sources allows us to achieve astrometric resolution at the square degree level [10].Although the two approaches essentially target the same signal in the limit of many detections, our method has a higher resolution and is more sensitive than background analyses.A first, straightforward extension of this work is to refine the test of BBH isotropy as more GW events are discovered.Tests of specific theoretical models of anisotropic distributions and cross-correlations with astrophysical populations in the EM domain are two additional applications.Our approach can also be directly extended to include information about the source distances.Statistical associations between the observed GW populations and other extragalactic populations may be within reach of current and next-generation GW detectors.This method will provide a means to rapidly detect and quantify any such associations.comments.We also thank Sylvia Biscoveanu for discussions regarding PE using bilby pipe, Derek Davis for discussion on detection statistics, Stuart Anderson for helping provide computing resources, and Reed Essick for useful comments on the manuscript.This manuscript was assigned LIGO-Document number P2300106.

Figure 1 :
Figure 1: Combined sky localization map of the O3 BBH events considered in the analysis.The sky localization of each event is generated with Bayestar [18] from the PE posterior samples for the declination and the right ascension.The map is created with the Healpy package [28, 29].

Figure 2 :
Figure 2: Combined sky localization map of the synthetic BBH events that are used to build the fiducial power spectrum.Their isotropic distribution in the sky is shown by the map.

Figure 3 :
Figure 3: The observed power spectrum of the O3 BBHevents considered in the analysis (red curve) and the fiducial power spectrum obtained from the 100 synthetic sets under the isotropic hypothesis (black curve).The gray-filled regions denote 1 -3σ deviations from the mean.

Figure 4 :
Figure 4: The cumulative distribution of observed p-values for the C l .The black solid line indicates the expected distribution under the isotropic hypothesis.The gray-filled regions correspond to 1 -3σ deviations from the expected distribution.

Figure 5 :
Figure 5: The observed correlation function of the O3 BBH events (red curve) and the fiducial correlation function under the isotropic hypothesis (black curve).The gray-filled regions denote 1 -3σ deviations from the mean.