Hierarchical test of general relativity with gravitational waves

TL;DR

This paper introduces a hierarchical method for testing general relativity using multiple gravitational wave detections, which models deviations as a population distribution rather than assuming they are common or unrelated, improving robustness and sensitivity.

Contribution

The paper presents a novel hierarchical framework that models deviation parameters as a population distribution, extending and improving upon existing methods for testing gravity with gravitational waves.

Findings

Population parameters are consistent with general relativity.

The method is robust to measurement uncertainties.

Current data shows no deviation from general relativity.

Abstract

We propose a hierarchical approach to testing general relativity with multiple gravitational wave detections. Unlike existing strategies, our method does not assume that parameters quantifying deviations from general relativity are either common or completely unrelated across all sources. We instead assume that these parameters follow some underlying distribution, which we parametrize and constrain. This can be then compared to the distribution expected from general relativity, i.e. no deviation in any of the events. We demonstrate that our method is robust to measurement uncertainties and can be applied to theories of gravity where the parameters beyond general relativity are related to each other, as generally expected. Our method contains the two extremes of common and unrelated parameters as limiting cases. We apply the hierarchical model to the population of 10 binary black hole systems so far detected by LIGO and Virgo. We do this for a parametrized test of gravitational wave generation, by modeling the population distribution of beyond-general-relativity parameters with a Gaussian distribution. We compute the mean and the variance of the population and show that both are consistent with general relativity for all parameters we consider. In the best case, we find that the population properties of the existing binary signals are consistent with general relativity at the ~1% level. This hierarchical approach subsumes and extends existing methodologies, and is more robust at revealing potential subtle deviations from general relativity with increasing number of detections.

Equations

$$p(x\mid \vec{d})\propto p(\vec{d}\mid x)p(x),$$

$$p(\delta {\hat{p}}_i^{(j)}\mid {\vec{d}}^{(j)})$$

$$\propto p({\vec{d}}^{(j)}\mid \delta {\hat{p}}_i^{(j)}),$$

$$p({\mu }_i,{\sigma }_i\mid \{{\vec{d}}^{(j)}{\}}_{j=1}^N)$$

$$\propto p({\mu }_i,{\sigma }_i)\prod _{j=1}^{N}p({\vec{d}}^{(j)}\mid {\mu }_i,{\sigma }_i),$$

$$p({\vec{d}}^{(j)}\mid {\mu }_i,{\sigma }_i)=\int$$

$$p(\delta {\hat{p}}_i^{(j)}\mid {\mu }_i,{\sigma }_i)\mathrm{d}[\delta {\hat{p}}_i^{(j)}],$$

$$p(\delta {\hat{p}}_i\mid \{{\vec{d}}^{(j)}{\}}_{j=1}^N)=\int$$

$$p({\mu }_i,{\sigma }_i\mid \{{\vec{d}}^{(j)}{\}}_{j=1}^N)\mathrm{d}{\mu }_i\mathrm{d}{\sigma }_i,$$

Full text

I Introduction

The ever-increasing number of binary coalescences Abbott et al. (2018a) detected by LIGO Aasi et al. (2015) and Virgo Acernese et al. (2015) has opened up avenues for rich new tests of general relativity (GR) Abbott et al. (2016a); Yunes et al. (2016); Abbott et al. (2016b, 2017, 2018b, 2019). This includes precision probes of strong-field orbital dynamics, the nature of the remnant object, and the properties of gravitational-wave (GW) propagation Abbott et al. (2016a, 2019). With the new data, however, comes the problem of properly interpreting constraints in a way that does not apply only to specific modified theories of gravity and that is not biased by hidden assumptions Abbott et al. (2019); Zimmerman et al. (2019).

In particular, there is an outstanding challenge to adequately combine information from different GW observations into a single statement about agreement with GR. Existing approaches Li et al. (2012a, b); Agathos et al. (2014); Meidam et al. (2014); Del Pozzo et al. (2011); Ghosh et al. (2016); Abbott et al. (2016b); Meidam et al. (2018); Ghosh et al. (2018); Brito et al. (2018) rely on strong assumptions about the space of potential GR deviations and their effect on the observable events, rendering them too restrictive Zimmerman et al. (2019). As a result, we might soon have a wealth of measurements from different techniques and events but no cohesive picture that brings them together.

