Coincident detection significance in multimessenger astronomy

Given two detections $a$ and $b$ in different data sets $D_{a}$ and $D_{b}$, we would like to assess the hypothesis that they originate from a common source. In general, the two detections will be described by different physical signal models $\mathcal{H}^{\textrm{S}}_{a}$ and $\mathcal{H}^{\textrm{S}}_{b}$, respectively. Each signal model will imply a likelihood, a set of parameters and an associated prior for those parameters. To assess whether they originate from a common source, the models must share a common set of parameters $\theta\in\Theta$.

We’ll use notation where $\mathcal{H}(\theta)\equiv[\mathcal{H}\;\textrm{and}\;\theta]$ denotes a hypothesis $\mathcal{H}$ with a particular choice for the parameters $\theta$, while $\mathcal{H}$ by itself denotes a hypothesis with unknown parameters, i.e. “for any choice of parameters $\theta$”. We can formally write this as $\mathcal{H}\equiv[\mathcal{H}(\theta)\text{ for any }\theta]$.

Then, we define the common-source hypothesis:

We also define $\mathcal{H}^{\textrm{N}}_{a/b}$ as the noise hypotheses (by which we mean any non-signal) for each data set. Then we can define any alternative hypothesis for which the observed detections in $a$ and $b$ are unrelated:

where $X$, $Y$ $\in\{\textrm{N},\textrm{S}\}$. We write this in a general form, but note that the noise hypothesis will not have any common model parameters. In total, there are four possible realizations of $\mathcal{H}^{XY}$, which we consider in detail below. However, $\mathcal{H}^{\textrm{S}\textrm{S}}$ is of particular interest in this work, being two unrelated signals from distinct sources.

These hypotheses imply priors on $\theta$ which in general differ from those implied by $\mathcal{H}^{\textrm{S}}_{a/b}$ individually: if a common source can only be detected in some subset of $\theta$, then $\mathcal{H}^{\textrm{C}}$ can only have prior support restricted to this subset. If this is not true, we identify the special case

The probability of the common-source hypothesis is given by

In this work, we will calculate the odds between $\mathcal{H}^{\textrm{C}}$ and different choices of $\mathcal{H}^{XY}$

where

is the Bayes factor and $P(\mathcal{H}^{\textrm{C}})/P(\mathcal{H}^{XY})$ is the prior odds. In Sec. 2.2 we discuss the calculation of the Bayes factor in general. The prior odds will depend on the context, but in Sec. 3.1.3 we calculate the prior odds modeling $\mathcal{H}^{\textrm{C}}$ and $\mathcal{H}^{\textrm{S}\textrm{S}}$ as realizations of a Poisson point process.

If both data sets contain a signal from the same event, then they are not independent: $P(D_{a},D_{b}|\mathcal{H}^{\textrm{C}})\neq P(D_{a}|\mathcal{H}^{\textrm{C}})P(D_{b}|\mathcal{H}^{\textrm{C}})$. Instead we must compute

where the domain of the integral in the second line is restricted to the prior support of $\mathcal{H}^{\textrm{C}}$, namely

The need for this restriction arises because assuming that $(\theta=\theta^{\prime})$ and $\mathcal{H}^{\textrm{C}}$ are both true would be a contradiction if $P(\theta^{\prime}|\mathcal{H}^{\textrm{C}})=0$, and so $P(D|\theta^{\prime},\mathcal{H}^{\textrm{C}})$ would be undefined. Rearranging the likelihood in the integrand

where in the second step we have used that the likelihoods conditional on $\theta$ can be separated for the two data sets, provided that $\theta$ is the set of all model parameters common between the two likelihoods. In the last step, we again used that $P(\theta|\mathcal{H}^{\textrm{C}})>0$ within the integration interval. A subtle point is that $P(\theta|D_{a/b},\mathcal{H}^{\textrm{C}})$ is the posterior distribution for the common model parameters (given either $D_{a/b}$) marginalized over all other model parameters and using the prior implied by $\mathcal{H}^{\textrm{C}}$.

Substituting Eq. (9) into Eq. (7)

where the posterior overlap integral

quantifies the agreement between the posterior distributions of $\theta$ derived independently. In this integral, the prior has the effect of setting a scale against which the degree of overlap can be compared.

Eq. (10)-(11) demonstrate how probabilities from separate data sets combine when each provides independent inferences about a common model parameter.

Returning to the Bayes factor, by our definition of the alternative hypothesis

So from Eq. (6),

We now specify three particular cases of interest for the alternative hypothesis. First, consider $\mathcal{H}^{\textrm{N}\textrm{N}}$: both $a$ and $b$ are caused by noise. Then Eq. (13) specializes to

