Supplement: The Rate of Binary Black Hole Mergers Inferred from Advanced LIGO Observations Surrounding GW150914

Search Pipelines

Both the pycbc and gstlal pipelines are based on matched filtering against a bank of template waveforms. See Abbott et al. (2016c) for a detailed description of the pipelines in operation around the time of GW150914; here we provide an abbreviated description.

In the pycbc pipeline, the single-detector signal-to-noise ratio (SNR) is re-weighted by a chi-squared factor (Allen, 2005) to account for template-data mismatch (Babak et al., 2013); the re-weighted single-detector SNR are combined in quadrature to produce a detection statistic for search triggers.

The gstlal pipeline’s detection statistic, however, is based on a likelihood ratio (Cannon et al., 2013, 2015) constructed from the single-detector SNR and a signal-consistency statistic. An analytic estimate of the distribution of astrophysical signals in multiple-detector SNR and signal consistency statistic space is compared to a measured distribution of single-detector triggers without a coincident counterpart in the other detector to form a multiple-detector likelihood ratio.

Both pipelines rely on an empirical estimate of the search background, making the assumption that triggers of terrestrial origin occur independently in the two detectors. The background estimate is built from observations of single-detector triggers over a long time (gstlal) or through searching over a data stream with one detector’s output shifted in time relative to the other’s by an interval that is longer than the light travel time between detectors, ensuring that no coincident astrophysical signals remain in the data (pycbc). For both pipelines it is not possible to produce an instantaneous background estimate at a particular time; this drives our choice of likelihood function as described in Section 2.

but they deviate from this at smaller SNR due to threshold effects in the search.

For the pycbc pipeline, a detection statistic $x^{\prime}\geq{10.1}{}$ corresponds to an estimated search false alarm rate (FAR) of one per century.

Derivation of Poisson Mixture Model Likelihood

In this section we derive the likelihood function in Eq. (3) of the Letter. Consider first a search of the type described in Section 1 over $N_{T}$ intervals of time of width $\delta_{i}$, $\left\{i=1,\ldots,N_{T}\right\}$. Triggers above some fixed threshold occur with an instantaneous rate in time and detection statistic $x$ given by the sum of the terrestrial and astrophysical rates:

where $R_{0}(t)$ is the instantaneous rate (number per unit time) of terrestrial triggers, $R_{1}(t)$ is the instantaneous rate density (number per unit time per unit comoving volume) of astrophysical triggers, $p_{0}$ is the instantaneous density in detection statistic of terrestrial triggers, $p_{1}$ is the instantaneous density in detection statistic of astrophysical triggers, and $V(t)$ is the instantaneous sensitive comoving redshifted volume (Abbott et al., 2016a, see also Eq. (15) of the Letter) of the detectors to the assumed source population. The astrophysical rate $R_{1}$ is to any reasonable approximation constant over our observations so we will drop the time dependence of this term from here on.The astrophysical rate can, in principle, also depend on redshift, but in this paper we assume that the BBH coalescence rate is constant in the comoving frame. Note that $R_{0}$ and $R_{1}$ have different units in this expression; the former is a rate (per time), while the latter is a rate density (per time-volume). The density $p_{1}$ is independent of source parameters as described in Section 3. Let

If the search intervals $\delta_{i}$ are sufficiently short, they will contain at most one trigger and the time-dependent terms in Eq. (2) will be approximately constant. Then the likelihood for a set of times and detection statistics of triggers, $\left\{(t_{j},x_{j})|j=1,\ldots,M\right\}$, is a product over intervals containing a trigger (indexed by $j$) and intervals that do not contain a trigger (indexed by $k$) of the corresponding Poisson likelihoods

(cf. Farr et al. (2015, Eq. (21)) or Loredo & Wasserman (1995, Eq. (2.8))).There is a typo in Eq. (2.8) of Loredo & Wasserman (1995). The second term in the final bracket is missing a factor of $\delta t$. Now let the width of the observation intervals $\delta_{i}$ go to zero uniformly as the number of intervals goes to infinity. Then the products of exponentials in Eq. (4) become an exponential of an integral, and we have

