On combining information from multiple gravitational wave sources

TL;DR

This paper critically examines methods for combining multiple gravitational wave observations to test general relativity, highlighting their assumptions and limitations, and proposing a hierarchical inference framework for more robust analysis.

Contribution

It reveals the assumptions behind common combination methods and introduces a hierarchical inference approach to improve tests of gravity theories using multiple sources.

Findings

Likelihood multiplication assumes identical beyond-GR parameters across events.

Bayes factor multiplication assumes independent new parameters per event.

Hierarchical inference offers a more flexible framework for combining data.

Abstract

In the coming years, advanced gravitational wave detectors will observe signals from a large number of compact binary coalescences. The majority of these signals will be relatively weak, making the precision measurement of subtle effects, such as deviations from general relativity, challenging in the individual events. However, many weak observations can be combined into precise inferences, if information from the individual signals is combined in an appropriate way. In this study we revisit common methods for combining multiple gravitational wave observations to test general relativity, namely (i) multiplying the individual likelihoods of beyond-general-relativity parameters and (ii) multiplying the Bayes Factor in favor of general relativity from each event. We discuss both methods and show that they make stringent assumptions about the modified theory of gravity they test. In particular, the former assumes that all events share the same beyond-general-relativity parameter, while the latter assumes that the theory of gravity has a new unrelated parameter for each detection. We show that each method can fail to detect deviations from general relativity when the modified theory being tested violates these assumptions. We argue that these two methods are the extreme limits of a more generic framework of hierarchical inference on hyperparameters that characterize the underlying distribution of single-event parameters. We illustrate our conclusions first using a simple model of Gaussian likelihoods, and also by applying parameter estimation techniques to a simulated dataset of gravitational waveforms in a model where the graviton is massive. We argue that combining information from multiple sources requires explicit assumptions that make the results inherently model-dependent.

Figures

Equations

$$\mathcal{B}$$

$$p(\mathbf{\lambda }|D,H_2)$$

$$=$$

$$\mathcal{B}=\frac{p(\mathbf{\lambda }=0|D,H_2)}{p(\mathbf{\lambda }=0|H_2)}.$$

$$p(D|\lambda )$$

$$\mathcal{B}$$

$$p({\lambda }_1,{\lambda }_2|D)$$

$$=\frac{p({\lambda }_1,{\lambda }_2)p(D_1|{\lambda }_1)p(D_2|{\lambda }_2)}{\int p({\lambda }_1^{\prime },{\lambda }_2^{\prime })p(D_1|{\lambda }_1^{\prime })p(D_2|{\lambda }_2^{\prime })d{\lambda }_1^{\prime }d{\lambda }_2^{\prime }}$$

$$=\frac{p({\lambda }_1,{\lambda }_2)p(D_1|{\lambda }_1)p(D_2|{\lambda }_2)}{\int p({\lambda }_1^{\prime })p({\lambda }_2^{\prime }|{\lambda }_1^{\prime })p(D_1|{\lambda }_1^{\prime })p(D_2|{\lambda }_2^{\prime })d{\lambda }_1^{\prime }d{\lambda }_2^{\prime }},$$

$$p({\lambda }_2|{\lambda }_1)=\delta ({\lambda }_1-{\lambda }_2).$$

$$p(\lambda |D)$$

$$=\int \frac{p({\lambda }_1,{\lambda }_2)p(D_1|{\lambda }_1)p(D_2|{\lambda }_2)}{\int p({\lambda }_1^{\prime })p({\lambda }_2^{\prime }|{\lambda }_1^{\prime })p(D_1|{\lambda }_1^{\prime })p(D_2|{\lambda }_2^{\prime })d{\lambda }_1^{\prime }d{\lambda }_2^{\prime }}d{\lambda }_2$$

$$=\frac{p(\lambda )p(D_1|\lambda )p(D_2|\lambda )}{\int p({\lambda }^{\prime })p(D_1|{\lambda }^{\prime })p(D_2|{\lambda }^{\prime })d{\lambda }^{\prime }}.$$

$$p(D_i|\lambda )$$

$${\sigma }_f^2$$

$$\mathcal{B}$$

$$p(\lambda |D)\propto p(\lambda )\prod _{i}p(D_i,\lambda ).$$

$${\sigma }^2$$

