The directional isotropy of LIGO-Virgo binaries

We demonstrate how to constrain the degree of absolute alignment of the total angular momenta of LIGO-Virgo binary black holes, looking for a special direction in space that would break isotropy. We also allow for inhomogeneities in the distribution of black holes over the sky. Making use of dipolar models for the spatial distribution and orientation of the sources, we analyze 57 signals with false-alarm rates<1/yr from the third LIGO-Virgo observing run. Accounting for selection biases, we find the population of LIGO-Virgo black holes to be fully consistent with both homogeneity and isotropy. We additionally find the data to constrain some directions of alignment more than others, and produce posteriors for the directions of total angular momentum of all binaries in our set. All code and data are made publicly available in https://github.com/maxisi/gwisotropy/.


I. INTRODUCTION
With the increasing number of binary black holes (BBHs) detected by LIGO [1] and Virgo [2], it has become possible to study the distribution of such gravitational wave (GW) sources over time and space [3][4][5][6][7][8][9].Since BBHs can be detected up to nonnegligible redshifts (currently, 1), we expect their distribution at large scales to reflect the homogeneity and isotropy that characterize the universe cosmologicallya departure from that expectation would reveal a major shortcoming in our understanding of the detection biases affecting the LIGO-Virgo instruments, or, more tantalizingly, point to fundamentally new physics or astrophysics.
The homogeneity [6][7][8][9] and isotropy [10] of BBHs have been studied before under different frameworks.In this work, we reconsider the problem from a new point of view by quantifying the degree of alignment of the BBH orbits; in other words, we ask the question: could the total angular momentum vectors of LIGO-Virgo BBHs be preferentially aligned with a special direction in space?We consider this possibility as we simultaneously search for angular inhomogeneities in the spatial distribution of sources, thus constraining the existence of special directions controlling both the alignment and location of BBHs.
Unlike previous studies, we look for a breaking of isotropy in the angular momenta through the existence of a special direction in space with reference to some absolute frame, like the cosmic microwave background or far away stars (Fig. 1, second panel), and not with respect to Earth.The discovery of such a special direction could reveal the presence of a vector field breaking Lorentz symmetry.This differs from previous works like Ref. [10], which checked for anomalies in the alignment of sources with respect to Earth, as reflected in the distribution of BBH inclinations relative to the line of sight (Fig. 1, third panel).Such studies are not sensitive to the kind of overall alignment in absolute space that we look for here.
In Sec.II we describe the population model we use to constrain anisotropies, as well as our assumptions about the astrophysical distribution of BBH properties like masses and redshifts, and our treatment of selection biases.In Sec.III we outline the data products used in our analysis.We present our results in Sec.IV, including constraints on the degree of orientation and location inisotropies, as well as maps of possible preferred directions.In Sec.V, we summarize validation tests of our infrastructure to contextualize our measurement.We conclude in Sec.VI.Appendices show how to obtain posteriors on the angular-momentum direction, discuss hyperparameter prior choices and display posteriors for the angular-momentum direction for all events in our set.