is the expected number of triggers of both types in the total observation time $T$.

If we ignore the trigger time, then the appropriate likelihood to use is a marginalization of Eq. (5) over the $t_{j}$. Let

If we assume that $R_{1}$ is constant in (comoving) time, and measure $p_{1}(x)$ by accumulating recovered injections throughout the run as we have done, then this expression reduces to the likelihood in Eq. (3) of the Letter. A similar argument with an additional population of triggers produces Eq. (10) of the Letter.

The procedure for estimating $p_{0}(x)$ in the pycbc pipeline also provides an estimate of the mean number of background events per experiment $\Lambda_{0}$ (Abbott et al., 2016c). The procedure for estimating $p_{0}$ used in the gstlal pipeline, however, does not naturally provide an estimate of $\Lambda_{0}$; instead gstlal estimates $\Lambda_{0}$ by fitting the observed number of triggers to a Poisson distribution. We have chosen to leave $\Lambda_{0}$ as a free parameter in our canonical analysis with a broad prior and infer it from the observed data, rather than using the pycbc background estimate to constrain the prior, which would result in a much narrower posterior on $\Lambda_{0}$. This is equivalent to the gstlal procedure for $\Lambda_{0}$ estimation in the absence of signals; the presence of a small number of signals in our data here do not substantially change the $\Lambda_{0}$ estimate due to the overwhelming number of background triggers in the data set.

Using a broad prior on $\Lambda_{0}$ is conservative in the sense that it will broaden the posterior on $\Lambda_{1}$ from which we infer rates. However, because there are so many more triggers in both searches of terrestrial origin than astrophysical there is little correlation between $\Lambda_{0}$ and $\Lambda_{1}$, and so there is little difference between the posterior we obtain on $\Lambda_{1}$ and the posterior we would have obtained had we implemented the tight prior on $\Lambda_{0}$. Figure 2 shows the two-dimensional posterior we obtain from Eq. (5) of the Letter on $\Lambda_{0}$ and $\Lambda_{1}$.

We have checked that using a $\delta$-function prior

in the pycbc analysis that is the result of the pipeline $\Lambda_{0}$ estimate from timeslidesWhile the statistical uncertainty on the pipeline $\Lambda_{0}$ estimate is not precisely zero, $\sigma_{\Lambda_{0}}/\Lambda_{0}\lesssim 10^{-3}$, it is so small that a $\delta$-function prior is appropriate. (Abbott et al., 2016c) and using a looser prior that is the result of a gstlal estimate on a single set of time-slid data produces no meaningful change in our results. Figure 3 shows our canonical rate posterior inferred with the pycbc $\Lambda_{0}$ prior in Eq. (11) and our canonical broad prior.

Universal Astrophysical Trigger Distribution

This result follows from the fact that the expected value of the SNR in a matched-filter search for compact binary coalescence (CBC) signals scales inversely with transverse comoving distance (Hogg, 1999):

where $A$ is an amplitude factor that depends on the intrinsic properties (source-frame masses and spins) of the source, the detector sensitivity expressed as a noise power spectral density $S(f)$ as a function of observer frequency and redshift $z$, and $B$ is an angular factor depending on the location of the source in the sky and the relative orientations of binary orbit and detector. The redshift enters $A$ only through shifting the source waveform to lower frequency at higher redshift, changing $A$ because the sensitivity varies with observer frequency $f$. For the redshifts to which we are sensitive to BBH in this observation period this effect on $A$ is small.

If we assume that the distribution of source parameters is constant over the range of distances to which we are sensitive, and ignore the small redshift-dependent sensitivity correction mentioned above, then the distribution of SNR will be governed entirely by the distribution of distances of the sources, which, in the local universe is approximately

yielding the distribution of SNR given in Eq. (12).