$$\ln C$$

$$\mathcal{B}$$

$$p({\lambda }_2|{\lambda }_1)=p({\lambda }_2),$$

$$p({\lambda }_1,{\lambda }_2|D)$$

$$=\frac{p({\lambda }_1)p({\lambda }_2)p(D_1|{\lambda }_1)p(D_2|{\lambda }_2)}{\int p({\lambda }_1^{\prime })p({\lambda }_2^{\prime })p(D_1|{\lambda }_1^{\prime })p(D_2|{\lambda }_2^{\prime })d{\lambda }_1^{\prime }d{\lambda }_2^{\prime }}$$

$$=\frac{p({\lambda }_1)p(D_1|{\lambda }_1)}{\int p({\lambda }_1^{\prime })p(D_1|{\lambda }_1^{\prime })d{\lambda }_1^{\prime }}\frac{p({\lambda }_2)p(D_2|{\lambda }_2)}{\int p({\lambda }_2^{\prime })p(D_2|{\lambda }_2^{\prime })d{\lambda }_2^{\prime }}.$$

$$\mathcal{B}$$

$$p({\lambda }_i|M,\alpha )$$

$$p({\lambda }_i)p(D_i|{\lambda }_i)$$

$${\tilde{\sigma }}_i^2$$

$$\ln {\tilde{C}}_i$$

$$\int p({\lambda }_i)p(D_i|{\lambda }_i)d{\lambda }_i=\frac{1}{{\sigma }_i}\frac{{\tilde{\sigma }}_i{\tilde{C}}_i}{\sqrt{2\pi }\alpha }.$$

$$\mathcal{B}$$

Full text

I Introduction

Two years after the first direct detection of gravitational waves (GWs) from the coalescence of two black holes (BHs), GW astronomy has transitioned into a booming field with eleven confirmed detections to date Abbott et al. (2018a). The rate of detections is only expected to increase in the coming years as the GW detector network becomes more sensitive and grows in number Abbott et al. (2018b). This has lead to increasing interest in the problem of combining information from multiple GW sources to strengthen inferences that might be inconclusive from single events. Examples include tests of general relativity (GR) Abbott et al. (2016), the measurement of the Hubble constant Chen et al. (2018); Feeney et al. (2018), as well as inference about population properties Abbott et al. (2018a, c).

Within the framework of Bayesian Inference, combining information from multiple sources is usually formulated in two complementary ways. The first involves computing the posterior probability distribution of parameters that characterize the entire population, such as the coupling constant of a beyond-GR theory or the shape of the mass distribution of stellar BHs. The second amounts to computing the Bayes Factor (BF), the ratio of the evidences of two competing models that describe the population.

In the context of testing GR with GWs, two approaches are commonly used. The first involves computing the BF in favor of GR compared to beyond-GR theories for each individual signal and then multiplying the individual BFs Li et al. (2012a, b); Agathos et al. (2014); Meidam et al. (2014). The second approach involves computing the posterior probability distribution of some parameter that quantifies the deviation from GR and then multiplying the individual likelihoods together in order to draw stronger inferences on this parameter Del Pozzo et al. (2011); Meidam et al. (2014); Ghosh et al. (2016); Abbott et al. (2016); Ghosh et al. (2018); Brito et al. (2018). This parameter might be a fundamental constant of the theory that is common between all binary systems, such as the Compton wavelength of the graviton, or a parameter that is expected to depend on the specific properties of each binary system, such as the modifications to the post-Newtonian expansion introduced within the parametrized post-Einsteinian (ppE) framework Yunes and Pretorius (2009); Chatziioannou et al. (2012).

In this study, we revisit both approaches in order to discuss the underlying assumptions about the properties of the beyond-GR theory they test. Unsurprisingly, we find that multiplying the individual likelihoods of a beyond-GR parameter is appropriate only for theories where this parameter is identical for every member of the population. Meanwhile, multiplying the BFs is appropriate only when each the beyond-GR parameter is different for every new system observed. Through a simple model of Gaussian likelihoods we show that both approaches lead to faulty conclusions if applied to a theory of gravity that does not follow their respective assumptions. In particular, we may fail to detect deviations from GR when using incorrect assumptions about the underlying modification to GR.