A. Isotropy modeling
In order to study the spatial distribution and orientation of BBHs, we must analyze the collection of detections under a hierarchical population model that allows for variable degrees of spatial and directional correlations (see, e.g., Ref. [9]).This requires modeling the distribution of source locations, N , and total-angular-momentum directions, Ĵ.For modeling purposes, the N and Ĵ vectors must be expressed through their Cartesian compo-FIG.1. BBH orientation models.By default, we expect BBH angular momenta (arrows) to be oriented randomly with respect to Earth (circle) or each other, reflecting isotropy (first panel).In this study, we consider the possibility that BBH orbits follow a special direction in space, the extreme of which is full alignment (second panel).Previous studies, like Ref. [10], have considered models in which binaries are (or are perceived to be) aligned anomalously with respect to Earth, e.g., pointing preferentially towards it (third panel).The first two panels both have a distribution of inclinations that looks isotropic to analyses like Ref. [10].
nents in some absolute reference frame.
We choose to work in geocentric, equatorial celestial coordinates, with the vernal equinox as the x axis.In that frame, the sky-location vector is, for each BBH, N = (cos α cos δ, sin α cos δ, sin δ) , for the right ascension α and declination δ; meanwhile, the orientation vector is Ĵ ≡ J/| J| where for the orbital angular momentum L and individual (dimensionful) black hole (BH) spin angular momenta S 1/2 .The Cartesian components of Ĵ can be computed in the above frame as a function of α, δ, and the polarization angle ψ [11], as well as the component masses m 1/2 , the dimensionless spin vectors χ 1/2 and the orbital phase φ ref at some reference frequency f ref ; we outline this calculation in Appendix A and release relevant code in [12].We prefer to work with Ĵ rather than L because the former is conserved over the coalescence to a high degree, even for precessing systems [13].
Our goal is to quantify the degree of isotropy in N and Ĵ.To this end, we model each of those vectors as drawn from an isotropic distribution with a dipolar correction of variable magnitude.This ad hoc modification may be generally seen as the first term in a harmonic expansion around isotropy, and is specifically well suited to capture the existence of a preferred direction in space.The dipolar components in the location and orientation distributions are defined by dipole vectors v N/J whose magnitudes control the degree of deviation from isotropy for N and Ĵ respectively.
Concretely, we implement the following population-level likelihood for each event (indexed by i) for some special-direction vectors v N/J , with 0 ≤ v N/J < 1, to be inferred as hyperparameters from the collection of detections.The fully isotropic case is recovered for v N = v J = 0. On the other hand, setting v J = 0 alone allows for a nonhomogenous distribution of sources in the sky while assuming isotropic source orientations; this reduces to the "dipole" model in [9].Modeling the likelihood as in Eq. (3), with v N/J < 1, ensures that the likelihood itself remains positive everywhere.
It is convenient to rephrase the constraint that 0 ≤ v N/J < 1 by re-parameterizing the population in Eq. (3), through two corresponding, auxiliary vectors u N/J defined such that In this way, we ensure 0 ≤ v N/J ≤ 1 for −∞ < u N/J < ∞; u N/J are unconstrained.This "decompactifying" transformation is more effective for sampling purposes than enforcing a sharp constraint on the magnitudes of v N/J .For small v N/J 1, v N/J u N/J to second order in v N/J .As a prior, we choose a three-dimensional Gaussian distribution on the components of each u N and u J , with zero mean and standard deviation σ = 0.4.This choice of prior is designed to be fairly uninformative about v N/J (i.e., lacking a strong gradient) while still peaking at v N/J = 0; see Appendix B for a description of this feature.

B. Reweighting to an astrophysical population
Besides the location and orientation modeling described above, we need to ensure that our assumptions about the parameters of each individual BBH, like masses and redshift, are astrophysically sensible.To that end, we assume a redshift distribution that follows the Madau-Dickinson star formation rate [14] in the comoving frame.For the masses, we assume a prior distribution inversely proportional to the heaviest component mass (∝ 1/m 1 ) and uniform in the mass ratio (constant in q = m 2 /m 1 ), and restrict our sample to BHs with (posterior-median) masses in the range 5 M < m 2 ≤ m 1 < 150 M ; within that range, this choice is a simple approximation to the measurement in [3].We assume the population of component spins are isotropically oriented with respect to the orbital angular momentum, with a uniform distribution over their dimensionless magnitude. 1A future analysis may fit the astrophysical distribution of these parameters jointly with the orientation and location of the binaries.

C. Selection biases
Since we are attempting to model the intrinsic distribution of all BBH sources, not merely those that were detected, we must account for the difference in LIGO-Virgo's sensitivity to various sources.This is true for both intrinsic parameters, like BH masses and spin magnitudes, as well as the location and orientation parameters in which we are interested for this work (namely, N and Ĵ).With knowledge of the instruments' sensitivity over parameter space, we use the measured selection function to obtain a measurement of the intrinsic distribution of parameters out of the distribution of detected sources [15,16].
Evaluating the detectors' sensitivity over parameter space requires large simulation campaigns that quantify the end-to-end performance of LIGO-Virgo detection pipelines by injecting and recovering synthetic signals.As in [9], we take advantage of the BBH dataset in [17] for this purpose. 2Since this injection campaign only covered LIGO-Virgo's third observing run, we only consider events detected during that period; together with the mass constraints cited above, this means there are 57 BBH events to be included in our analysis (listed in Table I).