where $\mathcal{B}_{\textrm{C}/\textrm{N}}$, in analogy with Eq. (6), is the common-source against noise Bayes factor. In the special case of Eq. (3), it can be shown that $\mathcal{B}_{\textrm{C}/\textrm{N}}(D_{a/b})=\mathcal{B}_{\textrm{S}/\textrm{N}}(D_{a/b})$, i.e. the independent signal against noise Bayes factor for each data stream.

This agrees with our intuition: if both signals are strong compared to the background noise and there is a good overlap of their common model parameters (quantified by the integral), we believe they originate from a common event. This is a powerful result as one can compute the joint Bayes factor from the common-source against noise Bayes factor for each detection individually, and the posterior overlap integral of $\theta$. Eq. (14) has analogous applications to the Fisher combined probability test used in Aartsen et al. (2014).

Second, consider $\mathcal{H}^{\textrm{S}\textrm{N}}$: $a$ was due to a signal, but $b$ was due to noise. For this case, Eq. (13) gives

For us to believe that detection $b$ is a real signal and originates from the same source as $a$, we require that the product of the Bayes factor for common-source against noise in $b$ and the posterior overlap be large. The case $\mathcal{B}_{\textrm{C}/\textrm{N}\textrm{S}}$ is analogous and the same special cases apply as mentioned previously.

Finally, consider $\mathcal{H}^{\textrm{S}\textrm{S}}$, the distinct-source hypothesis: both $a$ and $b$ are of the same nature as in the common-source hypothesis $\mathcal{H}^{\textrm{C}}$, but they are physically distinct (i.e. they belong to unrelated sources with different parameters $\theta_{a}\not=\theta_{b}$). Then,

This equation and the posterior overlap integral of Eq. (11) are the main results of this paper. This provides a simple and intuitive way to assess whether two detections originate from the same event, based on the posterior overlap of their common model parameters.

In the special case of Eq. (3), the prefactor to the posterior overlap integral is unity, such that

On the other hand, when Eq. (3) does not apply, the prefactor plays an important role in quantifying how the restricted prior implied by $\mathcal{H}^{\textrm{C}}$ affects the Bayes factor.

A similar result to Eq. (17) was obtained independently by Haris et al. (2017) in the context of strongly lensed gravitational wave signals from binary black hole mergers.

When calculating $\mathcal{I}_{\theta}$, it is often convenient to factorize the posterior overlap integral, e.g. $\mathcal{I}_{\theta}=\mathcal{I}_{\phi}\mathcal{I}_{\psi}$ where $\phi\subsetneq\theta$ and $\psi=\theta\setminus\phi$. This factorization can only be performed, however, if $P(\phi|\psi,D_{A/B},\mathcal{H}^{\textrm{S}})=P(\phi|D_{A/B},\mathcal{H}^{\textrm{S}})$. There are situations in which this is the case, for example if the joint posterior distribution is an uncorrelated multivariate normal distribution. But generally, this will not be the case and the posterior over the full common parameter space must be used. There are however cases where, under certain assumptions, the integral can be approximately factorized. We will explore one such setting in the next Section.

To calculate the Bayes factor, Eq. (17), we first write down the posterior overlap integral over the conditional joint distribution of the spatial and temporal parameters

We will now show that this can be factorized into a spatial and temporal overlap under the following assumptions. First, that the prior itself factors, $P(\boldsymbol{\Omega},t_{\textrm{c}}|\mathcal{H}^{\textrm{S}})=P(\boldsymbol{\Omega}|\mathcal{H}^{\textrm{S}})P(t_{\textrm{c}}|\mathcal{H}^{\textrm{S}})$. Second, that $t_{\textrm{c}}$ inferred from the GW data is exactly determined, i.e.

where the “hat” indicates the observed value. Then, Eq. (18) can be factorized as $\mathcal{I}_{\boldsymbol{\Omega},t_{\textrm{c}}}=\mathcal{I}_{t_{\textrm{c}}}\mathcal{I}_{\boldsymbol{\Omega}}$, where

and

This factorization is exact under the two assumptions made. However, the coalescence time is typically known with a nonzero uncertainty. For this case, taking $\widehat{t_{\textrm{c}}}$ to be a point estimate (the mean for example), the factorization is approximate, but applicable provided that over the uncertainty in $t_{\textrm{c}}$, $P(t_{\textrm{c}}|D_{\textrm{EM}},\mathcal{H}^{\textrm{S}}_{\textrm{EM}})$, $P(\boldsymbol{\Omega}|t_{\textrm{c}},D_{\textrm{GW}},\mathcal{H}^{\textrm{S}}_{\textrm{GW}})$, and $P(\boldsymbol{\Omega}|t_{\textrm{c}},D_{\textrm{EM}},\mathcal{H}^{\textrm{S}}_{\textrm{EM}})$ do not vary substantially. In Sec. 3.1.1 and Sec. 3.1.2, we will provide approximations for Eq. (20) and (21) under some reasonable assumptions and illustrate some of the subtleties in their calculation.