The two assumptions above, namely that parameters are equal or completely uncorrelated are arguably extreme assumptions that generic theories of gravity are not required to obey. Instead it is more reasonable to expect the middle ground of parameters that are related to each other through an underlying population distribution of the beyond-GR parameters. This scenario is akin to, for example, the expectation that the masses of stellar-mass BHs in binaries are drawn from an underlying population mass distribution. We show that hierarchical analyses where one parametrizes and infers the underlying distribution can be interpreted as an intermediate case that encompasses the two extremes of multiplying likelihoods or BFs.

Our results raise the issue of the model-dependency of combining information from multiple signals, something particularly critical for tests of GR with GWs. We argue that even though model-independent tests are feasible for single sources, this notion of model-independent tests does not carry over to multiple detections. The rest of the paper presents the details of our argument. Section II studies the assumptions inherently made when combining inferences and illustrates these ideas using a simple model of Gaussian likelihoods. Section III presents a concrete example in terms of measuring the Compton wavelength of the graviton. Section IV summarizes our main conclusions.

II Assumptions when combining inferences

In this section we outline assumptions made when using two common methods for combining inferences from multiple detections in order to test deviations from relativity. Each of these is an extreme case of a more general framework, where the parameters controlling the deviation of a set of signals from the predictions of GR are drawn from some common (and generally unknown) distribution. In order to better illustrate these ideas and how constraints of modified theories scale with the number of signals analyzed, we use a simple model of Gaussian likelihoods. We then illustrate how the two common methods of combining inferences can fail to detect deviations from relativity when synthetic signals violate the underlying assumptions of each method.

II.1 Model comparison

Consider two competing hypotheses $H_{1}$ and $H_{2}$ that we wish to test with data $D$. The BF in favor of $H_{1}$ is given by

[TABLE]

where $p(D|H_{i})=\int p(\bm{\theta}_{i})p(D|\bm{\theta}_{i},H_{i})d\bm{\theta}_{i}$ is the evidence for $H_{i}$, $\bm{\theta}_{i}$ are the model parameters, $p(\bm{\theta}_{i})$ is the prior on the parameters, and $p(D|\bm{\theta}_{i},H_{i})$ is the likelihood.

We assume that $H_{1}$ and $H_{2}$ are nested models, i.e. $\bm{\theta}_{2}=\{\bm{\theta}_{1},\bm{\lambda}\}$ and $p(D|\bm{\theta}_{1},H_{1})=p(D|\bm{\theta}_{1},\bm{\lambda}=0,H_{2})$. The posterior for the extra parameters $\bm{\lambda}$ is

[TABLE]

and using the Savage-Dickey density ratio Dickey (1971) the BF simplifies to

[TABLE]

From here on we suppress the explicit dependence on $H_{1}$ and $H_{2}$.

As an example, assume that $H_{2}$ has only one additional parameter, $\lambda$, for which the likelihood is Gaussian with a median $\mu$ and standard deviation $\sigma$,

[TABLE]

and that the prior is uniform over the width $\Delta\lambda$, $p(\lambda)=1/\Delta\lambda$. The BF is then

[TABLE]

The last term, which is greater than $1$, is the usual Occam’s penalty for the fact that $H_{2}$ has additional parameters and is penalized as compared to $H_{1}$.

Consider now two independent measurements, for example two GW events with $D=\{D_{1},D_{2}\}$, and that the independent events have a parameters $\lambda_{1}$ and $\lambda_{2}$, respectively. The posterior for the deviation parameters after both measurements is

[TABLE]

where to obtain the second line we have factorized the likelihood by assuming that each measurement is independent, and to obtain the third line we have used the identity $p(\lambda_{1},\lambda_{2})=p(\lambda_{1})p(\lambda_{2}|\lambda_{1})$.

Crucially, before we can proceed, we need to specify $p(\lambda_{2}|\lambda_{1})$, which quantifies the relation between $\lambda_{1}$ and $\lambda_{2}$. Notice that $\lambda_{1}$ and $\lambda_{2}$ are parameters describing different systems, but measurement of one can influence our prior belief on the value of the other. This is similar to, for example, the fact that measurement of the mass of one binary affects our prior expectation for the mass of a future binary as the masses are linked through a common astrophysical mass distribution Abbott et al. (2018c). We consider two extreme cases which are commonly used in the literature: one where each GW signal shares a common parameter $\lambda$, and one where the two additional parameters $\lambda_{i}$ are unrelated. We also discuss the most general case, where the properties of each individual signal are related by some common, underlying distribution.