Both the pycbc and gstlal pipelines use goodness-of-fit statistics in addition to SNR and employ a more complicated system of thresholds than this simple model, but the empirical distribution of detection statistics remains, to an approximation suitable for our purposes, independent of the source parameters. Figure 4 shows the distribution of recovered detection statistics for the various injection campaigns with varying source distribution used to estimate sensitive time-volumes in the pycbc pipeline. In each injection campaign $\mathcal{O}(1000)$ signals were recovered. For loud signals, the detection statistic is proportional to SNR in this pipeline, and the distribution is not sensitive to the complicated thresholding in the pipeline, so we recover Eq. (12); for quiet signals the interaction of various single-detector thresolds in the pipeline causes the distribution to deviate from this analytic approximation, but it remains independent of the distribution of sources. Note that the empirical distribution of detection statistics, not the analytic one, forms the basis for $p_{1}$, the foreground distribution used in this rate estimation work.

To quantify the deviations from universality, we have preformed two-sample Kolmogorov-Smirnov (KS) tests between all six pairings of the sets of detections statistics recovered in the injection campaigns described in Sections 2 and 3 of the Letter and featured in Figure 4. The most extreme KS $p$-value occurred with the comparison between the injection set with BBH masses drawn flat in $\log m$ and the one with masses drawn from a power law (both described in Section 3 of the Letter); this test gave a $p$-value of $0.013$. Given that we have performed six identical comparisons we cannot reject the null hypothesis that the empirical distributions used for rate estimation from the pycbc pipeline are identical even at the relatively weak significance $\alpha=0.05$. Certainly any differences in detection statistic distribution attributable to the BBH population are far too small to matter with the few astrophysical signals in our data set (compared with $\mathcal{O}(1000)$ recovered injections in each campaign).

Because the distribution of detection statistics is, to a very good approximation, universal, we cannot learn anything about the source population from the detection statistic alone; we must instead resort to parameter estimation (PE) followup (Veitch et al., 2015; Abbott et al., 2016e) of triggers to determine their parameters. The parameters of the waveform template that produced the trigger can be used to guess the parameters of the source that generated that trigger, but the bias and uncertainty in this estimate are very large compared to the PE estimate. We therefore ignore the parameters of the waveform template that generated the trigger in the assignment of triggers to BBH classes.

Count Posterior

We impose a prior on the $\Lambda$ parameters of:

The posterior on expected counts is proportional to the product of the likelihood from Eq. (3) of the Letter and the prior from Eq. (15):

For estimation of the Poisson rate parameter in a simple Poisson model, the Jeffreys prior is $1/\sqrt{\Lambda}$. With this prior, the posterior mean on $\Lambda$ is $N+1/2$ for $N$ observed counts. With a prior proportional to $1/\Lambda$ the mean is $N$ for $N>0$, but the posterior is improper when $N=0$. For a flat prior, the mean is $N+1$. Though the behaviour of the mean is not identical with our mixture model posterior, it is similar; because we find $\left\langle\Lambda_{1}\right\rangle\gg 1/2$, the choice of prior among these three reasonable options has little influence on our results here.

For the pycbc data set we find the posterior median and $90\%$ credible range $\Lambda_{1}={{3.2}^{+{4.9}}_{-{2.4}}}$ above our threshold. For the gstlal set we find the posterior median and 90% credible range $\Lambda_{1}={{4.8}^{+{7.9}}_{-{3.8}}}{}$. Though we have only one event (GW150914) at exceptionally high significance, and one other at marginal significance (LVT151012), the counting analysis shows these to be consistent with the possible presence of several more events of astrophysical origin at lower detection statistic in both pipelines.

The thresholds applied to the pycbc and gstlal triggers for this analysis are not equivalent to each other in terms of either SNR or false alarm rate; instead, both thresholds have been chosen so that the rate of triggers of terrestrial origin ($\Lambda_{0}p_{0}$) dominates near threshold. Since the threshold is set at different values for each pipeline, we do not expect the counts to be the same between pipelines.