In this paper, we present a flexible and robust solution to this problem by framing it in the language of hierarchical inference. The result is an easy-to-interpret null-test of GR that can incorporate multiple measurements from different events, without strong restrictions to specific theories of gravity or subclasses of events, and without the need to explicitly weigh events based on their significance. We demonstrate that our method can produce strong combined constraints on deviations from GR. If deviations are present, it can detect them even if they affect our measurements nontrivially, e.g. by altering waveforms in ways that depend on the properties of each source. We apply our method to GW detections from the GWTC–1 catalog of compact binaries Abbott et al. (2018a, 2019), using publicly available posterior samples for parameters controlling waveform deviations from the GR prediction LIGO Scientific Collaboration and Virgo Collaboration (2019). We obtain joint constraints on deviations from GR that apply to generic theories of gravity, and find the data to be in agreement with Einstein’s theory up to the ${\sim}1\%$ level.

II Method

Waveform models for quasicircular compact binaries so far exist only within GR. They are generically parametrized by $15$ parameters that describe the signal observed by an interferometric detector: component masses, component spins, location, orientation, and phase. To date no parametrized waveform model exists that describes the inspiral, merger, and ringdown of generic binaries in any beyond-GR theory. For this reason, and guided by the desire for model-idependent tests that do not conform to specific theories, many studes are based on parametrized deviations away from the GR waveform. In these tests, new degrees of freedom $\delta\hat{p}_{i}$ are introduced, with $\delta\hat{p}_{i}=0$ corresponding to GR. These parameters are introduced to modify different aspects of the waveform’s frequency or amplitude evolution and, together with the $15$ GR parameters, define a generalized waveform model. For more details on these parametric tests, see the Supplemental Material and Abbott et al. (2019) and references therein.

Unless they are somehow fixed to a constant by the true theory of gravity, we should generally expect the $\delta\hat{p}_{i}$’s to vary across different GW events. For instance, the GW deformation could depend on the binary mass ratio or other properties of the system, and different combinations of $\delta\hat{p}_{i}$’s could come into play under different circumstances. Without assuming a theory of gravity, it is not possible to constrain the functional form of the $\delta\hat{p}_{i}$’s, making it difficult to combine measurements from different events Abbott et al. (2019); Zimmerman et al. (2019). To tackle this problem, we follow Zimmerman et al. (2019) and employ a hierarchical formalism wherein we assume that the true value of the beyond-GR parameters for each of the events is drawn from some common unknown distribution Loredo (2004). If there are $P$ parameters measured for $N$ events, this amounts to $P\times N$ random variables, which we denote $\delta\hat{p}_{i}^{(j)}$ for $i=1,\dots,P$ and $j=1,\dots,N$. Then each set of $N$ variables corresponding to a given $\delta\hat{p}_{i}$ should follow a shared distribution, implicitly determined by the underlying theory of gravity and the source population properties. The goal of the hierarchical approach (vividly named “extreme deconvolution” by some (Bovy et al., 2009)) is to infer the properties of the underlying distributions based on imperfect measurements from a population of events.

The first step is to select a functional form for the distribution of $\delta\hat{p}_{i}$, which is in principle nontrivial. Given the small number of detections, here we only attempt to measure its mean $\mu_{i}$ and standard deviation $\sigma_{i}$. Higher moments, such as the skewness, could become measurable with an increasing number of detections. In our case, and under a minimum-information assumption, we can thus model the population distribution with a Gaussian, i.e. we will take the population distribution to be $\delta\hat{p}_{i}\sim\mathcal{N}(\mu_{i},\,\sigma_{i})$. A more complex function could be chosen as needed, with little impact on the method. This potentially includes explicitly considering correlations among different $\delta\hat{p}_{i}$’s, although we demonstrate below that this is not strictly necessary.

With the above choice of likelihood and appropriate values of $\sigma_{i}$, our method reduces to traditional nonhierarchical approaches for combining events Zimmerman et al. (2019). Setting $\sigma_{i}=0$ amounts to assuming that all systems share the same beyond-GR parameter $\delta\hat{p}^{(j)}_{i}=\mu_{i}$. The results are equivalent to multiplying the likelihood functions of the $\delta\hat{p}^{(j)}_{i}$ for all detections $j$. On the opposite extreme, letting $\sigma_{i}\rightarrow\infty$, the $\delta\hat{p}^{(j)}_{i}$ are drawn from an effectively flat distribution and, as a result, measurement of one does not inform the others. This corresponds to testing a theory of gravity in which each system is described by its own fundamental constant Zimmerman et al. (2019). The results are equivalent to multiplying the Bayes factors from individual detections (assuming that a flat prior is imposed on each beyond-GR parameter). However, both these assumptions can lead to incorrect conclusions if they do not apply to the true theory of gravity Abbott et al. (2019); Zimmerman et al. (2019).