III. DATA
Our analysis starts from posterior samples for individual events reported by LIGO-Virgo in [18,19] and 1 This implies that our modeling of Ĵ following Eq.( 3) can be reinterpreted as a nontrivial modeling of L through Eq. ( 2). 2 Specifically, the endo3 bbhpop-LIGO-T2100113-v12 injection set.publicly released in [20,21] through the Gravitational Wave Open Science Center [22][23][24].Specifically, we make use of results obtained with the IMRPhenomXPHM waveform [25][26][27][28] that have been already reweighted to a distance prior uniform in comoving volume.The single-event inference was carried out by the LIGO-Virgo collaborations using the Bilby parameter estimation pipeline [29,30], as detailed in [18,19].We reweight those samples following the astrophysical population described above (see, e.g., [31,32] for methodology) and then use them to produce distributions for the components of N and Ĵ, which we use as input for our hierarchical analysis based on Eq. (3).
The six-dimensional posteriors for the components of Ĵ and N for each event are the primary input for our hierarchical isotropy analysis.We show the posteriors for Ĵ in Fig. 8 as skymaps for all events in our set, produced using standard LIGO-Virgo tools for representing probability densities over the sky [33][34][35]; we make these posteriors, resulting from the calculation detailed in Appendix A, available in our data release [12].The equivalent figures for N are nothing but the sky-localization maps already made available by LIGO-Virgo [20,21].

IV. RESULTS
We showcase the full result of our analysis in Fig. 2, which represents the simultaneous measurement of location and orientation anisotropies through the sixdimensional posterior on the components of v N and v J .The result is fully consistent with both kinds of isotropy, with v N = 0 and v J = 0 supported with high credibility, falling close to the peak of the marginal distributions on Isotropy measurement.Result of the simultaneous measurement of location and orientation isotropy through the model in Eq. ( 3), as represented by the posterior distribution on the dipole vectors v N/J (corner plot), and the corresponding projections over the sky (Mollweide insets).The six-dimensional posterior distribution is represented through credible levels over two-dimensional slices (blue contours, spaced at intervals corresponding to 10% increments in probability mass, with the outer contour enclosing 90% of the probability), and one-dimensional marginals (diagonal).The upper-left and lower-right subcorners encode constraints on the individual components of each vN and vJ respectively (highlighted with vertical and horizontal lines in the margin), while the other panels encode potential correlations between the location and orientation inisotropies.The measurements for v N/J can be projected into distribution over the sky as in the top-right insets, which show the allowed dipole orientations for vN ≡ vN /| vN | (top) or vJ ≡ vJ /| vJ | (bottom), with lighter colors encoding higher probability density over the celestial sphere [33][34][35]; inhomogeneities in these sky-maps do not constitute evidence for anisotropies.Isotropy is recovered for vN = vJ = 0 (dotted lines), which is well supported by the 6D posterior.3. 3D distributions.Three-dimensional representation of the v J/N measurement in Fig. 2 (first two panels), in comparison to the prior (last panel).Each point is drawn from the corresponding three-dimensional distribution, with color proportional to the probability density (lighter colors for higher density).The origin, representing isotropy, is well favored in all cases (intersection of gray dashed lines).The posteriors are tighter with respect to the prior, as is also seen in Fig. 4.
the 1% and 49% quantiles respectively; 3 this feature is also reflected in the fact that the origin is well supported in all the panels of the corner plot in Fig. 2. The result for v N is consistent with previous studies [9], which did not find evidence against isotropy in the location of LIGO-Virgo sources.
To the extent that there is any support for nonzero dipolar contributions to the location or orientation densities (namely, for | v J/N | > 0), their possible directions in the sky are represented by the insets on the top right of Fig. 2: a higher density indicates a potentially allowed orientation for the v N (top) or v J (bottom) dipole vectors.These skymaps reveal that the data are not in conflict with the existence of a weak dipole along the vernal equinox in the celestial equatorial plane (the direction of the x axis in our Cartesian coordinate system), as implied by the marginal on v J,x in Fig. 2; this dipole is allowed, but not required, by the data, since the posterior is fully consistent with v J = 0.
Indeed, although inhomogeneities appear in these Mollweide projections, this should not be interpreted as evidence for anisotropies: the density of points in those maps only encodes permissible directions for the dipole, without implications for its magnitude.In fact, inhomogeneities will appear in such plots any time the posterior does not happen to peak exactly at the three-dimensional origin of v J/N , as we expect to be commonly the case even if v J/N = 0 is the underlying truth. 4 In Fig. 3, we provide an additional representation 3 These three-dimensional quantiles correspond to the fraction of v J/N samples with higher probability density than the origin. 4With a finite number of events in the catalog, even if v J/N is zero in truth, the maximum-likelihood estimate will not, in general, of this posterior projected onto the three-dimensional spaces of v J and v N ; we also present the prior for comparison.Large dipolar contributions are disfavored by the posterior in all cases, and the posterior distributions are more concentrated than the prior, indicating that the data are informative.The posterior standard deviations for the three Cartesian components of v J and v N are smaller by {15%, 10%, 9%} and {25%, 20%, 15%} respectively with respect to the prior.The stronger tightening be zero.The posterior will therefore peak away from zero (but be fully consistent with zero).In this situation inhomogeneity can appear as in Fig. 2.
in the v N distribution is a consequence of N being better constrained than Ĵ in individual events.Figure 3 again makes it apparent that data disfavor certain directions for the v J dipole (towards the positive-x quadrants, away from the vernal equinox).We may translate the result in Fig. 2 into a posterior on the magnitudes | v J/N | of the dipole components, as we do in Fig. 4.However, the one-dimensional posteriors on these quantities is heavily dominated by the dimensionality of the problem, which results in a Jacobian disfavoring small values of | v J/N | due to the limited phase space near v J = 0 or v N = 0.5 Therefore, even though our three-dimensional prior does not treat v J/N = 0 as a special point (Appendix B), the effective prior induced on the magnitudes heavily disfavors the origin in v J/N due to the reduced available prior volume (gray curve in Fig. 4); that explains why the posteriors in Fig. 4 themselves appear to disfavor | v J/N | = 0.With that in mind, the influence of the data can be seen in the leftward shift of the | v J/N | distributions with respect to the prior; this is effect is more pronounced for v N , which is a consequence of the fact that N is generally better measured than J.Although the shift in | v J | is slight, the data are informative about v J -constraining some of its possible orientations (Figs. 2 and 3), if not its overall magnitude.