II.2 Case 1: Common parameter and multiplying likelihoods

In the case where we have identified a common parameter for all GWs then $\lambda_{1}=\lambda_{2}\equiv\lambda$ and

[TABLE]

Now in order to make sense of the posterior, we marginalize over $\lambda_{2}$, then simplify the result as

[TABLE]

This equation shows that this case is equivalent to multiplying the individual likelihoods for $\lambda$ obtained by each event, then multiplying by one instance of the common prior and normalizing the resulting probability distribution. In other words, this is the standard situation of measuring a common parameter using multiple, independent trials, where the posterior distribution from the previous measurement can be used as a prior for the next measurement.

To clearly illustrate the result we specialize to uniform priors and Gaussian likelihoods, so that

[TABLE]

The product of the two Gaussian likelihoods is proportional to another Gaussian with

[TABLE]

With this the BF is

[TABLE]

but now $\sigma_{f}$ is smaller than the individual standard deviations, and $\mu_{f}$ is a weighted average of the maximum likelihood value for $\lambda$ coming from each measurement. The Occam’s penalty factor is modified from the case of a single measurement only by the change in $\sigma$.

From here is it straightforward to generalize to $N$ events. For each of the $N$ events, we measure the deviation parameter $\lambda_{i}$. Since the events are independent, the overall likelihood again factorizes, and only the prior $p(\lambda_{1},\dots,\lambda_{N})$ links the inferences. Repeating the proceeding steps for multiple common parameters, where all $\lambda_{i}=\lambda$, we see that the posterior is proportional to the product of the individual likelihoods and a single instance of the prior,

[TABLE]

To show how this influences the BF, we again assume the likelihoods are proportional to Gaussians, and note that for the multiplication of $N$ Gaussians, the result is again proportional to a Gaussian with parameters

[TABLE]

Now the BF is given by Eq. (11) with the above values of $\sigma$ and $\mu$ (the factors of $C$ cancel in the numerator and denominator). To get a clearer picture of how this scales with $N$, we assume all events have the same $\mu_{i}$ and $\sigma_{i}$, in which case

[TABLE]

In the case where $H_{1}$ is correct and so $\mu_{i}=0$, our BF grows with $\sqrt{N}$. In reality even in the case where GR is correct and $\mu_{i}=0$, the specific noise realization of the detector noise will cause variations in $\mu_{i}$ with a standard deviation of $\sigma_{i}$. Below we numerically show that even if we include this variation in our results the $\sqrt{N}$ scaling is preserved on average.

II.3 Case 2: Unrelated parameters and multiplying BFs

In the opposite extreme, we assume that all systems have their own additional parameter, so that we learn nothing about $\lambda_{2}$ from measuring $\lambda_{1}$. Returning briefly to the case where we have only two events, in Eq. (II.1) we set

[TABLE]

and then Eq. (II.1) becomes

[TABLE]

From this posterior and using Eq. (3) we see that the total BF in this case reduces to the product of the individual BFs. Generalizing to $N$ measurements and under the assumption that all $\lambda_{i}$ are unrelated, the total BF is given by the product of the individual BFs, a simplification sometimes used in the literature Li et al. (2012a, b); Agathos et al. (2014); Meidam et al. (2014).

If we specialize once again to Gaussian likelihoods, we find

[TABLE]

where $C$ is the constant prefactor of the multiplied Gaussians, see Eq. (14). We see when comparing to Eq. (11) that $H_{2}$ is penalized by two Occam’s factors, one for each $\lambda_{i}$, rather than the single factor for the extra parameter $\lambda$. This makes sense, since now $H_{2}$ includes two new independent parameters.

Meanwhile, for $N$ measurements where the $\lambda_{i}$ are all unrelated to each other we have that the product of BFs is

[TABLE]

We see that we have $N$ Occam’s factors which penalize $H_{2}$. Again to get a concrete picture we let all the likelihoods be equal, in which case $C=1$, and all prior are ranges equal, so that

[TABLE]

In the case where $H_{1}$ is the true theory, $\mu_{i}=0$ and this BF grows exponentially with $N$.