In its general form, our hierarchical method is not limited by those assumptions and provides a robust way of detecting a deviation from GR even when the non-GR parameters are not trivially related to each other. If GR is correct, then both hyperparameters, $\mu_{i}$ and $\sigma_{i}$, are expected to be consistent with zero. If we find a nonzero $\mu_{i}$, this is an obvious deviation from GR or a systematic error in the analysis of one or several of the events under consideration. Alternately, the true $\delta\hat{p}^{(j)}_{i}$ could be symmetrically distributed around $\mu_{i}=0$. In this case, the inferred $\mu_{i}$ will be consistent with [math], but the $\sigma_{i}$ posterior will peak away from zero, signaling that the scatter in $\delta\hat{p}^{(j)}_{i}$ is larger than expected from statistical measurement errors, again revealing a beyond-GR effect (or modeling error).

III Simulation: GR is right

Given the long history of GR’s experimental success Will (2014), it is unavoidable to imagine that GW observations may also fail to reveal any shortcomings of the theory. Accordingly, we begin by demonstrating our method on simulated signals that obey GR. For simplicity, we take the measurement of each beyond-GR parameter to be summarized by a Gaussian likelihood with mean $\tilde{\mu}^{(j)}_{i}$ and standard deviation $\tilde{\sigma}^{(j)}_{i}$, i.e. $p(\mathrm{data}^{(j)}\mid\delta\hat{p}_{i}^{(j)})=\mathcal{N}(\tilde{\mu}^{(j)}_{i},\tilde{\sigma}^{(j)}_{i})$. Such a likelihood is hardly realistic, especially for weak signals, but it suffices to illustrate our method and its scaling with the number of detections. Note that $\tilde{\mu}^{(j)}_{i}$ and $\tilde{\sigma}^{(j)}_{i}$ describe the idealized measurement of parameter $\delta\hat{p}_{i}$ in the $j$th event, while $\mu_{i}$ and $\sigma_{i}$ define the distribution of true values of $\delta\hat{p}_{i}$ across events.

We simulate a population of $N$ observations as follows: first, we assign a random signal-to-noise ratio (SNR) to each event $j$ with the expected probability $\mathrm{SNR}^{(j)}\sim 1/\rm{SNR}^{4}$ Chen and Holz (2014); then, for each $\delta\hat{p}_{i}^{(j)}$, we assign a value of $\tilde{\sigma}^{(j)}_{i}$ proportional to $1/\mathrm{SNR}^{(j)}$; finally, we choose a value of $\tilde{\mu}^{(j)}_{i}$ consistent with $\tilde{\sigma}^{(j)}_{i}$ by drawing it from ${\cal{N}}(0,\tilde{\sigma}^{(j)}_{i})$, mimicking the expected scatter due to noise in the detector. For concreteness, we consider only three non-GR parameters $\delta\hat{\varphi}_{i}$, $i=0,1,2$. These are defined as in Abbott et al. (2019) and are related to the parametrized post-Einsteinian (ppE) framework of Yunes and Pretorius (2009), as discussed in the supplement. We set the overall scale of the $\tilde{\sigma}^{(j)}_{i}$’s based on the uncertainty of measurements from GW150914 data, namely 68%-level widths of 0.06, 0.3 and 0.2 for $\delta\hat{\varphi}_{0}$, $\delta\hat{\varphi}_{1}$ and $\delta\hat{\varphi}_{2}$ respectively LIGO Scientific Collaboration and Virgo Collaboration (2019).

Figure 1 shows the projected constraints on $\mu_{i}$ (top) and $\sigma_{i}$ (bottom) for the ppE-like coefficients $\delta\hat{\varphi}_{0}$, $\delta\hat{\varphi}_{1}$ and $\delta\hat{\varphi}_{2}$ as the number of detections grows. Colored bands represent the $1\sigma$ variation over 200 simulated populations. The dashed line is proportional to $1/\sqrt{N}$ and demonstrates that bounds scale with the number of detections as expected. Our method improves with increasing number of signals at a rate similar to the simple approach of multiplying the likelihoods, in spite of the presence of an additional parameter, $\sigma_{i}$. This is because $\mu_{i}$ and $\sigma_{i}$ are uncorrelated, so we can safely add $\sigma_{i}$ to our model without affecting the $1/\sqrt{N}$ scaling of $\mu_{i}$, and vice versa.