We note that a similar result to Eq. (21) was previously derived in Urban (2016); in particular, Eq. (3.6) of that work is equivalent to Eq. (21) assuming an isotropic prior. Then the resulting joint likelihood ratio is defined using the alternative hypothesis that went into Eq. (15).

We now apply the example calculated in Sec. 3.1 to GW170817 and GRB 170817A, the result of which can be compared with Abbott et al. (2017a). We note that, the calculation presented here could be improved by using the full joint distribution without making assumptions that allow the result to be factorized, and including other pertinent model parameters such as the luminosity distance (for Fermi-GBM, this may be as simple as estimating the range of conceivable values).

The sky localization for the BNS inspiral and short GRB can be seen in Fig. 1 of Abbott et al. (2017a). Using the published localization FITS files (Singer, 2017; Goldstein et al., 2017) and a uniform prior distribution on the whole sky, Eq. (21) yields $\mathcal{I}_{\boldsymbol{\Omega}}=32.4$. The spatial overlap alone provides moderate support for the common-event model, the main limitation being the uncertainty on the localization of GRB 170817A.

If we did not have the actual FITS files, we could still use the published 90% confidence intervals for the sky localization of GW170817 and GRB 170817A and take the localization posterior distributions to be uniform within those intervals. The intervals cover respectively $28\text{\,}\textrm{deg}\mathrm{{}^{2}}$ (Abbott et al., 2017c) and $1100\text{\,}\textrm{deg}\mathrm{{}^{2}}$ (Goldstein et al., 2017) and the GW170817 interval is entirely contained within GRB 170817A. A straightforward application of Eq. (29) then yields an approximate spatial Bayes factor of $\mathcal{I}_{\boldsymbol{\Omega}}=37.5$, which is close to the exact value. Since we do have the full posteriors, however, we can repeat this calculation with different confidence levels. We find that $\mathcal{I}_{\boldsymbol{\Omega}}$ can be biased by a factor of a few in both directions, depending on what confidence level is used; the 90% interval just happens to produce a particularly close number.

In calculating the odds from Eq. (31), there are large uncertainties on the three rates. However, the rate of short GRB detections by Fermi-GBM is well known and must satisfy $R_{\textrm{\emph{Fermi}}}\approx R_{\textrm{EM}}+R_{\textrm{GW,EM}}$. One could model the uncertainties to produce beaming and volume corrected estimates for these rates (see e.g. Fong et al. (2015); Siellez et al. (2016)). However, for a simple estimate we assume that $R_{\textrm{GW,EM}}$ and $R_{\textrm{GW}}$ are of similar magnitude, and take $R_{\textrm{EM}}$ to be approximately $R_{\textrm{\emph{Fermi}}}=0.124$ per-day (Abbott et al., 2017a). Then using $[\Delta t^{\mathrm{min}},\Delta t^{\mathrm{max}}]=[-1,5]$ s (as used in the Abadie et al. (2012) search) in Eq. (31), the odds, including the spatial overlap, are $\mathcal{O}_{\textrm{C}/\textrm{S}\textrm{S}}(D_{\textrm{GW}},D_{\textrm{EM}})\gtrsim 10^{6}$: the odds provide decisive evidence that the two detections originate from the same event.

These numbers are consistent with the p-values estimated in (Abbott et al., 2017a): the time overlap dominates, the spatial part is small but supports the hypothesis, and the overall factor is highly significant (the total p-value was found to be $5\times 10^{-8}$ (Abbott et al., 2017a)).

There are parallels that can be drawn between the odds calculated in Sec. 3.1 and the p-value approach of (Abbott et al., 2017a); namely, the spatial overlap Eq. (21) with the $\mathcal{S}$ statistic and the form of the temporal overlap (i.e. inversely proportional to the background rate).

However, the two methods are not equivalent and the numerical values themselves cannot be directly compared, as they answer different questions. The odds is exactly our relative degrees of belief for the common- vs. distinct-source hypotheses, given the assumptions made in the calculation; the p-value tests whether the data is consistent or not with the null (distinct-source) hypothesis, and is typically interpreted as the rate at which a rule for deciding the significance of a joint detection leads to false positives (Finn et al., 1999).

Finally, a more practical difference is that while Eq. (21) can be directly interpreted as a Bayes factor for the spatial overlap, interpreting the $\mathcal{S}$ statistic requires numerical calculation of the background by randomly rotating a set of observed short GRB sky localizations.