On combining information from multiple gravitational wave sources
Aaron Zimmerman Carl-Johan Haster Katerina Chatziioannou














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.


On combining information from multiple gravitational wave sources
    Katerina Chatziioannou
    December 30, 2023
=================================================================



§ 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 <cit.>. 
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 <cit.>. 
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) <cit.>, the measurement of the Hubble constant <cit.>, as well as inference about population properties <cit.>.

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 <cit.>. 
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 <cit.>. 
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 <cit.>.

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 <ref> studies the assumptions inherently made when combining inferences and illustrates these ideas using a simple model of Gaussian likelihoods. 
Section <ref> presents a concrete example in terms of measuring the Compton wavelength of the graviton. 
Section <ref> summarizes our main conclusions.








§ 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.





 §.§ 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


    = p(D | H_1)/p(D|H_2) ,


where p(D|H_i)=∫ p(θ_i) p(D|θ_i, H_i) d θ_i is the evidence for H_i, θ_i are the model parameters, p(θ_i) is the prior on the parameters, and p(D|θ_i, H_i) is the likelihood.

We assume that H_1 and H_2 are nested models, i.e. θ_2= {θ_1,λ} and p(D|θ_1, H_1)=p(D|θ_1,λ=0,H_2). 
The posterior for the extra parameters λ is


    p(λ|D,H_2)   =p(λ| H_2)p(D|λ, H_2)/p(D|H_2)
    
    =   p(D| λ,H_2)/∫ p(λ',H_2) p(D| λ',H_2) d λ' .


and using the Savage-Dickey density ratio <cit.> the BF simplifies to 


    =p(λ=0|D,H_2)/p(λ=0| H_2) .


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, λ, for which the likelihood is Gaussian with a median μ and standard deviation σ,


    p(D|λ)     = c e^-(λ - μ)^2/(2σ^2) ,