IV Simulation: GR is wrong

We now turn to the tantalizing scenario that GR disagrees with experiment. In such a case, we should generally expect the deviation from GR to manifest itself in multiple $\delta\hat{p}_{i}$’s, even if it intrinsically occurs at a specific post-Newtonian (PN) order Sampson et al. (2013); Meidam et al. (2018). This is because the phenomenological effect of modifications at different PN orders are not necessarily orthogonal, introducing degeneracies in our measurement. Consequently, a deviation from GR affecting a given $\delta\hat{p}_{i}$ could be measured through the $\mu_{i}$ and $\sigma_{i}$ of multiple parameters, not just the one that is actually modified by the theory.

To demonstrate this effect, we construct a simple mock alternative theory of gravity that differs from GR at the 1PN order, affecting all binaries equally. This intrinsic waveform correction is independent of source parameters, making it amenable to multiplication of the individual parameter likelihoods. Generally, of course, this is not the case Cornish et al. (2011); Yunes et al. (2016). Even with this simplifying assumption, the measured $\delta\hat{p}_{i}$’s may vary in a nontrivial way with source properties as signals with different frequency contents may be affected by the same deviation differently.

Following Meidam et al. (2018), we assume that the measured non-GR parameters $\delta\hat{p}_{i}$ depend nontrivially on the true values $\delta\hat{p}^{\text{true}}_{i}$. Generally, such relation could always be expressed via some measurement matrix $M$, such that $\delta\hat{p}_{i}=M\delta\hat{p}^{\text{true}}_{i}$, where the components of $M$ could depend on the specific properties of each system. For our example, we again consider the three ppE-like parameters $\delta\hat{p}_{i}=(\delta\hat{\varphi}_{0},\,\delta\hat{\varphi}_{1},\,\delta\hat{\varphi}_{2})$ and we imagine $\delta\hat{p}^{\text{true}}_{i}=(0,\,0,\,0.1)$, i.e. the only parameter in which the modified theory deviates from GR is $\delta\hat{\varphi}_{2}$. As an illustration, we arbitrarily pick a matrix $M$ that yields $\delta\hat{p}_{i}=(1.1-2q,\,0,\,0.1)$, where $q$ is the mass ratio of the system. This is inspired by the degeneracy between high and low-order PN corrections demonstrated in Meidam et al. (2018). Quantitative results will be highly dependent on the true measurement matrix, though we only wish to demonstrate the qualitative effect here.

We simulate a population of observations by drawing $q$ uniformly from $[0.1,1]$, and using those values to produce the measured parameters $\delta\hat{p}_{i}$. To simulate the corresponding posteriors, we draw the event SNRs and add a scatter due to noise as in the previous section. As a result of the nontrivial dependence on $q$, the resulting population of each $\delta\hat{p}_{i}$ is not normally distributed. In spite of this, we demonstrate that our simple Gaussian model can detect the deviation from GR.

Figure 2 shows the posteriors for $\mu_{i}$ and $\sigma_{i}$ for a population of 100 events. As expected, we find that the posterior for $\mu_{2}$ peaks at the injected value of $0.1$ and excludes GR at the $96$% credible level. Additionally, we find that $\sigma_{0}$ is not consistent with GR at the ${\gtrsim}99.99\%$ credible level. This means that the scatter in $\delta\hat{\varphi}^{(j)}_{0}$ is too large to be accounted for by statistical noise. Indeed, part of the scatter in $\delta\hat{\varphi}^{(j)}_{0}$ is caused by the deviation from GR. This illustrates that, even if we did not take $\delta\hat{\varphi}_{2}$ into account, we would have detected this deviation from GR solely through the lower PN order coefficient. Additionally, the $\sigma_{0}$ posterior is farther from GR than the $\mu_{2}$ one, suggesting that this deviation could be detected first with a lower PN-order parameter. We emphasize that these results are illustrative only: the properties of the posteriors in a real analysis will depend on the nature of the true measurement matrix, which is generally unknown.

V Real Events