We conclude that this method leads in general to larger BFs in favor of GR and tighter constraints compared to multiplying the likelihoods. However, this method assumes that the theory of gravity has $N$ coupling constants with respect to GR, one for every signal detected.

II.4 General case: Deviation parameters drawn from a common distribution

So far we have described two extreme cases for computing joint inference for nested hypotheses. In each case, certain assumptions were used to simplify the conditional probability on the deviation parameters $p(\lambda_{1},\dots,\lambda_{N})$. The most general case would allow for a nontrivial distribution for $p(\lambda_{1},\dots,\lambda_{N})$, though not all distributions are physically realistic for GW signals. The various datasets are still assumed to be independent of each other, but now the parameters $\lambda_{i}$ are drawn from a common distribution, similar for example to the expectation that BBH masses are drawn from a common mass distribution. We can then assume that the prior distributions for $\lambda_{1}$ and $\lambda_{2}$ can be related by a common set of additional hyperparameters Loredo (2004).

As a concrete example for an underlying distribution, we return to our Gaussian theme and consider $\lambda_{i}$ drawn from a Gaussian distribution centered on $M$, with standard deviation $\alpha$,

[TABLE]

With this, the posteriors of the individual events will be modified Gaussians with widths $\tilde{\sigma}_{i}$ and means $\tilde{\mu}_{i}$,

[TABLE]

where

[TABLE]

and the overall normalization factor is

[TABLE]

With these definitions, we find that the evidences of the individual posteriors are

[TABLE]

The combined BF can be shown to be similar to Eq. (18), with the numerator unchanged and the denominator modified as above. The result is

[TABLE]

Note that as we hold $\alpha$ fixed and vary $\sigma_{i}$, we always have $\alpha\geq\tilde{\sigma}_{i}$.

The above analysis gives the BF in the case where $H_{2}$ is associated with a particular choice of $M$ and $\alpha$. If we are interested in the BF for a hypothesis where the $\lambda_{i}$ are drawn from some Gaussian with unknown $M$ and $\alpha$, then we should marginalize over these hyperparameters, using some prior $p(M,\alpha)$. This leads to

[TABLE]

For example, our first case where all the $\lambda_{i}$ are equal corresponds to drawing the $\lambda_{i}$ from a Gaussian with zero width $\alpha=0$, centered on some unknown $\lambda=M$. To accomplish this we set $p(M,\alpha)=p(M)\delta(\alpha)$. Assuming a flat hyperprior on $M$, $p(M)=1/(\Delta M)$, we integrate in order to recover the BF of Eq. (11), with the substitution $\Delta\lambda\to\Delta M$.

The second case, where the $\lambda_{i}$ are treated as unrelated parameters, is given by taking $\alpha$ large and enforcing a cutoff on the prior range. In this case, we fix $M$ to some value within the prior range and let $\alpha\to\infty$, which we accomplish by setting $p(M,\alpha)$ to the appropriate delta functions and integrating over them. Then the BF reduces to Eq. (18).

The framework above, where each $\lambda_{i}$ are drawn from some common distribution controlled by a set of hyperparameters, is the most general framework for systems which do not interact. We have shown how each of the two common methods for combining multiple GW events arise out of simple limits from this framework.

II.5 Combining synthetic observations from multiple events

The expressions for the BF presented in Eqs. (15) and (19) were obtained under specific assumptions about the parameters $\lambda_{i}$ that describe the deviation from GR. In order to study the implications of indiscriminately applying these formulas to situations and datasets that violate these assumptions, we return to our toy model of Gaussian likelihoods and numerically compute the BFs by either multiplying the likelihoods, Eq. (15), or by multiplying the BFs of the individual measurements together, Eq. (19). To simulate a mock population of measurements for $\lambda_{i}$ we draw the mean $\mu_{i}$ of the Gaussian likelihood from various distributions and for simplicity set the standard deviation to $\sigma_{i}=1$. We also assume a flat prior on $\lambda_{i}$.

We focus on $5$ example cases for the distribution of $\mu_{i}$:

$\mu_{i}=0+{\cal{N}}(0,\sigma_{i}=1)$: This corresponds to the case that GR is correct, but the posteriors have a scatter around the true value of $\mu_{i}=0$ due to the detector noise realization. Here ${\cal N}(\mu,\sigma)$ indicates a normal distribution of mean $\mu$ and standard deviation $\sigma$. 2. 2.

