key: cord-0659961-0fgxqzhp
authors: Payne, Ethan; Banagiri, Sharan; Lasky, Paul; Thrane, Eric
title: Searching for anisotropy in the distribution of binary black hole mergers
date: 2020-06-22
journal: nan
DOI: nan
sha: e13e137051b850711e2c8b6a4ac22c1d8ac82ec7
doc_id: 659961
cord_uid: 0fgxqzhp

The standard model of cosmology is underpinned by the assumption of the statistical isotropy of the Universe. Observations of the cosmic microwave background, galaxy distributions, and supernovae, among other media, support the assumption of isotropy at scales $gtrsim 100$,Mpc. The recent detections of gravitational waves from merging stellar-mass binary black holes provide a new probe of anisotropy; complementary and independent of all other probes of the matter distribution in the Universe. We present an analysis using a spherical harmonic model to determine the level of anisotropy in the first LIGO/Virgo transient catalog. We find that the ten binary black hole mergers within the first transient catalog are consistent with an isotropic distribution. We carry out a study of simulated events to assess the prospects for future probes of anisotropy. Within a single year of operation, third-generation gravitational-wave observatories will probe anisotropies with an angular scale of $sim36^circ$ at the level of $lesssim0.1%$.

Two key assumptions in the standard ΛCDM model of modern cosmology are the isotropy and homogeneity of the Universe [1] [2] [3] . Since ΛCDM cosmology underlies our current understanding of the Universe, rigorous tests are required to validate this fundamental assumption. The Universe has been observed to be isotropic at scales 100 Mpc, with smaller-scale deviations considered to be statistically isotropic on an overall homogeneous structure produced by phenomena such as baryon acoustic oscillations or gravitational interactions [4] [5] [6] . Verifying the lack of large scale anisotropy in the Universe is important for the validity of the ΛCDM cosmology.