We now apply our hierarchical model to the confident binary black hole detections presented in GWTC–1 Abbott et al. (2018a). As a starting point, we use posterior samples for all $\delta\hat{p}_{i}$ parameters from Abbott et al. (2019); LIGO Scientific Collaboration and Virgo Collaboration (2019), obtained with the IMRPhenomPv2 waveform model Hannam et al. (2014); Khan et al. (2016). This study did not perform both sets of tests on all detected signals, but rather imposed certain thresholds on the SNR of the signals to determine whether to look for deviations in the inspiral or postinspiral regime, or both. As a result, $5$ signals where analyzed for inspiral deviations and $9$ for postinspiral ones. See Abbott et al. (2019) for details.

As demonstrated in the previous sections, for each ppE-like coefficient $\delta\hat{p}_{i}$, we obtain a posterior distribution for the corresponding hyperparameters $\mu_{i}$ and $\sigma_{i}$. We find that the population of the analyzed BBHs is consistent with GR both in terms of $\mu_{i}$ and $\sigma_{i}$ for all beyond-GR parameters. All $\mu_{i}$ posteriors are consistent with [math] at the 0.5$\sigma$ level or better, while all $\sigma_{i}$ posteriors peak at [math]. This is a novel test, sensitive to generic beyond-GR effects in a population of detections. With more signals we expect this analysis to either result in improved bounds (Fig. 1), or to reveal deviations from the ensemble properties expected from binaries in GR (Fig. 2).

From the hyperparameter posteriors, we also compute the inferred population distributions for the $\delta\hat{p}_{i}$’s, formally defined in Eq. (B) of the supplemental material and plotted in Fig. 3. These distributions, ${\rm d}N/{\rm d}(\delta\hat{p}_{i})$, represent our best knowledge of the population from which the allowed deviations from GR, $\delta\hat{p}_{i}^{(j)}$, were drawn for each binary black hole signal. Because all of these distributions contain zero with high probability, their width indicates the level to which we can constrain our measurements of the ppE-like coefficients to agree with GR. In the best case, for $\delta\hat{\varphi}_{-2}$, we find consistency with GR at the ${\sim}1\%$ level with $68\%$ credibility. In the future, and assuming GR is correct, we should find that the distributions tend to a $\delta$-function at $\delta\hat{p}_{i}=0$ as we accumulate more observations (cf. Fig. 1). We emphasize that, unlike Fig. 3 in Abbott et al. (2019), Fig. 3 here does not show the inferred posterior on the ppE-like parameters assuming all signals share the same value. Instead, Fig. 3 summarizes our inference for the distributions from which the potentially unequal ppE-like parameters of each signal were drawn 111This is not the same as the most likely such distribution, which by construction would be a Gaussian with mean and variance given by the peak of the posterior on $\mu_{i}$ and $\sigma_{i}$. Rather, Fig. 3 marginalizes over $\mu_{i}$ and $\sigma_{i}$, which is why the ${\rm d}N/{\rm d}(\delta\hat{p}_{i})$ distributions in Fig. 3 are not Gaussians. See supplement for details..

These results are subject to the thresholds imposed in Abbott et al. (2019) that determine which GW events are subject to each test. They would thus be vulnerable to the same potential selection effects. This includes the requirement that signals be sufficiently loud and akin to GR, such that they are detectable by matched-filtering procedures looking for GR signals. Reference Abbott et al. (2019) argues that both types of potential selection effects are partially mitigated by the fact that more generic searches are also employed alongside matched-filter ones. With this caveat in mind, we find no evidence of any deviation from GR.

VI Conclusions

We use a hierarchical approach to test GR with GWs by assuming that beyond-GR parameters in each event are drawn from a common underlying distribution. This approach is flexible and powerful, as it can encompass generic population distributions even if the chosen parametrization is inaccurate. It can trivially incorporate future detections and can be applied to different kinds of tests of GR, including searches for modified dispersion relations Mirshekari et al. (2012); Abbott et al. (2017) or inspiral-merger-ringdown consistency checks Ghosh et al. (2016, 2018). We apply this method to the current 10 confident binary black hole detections Abbott et al. (2018a), measuring posterior distributions for the mean and standard deviation of the population of ppE-like parameters $\delta\hat{p}_{i}$ LIGO Scientific Collaboration and Virgo Collaboration (2019). This is a conceptually new test that examines the ensemble properties of GW signals rather than the properties of individual events; we find the set of measurements to be consistent with GR (Fig. 3).