V. VALIDATION
In this section, we validate our setup by studying simulated datasets in which Ĵ and N are isotropically distributed (V A).We also revisit our assumptions about the astrophysical distribution of BBH parameters, described in Sec.II B, and show that they are robust.

A. Injections from selection set
We validate our infrastructure on the set of injections used to evaluate the selection function of the instruments, which was drawn from an intrinsically isotropic distribution [17].This provides an end-to-end test of our setup, including the computation of location and orientation vectors, the selection function and the inference process.
Concretely, we simulate catalogs of detections drawn from the injection set used to evaluate the selection function as described in Sec.II C [17].For simplicity, we treat the injection parameters as a single sample of a fictitious posterior: the input to our hierarchical analysis is one sample per synthetic event, drawn from the distribution in Sec.II C. At each iteration, we double the size of the simulated catalog.Since the injection distribution was constructed to be isotropic [17], we expect this experiment to indicate v J/N = 0, with certainty growing with catalog size.We show the result in Fig. 5 for v J (the result for v N is similar).As expected, v J/N = 0 is supported with increasing certainty as the synthetic catalog grows.This is what we expect if our infrastructure is working as designed.

B. Astrophysical distributions
As described, in Sec.II B, analysis above (Figs.2-4) hinged on simplified assumptions about the astrophysical distribution of BH masses, spins and redshifts.To evaluate the impact of this simplification, we repeat our analysis but now reweighting to a different astrophysical population based on the measurements in [3].Specifically, we make use of the highest-probability instantiation of the Power Law + Peak parametric mass model and the vN,y Effect of astrophysical population.We repeat the isotropy analysis with a different assumption about the underlying distribution of BBH parameters based on the measurement in [3].This yields the result shown in black, as opposed to the result in Fig. 2, which is reproduced here in blue for comparison.The black contours enclose 90%, 50% and 10% of the marginal probability mass.
Default spin model, with parameters obtained from the posterior samples released in [36].This model treats the astrophysical distribution of BH masses as a power law plus a Gaussian peak, with density evolving over comoving volume as a power law; the spins are isotropic with a possible Gaussian overdensity in alignment with the orbital angular momentum, and with magnitudes following a Beta distribution (see [3] for details).To implement this new astrophysical prior, both individual detections and the selection injections are reweighted accordingly.
The result of assuming this different astrophysical distribution is shown in Fig. 6, compared to the main result above.Although the new posterior differs slightly from our primary one in Fig. 2, as we might expect given the difference in models, the change is quite limited and does not qualitatively impact the discussion above.The discrepancy is somewhat more pronounced for the components of v J , as we might expect from the fact that the reconstruction of Ĵ must factor our inference on the masses and three dimensional spins of each BBH.In the future, a more comprehensive analysis might simultaneously measure v J/N and the distribution of astrophysical properties.