$\mu_{i}=1+{\cal{N}}(0,\sigma_{i}=1)$: In this case GR is not correct, but the true theory of gravity predicts the same beyond-GR parameter for each event, $\mu_{i}=1$, plus a scatter due to detector noise. This corresponds to the situation discussed in Sec. II.2. 3. 3.

$\mu_{i}\in[-1,1]+{\cal{N}}(0,\sigma_{i}=1)$: Here GR is not correct and the true theory of gravity predicts beyond-GR parameters that are unrelated to each other, plus a scatter due to detector noise. This theory breaks the assumptions of Sec. II.2. 4. 4.

$\mu_{i}\in[-4,4]+{\cal{N}}(0,\sigma_{i}=1)$: This is the same as 3 but here we draw $\mu_{i}$ from a broader distribution so that we clearly correspond to the situation discussed in Sec. II.3. 5. 5.

$\mu_{i}=0.1+{\cal{N}}(0,\sigma_{i}=1)$: This is the same as 2, but with a smaller deviation from GR as compared to the assumed level of the noise in the detectors. This case breaks the assumptions of Sec. II.3.

We then compute the resulting BFs for each scenario of Secs. II.2 and II.2 numerically. The results are in Fig. 1 whose top panel shows the BF as a function of the number of events computed according to Eq. (15), while the bottom panel corresponds to the case of Eq. (19). On both panels the black lines show the expected BF scaling for GR and non-GR, in the case of zero-noise realizations. The green dots refer to the case where GR is the correct theory of gravity but also accounts for the specific noise realization in the detectors. In both cases we find that even with noise the combined BFs favor GR as expected, although the exact scaling is affected compared to a set of zero-noise realizations.

In both panels the purple dots correspond to the BFs computed given a non-GR theory of gravity that obeys the assumptions under which each individual BF expression was derived. Specifically, in the top panel we show the second injection case above, in which all GW systems share a common nonzero parameter, while in the bottom panel we show the fourth injection case which corresponds to the scenario where each system has a randomly chosen parameter in the range $[-4,4]$. Both these theories lead to sharply decreasing BFs correctly signaling a deviation from GR.

Finally, the pink dots in both plots correspond to non-GR theories of gravity that break the assumptions under which the BF expression was derived. In the top panel, we show the results from the third case, while in the bottom panel we show the fifth case. Both lead to BFs in favor of GR, and as a consequence we would conclude that GR is correct and fail to detect the deviation. Specifically, if the correct theory of gravity predicted parameters distributed symmetrically around [math] and we analyze these systems assuming they have the same parameter (which is equivalent to multiplying the individual likelihoods), we would erroneously build confidence in favor of the wrong theory of gravity, as shown on the top panel. We encounter a similar situation in the right panel, where the pink dots are in fact more in favor of GR than the case where GR is correct and the data is noisy.

The above considerations lead us to a seemingly obvious conclusion: if we perform an analysis under certain assumptions which in reality are violated, we are in peril of obtaining biased answers. However in the case of testing GR this statement has additional implications. In Sec. II we argued that certain assumption needs to be made when combining many events, in the form of selecting an expression for the term $p(\lambda_{2}|\lambda_{1})$. As a consequence, it is not possible to test all possible theories of gravity with a common analysis, and so there can be no truly model-independent test of GR using multiple GW detections.

III Full GW analysis: Measuring the graviton mass

As a concrete example of how a joint analysis, calling on information from several independent observations, can be either powerful or misleading we present the case of a theory where the graviton has an unknown but finite mass. The emitted GW will then be dispersed as it propagates through the universe. This will manifest as a modification of the GW phase relative to the GR prediction where the graviton is assumed to be massless.

In the case of a nonzero graviton mass, the waveform will be modified starting at the 1PN order111In the post-Newtonian framework, an $N$PN order term is a factor of $(u/c)^{2N}$ over the leading order term, where $u$ is some characteristic velocity of the system and $c$ is the speed of light. Will (1998) and the GW phase will get an additional leading-order term of

[TABLE]