Parametrized tests, such as the ones studied here, are powerful probes of beyond-GR effects. Yet, it has long been appreciated that their interpretation demands caution: correlations between parameters make it necessary to have a consistent model to characterize a detected deviation. Our method provides a framework to execute a null test of GR with several detections, largely without the need for specific models of potential deviations. It improves with increasing number of signals at a rate similar to simpler approaches (Fig. 1). Furthermore, hierarchical methods could exploit degeneracies in our measurements to detect otherwise inaccessible deviations from GR, e.g. because they intrinsically occur at a higher PN order than can be directly probed (Fig. 2).

The framework presented here is not restricted to tests of GR with GWs, but can be generalized to include information from other observations. For example, the measured likelihood for $\delta\hat{\varphi}_{-2}$ from GWs could be combined with corresponding constraints from binary pulsar measurements. Our hierarchical method not only unifies the signals seen by ground-based detectors, but also offers a way to consider multiple tests of GR simultaneously.

Acknowledgements.

VII Acknowledgments

We thank Aaron Zimmerman and Carl-Johan Haster for useful discussions. We thank Nathan Johnson-McDaniel for comments on the draft. Samples from the $\mu_{i}$ and $\sigma_{i}$ posteriors were drawn with stan Carpenter et al. (2017), and plots were produced with matplotlib Hunter (2007). M.I. is supported by NASA through the NASA Hubble Fellowship grant #HST–HF2–51410.001–A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5–26555. The Flatiron Institute is supported by the Simons Foundation. This paper carries LIGO document number LIGO–P1900109.

Appendix A APPENDIX: TECHNICAL DETAILS

Appendix B Hierarchical inference

Our goal is to estimate the posterior probability density representing the measurement of some parameters of interest, conditional on some data or observation. In order to do this, it is necessary to write down a prior probability density function that encodes knowledge about the possible values of the parameters before obtaining the data. We use Bayes’ theorem to obtain the posterior from the likelihood, which is the probability of the data conditional on the values of the parameters. For a parameter $x$ and set of data $\vec{d}$, Bayes’ theorem is just

[TABLE]

relating the posterior (left) to the likelihood times the prior (right), with a proportionality constant that ensures unitarity.

The specific functional form of the prior is determined according to different expectations for the parameters. For example, in situations of high ignorance about the expected values, common choices for the prior are a uniform distribution across some broad range or a least-informative (aka, “Jeffrey’s”) prior. Alternatively, the prior may take the form of a multivariate normal distribution with some mean vector and covariance matrix. The mean informs the analysis of the likely location of the parameter values in the $P$-dimensional parameter space, and the covariance sets the scale of the corresponding uncertainty. Just as with the boundaries of the flat prior, the mean and covariance matrix can be set arbitrarily.

An extension to the above fixed priors is the case where we allow this distribution to vary in some controlled way. One such example is the case where we wish to marginalize over assumptions about the provenance of a parameter. Alternatively, as in the main text, besides the values of the parameters for any given observation, we might also be interested in the distribution of said parameter in a population of observed data. Both types of analyses make use of hierarchical models, with multiple nested levels of inference with corresponding priors. Priors on parameters that control distributions of other parameters (e.g. the mean and the covariance of the example priors above) are often called hyperpriors.

In this paper, we perform two inference levels: the first corresponds to the measurement of ppE-like parameters from individual events, and the second is the measurement of the properties of the population of ppE-like parameters. Concretely, we begin from the posteriors on $P$ parameters $\delta\hat{p}_{i}^{(j)}$ from the data $\vec{d}^{(j)}$ for the $N$ events. The posterior obtained from event $j$ for parameter $i$ takes the form:

[TABLE]

where the last line is only valid under the assumption of a uniform prior on $\delta\hat{p}_{i}^{(j)}$, as was applied in the original measurements by LIGO and Virgo Abbott et al. (2019); LIGO Scientific Collaboration and Virgo Collaboration (2019).

We then assume that the values for each event are drawn from a normal distribution with mean $\mu_{i}$ and standard deviation $\sigma_{i}$, characteristic of each parameter. Since the goal is to measure those quantities from the data from all events, for each $i$ we want the posterior

[TABLE]

where we took advantage of the fact that events are statistically independent to factorize the likelihood in the last line. Each of those factors may be written explicitly in terms of the parameters $\delta\hat{p}_{i}^{(j)}$,

[TABLE]