The estimated astrophysical and terrestrial trigger rate densities (Eq. (1) of the Letter) for pycbc are plotted in Figure 5. We select triggers from a subset of the search parameter space (i.e. our bank of template waveforms) that contains GW150914 as well as the mass range considered for possible alternative populations of BBH binaries in Section 3 of the Letter. There are $M^{\prime}={270}{}$ two-detector coincident triggers in this range in the pycbc search (Abbott et al., 2016c). Figure 5 also shows an estimate of the density of triggers that comprise our data set which agrees well with our inference of the trigger rate.

Based on the probability of astrophysical origin inferred for LVT151012 from the two-component mixture model in Eq. (16) and shown in Figure 6, we introduce a third class of signals and use a three-component mixture model with expected counts $\Lambda_{0}$ (terrestrial), $\Lambda_{1}$ (GW150914-like), and $\Lambda_{2}$ (LVT151012-like) to infer rates in Sections 2.1 of the Letter and 2.2 of the Letter.

We use the Stan and emcee Markov-Chain Monte Carlo samplers (Foreman-Mackey et al., 2013; Stan Development Team, 2015b, a) to draw samples from the posterior in Eq. (5) of the Letter for the two pipelines. We have assessed the convergence and mixing of our chains using empirical estimates of the autocorrelation length in each parameter (Sokal, 1996), the Gelman-Rubin $R$ convergence statistic (Gelman & Rubin, 1992), and through visual inspection of chain plots. By all measures, the chains appear well-converged to the posterior distribution.

Table 1 contains the full results on expected counts and associated sensitive time-volumes for both pipelines.

Calibration Uncertainty

The SNR reported by our searches are quadratically sensitive to calibration errors because they are maximized over arrival time, waveform phase, and a template bank of waveforms (Allen, 1996; Brown & LIGO Scientific Collaboration, 2004). Abbott et al. (2016c) demonstrates that the other search pipeline outputs are also not affected to a significant degree by the calibration uncertainty present during our observing run. Therefore, we ignore effects of calibration on the pipeline detection statistics $x$ and $x^{\prime}$ we use here to estimate rates from the pycbc and gstlal pipelines.

The amplitude calibration uncertainty in the detector results at leading order in a corresponding uncertainty between the luminosity distances of sources measured from real detector outputs (Abbott et al., 2016e) and the luminosity distances used to produce injected waveforms used to estimate sensitive time-volumes in this work. A 10% uncertainty in $d_{L}$ at these redshifts corresponds to an approximately 30% uncertainty in volume. We model this uncertainty by treating $\left\langle VT\right\rangle$ as a parameter in our analysis, and imposing a log-normal prior:

where $\mu$ is the Monte-Carlo estimate of sensitive time-volume produced from the injection campaigns described in Section 2.2 of the Letter and

Since the likelihood in Eqs. (3) of the Letter or (10) of the Letter does not constrain $\left\langle VT\right\rangle$ independently of $R$, sampling over $\left\langle VT\right\rangle$ at the same time as $\Lambda$ and $R$ has the effect of convolving the log-normal distribution of $\left\langle VT\right\rangle$ with the posterior on $\Lambda$ in the inference of $R$. In spite of the 30% relative uncertainty in $\left\langle VT\right\rangle$ from calibration uncertainty, the counting uncertainty on $R$ from the small number of detected events dominates the width of the posterior on $R$.

Analytic Sensitivity Estimate

Generate inspiral–merger–ringdown waveforms in a single detector at various redshifts from the source distribution $s(\theta)$ with random orientations and sky positions.

Using the high-sensitivity early Advanced LIGO noise power spectral density from Abbott et al. (2016f), compute the SNR in a single detector.

Consider a signal found if the SNR is greater than $8$.

estimated from this procedure for systems with various parameters superimposed on the Monte-Carlo estimates from the injection campaign described above.

References