where ${\cal{M}}$ is the detector frame chirp mass of the system, $\lambda_{G}$ is the graviton Compton wavelength (with GR assuming $\lambda_{G}=\infty$), $z$ is the redshift, and $D$ is

[TABLE]

with $D_{L}$ being the luminosity distance. The frequency dependence of this term is encoded in $u=\pi{\cal{M}}f$.

In the ppE framework Yunes and Pretorius (2009), the 1PN ppE parameter $\beta_{1}$ is nonzero and equal to $\Psi_{G}$, while all other ppE terms $\beta_{i}$ vanish. Using the redefinition of Agathos et al. (2014) where the ppE deviations are expressed as relative phase shifts $d\chi_{i}$, compared to the GR PN term, the only nonzero modification is

[TABLE]

Restricting to this specific beyond-GR model allows us to study combining information under different assumptions, as the theory is both well understood and easily quantified through a universal additional parameter $\lambda_{G}$.

III.1 Simulated BBH systems

We simulate $18$ binary black hole (BBH) systems with masses drawn from a distribution between $4-45M_{\odot}$, isotropic spin directions, and dimensionless spin magnitudes distributed uniformly between $0-0.8$. The emitted gravitational waves are described by the IMRPhenomPv2 waveform family Hannam et al. (2014); Husa et al. (2016); Khan et al. (2016), with all systems assuming a common graviton Compton wavelength of $\lambda_{G}^{\mathrm{true}}=10^{12}$ km. Notice that while this value of $\lambda_{G}$ has been ruled out by recent GW observations Abbott et al. (2017), here we are primarily interested in an illustration of how a deviation of GR could be extracted from multiple signals, and not to directly show how well the graviton mass can be constrained.

The BBHs are distributed uniformly across the sky with distances drawn from a distribution uniform in co-moving volume, assuming the cosmology defined in Ade et al. (2016). We inject those systems in a noise-free network of 3 advanced LIGO detectors Aasi et al. (2015) operating at design sensitivity aLI , including the currently-operational LIGO-Hanford and LIGO-Livingston detectors, and the under-construction detector LIGO-India. All of the 18 systems considered here have a network signal-to-noise ratio above 12.

We perform parameter estimation on those simulated signals with the publicly-available software library LALInference Veitch et al. (2015); LIGO Scientific Collaboration under the same prior assumptions as the BBH analyses presented in Abbott et al. (2018a), with the addition of a ppE-type deviation from GR at the 1PN order parametrized as in Eq. (30), with a prior defined as uniform across $-30\geq d\chi_{2}\geq 30$.

III.2 General ppE test

Following Eqs. (28) and (30), the specific true value of the parameter $d\chi_{2}$ differs from source to source as their respective masses and distances vary. This means that a joint analysis assuming a common value for $d\chi_{2}$ for all the BBHs will result in incorrect conclusions when applied to the massive graviton model. We present such an analysis here to illustrate what can occur when ppE parameters are assumed to be common among sets of GW detections.222Note that published constraints on the graviton mass using multiple GW events do correctly combine inferences on the common parameter $\lambda_{G}$ Abbott et al. (2017), contrary to our illustration in this section.

In the case where the individual posteriors for $d\chi_{2}$ from each event are broad and overlap each other, assuming a common $d\chi_{2}$ and combining likelihoods as in Eq. (II.2) is expected to lead to uninformative joint posteriors. We would expect to recover the prior on $d\chi_{2}$, and combining BFs as in Eq. (15) would therefore be uninformative. Here we face another situation, as illustrated in Fig. 2, where we show the posterior distributions for $d\chi_{2}$ for all our $18$ events. Because our chosen $\lambda_{G}$ provides a relatively strong deviation in the waveforms, each single-event analysis can lead to the conclusion that a deviation from GR is present. None of the posteriors is consistent with the GR prediction of $d\chi_{2}=0$. However, many of the $d\chi_{2}$ posteriors shown in Fig. 2 are disjoint, so the joint posterior distribution would have no support over its prior range. In this case it is not possible to make a statement about the nature of the observed deviation, or combine the observations in order to build further evidence of the deviation, without modeling the observed $d\chi_{2}$ distributions as arising from an underlying population.

Hence, treating the general $d\chi_{2}$ parameter as universal, and as described in Sec. II.5, combining their likelihoods accordingly, will not provide any relevant information about the validity of GR from the ensemble of observations.