Evidence for large-scale isotropy of the Universe has been presented through numerous observations of various sources [6] [7] [8] [9] [10] [11] [12] [13] [14] . The most stringent measurements on the anisotropy of the Universe come from the cosmic microwave background, which generally show small-scale statistical deviations on the order of 10 −4 -10 −5 [13, 14] . Meanwhile, multiple studies hint at the existence of deviations from isotropy at large angular scales in the cosmic microwave background [15, 16] , supernovae [17, 18] and galaxy [19] distributions, as well as large bulk flows [20, 21] . While these inferences are speculative [e.g., [22] [23] [24] , they can be supported or contradicted by independent measures of the Universe's isotropy using gravitational-wave observations.

Studies of the gravitational-wave stochastic background have placed limits of anisotropy from unresolved sources [25] [26] [27] [28] [29] [30] , however observations of gravitational waves from resolved binary black hole (BBH) mergers present another tool to probe anisotropy. Prior to the third observing run, Advanced LIGO (aLIGO) [31] and * Electronic address: ethan.payne@ligo.org Advanced Virgo (aVirgo) [32] released the details for ten binary black holes over the first (O1) and second (O2) observing runs [33] . These observations are collated in the first gravitational-wave transient catalog, GWTC-1. With the third observing run (O3) complete, the total number of gravitational-wave BBH merger candidates has increased to more than 50 [34] . Furthermore, the addition of aVirgo for the entirety of O3 has already resulted in many well-localized sources (e.g. [35] ), providing further motivation for utilizing BBH mergers to study the anisotropy of the Universe.

In this paper, we use LIGO/Virgo data to probe anisotropies in the distribution of BBH mergers, taking care to handle selection bias associated with the detection of gravitational-wave sources. We explore the future prospects of anisotropy measurements with gravitational waves. In Ref. [36] , a HEALPix [37] basis was adopted to parameterize the anisotropy of binary black hole mergers. In Ref. [38] , a two-dimensional correlation function was implemented (see also [39] ). Both analyses find results consistent with isotropy. In contrast, our method utilizes a spherical harmonic basis to define a probability distribution, providing results in terms of typical spherical harmonic functions. Our results are qualitatively similar to those of [36, 38] . However, there are some potential advantages to the spherical harmonic approach: the method lends itself straightforwardly to comparison with results from the cosmic microwave background, and it could be argued that at least some plausible deviations from anisotropy are more clearly visible in the spherical harmonic basis than the pixel basis, which emphasizes hot spots.

The remainder of the manuscript is structured as follows. In Sec. II, we outline a method of hierarchically analysing binary black hole mergers to determine the level of anisotropy, including a parameterization of the distribution of BBH merger events in spherical harmonics and a discussion of the observational selection biases. The analysis of the previously observed BBH events in GWTC-1 is presented in Sec. III. In Sec. IV, we look to the near future by demonstrating the recovery of an isotropic and anisotropic universe with simulated gravitational-wave events. We conclude with an investigation into the future of using gravitational waves as an indicator for anisotropy in the matter distribution of the Universe and the associated implications.

In order to test the assumption of an isotropic universe with gravitational-wave observations, we employ a spherical harmonic basis. We do not employ a linear spherical harmonic basis expansion as this can result in a negative probability density. Instead, we employ a non-linear spherical harmonic basis to guarantee positivity. We write the probability density function for mergers π(Ω) as

Here, Ω is the sky position unit vector, Y m (Ω) are spherical harmonics, and {b m } is the set of coefficients, up to = max , that parameterize the anisotropy of the distribution. The total number of parameters for a given model is max ( max + 2), where we set the monopole term to b 00 = 1 without loss of generality. Ideally, the most complete analysis of anisotropy would use max → ∞. However, the quadratic scaling of the number of parameters with max computationally limits the analysis to a subset of spherical harmonics. In this study, we restrict our analysis to max ≤ 5.

Note that π(Ω) is given not as a linear combination of spherical harmonic coefficients, but as the square of this sum. Coupled with,

implying the b m are complex for m = 0 and that the b 0 are real, we ensure that the probability density function is positive definite. Since the m < 0 modes are uniquely determined from the m > 0 modes, we consider the set of parameters (which we sample over to be) {b m } with ≥ 1, m ≥ 0, but all summations occur over the full range of allowed m's.

The orthogonality of spherical harmonics allows us to express π(Ω) as a linear combination of spherical harmonics

The a m 's are related to the b m by Clebsch-Gordan coefficients:

where

Here, C LM m, m is the set of well-known Clebsch Gordan coefficients [40, 41] . The reconstructed values of a m 's extend up to 2 max . Throughout the manuscript, max refers to the 's of the b m 's.

The relation of the b m to the a m can be further simplified in the limit of small deviations from isotropy, b m /b 00 1, to

where the monopole is always a 00 = 1/ √ 4π. Equation (7) is not used for any calculations in this manuscript, but it is useful for understanding the relationship between the b m and a m . With this model, we aim to measure the posterior distributions of the a m 's, where non-zero values imply the existence of an anisotropic BBH merger distribution.

In order to fit the b m 's, we employ hierarchical inference [42, 43] . We marginalize over Ω to obtain "hyperlikelihoods" for the data given the b m . The hyperlikelihood for a set of data from N events, {d i }, is

Here, {b m } are the spherical harmonic coefficients for the anisotropic model in Eq. (1), up to = max , and p det ({b m }|N ) is the probability of detection given a set of hyper-parameters. The isotropic probability distribution is defined as π(Ω|iso), defined as b m = 0 for ≥ 1, Z iso (d i ) is the ith event's evidence assuming the isotropic sky location probability distribution, and n i is the total number of samples in the ith event's posterior distribution. We note that this implementation does not account for false positive detections of BBH mergers. Procedures to include the possibility of falsely identified sources rely heavily on the probability that a given event was astrophysical [44] . These probabilities vary between pipelines [33] , and therefore we have assumed all detections are astrophysical in origin for simplicity. The detection probability term depends on the selection bias of the interferometers and is discussed in depth below. We use the samples and Bayesian evidence from astrophysical inference of gravitational-wave events. By default, the evidence for event i is calculated assuming an isotropic distribution, Z iso (d i ). The total evidence for isotropy, Z iso , is simply the product of all individual event evidences. The hyper-posterior is

where Z ani is the evidence for the anisotropic model,

and π({b m }) is our hyper-prior.

To present the predicted BBH merger sky position probability density function, we calculate the posterior predictive distribution (PPD)

The PPD predicts the angular distribution given the previously observed events and can be approximated numerically from the hyper-posterior samples. We note that although the PPD allows visualisation of the average posterior, it does not capture the full details of the posterior, hiding the full range of possible distributions. We use the Bayes factor

to determine if one model is preferred relative to another. We adopt a threshold of ln BF ani iso = 8 for a statistically significant signature of anisotropy.

Gravitational-wave observatories do not observe the sky with an isotropic sensitivity. The detector response is characterized by an "antenna factor," which varies depending on the sky position of the source relative to the L-shaped geometry of the interferometer [45] . During a sidereal day, this antenna pattern function is swept about right ascension. Additionally, gravitational-wave observatories are more likely to be operating during the night due to reduced anthropogenic activity [46] . Furthermore, even when the observatories are operating during the day, the noise power spectral density is on average higher than at night.

In order to account for selection effects, we first calculate the detection probability as a function of the sky position of the source, denoted p det (Ω). We calculate the detection probability by determining the fraction of binaries detected at a particular distance, time, and sky position [47] ,

(13) Here, V tot is the total volume of our population model

The variable, V c is the comoving volume, f (z, Ω, t) is the selection function determining the fraction of binaries observed at a given sky position, time, and redshift z.

Finally the variables t 0 and T correspond to the start of the observing period and the duration respectively. We marginalize over time within the network to ensure that the detection probability accounts for all possible times during observing runs when a BBH merger event could be detected. The selection function is calculated by simulating gravitational-wave events from an astrophysically motivated population distribution, π pop (θ), and determining the fraction of detected events from their SNR. For each simulated event, we calculate the signal-to-noise ratio (SNR) for individual observatories ρ ifo and for the entire network ρ net . The expectation value of the matched-filter SNR within a single interferometer for a particular event is [48] ,

(15) This is often referred to as the "optimal" signal-to-noise ratio. Here, F +,× denote the antenna pattern function for plus and cross polarizations, h +,× are the associated gravitational-wave strains for each polarization, and S n (f ) denotes the noise power spectral density (PSD) of the interferometer. For a network with N ifo interferometers, the network signal-to-noise ratio is

We adopt the criteria that a signal is detected if ρ ifo > 8 or ρ net > 12.

For the population distribution, we assume a power law distribution for the source frame masses of the binary systems as given in Ref. [49] with a primary mass power law of π pop (m 1 ) ∝ m −1.6 1 , and a minimum and maximum mass of 7.9 M and 42.0 M , respectively. The mass ratio distribution is assumed to be π pop (q) ∝ q 6.7 . We assume a uniform prior in dimensionless spin magnitude between (0,0.9) with isotropic spin orientations. We employ standard priors for extrinsic parameters with a uniform-in-source-frame volume prior for distance.

We use Monte Carlo integration to calculate the integral in Eq. (13) , which yields p det (Ω). We generate 10 4 realizations of the noise power spectral density for observatories in our network. Each realization is assigned a random time from the O1/O2 observing runs. Using this set of PSDs, we calculate the detection probability, p det (Ω, t i ), as a function of sky position at each time t i . The calculation is carried out using a 3072-pixel HEALPix [37] grid. The fraction of events from the BBH injections that exceed the signal-to-noise ratio detection threshold at each sky position determines p det (Ω, t i ). All realizations of the detection probability at different times are then averaged to determine p det (Ω). Figure 1 presents the detection probabilities (a) for a particular instant in time, (b) marginalized over the first observing run, and (c) marginalized over both the first and second observing runs. Then detection probability at any given time, as seen in Fig. 1(a) , is dominated by the antenna pattern functions of the interferometers. The marginalized detection probability skymaps in Fig.  1(b) and 1(c) , present more non-trivial behaviour due to the interplay of the antenna factors with the timedependence of the detectors' performance. In general, the probability of detection is smeared over the sky over extended periods of times, suppressing the lobe structure due to the antenna pattern function. Within our analysis of GWTC-1, only the combined detection probability is utilized. However the O1 detection probability is presented to highlight the strong selection bias present in O1 analysis. This was primarily due to both LIGO interferometers preferring a mid-declination due to their antenna pattern function maxima, and the consistency of a diurnal cycle during the first observing run [46] . These features are less clear in the combined detection probability due to the duration of the second observing run and improved duty cycle. To determine the detection probability as a function of the hyper-parameters, we compute

where N is the number of events observed [55] . We implement the method outlined above with the Bayesian inference library, Bilby [50] , using the nested sampling algorithm, Dynesty [51] . To implement the spherical harmonic model in the nested sampler, we sample in b m /b 00 for > 0 without loss of generality. To first order, b m /b 00 ≈ a m /2a 00 , implying sampling in b m /b 00 directly presents the degree of anisotropy in the distribution. For each spherical harmonic coefficient, we define a uniform hyper-prior on the magnitude from 0 to 1/ √ 2. For m = 0 modes, this applies to the coefficient b m /b 00 , sampling uniformly over [−1/ √ 2, 1/ √ 2]. For m > 0 modes, we sample in the magnitude and phase of b m /b 00 using the uniform magnitude prior and a uniform phase prior. The prior distributions are constructed to fully span the a m parameter range for ≤ max , while limiting the magnitude of a m 's for > max . This significantly reduces multi-modality of the likelihood, ensuring improved convergence of the nested sampling algorithm.

We apply the above methodology using max = 1 through 5 to data from GWTC-1 [33] , utilizing the samples from [52] and reweighting to use the population distribution as our prior. We down-sample to 3000 individual samples for each event to make the calculation more tractable. We apply the detection probability skymap calculated for the combination of O1 and O2; see Fig.  1 (c). The Bayes factors between anisotropy and isotropy (b m = 0 for all > 0) for the different models are reported in Table I .

We find the isotropic model is preferred over the anisotropic models with the Bayes factors ranging from ln BF = −0.45 for max = 1 to −0.12 for max = 5. We present the posterior distributions of a m /a 00 for ≤ 2 for the max = 2 anisotropy model in Fig. 2 . The real and imaginary components are shown for the m = 0 spherical harmonic coefficients. The posterior distributions for the other models are presented in Appendix A. The posterior distribution for all anisotropy models differ from the prior distribution, indicating that the ten events in GWTC-1 do provide some information about the overall sky distribution of events. However, the posterior distributions are still consistent with an isotropic universe (a m = 0 for ≥ 1). There this support for all a m = 0 present to at least the 95% level. We also present the posterior predictive distribution for the max = 2 anisotropic analysis in Fig. 2 . Although there is no preference for the max = 2 anisotropic model, the posterior predictive distribution illustrates what the anisotropy may look like with the anisotropic model if it was detected in GWTC-1. These results are qualitatively similar to Ref. [36] , where similar preferences for BBH sky positions are observed using a HEALPix [37] based model with 12 individual pixels, and three rotation angles. They report mild support for isotropy.

In this section, we assess how future observations will be able to probe anisotropy of compact binary mergers. We consider a network consisting of two aLIGO observatories and one aVirgo observatory, all operating at design sensitivity. We assume this network operates with a 100% duty cycle. The detection probability skymap is calculated following an identical procedure to the analysis of GWTC-1, except there is no downtime and we assume stationary Gaussian noise. We simulate gravitational-wave signals from binary black hole mergers with parameters distributed according to the aforementioned model over the course of one year of three-detector coincident operation.

We create simulated data for three angular distributions: isotropic, dipolar (a 11 /a 00 = 0.51), and quadrupolar (a 22 /a 00 = 0.51) distributions with all other a m = 0 for = 0. In each case, we draw 269 random binary black hole events. For the remainder of the manuscript, all analyses focus on the max = 5 model, as it is not currently computationally practical to include more than the 35 parameters required for max = 5. We recover the injected sky distribution, shown in Fig. 3 . The posterior predictive distribution in the right column closely matches the true distribution, with smaller angular resolution deviations from the true population distribution. These deviations are not significant. We report the Bayes factors for each of these scenarios in Table II .

In order to study the detectability of anisotropy with a large numbers of events, we simulate thousands of events with perfectly reconstructed sky position, from 100 different realizations of a universe. By analysing this idealized injection set, we can measure the "cosmic variance" in our measurements of the b m -the uncertainty arising from the finite number of detections. The total uncertainty, of course, includes contributions from the uncertain sky localization as well as cosmic variance. Thus, this study yields an optimistic view of what might be possible in the future.

In Fig. 4 , we plot ln BF ani iso as a function of the number of detections N . The shaded regions indicate 1σ confidence intervals. Since an isotropic model is a subset of the general anisotropic model, as more events are added, the Bayes factor can provide support for the model with a more compact parameter space, until enough evidence mounts for the anisotropic model. Therefore, we cannot use a method such as this to "prove" isotropy; rather we may provide an increasingly stringent limit on anisotropy. In Fig. 5 we plot the 99.7% upper limit on a m /a 00 found in each analysis as a function of the number of events N , (a m /a 00 ) UL . For a low number of events, (a m /a 00 ) UL slowly decreases with the number of events before reaching a consistent power law relation in the high event region. We find that (a m /a 00 ) UL scales as ∼ N −1/2 at high number of detected events (≥ 400 events).

Looking to the future, as the gravitational-wave network improves in sensitivity, the number of binary black hole events observed will significantly increase. Furthermore, the fraction of the those with improved sky localization will also improve [53] . This will allow for increasing improvements to the estimates of the anisotropy This results in notable features such as the a20 posterior peaking away from a20 = 0. These effects are suppressed with increasing max. While the corner plot does not meaningfully exclude zero, the posterior predictive distribution does present regions of higher probability as Ref. [36] has also observed. However, isotropy is not excluded by this result and the Bayes factor for max = 2 anisotropy does not prefer either model. The white contours in the posterior predictive distribution correspond to the 95% confidence intervals for the locations of all BBH merger events in GWTC-1.

of the binary black hole gravitational-wave sky location distribution. Given that upper limit in the estimate of anisotropy scales as ∼ N −1/2 (Fig. 5 ) in the regime of many events, we can estimate that within one year of operation of third-generation detectors, such as Cosmic Explorer and the Einstein Telescope [53, 54] , gravitational waves will be able to probe anisotropies at the ∼ 0.1% level. Coupled with the precise location of sources, this will enable a precise measurement of the anisotropy of the stellar mass binary hole distribution in the universe. We stress that the inference of (an)isotropy complements (but in no way replaces) the measurements of anisotropy from the cosmic microwave background, which is associated with a very different time in the history of the Universe. Gravitational-wave observations will provide an additional method to probe the large-scale distribution of stellar mass binary black hole mergers within the Universe.

This analysis provides an approach way to search for the anisotropy of the stellar mass binary black hole sky distribution. The results presented in Ref. [36] qualitatively agree for the predicted population distribution for the sky position of stellar mass binary black holes. Our results demonstrate similar probability overdensities with an independent method and parameterization of the anisotropic distribution. We also demonstrate the study of both simulated binary black hole mergers and point estimate sources to determine an approximate scaling relation of the inferred maximum anisotropy as a function of the number of detected BBH events. FIG. 3: Injected binary black hole merger sky distribution (left) and the associated recovered posterior predictive distribution (right), for (a) an isotropic universe, (b) a universe with an a11/a00 = 0.51 dipolar anisotropy, and (c) a universe with an a22/a00 = 0.51 quadrupolar anisotropy. The sky distribution is recovered from 269 events drawn from the population prior. All events are drawn such that they exceed the detection threshold of the three detector network. We determine the detection probability as a function of signal position to remove this selection bias. We recover an approximately similar sky distribution to the injected distribution, demonstrating the accurate recovery of anisotropy. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. Computing was performed on LIGO Laboratory computing cluster at Cali-fornia Institute of Technology. We would like to thank all of the essential workers who put their health at risk during the COVID-19 pandemic, without whom we would not have been able to complete this work. The bounded regions correspond to the 68% confidence intervals for the log Bayes factors recovered from point source estimates. We observe the ln BF ani iso trend for the simulated binary injections to approximately follow the point source estimates. Since isotropy is a compact subset of anisotropic models, the Bayes factor will prefer isotropy until enough events are observed. This is observed in the behavior of the a m/a00 = 0.10 results.

for all spherical harmonics up to max . Inferred maximum 99.7% upper limits on the anisotropy present in the simulated point source analyses given a number of events. The maximum anisotropy is presented for five different scenarios analysed for 100 different realizations with max = 5 anisotropic models. In the regime of a high number of detected events, the maximum anisotropy follows a ∼ N −1/2 relation. As the number of events increases, the maximum anisotropy measured plateaus towards the true anisotropy. 

Large-scale background temperature and mass fluctuations due to scale-invariant primeval perturbations

The LIGO Scientific Collaboration and the Virgo Collaboration)

Gruppentheorie und ihre Anwendung auf die Quantenmechanik der Atomspektren

Bayesian data analysis

6: Posterior distribution and posterior predictive distribution for GWTC-1 using a max = 1 anisotropy model

This implicitly assumes a uniform-in-log prior on the merger rate [43]. This is consistent with the population inference of GWTC-1 undertaken in Ref

We thank Isobel Romero-Shaw for sharing posterior samples for GWTC-1 from the Bilby GWTC-1 catalog paper [52] . The authors appreciate the helpful discussions with Colm Talbot and Vuk Mandic. We thank Marco Cavaglia and Aditya Vijaykumar for comments on the manuscript. This work is supported through Aus-