obtaining an expression in terms of the individual likelihoods, Eq. (B). These likelihoods, $p(\vec{d}^{(j)}\mid\delta\hat{p}_{i}^{(j)})$—or, rather, the corresponding posteriors with a suitable choice of prior—are what is computed in a regular (nonhierarchical) analysis with the goal of measuring $\delta\hat{p}_{i}^{(j)}$, as is the case in Abbott et al. (2019).

The last two equations allow us, then, to measure the population mean and standard deviation for the ppE-like parameters starting from the posterior samples released by the LIGO and Virgo collaborations LIGO Scientific Collaboration and Virgo Collaboration (2019), without re-analyzing the raw data (Fig. 3). Specifically, we obtain posteriors on the population parameters, Eq. (B), by sampling over $\mu_{i}$ and $\sigma_{i}$ with a Hamiltonian Monte Carlo algorithm implemented in the stan package Carpenter et al. (2017). To simplify the integration over the $\delta\hat{p}_{i}$’s in Eq. (B), we internally represent the LIGO-Virgo posteriors via a Gaussian-kernel density fit produced from the samples in LIGO Scientific Collaboration and Virgo Collaboration (2019).

From a posterior on the hyperparameters, $p(\mu_{i},\sigma_{i}\mid\{\vec{d}^{(j)}\}_{j=1}^{N})$, we may also infer the shape of the population distributions themselves. This can be done by marginalizing over $\mu_{i}$ and $\sigma_{i}$,

[TABLE]

where $p(\delta\hat{\varphi}_{i}\mid\mu_{i},\sigma_{i})={\cal N}(\mu_{i},\sigma_{i})$ by construction. In the main text we label the inferred distribution of Eq. (B) as “${\rm d}N/{\rm d}(\delta\hat{p}_{i})$,” since it can be interpreted as the expected fractional number of events ${\rm d}N$ with a value of the ppE-like coefficient $\delta\hat{p}_{i}^{(j)}$ within $[\delta\hat{p}_{i},\delta\hat{p}_{i}+{\rm d}(\delta\hat{p}_{i})]$ .

Appendix C Parametrized tests of GR

In this study, we consider parametric deformations to the GW signal described by some non-GR quantities $\delta\hat{p}_{i}$ indexed by $i$, with $\delta\hat{p}_{i}=0$ corresponding to GR. These new parameters are introduced on top of the 15 usual parameters that describe a GW within GR, with the goal of providing a model-independent framework with which to test the theory. The parameters are generally chosen such that they modify a specific aspect of the waveform, and hence a non zero value would signal a specific type of violation of GR. The parametric tests we will consider here are introduced to test the waveform in the different regimes of a binary coalescence: the inspiral, the merger, and the ringdown.

The inspiral phase is deformed through the post-Einsteinian (ppE) inspiral parameters $\delta\hat{p}_{i}=\delta\hat{\varphi}_{i}$ Yunes and Pretorius (2009). We define these $\delta\hat{\varphi}_{i}$’s as in Abbott et al. (2019), a choice that differs slightly from that of Yunes and Pretorius (2009) but which is analytically equivalent Yunes et al. (2016). The usual inspiral phase within GR is usually expressed as an expansion in small velocities, referred to as the post-Newtonian (PN) expansion. Each coefficient in the expansion depends on the system parameters and measurement of the GW phase evolution amounts to measuring said parameters. The ppE parameters are introduced as relative (or absolute, in some cases where the GR term vanishes) deformations of the normal PN coefficients such that $\delta\hat{\varphi}_{i}$’s encodes a correction at the $i/2$ PN order. For example, $\delta\hat{\varphi}_{-2}$ corresponds to a -1PN correction, associated with dipole radiation.

The merger-ringdown and the intermediate regimes are deformed through a different set of post-inspiral parameters, denoted $\delta\hat{\alpha}_{i}$ and $\delta\hat{\beta}_{i}$ respectively Meidam et al. (2018). These parameters encode modifications to the analytic description of those post-inspiral stages as implemented in the IMRPhenomPv2 waveform model Hannam et al. (2014); Khan et al. (2016). In particular, the $\delta\hat{\alpha}_{i}$’s control the merger-ringdown coefficients $\alpha_{i}$, which are obtained both from phenomenological fits and black-hole perturbation theory Khan et al. (2016). The $\delta\hat{\beta}_{i}$’s control deviations from the NR-calibrated phenomenological coefficients $\beta_{i}$ of the intermediate stage.