VI. CONCLUSION
We have demonstrated a measurement constraining the potential alignment of the total angular momenta of BBHs detected by LIGO and Virgo, using 57 detections from their third observing run and duly accounting for selection effects.In addition to alignment of momenta, we simultaneously looked for inhomogeneities in the distribution of sources over the sky.We found no evidence against isotropy in either the orientation or, consistent with previous works, the location of LIGO-Virgo BBHs.Additionally, we determined that the GWTC-3 data disfavor certain orientations of the potential preferred alignment of angular momenta more than others.
Future measurements will improve as the LIGO, Virgo and KAGRA detectors grow in sensitivity, resulting in many more detections at higher signal-to-noise ratios, which will result in more precise isotropy constraints.The advent of next generation detectors, like Cosmic Explorer [37][38][39] or the Einstein Telescope [40], will enable the exploration of higher order anisotropies and other interesting effects, like correlations with redshift or, potentially, local correlations in the directions N and Ĵ on the sky.Prior.Two-dimensional slice of the Gaussian prior on the components of v J/N for σ = 0.4, as used for the analysis in the main paper.With this choice, the origin is not disfavored by the three-dimensional prior probability density.on u N/J induces a prior on v N/J that is The derivative of the density with respect to v N/J vanishes at the origin, as it must by symmetry.But the second derivative does not, and is Thus, such a Gaussian prior will have a maximum at v N/J = 0 only if σ < 1/ √ 5 0.45; otherwise it places too much prior mass on large u N/J which generates a ring-shaped maximum in the prior for some v N/J > 0 and a minimum at the origin.With the desire to be uninformative about the typical scale of the components of v N/J while keeping the maximum prior density at the origin, we choose σ = 0.4 < 1/ √ 5 for the analysis in this paper.A two-dimensional slice of the prior on the components of v N/J for this choice is illustrated in Fig. 7.
Appendix C: Direction of total angular momentum for individual events In Fig. 8, we present posteriors on the direction of the total angular momentum Ĵ for all 57 events in our set.We produced these posteriors by applying the calculation described in Appendix A to the samples released by LIGO-Virgo [20,21], after reweighting as outlined in Sec.II B.
The skymaps in the figure provide a Molleweide projection of probability density over the celestial sphere, in the standard equatorial, geocentric coordinates used in the main text.A darker color corresponds to a direction in the sky with more probability density.From these maps, it is clear that not all events are equally informative about Ĵ; better constrained events will tend to dominate our hierarchical measurement.
The skymaps were produced using the ligo.skymappackage [33][34][35], a set of standard LIGO-Virgo tools for the processing of probability densities over the sphere.We make these figures, corresponding Flexible Image Transport System (FITS) files, and code used to generate them, available in our data release [12].FIG. 8. Measurements of the total angular momentum direction, Ĵ, for the events in our set, in a Mollweide projection of Earth-centric Celestial coordinates; darker color represents higher probability density for that direction in space (cont.).
FIG.3.3D distributions.Three-dimensional representation of the v J/N measurement in Fig.2(first two panels), in comparison to the prior (last panel).Each point is drawn from the corresponding three-dimensional distribution, with color proportional to the probability density (lighter colors for higher density).The origin, representing isotropy, is well favored in all cases (intersection of gray dashed lines).The posteriors are tighter with respect to the prior, as is also seen in Fig.4.

FIG. 4 .
FIG. 4. Posterior on dipole magnitudes.The measurement of Fig. 2 translated into one-dimensional posteriors on the | vJ | (blue) and | vN | (orange) magnitudes.These are shifted rightward with respect to the implied prior (gray), which itself heavily disfavors | v J/N | = 0 due to the reduced phase space near the origin in a three-dimensional space.

FIG. 5 .
FIG.5.Validation on synthetic catalogs drawn from selection injection set.Result of hierarchical analyses on synthetic catalogs of increasing size N (color), obtained by taking draws from the injection set used to evaluate the selection function (Sec.II C) and using them as single samples from the Ĵ and N posterior of synthetic events.We show the posterior on the components of vJ (main corner), and the implied posterior on the magnitude | vJ | (upper right).A lighter color corresponds to a larger catalog: starting with N = 4 events for the darkest color and progressively doubling the catalog size 11 times to reach N = 8192 events for the lightest color.The hierarchical analysis measures vJ = 0 with growing precision as the size of the catalog increases.
FIG. 7.Prior.Two-dimensional slice of the Gaussian prior on the components of v J/N for σ = 0.4, as used for the analysis in the main paper.With this choice, the origin is not disfavored by the three-dimensional prior probability density.

FIG. 8 .
FIG.8.Measurements of the total angular momentum direction, Ĵ, for the events in our set, in a Mollweide projection of Earth-centric Celestial coordinates; darker color represents higher probability density for that direction in space.

TABLE I .
Events considered.