and that the prior is uniform over the width Δλ, p(λ) = 1/Δλ.
The BF is then

    = 1/√(2 π)e^-μ^2/(2 σ^2)Δλ/σ .

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 λ_1 and λ_2, respectively. 
The posterior for the deviation parameters after both measurements is


    p(λ_1,λ_2|D)     = p(λ_1,λ_2) p(D| λ_1,λ_2)/∫ p(λ_1',λ_2') p(D|λ_1',λ_2') d λ_1' d λ_2'
       = p(λ_1,λ_2) p(D_1| λ_1)p(D_2| λ_2)/∫ p(λ_1',λ_2') p(D_1|λ_1')p(D_2|λ_2') d λ_1' d λ_2'
       = p(λ_1,λ_2) p(D_1| λ_1 )p(D_2| λ_2)/∫ p(λ_1')p(λ_2'|λ_1') p(D_1|λ_1')p(D_2|λ_2') d λ_1' d λ_2' ,


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(λ_1,λ_2)=p(λ_1)p(λ_2|λ_1).

Crucially, before we can proceed, we need to specify p(λ_2|λ_1), which quantifies the relation between λ_1 and λ_2. 
Notice that λ_1 and λ_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 <cit.>.
We consider two extreme cases which are commonly used in the literature: one where each GW signal shares a common parameter λ, and one where the two additional parameters λ_i are unrelated.
We also discuss the most general case, where the properties of each individual signal are related by some common, underlying distribution.





 §.§ Case 1: Common parameter and multiplying likelihoods



In the case where we have identified a common parameter for all GWs then λ_1=λ_2≡λ and


    p(λ_2|λ_1) = δ(λ_1-λ_2).


Now in order to make sense of the posterior, we marginalize over λ_2, then simplify the result as


    p(λ|D)   =p(λ_1|D)=∫ p(λ_1,λ_2|D)d λ_2 
       =∫p(λ_1,λ_2) p(D_1| λ_1)p(D_2| λ_2)/∫ p(λ_1')p(λ_2'|λ_1') p(D_1|λ_1')p(D_2|λ_2') d λ_1' d λ_2'd λ_2
        = p(λ)p(D_1| λ)p(D_2| λ)/∫ p(λ') p(D_1|λ')p(D_2|λ') d λ'  .


This equation shows that this case is equivalent to multiplying the individual likelihoods for λ 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


    p(D_i|λ)     = c_i e^-(λ - μ_i)^2/(2σ_i^2) .


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


    σ_f^2     = σ_1^2 σ_2^2/σ_1^2 + σ_2^2 ,    μ_f     = μ_1 σ_2^2 + μ_2 σ_1^2/σ_1^2 + σ_2^2 .


With this the BF is


    = e^-μ_f^2/(2 σ_f^2)Δλ/√(2 π)σ_f ,


but now σ_f is smaller than the individual standard deviations, and μ_f is a weighted average of the maximum likelihood value for λ coming from each measurement.
The Occam's penalty factor is modified from the case of a single measurement only by the change in σ. 

From here is it straightforward to generalize to N events. 
For each of the N events, we measure the deviation parameter λ_i.
Since the events are independent, the overall likelihood again factorizes, and only the prior p(λ_1, …, λ_N) links the inferences.
Repeating the proceeding steps for multiple common parameters, where all λ_i = λ, we see that the posterior is proportional to the product of the individual likelihoods and a single instance of the prior,

    p(λ|D) ∝ p(λ) ∏_i p(D_i,λ)  .


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

    σ^2     = (∑_i^N 1/σ_i^2)^-1 ,    μ    = σ^2 ∑_i^N μ_i/σ_i^2 ,
    
    ln C     = -1/2∑_i^N μ_i^2/σ_i^2 + μ^2/2σ^2 .

Now the BF is given by Eq. (<ref>) with the above values of σ and μ (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 μ_i and σ_i, in which case


    = √(N) e^-N μ_i^2/(2 σ_i^2)Δλ/√(2 π)σ_i .


In the case where H_1 is correct and so μ_i =0, our BF grows with √(N). 
In reality even in the case where GR is correct and μ_i =0, the specific noise realization of the detector noise will cause variations in μ_i with a standard deviation of σ_i. 
Below we numerically show that even if we include this variation in our results the √(N) scaling is preserved on average. 




 §.§ 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 λ_2 from measuring λ_1.
Returning briefly to the case where we have only two events, in Eq. (<ref>) we set


    p(λ_2|λ_1) = p(λ_2),


and then Eq. (<ref>) becomes


    p(λ_1,λ_2|D)     = p(λ_1,λ_2)p(D_1| λ_1)p(D_2| λ_2)/∫ p(λ_1')p(λ_2'|λ_1') p(D_1|λ_1')p(D_2|λ_2') d λ_1' d λ_2'
        = p(λ_1) p(λ_2) p(D_1| λ_1)p(D_2| λ_2)/∫ p(λ_1') p(λ_2') p(D_1|λ_1')p(D_2|λ_2') d λ_1' d λ_2'
        = p(λ_1) p(D_1| λ_1)/∫ p(λ_1') p(D_1|λ_1') d λ_1' p(λ_2)p(D_2| λ_2)/∫ p(λ_2') p(D_2|λ_2') d λ_2' .


From this posterior and using Eq. (<ref>) 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 λ_i are unrelated, the total BF is given by the product of the individual BFs, a simplification sometimes used in the literature <cit.>.

If we specialize once again to Gaussian likelihoods, we find 


    = C e^-μ_f^2/(2 σ_f^2)Δλ_1/√(2 π)σ_1Δλ_2/√(2 π)σ_2 ,


where C is the constant prefactor of the multiplied Gaussians, see Eq. (<ref>).
We see when comparing to Eq. (<ref>) that H_2 is penalized by two Occam's factors, one for each λ_i, rather than the single factor for the extra parameter λ. 
This makes sense, since now H_2 includes two new independent parameters.

Meanwhile, for N measurements where the λ_i are all unrelated to each other we have that the product of BFs is


    = C e^- μ^2/(2 σ^2)∏_i Δλ_i/√(2 πσ_i) .


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

    =  e^-N μ_i^2/(2 σ_i^2)( Δλ/√(2π)σ_i)^N  .
 
In the case where H_1 is the true theory, μ_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. 




 §.§ 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(λ_1, …, λ_N).
The most general case would allow for a nontrivial distribution for p(λ_1, …, λ_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 λ_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 λ_1 and λ_2 can be related by a common set of additional hyperparameters <cit.>. 

As a concrete example for an underlying distribution, we return to our Gaussian theme and consider λ_i drawn from a Gaussian distribution centered on M, with standard deviation α,


    p(λ_i|M,α)     = 1/√(2π)αe^-(λ_i - M)^2/(2α^2) .


With this, the posteriors of the individual events will be modified Gaussians with widths σ̃_i and means μ̃_i,


    p(λ_i) p(D_i|λ_i)     = 1/2 πσ_i α_iC̃ e^-(λ_i - μ̃_i)^2/(2 σ̃_i^2) ,


where


    σ̃_i^2    = σ_i^2 α^2/σ_i^2 + α^2 ,    μ̃_i     = μ_i α^2 + M σ_i^2/σ_i^2 + α^2 ,


and the overall normalization factor is


    lnC̃_i     = 1/2(μ̃_i^2/σ̃_i^2 - μ_i^2/σ_i^2 - M^2/α^2)   = - (μ_i - M)^2/2(σ_i^2 + α^2) .


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


    ∫ p(λ_i) p(D_i|λ_i) d λ_i = 1/σ_iσ̃_i C̃_i/√(2 π)α .


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


    = C/C̃_1 C̃_2 e^-μ_f^2/(2 σ_f^2)α/σ̃_1α/σ̃_2 .


Note that as we hold α fixed and vary σ_i, we always have α≥σ̃_i.

The above analysis gives the BF in the case where H_2 is associated with a particular choice of M and α.
If we are interested in the BF for a hypothesis where the λ_i are drawn from some Gaussian with unknown M and α, then we should marginalize over these hyperparameters, using some prior p(M, α).
This leads to


    = Ce^-μ_f^2/(2 σ_f^2)( ∫C̃_1 σ̃_1/αC̃_2 σ̃_2/α 
    p(M, α) dM dα)^-1 .



For example, our first case where all the λ_i are equal corresponds to drawing the λ_i from a Gaussian with zero width α=0, centered on some unknown λ = M. 
To accomplish this we set p(M,α) = p(M) δ(α).
Assuming a flat hyperprior on M, p(M) = 1/(Δ M), we integrate in order to recover the BF of Eq. (<ref>), with the substitution Δλ→Δ M.

The second case, where the λ_i are treated as unrelated parameters, is given by taking α large and enforcing a cutoff on the prior range.
In this case, we fix M to some value within the prior range and let α→∞, which we accomplish by setting p(M, α) to the appropriate delta functions and integrating over them.
Then the BF reduces to Eq. (<ref>).

The framework above, where each λ_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.

 


 §.§ Combining synthetic observations from multiple events


 
The expressions for the BF presented in Eqs. (<ref>) and (<ref>) were obtained under specific assumptions about the parameters λ_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. (<ref>), or by multiplying the BFs of the individual measurements together, Eq. (<ref>).
To simulate a mock population of measurements for λ_i we draw the mean μ_i of the Gaussian likelihood from various distributions and for simplicity set the standard deviation to σ_i=1. 
We also assume a flat prior on λ_i.
 
We focus on 5 example cases for the distribution of μ_i: 



  * μ_i=0+ N(0,σ_i=1): This corresponds to the case that GR is correct, but the posteriors have 
a scatter around the true value of μ_i=0 due to the detector noise realization. Here N(μ, σ) indicates a normal distribution of mean μ and standard deviation σ.

  * μ_i=1+ N(0,σ_i=1): In this case GR is not correct, but the true theory of gravity predicts
the same beyond-GR parameter for each event, μ_i=1, plus a scatter due to detector noise. This corresponds
to the situation discussed in Sec. <ref>.

  * μ_i∈ [-1,1]+ N(0,σ_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. <ref>.

  * μ_i ∈ [-4,4] + N(0,σ_i=1): This is the same as 3 but here we draw μ_i from a broader distribution so that we clearly correspond to the situation discussed in Sec. <ref>.

  * μ_i =0.1 + N(0,σ_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. <ref>.



 



We then compute the resulting BFs for each scenario of Secs. <ref> and <ref> numerically. 
The results are in Fig. <ref> whose top panel shows the BF as a function of the number of events computed according to Eq. (<ref>), while the bottom panel corresponds to the case of Eq. (<ref>). 
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 0 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. <ref> we argued that certain assumption needs to be made when combining many events, in the form of selecting an expression for the term p(λ_2|λ_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.

 
 


§ 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 order[In the post-Newtonian framework, an NPN 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.] <cit.> and the GW phase will get an additional leading-order term of


    Ψ_g=-π^2 D M/λ_G^2 (1+z)u^-1,


where M is the detector frame chirp mass of the system, λ_G is the graviton Compton wavelength (with GR assuming λ_G=∞), z is the redshift, and D is


    D=D_L 1+(2+z)(1+z+√(1+z))/5(1+z)^2,


with D_L being the luminosity distance. 
The frequency dependence of this term is encoded in u=πMf.

In the ppE framework <cit.>, the 1PN ppE parameter β_1 is nonzero and equal to Ψ_G, while all other ppE terms β_i vanish. 
Using the redefinition of <cit.> where the ppE deviations are expressed as relative phase shifts dχ_i, compared to the GR PN term, the only nonzero modification is


    dχ_2=Ψ_G/5/96( 743/336+11/4η)η^-2/5u^-1.



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 λ_G.


  
 

 §.§ Simulated BBH systems

  
We simulate 18 binary black hole (BBH) systems with masses drawn from a distribution between 4-45 M_⊙, isotropic spin directions, and dimensionless spin magnitudes distributed uniformly between 0-0.8.
The emitted gravitational waves are described by the IMRPhenomPv2 waveform family <cit.>, with all systems assuming a common graviton Compton wavelength of λ_G^true=10^12 km. 
Notice that while this value of λ_G has been ruled out by recent GW observations <cit.>, 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 <cit.>.
We inject those systems in a noise-free network of 3 advanced LIGO detectors <cit.> operating at design sensitivity <cit.>, 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 <cit.> under the same prior assumptions as the BBH analyses presented in <cit.>,  with the addition of a ppE-type deviation from GR at the 1PN order parametrized as in Eq. (<ref>), with a prior defined as uniform across -30 ≥ dχ_2 ≥ 30.
 

 

 §.§ General ppE test


Following Eqs. (<ref>) and (<ref>), the specific true value of the parameter dχ_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χ_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.[Note that published constraints on the graviton mass using multiple GW events do correctly combine inferences on the common parameter λ_G <cit.>, contrary to our illustration in this section.]

In the case where the individual posteriors for dχ_2 from each event are broad and overlap each other, assuming a common dχ_2 and combining likelihoods as in Eq. (<ref>) is expected to lead to uninformative joint posteriors.
We would expect to recover the prior on dχ_2, and combining BFs as in Eq. (<ref>) would therefore be uninformative.
Here we face another situation, as illustrated in Fig. <ref>, where we show the posterior distributions for dχ_2 for all our 18 events. 
Because our chosen λ_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χ_2 = 0.
However, many of the dχ_2 posteriors shown in Fig. <ref> 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χ_2 distributions as arising from an underlying population.


 
Hence, treating the general dχ_2 parameter as universal, and as described in Sec. <ref>, combining their likelihoods accordingly, will not provide any relevant information about the validity of GR from the ensemble of observations.
 
  
 

 §.§ Graviton mass analysis


If we instead make the initial, and very strong, assumption that the observed deviations shown in Fig. <ref> 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 λ_G.
For each sample in the posteriors, we compute λ_G from dχ_2, D_L, our assumed cosmology, and the masses using Eqs. (<ref>)–(<ref>) in order to arrive at a posterior for λ_G for each event.
Now the assumptions about creating a joint posterior distribution described in Sec. <ref> 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. <ref>.

The joint measurement of λ_G from these events is 1.0115^+0.0086_-0.0118× 10^12 km given as the maximum posterior value and its associated 90% credible interval, which encompasses the true value of λ_G^true=10^12 km.
Fig. <ref> shows a slight bias towards larger values of λ_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) ∝ D_L^2. Our inferences on λ_G are then affected due to the fact that λ_G ∝√(D_L/dχ_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.


 
 
 

§ 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.
 
 
 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.