III.3 Graviton mass analysis

If we instead make the initial, and very strong, assumption that the observed deviations shown in Fig. 2 are generated by the presence of a massive graviton, we can reinterpret them as describing the fundamental constant associated with this specific extension of GR, in this case the graviton Compton wavelength $\lambda_{G}$. For each sample in the posteriors, we compute $\lambda_{G}$ from $d\chi_{2}$, $D_{L}$, our assumed cosmology, and the masses using Eqs. (28)–(30) in order to arrive at a posterior for $\lambda_{G}$ for each event. Now the assumptions about creating a joint posterior distribution described in Sec. II.2 are applicable, since the likelihood distributions from the individual analyses can be multiplied together to constructively build up the available information, as is shown in Fig. 3.

The joint measurement of $\lambda_{G}$ from these events is $1.0115^{+0.0086}_{-0.0118}\times 10^{12}$ km given as the maximum posterior value and its associated $90\%$ credible interval, which encompasses the true value of $\lambda_{G}^{\mathrm{true}}=10^{12}$ km. Fig. 3 shows a slight bias towards larger values of $\lambda_{G}$. This is caused by the fact that the distance posteriors are relatively uninformative for these signals, and largely follow the assumed prior distribution where $p(D_{L})\propto D_{L}^{2}$. Our inferences on $\lambda_{G}$ are then affected due to the fact that $\lambda_{G}\propto\sqrt{D_{L}/d\chi_{2}}$. For a larger set of observations, with a wider spread in both their source parameters and their SNR distribution, this bias is expected to diminish.

IV Conclusions

In this work, we revisit two common methods for combining information from multiple sources employed in the testing-GR literature and argue that they are the limiting cases of a more generic hierarchical framework. Moreover, we show that both methods make certain assumptions about the underlying theory of gravity being tested. In particular, multiplying the likelihood functions of beyond-GR parameters from individual binary merger events amounts to imposing that all signals share the same beyond-GR parameter. Meanwhile, computing the BF in favor of GR from each event and multiplying the BFs together makes the implicit assumption that each binary has an additional beyond-GR parameter, and these parameters are unrelated to each other.

We argue that both assumptions are expected to fail for generic theories of gravity. We present numerical results of BFs where we show that if the true theory of gravity does not obey the assumptions under which information was combined and the BF was calculated, we might fail to detect a deviation from GR.

Our results suggest that combining inferences from multiple sources cannot be performed in a completely model-independent way such that it applies to generic theories of gravity. Even in the generic hierarchical analysis, a model for the population distribution needs to be specified in some fashion. We present here the simple case where the beyond-GR parameters are distributed according to a normal distribution, though with enough detections multiple models for the common distribution can be explored and compared, while their hyperparameters are measured.

We have also applied these ideas to a more realistic study of signals generated assuming a theory where the graviton is massive, which we have analyzed using a full parameter estimation analysis. We recover the signal parameters in two ways: either we make inferences on the Compton wavelength of the graviton directly, in which case the deviation parameter is the same among all the signals, or we use a generic ppE parameterization where the beyond-GR parameter differs from signal to signal. Although for illustration we select a graviton mass where the deviations from GR are detectable in the individual signals, we argue that when using incorrect assumptions on how to combine these signals we would fail to detect a deviation from GR in the joint analysis. Meanwhile, using the correct method for combining the posteriors allows us to draw tight constraints on the graviton mass using the full dataset. This exercise highlights the needs to re-interpret the ppE parameters in terms of specific theories of gravity before combining measurements of them.

We emphasize that even though we frame our results in the context of testing for the true theory of gravity, our conclusions are applicable to any problem that involves combining information, for example in making cosmological measurements and in inferring populations of compact binaries, where these ideas have long been in use.

Acknowledgements.

We thank Will Farr, Max Isi, Walter Del Pozzo, and Salvatore Vitale for helpful discussions. We are grateful for computational resources provided by Cardiff University, and funded by an STFC grant supporting UK Involvement in the Operation of Advanced LIGO. C.-J. H. acknowledge support of the MIT physics department through the Solomon Buchsbaum Research Fund, the National Science Foundation, and the LIGO Laboratory. The Flatiron Institute is supported by the Simons Foundation. This is LIGO Document Number P1900098.