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

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.

The gstlal pipeline natively determines the functions $p_{0}(x)$ and $p_{1}(x)$ for its detection statistic $x$. For this analysis a threshold of $x_{\mathrm{min}}=5$ was applied, which is sufficiently low that the trigger density is dominated by terrestrial triggers near threshold. There were $M=15\,848$ triggers observed above this threshold in the 17 days of observation time analyzed by gstlal.

For pycbc, the quantity $x^{\prime}$ is the re-weighted SNR detection statistic.¹¹1When quoting pipeline-specific values we distinguish pycbc quantities with a prime. We set a threshold $x_{\mathrm{min}}^{\prime}={8}{}$, above which $M^{\prime}=270$ triggers remain in the search. We use a histogram of triggers collected from time-shifted data to estimate the terrestrial trigger density, $p_{0}\left(x^{\prime}\right)$, and a histogram of the recovered triggers from the injection sets described in Section 2.2 of the Letter to estimate the astrophysical trigger density, $p_{1}\left(x^{\prime}\right)$. These estimates are shown in Figure 1. The uncertainty in the distribution of triggers from this estimation procedure is much smaller than the uncertainty in overall rate from the finite number statistics (see, for example, Figure 5). The empirical estimate is necessary to properly account for the interaction of the various single- and double-interferometer thresholds in the pycbc search (Abbott et al., 2016c). At high SNR, where these thresholds are irrelevant, the astrophysical triggers follow an approximately flat-space volumetric density (see Section 3) of

$$p_{1}(x^{\prime})\simeq\frac{3x_{\mathrm{min}}^{\prime 3}}{x^{\prime 4}},$$

(1)

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.

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:

$$\dfrac{\mathrm{d}{N}}{\mathrm{d}{t\mathrm{d}x}}(t,x)=R_{0}(t)p_{0}(x;t)+R_{1}(t)V(t)p_{1}(x;t),$$

(2)

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.²²2The 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

$$\dfrac{\mathrm{d}{N}}{\mathrm{d}{t}}\equiv\int\mathrm{d}x\,\dfrac{\mathrm{d}{N}}{\mathrm{d}{t\mathrm{d}x}}=R_{0}(t)+R_{1}V(t).$$

(3)

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

$$\mathcal{L}=\left\{\prod_{j=1}^{M}\dfrac{\mathrm{d}{N}}{\mathrm{d}{t\mathrm{d}x}}\left(t_{j},x_{j}\right)\exp\left[-\delta_{j}\dfrac{\mathrm{d}{N}}{\mathrm{d}{t}}\left(t_{j}\right)\right]\right\}\\ \times\left\{\prod_{k=1}^{N_{T}-M}\exp\left[-\delta_{k}\dfrac{\mathrm{d}{N}}{\mathrm{d}{t}}\left(t_{k}\right)\right]\right\}$$

(4)

(cf. Farr et al. (2015, Eq. (21)) or Loredo & Wasserman (1995, Eq. (2.8))).³³3There 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

$$\mathcal{L}=\prod_{j=1}^{M}\left[\dfrac{\mathrm{d}{N}}{\mathrm{d}{t\mathrm{d}x}}\left(t_{j},x_{j}\right)\right]\exp\left[-N\right],$$

(5)

where

$$N=\int\mathrm{d}t\,\dfrac{\mathrm{d}{N}}{\mathrm{d}{t}}$$

(6)

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

As discussed in Section 1, in our search we observe that $R_{0}$ remains approximately constant and that $p_{0}$ retains its shape over the observation time discussed here; this assumption is used in our search background estimation procedure (Abbott et al., 2016c). The astrophysical distribution of triggers is universal (Section 3) and also time-independent. Finally, the detector sensitivity is observed to be stable over our ${16}~\mathrm{days}$ of coincident observations, so $V(t)\simeq\mathrm{const}$ (Abbott et al., 2016b). We therefore choose to simply ignore the time dimension in our trigger set. This generates an estimate of the rate that is sub-optimal (i.e. has larger uncertainty) but consistent with using the full data set to the extent that the detector sensitivity varies in time; since this variation is small, the loss of information about the rate will be correspondingly small. We do capture any variation in the sensitivity with time in our Monte-Carlo procedure for estimating $\left\langle VT\right\rangle{}$ that is described in Section 2.2 of the Letter.

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

$$\bar{\mathcal{L}}\equiv\int\left[\prod_{j}\mathrm{d}t_{j}\right]\,\mathcal{L}\\ =\prod_{j}\left[\Lambda_{0}p_{0}\left(x_{j}\right)+\Lambda_{1}p_{1}\left(x_{j}\right)\right]\exp\left[-\Lambda_{0}-\Lambda_{1}\right],$$

(7)

where

$$\Lambda_{0}p_{0}(x)=\int\mathrm{d}t\,R_{0}(t)p_{0}\left(x;t\right),$$

(8)

and

$$\Lambda_{1}p_{1}(x)=\int\mathrm{d}t\,R_{1}V(t)p_{1}\left(x;t\right),$$

(9)

with

$$\int\mathrm{d}x\,p_{0}(x)=\int\mathrm{d}x\,p_{1}(x)=1.$$

(10)

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.