A New Search Pipeline for Compact Binary Mergers: Results for Binary Black Holes in the First Observing Run of Advanced LIGO

I Introduction

The LIGO and Virgo observatories reported the detection of several gravitational wave (GW) events from compact binary coalescence in their First and Second Observing Runs (O1 and O2 respectively) Abbott et al. (2018). These detections required technically sophisticated analysis pipelines to reduce the strain data. This is because typical events are buried under the detector noise, and cannot be simply “seen” in raw data at current sensitivities. Hence, any search for signals in the data needs to properly and precisely model the detector noise.

The simplest model is that the detector noise is stationary and Gaussian in nature. Under these assumptions, the best method to detect signals is matched-filtering, which involves creating a bank of possible signals, constructing optimal filters (or templates) for the signals given the noise model, and running the templates over the data. The resulting scores are distributed according to known (chi-squared) distributions in the presence or absence of real signals Jaranowski and Królak (2012).

Unfortunately, both the assumptions underlying matched-filtering fail at some level: the noise statistics vary even on the timescales of the (putative) signals, and there are intermittent non-astrophysical artifacts which are clearly not produced by Gaussian random noise (“glitches”) Cabero et al. (2019), examples of such disturbances can be found in Ref. Zevin et al. (2017). These systematics pollute the distribution of the matched-filtering scores. Moreover, the templates describing different astrophysical signals have finite overlaps, and thus often trigger on the same underlying noise transients. Detectable real events lie in the tails of the score distribution, and hence it is crucial to properly correct for systematics in order to maximize the sensitivity to GW events, and to quote reliable false-alarm rates (FARs).

The official LVC catalog of GW events comprises candidates from two independent pipelines: PyCBC Usman et al. (2016) and GstLAL Messick et al. (2017). Additional analysis of the data was presented in Ref. Nitz et al. (2018). Each of these pipelines has developed solutions for the data complexities described above. In this paper, we describe a new and independent analysis pipeline that we have developed for analyzing the publicly available data from the first observing run of advanced LIGO Vallisneri et al. (2015). Our solutions and implementation choices were guided by the desire to attain, as much as possible, the ideal of the distributions in the Gaussian case, which are easily understood and interpreted.

First, we developed a method to construct template banks that enumerates not over physical waveforms, but over linear combinations of a complete set of basis functions for their phases. Correlations between templates have a uniform and isotropic metric in this space.

Second, when dealing with systematics, we use procedures with analytically tractable behavior in the case of Gaussian random noise, which enables us to set thresholds based on well-defined probabilities. We developed a simple method to empirically correct for the non-stationary nature of the detector noise (PSD drift). Under this procedure, segments of data with no apparent glitches produce trigger scores with perfect chi-squared distributions. At the first pass, we attempt to veto out residual “glitches” using a collection of simple tests (either at the signal-processing level or after triggering), while still using the matched-filtering scores as the ranking statistics to leave the Gaussian “floor” untouched. We also developed methods to condition masked data in a way that guarantees that the following matched filtering step would have zero response to the masked data segments.

Finally, we estimate the background of coincident triggers between the two detectors using time slides (akin to PyCBC). Our pipeline includes methods to use the information from background triggers to combine physical triggers from different detectors in a statistically optimal manner for distinguishing astrophysical events from noise transients.

Our paper is organized as follows: Section II provides an overview of the stages in the pipeline. Section III expands upon each of the stages while omitting derivations and precise details, which we present in accompanying papers Roulet et al. ; Zackay et al. ; Venumadhav et al. . In Section IV we present the results of our search for binary black hole mergers in O1.

II Pipeline stages

We construct our pipeline in several stages, which are organized as follows:

Construction of a template bank: We divide the mergers into banks with logarithmic spacing in the chirp mass, and analyze each bank separately. Section III.1 provides further details on the underlying method, and the properties of the resulting banks.

Analysis of single detector data: We first analyze the data streams from the Hanford (H1) and Livingston (L1) detectors separately, as follows:

We iteratively whiten the data stream, perform several tests to detect and remove bad data segments (“glitches”), and condition the remaining data to preserve astrophysical signals. Sections III.3 and III.4 describe this procedure.

We correct for the non-stationary nature of the noise (PSD drift), which if untreated, systematically pollutes the connection between the matched-filtering scores and probability. Section III.6 provides more details.

We generate matched-filtering overlaps for the waveforms in our banks with the whitened data stream, apply the PSD drift correction, and record triggers whose matched-filtering scores are above a chosen threshold (Section III.5).

Coincidence analysis between detectors: We analyze triggers that are coincident in H1 and L1. In Section III.7, we describe how we collect coincident triggers with combined incoherent score above a threshold, at both physical (candidates) and unphysical (background) time delays.

Refining on a fine grid: We refine the parameters of the candidates and the background on a finer grid around the triggers in order to account for template bank inefficiency, and allow room for more stringent signal quality vetoes.

Trigger quality vetoes: We apply vetoes on the triggers based on the signal quality, as well as the data quality. The vetoes have to be applied either at the single-detector level, to avoid biasing the calculation of the coincident background using time slides. Section III.9 lists the vetoes we applied to the triggers.

Estimating the significance of candidates: We use the set of background triggers to estimate the FAR for the candidates at physical lags between H1 and L1. We do this in two stages:

We first compute a ranking score that is purely a function of the incoherent scores of the triggers, under the assumption that the noise processes that produce the background are independent between detectors (Section III.10).

Section III.11 describes our coherent score, which adds all the information encapsulated in the phase, amplitude, relative sensitivity and arrival time differences between the detectors to create our final candidate ranking statistic.

Section III.12 describes how we construct an estimate for the probability of a coincident event being of astrophysical origin given an astrophysical event rate.

III Concise description of the pipeline stages

We perform our search by matching the strain data to a discrete set of waveform templates that sufficiently closely resemble any gravitational wave signal within our target parameter space. We target our search at coalescing binary black holes (BBH), defined here as compact binary objects with individual masses between $3\,M_{\odot}$ and $100\,M_{\odot}$ and with aligned spins. We allow spin magnitudes up to $\absolutevalue{\chi_{1,2}}<0.85$. We restrict the mass ratios to be $q^{-1}<18$.

As described in Ref. Roulet et al. , we construct five BBH template banks (BBH 0-4) that together span this target parameter space, and conduct a separate search within each of them. The banks are defined by regions in the plane of component masses, as shown in Fig. 1. We place the bounds between adjacent banks at $\mathcal{M}=\{5,10,20,40\}\,M_{\odot}$, where $\mathcal{M}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ is the chirp mass and $m_{1,2}$ are the individual masses. We find several motivations for dividing the search. The low-mass banks have many more templates than the heavier banks, and thus they inherently have a larger look-elsewhere penalty. Dividing the search prevents this from strongly affecting the sensitivity of the high-mass searches: in this way, on astrophysical grounds we might expect roughly comparable numbers of signals in each bank, regardless of the largely different number of templates they have. Moreover, this splitting enables us to discriminate between the different types of background events that each search is subject to. The different duration of the signals in each bank will require us to use different thresholds when masking bad data segments (see Section III.3). The prevalence of non-Gaussian glitches will be different in each bank and thus the score we assign to events with the same signal-to-noise ratio (SNR) is different in each bank (see Section III.10). Table 1 summarizes the template bank parameter ranges and sizes.

The template bank needs to be effectual, that is, to guarantee a sufficiently high match between a GW waveform and at least one template in the bank. We define the inner product between waveforms $h_{i},h_{j}$

where $S_{n}(f)$ is the one-sided noise power spectral density (PSD) of the detector and a tilde indicates a Fourier transform into the frequency domain. It is used to define the match

throughout this section we assume that all waveforms are normalized to $(h\mid h)=1$. We assess the effectualness $E$ of each bank by computing the best match with ${10}^{4}$ random waveforms in its target parameter space. We apply the down-sampling and sinc-interpolation described in Section III.5 and the waveform optimization described in Section III.8 to the test waveforms, to properly simulate the search procedure. We report the effectualness of the banks in Table 1. When designing banks, we set the reference PSD to be the aLIGO_MID_LOW PSD (LIGO Scientific Collaboration, 2018), which is representative of O1.

In order to correct the PSD drift at manageable computational cost, our search pipeline requires that the frequency domain templates, of the form

share a common amplitude profile $\overline{A}(f)$ (see Section III.6) and differ only in the phase $\psi(f)$. In order to avoid excessive loss of effectualness due to this approximation, we split each bank into several subbanks, each of which is assigned a different $\overline{A}(f)$ profile. We use the method of “stochastic placement” to determine as many subbanks as needed to guarantee that every waveform within the target parameter range has an amplitude match,

with at least one subbank. The resultant divisions into subbanks are color-coded in Fig. 1.

The remaining task is to place templates in each subbank to efficiently capture the possible phase shapes $\psi(f)$. We achieve that with a geometric approach, where we use the mismatch between templates to define a mismatch distance, which quantifies the similarity between any two waveforms. We abandon the physical parameters as a description of the templates in favor of a new basis of coordinates $\bm{c}$, in which the mismatch distance induces an Euclidean metric. We then set up a regular grid in this space. Our templates take the form

where $\overline{\psi}(f)$ is the average phase, and $\{\psi_{\alpha}(f)\}$ are phase basis functions which are orthonormalized such that the mismatch distance satisfies

An input set of physical waveforms representing the target signals are used, first to define the subbanks and then to determine the appropriate phase basis functions. The input waveforms may be generated with any frequency-domain model; we use the IMRPhenomD approximant (Khan et al., 2016). The phase basis functions are found from a singular value decomposition of the input waveforms which identifies the minimal set of linear independent components that need to be kept. A small number of basis functions are enough to approximate all possible phases to sufficient accuracy. All banks require five linearly independent bases or fewer, with about half of them having only three or fewer. While the coefficient for the lowest order bases may vary over a range of several hundred units, the coefficients for the highest order bases vary within narrow ranges, sometimes by less than one unit.

III.2 Loading and preprocessing the data

The next step after loading the data is to estimate its PSD. We use Welch’s method Welch (1967), in which several overlapping chunks of data are windowed and their periodograms are averaged (we use the implementation in scipy.signal with a Hann window). We make our PSD estimation robust to bad data by (a) disregarding chunks that overlap with segments that were marked by LIGO’s quality flags, and (b) averaging using the median instead of the mean (see Appendix B of Ref. Allen et al. (2012)).

III.3 Identifying bad data segments

Advanced LIGO data contains intermittent loud disturbances that are not marked by the provided data quality flags. We need to flag and remove these segments to prevent them from polluting our search, while taking care to preserve astrophysical signals of interest. This is the fourth analysis of the data, and hence we assume that any new signals we find will have an integrated matched filter SNR $\rho<30$ in a single detector. This assumption allows us to bound the influence of a true signal on our procedure.

We devise several complementary tests to flag bad data segments. We design our tests to satisfy the following conditions:

The test statistics have analytically known distributions for Gaussian random noise.

The thresholds are set to values of the test statistics achieved by waveforms with single-detector $\rho=30$ in noiseless data. Signals at this SNR have a probability of $\simeq 0.5$ of triggering a single test in the presence of Gaussian random noise. We found empirically that signals satisfying $\rho\leq 20$ are almost always retained.

If the above thresholds are too low, they are adjusted so that a single test is triggered at most once per five files due to Gaussian random noise alone. This is important for template banks with long waveforms.

These conditions ensure that we are sensitive to gravitational waves while not over-flagging the data. It is important that the tests be done at the single-detector level in order to avoid biasing the calculation of the background using time slides.

Our tests trigger on the following anomalies: (a) outliers in the whitened data-stream, (b) sine-Gaussian transients in particular bands, (c) excess power localized to particular bands and timescales, and (d) excess power (summed over frequencies) on particular timescales. We picked timescales and frequency bands for the tests based on inspecting the spectrograms of the bad segments; Table 2 details the choices.

The data has spectral lines at which the PSD is several orders of magnitude higher than in the continuum. The power in these lines often significantly varies in a non-Gaussian manner within a single file. The lines do not contribute to the matched-filtering overlap, since the PSD is effectively infinite at their frequencies. Hence it is preferable that varying lines do not trigger our tests.

We detect sine-Gaussian artifacts in a given band by matched-filtering with a complex waveform that saturates the time-frequency uncertainty principle and contains most of its power in the band. We apply notch filters to the sine-Gaussian template to remove any overlap with spectral lines. We flag any outliers in the matched-filtering results above a threshold defined to satisfy the aforementioned conditions (see second paragraph of Sec. III.3), which is a procedure safe to any relevant events)

We detect excess power using a spectrogram (computed using the spectrogram function in scipy.signal with its default Tukey window). We sum the power in the frequency ranges of interest, disregarding frequency bins that overlap with varying lines. For Gaussian random noise, this sum has a chi-squared distribution. This is not achieved in practice unless correcting for the effects of PSD changes. We make the excess power statistic robust to the drifting of the PSD by comparing the instantaneous excess power with with a local moving-average power baseline.

The simplest test is to look for outliers in the whitened strain, since individual samples should be independent and normally distributed with unit variance. We flag segments of whitened data, with a safety margin in time, around outliers above a chosen threshold.

Whenever one or more of these tests fire, we excise the offending segments (which we refer to as “holes”) and inpaint the raw data within as described in Section III.4. In practice, we observe that the outlier test often does not catch all of the “bad” data, in which case the inpainted and whitened data contain further outliers. Hence, we iterate over the “identify bad segments, inpaint, whiten” cycle multiple ($<7$) times, increasing the safety margin in time by successively larger multiples of 0.1 s, until the process converges.

We treat any part of the data that was marked with any of the LIGO quality flags as if it contained large disturbances. After all the data quality tests done in this section, we are left with roughly 46 days of coincident on-time between the detectors, with slight changes from bank to bank, as all the test thresholds are waveform dependent.

III.4 Inpainting bad data segments

The matched-filtering score for a template $h$ with data $d$ with a noise covariance matrix $C$ is:

To deal with this problem, if we consider a fraction of the data of length $N_{d}$ in which we have masked $N_{h}$ samples, we filter the data with a filter $F$ and define a new score by:

where the matrix $W$ has one column of length $N_{d}$ for every sample that is masked with all the entries zero except for a one at the position of the masked sample and $M$ is the $N_{h}\times N_{h}$ matrix $M=W^{T}\,C^{-1}\,W$. The computationally expensive part of this filtering procedure is to invert the matrix $M$.

Figure 2 shows an example of a small section of the data containing a “glitch” artifact. We show the difference between ‘gating’ the bad data by applying a window function to it, and creating a hole and inpainting it with the algorithm we described. We can see that gating substantially changes the standard deviation of the samples in the hole and the few seconds surrounding it, which can potentially create spurious triggers, and can damage any real signals that happen to be in the data at the same time. In our method, the “blued” data is set to be identically zero inside the hole.

III.5 Matched filtering

Given the whitened, hole-filled data, we compute the overlaps with all templates in the template bank, and register the times and templates when the ${\rm SNR}^{2}$ is above a triggering threshold. The choice of the threshold was driven by the requirement to produce a manageable number of triggers per file, and was generally in the range $20<{\rm SNR}^{2}_{\rm thresh}<25$ for the various banks and subbanks.

III.6 Applying corrections due to the varying Power Spectral Density of the Noise

While at first sight it may seem impossible to both capture the width of the lines and track the fast variation in the PSD, we accomplish it by correcting the first order effect of PSD mis-estimation on time-scales that are as short as the PSD changes, to precision of $\sim 1\%$.

This correction is basically a local estimate of the standard deviation of the overlaps, and is derived (along with some other nice properties it has) in Ref. Zackay et al. . In Figure 4, we present a histogram of the distribution of the local variance estimates. Notice the large deviations from unity in both directions. We note that the tail reaches values as high as $1.5$; at such high values, there are visible disturbances in the spectrogram, sometimes referred to as glitches. However, at values in the range $[0.85,1.2]$, the data mostly behaves in a regular fashion, and there is no apparent sign something bad is going on in the spectrogram of the data. These changes can cause substantial loss of sensitivity in binary coalescence analyses that neglect this effectAfter this manuscript was made public, we were informed that fluctuations in the SNR integral (due to short-timescale variations in the PSD) at comparable levels were previously noted, but the mitigation steps were not incorporated into the search pipelines used in the catalog paper (Thomas Dent, private communication)..

To illustrate why correcting for these variance estimates is crucial for determining the exact significance of a candidate event, we point out that the most economic way of creating a (spurious) $\rho=8$ event is to wait for a lucky time where the PSD mis-estimation is large (say, 1.2), and then create a (genuine) $\rho=7.3$ fluctuation. In Figure 5, we see the tail of the trigger distribution is substantially inflated if the PSD drift is not corrected.

III.7 Coincidence Analysis of the two detectors

We then take each remaining trigger, and insert it into a dictionary according to the template key. This would allow us to immediately find all the times at which this template triggered. Using queries to the dictionary, we find all the pairs of triggers that belong to either the background or the foreground group, and pass the threshold $\rho_{\rm collect}$. This threshold depends on the bank via computing the Gaussian noise threshold for obtaining one significant event per O1, and then multiplied by the bank effectualness, to guarantee that every trigger that can acquire the one-per-O1 significance after optimization is included.

We now view the H1 component of all pairs of triggers and group them to groups of 0.1s. We use the less stringent version of the veto to vet the trigger with the highest SNR in each group, and upon failure discard the entire group (the logic here is that similar triggers are all passing or failing the veto together). We do the same for the L1 component of all remaining trigger pairs.

We then optimize every trigger by computing the overlaps with the data of every template in the sub-grid $\bm{c}$ values (see Sec. III.8). We further sinc interpolate with a long support to obtain further time resolution for the overlaps. We then choose the sub-grid template that maximizes the quadrature sum of the single detector SNRs. This trigger pair is now vetoed with the stringent veto. If a trigger pair passes all these, it is registered.

III.8 Refining triggers on a finer template grid

The template-bank is organized as a regular grid, which facilitates refinement in places of interest. This enables us to squeeze more sensitivity and imitate the strategy of a continuous template bank, which is more objective than an arbitrarily chosen grid. The effectualness achieved by the top 99.9% of injections with the template banks used for the search varies between 0.9 and 0.96. Refining the grid by a factor of two in each dimension would bring it to $>0.96$ in all cases, but would also substantially increase the number of waveforms in the bank (which in turn increases the computational complexity and memory requirements of our search). We therefore take the approach of refining every candidate and background trigger pair. Since we know the maximum amount of SNR increase that is possible for a real event, we refine all candidates that have a score that is high enough to have a chance of reaching a FAR of 1/O1 after refinement. We greatly speed up the candidate refinement by calculating the likelihood using the relative binning method Zackay et al. (2018) (using the original grid-point trigger as the reference waveform). Table 1 reports the improvement in effectualness achieved by this procedure for our banks.

III.9 Vetoing triggers

The matched-filtering score is the optimal statistic for detecting signals buried in Gaussian random noise. As emphasized in the previous sections, the LIGO strain data is not well-described by purely Gaussian random noise, and hence, the matched-filtering score may be triggered (i.e., pushed above the Gaussian-noise significance threshold) by either transient or prolonged disturbances in the detector. Our pipeline attempts to reject these candidates by identifying bad segments at the preprocessing level (Section III.3), or downweighting the scores by their large (empirically measured) variance (Section III.6). However, this is not enough to bring us down to the Gaussian detection limit, especially for the heavier black hole banks. Thus, we need additional vetoes at the final stage to reject glitches. We use vetoes that are based on either the quality of the neighboring data, as well as that of the signal.

Our most selective vetoes are based on signal quality, and check that the matched-filtering SNR builds up the right way with frequency. We perform the following tests:

We subtract the best-fit waveforms from the data and repeat the excess power tests of Section III.3, but with lower thresholds computed using waveforms with $\rho=3$ (and bounded to fire once per 10 files due to Gaussian noise). Moreover, when we see excess power in a particular band and at a particular time, we only reject candidates with power at the same time in their best-fit waveforms (in order to avoid vetoing candidates due to unrelated excess power).

We split the best-fit waveform into disjoint chunks, and check for consistency between their individual matched-filtering scores. This test is similar in philosophy to the chi-squared veto described in Ref. Allen (2005), but improves upon it by accounting for the mis-estimation of the PSD (which is an inevitable consequence of PSD drift) and by projecting out the effects of small mismatch with the template bank grid.

We empirically find triggers that systematically miss the low-frequency parts of the waveforms, or have large scores at intermediate frequencies. The check described above is agnostic to the way the matched-filtering scores in various chunks disagree, and hence is not the most selective test for these triggers. We reject these triggers by using “split-tests” that optimally contrast scores within two sets of chunks.

The final two tests are the most selective vetoes, and hence their thresholds must be set with care. Our method for constructing template banks enables us to set these thresholds in a rigorous and statistically well-defined manner to ensure a given worst-case false-positive probability, which, accounting for the inefficiency in the bank, is achieved with adversarial template mismatches. Hence we set the worst-case false-positive probability of $10^{-2}$ for each of these tests. The details of the tests, and the methods to set thresholds, are described in Ref. Venumadhav et al. . We note that all hardware injections that triggered passed the single-detector signal-based veto.

The data-quality vetoes are relatively simple in nature, and motivated by segments with excess power (as observed in spectrograms) that slip through the combination of the flagging procedure (of Sec. III.3) and PSD drift correction (of Sec. III.6). The tests are as follows:

Finally, we account for rare cases with significant PSD drifts on finer timescales than the ones used while triggering (described in Section III.6 and Ref. Zackay et al. ). When this PSD drift is statistically significant, we veto coincidence candidates (both at zero-lag and in timeslides) whose combined incoherent scores, after accounting for the finer PSD drift correction, are brought down below our collection threshold.

Figure 6 shows the cumulative effect of our vetoes on the score distribution of the triggers in the BBH 3 bank, which contains short waveforms of heavy binary black hole mergers. Also shown are the hardware injections present in the data stream and GW150914 which belongs to this bank’s chirp mass range. We note that the veto retained every hardware injection in this chirp mass domain that passed the flagging procedure of Section III.3. It is interesting to note that GW150914 does not stand out from the single detector trigger distribution before the application of the veto, and is clearly detected even without resorting to coincidence after it.

III.10 Incoherent Ranking

When constructing a statistic to rank events an important part is $P(\rho_{\rm H}^{2},\rho_{\rm L}^{2}\mid H_{0})$, the probability of obtaining a trigger with squared SNRs $(\rho^{2}_{\rm H},\rho^{2}_{\rm L})$ in each detector under the null hypothesis $H_{0}$. Under the assumption that the noise in both detectors is independent,

If the noise in each detector was Gaussian,

Under this assumption it is optimal to use $\rho_{\rm H}^{2}+\rho_{\rm L}^{2}$ to rank candidate events. Unfortunately this is an invalid assumption for two reasons: firstly, even for Gaussian noise, at high ${\rm SNR}$ the maximization over templates, phase and arrival time leads to

where the constant $c$ depends on the bank dimension. However, in practice this is a minor correction, the more substantial problem is the non-Gaussian tail of the noise, the so-called glitches. In the high-SNR limit $P(\rho\mid H_{0})$ is much larger than the Gaussian prediction.

The non-Gaussian tail in the $\rho$ distribution has an important consequence when combining the scores of multiple detectors. If we were simply to use $\rho_{\rm H}^{2}+\rho_{\rm L}^{2}$ as a score, we would be ranking coincidences in which the trigger in one of the detectors is coming from this non-Gaussian tail, as we would be misjudging its probability by many orders of magnitude.

To correct this problem we empirically determine $\log[P(\rho_{i}\mid H_{0})]$ for each detector. We do so by taking our triggers and ranking them according to decreasing $\rho_{i}$ for each detector $i$. We then model

which is a good approximation for distributions with exponential or polynomial tails. We denote

as a robust approximation of the optimal score. In principle, a parametric model for the probability density might outperform the rank estimate, but practical reasons as too few surviving glitches made such estimates prone to fine tuning. Moreover, at the high SNR parts of the distribution, single-detector glitches find background in many timeslides, which makes it problematic to estimate the uncertainty in any such procedure. For this reason, and to maintain simplicity, we chose to use the rank function as a proxy for the single detector trigger probability distribution function.

In Figure 8 we show the two-dimensional histogram of the background obtained by adding $20\,000$ unphysical time shifts between detectors to the O1 LIGO data (so as to recreate an equivalent of $20\,000$ O1 observing runs) for banks BBH 2 and BBH 3. In the left panels we show the distribution of background triggers using $\rho$ as the score. The tail of non-Gaussian glitches is clearly visible leading to an overproduction of triggers where the SNR in one detector is much larger than in the other. On the right panels we show the distribution of the same triggers but now using our rank score to bin them. The lines of constant probability are now straight. Our sub-threshold candidates in these banks are shown together with GW151012, which is a clear outlier, and with GW151216.

III.11 Coherent Score

In this section we further improve the statistic used to rank candidates by exploiting the information encapsulated in the relative phases, amplitudes and arrival times to the different detectors. We begin with the standard expression:

where $T$ is a template in the continuous template bank. Because the maximization procedure on $T$ is done incoherently, and prior to the application of all these terms, we will drop it from the notation. Note that in principle we should have maximized the full expression, but for practical reasons we decided to do the maximization prior to the coherent analysis. In favor of this approximation stands the fact that to linear order, the phase and time shifts are built to be orthogonal to the template identity Roulet et al. , so the template’s fine optimization is expected to preserve the $\phi$ and $\delta t$ of a candidate to high accuracy. We further develop this expression using Bayes rule (and using some basic independence arguments):

where $n_{i}$ is the momentary response of detector $i$ computed from the measured PSD, PSD drift correction and the ovelap of the waveform with holes using the data of detector $i$. $\Delta\phi$ is the difference between detectors in overlap phase of matched filtering the best-fit $T$ with the data. $\Delta t$ is the difference in arrival time of the maximum score between the detectors. $P(\rho_{\rm H}^{2},\rho_{\rm L}^{2}\mid H_{0})$ was computed using the ranking approximation detailed in Section III.10.

$P(\Delta\phi,\Delta t\mid H_{0})$ is taken to be the uniform distribution by symmetry. Here we note that in principle, $P(\rho_{i}\mid t,H_{0})$ can be non-uniform, if there are bad times where glitches conglomerate. Also, glitches could have a waveform model that prefers a particular phase for a particular template. We currently choose not to introduce these complications (other than the bad times veto applied in Sec. III.9).

P\big{(}\rho^{2}_{\rm H},\rho^{2}_{\rm L},\Delta\phi,\Delta t\,\big{|}\,n_{\rm H}/n_{\rm L},H_{1}\big{)} is measured by drawing samples that are uniformly distributed in volume out to a distance where the expected value of the SNR is four, calculating the detector response, and adding noise with the standard complex normal distribution. Out of these samples, we have created a binned histogram of the observed meaningful values $\Delta t,\Delta\phi,\rho_{\rm H}^{2},\rho_{\rm L}^{2}$; the probability of an observed configuration given the signal hypothesis is proportional to the histogram’s occupancy. The same number of samples is used for all values of $n_{\rm H}/n_{\rm L}$ so that the pipeline’s preference for detecting events with equal response between the detectors could be evaluated. This is very similar to the coherent score used in Nitz et al. (2017).

reflects the changes in sensitivity in the detector as a function of time. Including it allows to analyze different segments of data with very different sensitivities, including multiple runs together (say O1 and O2) while maintaining a consistent detection bar, down-weighting the significance of spurious events from less sensitive detector times. One important note is that once we include this term, the FAR does not have units of inverse time, but units of inverse volume time.

III.12 Determination of FAR

III.13 Determination of the probability of a source being of astrophysical origin

While the FAR is largely agnostic of the astrophysical rates (beyond the use of the model in constructing the detection statistic) and is objectively and accurately measurable through time-slides, it is hard to convert to an assessment of the astrophysical origin of a particular event. Such an assessment depends both on the exact (potentially multidimensional) noise probability density at the event’s location (contrast with the one dimensional cumulative probability density the FAR depends on) and the exact probability density given the astrophysical model, including the unknown rate (also as a function of physical parameters). Essentially, if all exact details in the model were known, the probability of an event being of astrophysical origin would be exactly computable, but in the presence of rate uncertainties, especially when considering the rate as a function of physical parameters, the determination of $p_{\rm astro}$ may be dominated by rate uncertainties and astrophysical prejudice. Nevertheless, the objectivity of $p_{\rm astro}$ to ranking functions and its immunity to the existence of the few last glitches that are left after our heavy vetoing are compelling, and we therefore proceed in computing it.

To do that, we strictly assume all templates inside a bank are equally probable (even though parameter dependant rate differences probably exist). We further assume that the background probability density is uniform in time and phase, an assumption we find is extremely good when the SNR value is in the region where the Gaussian noise is dominant.

We then compute the rate at which we observe such an event in coincidence between the two detectors:

where $T$ is the allowed physical time shift between the detector, and $P(\rho^{2}_{\rm H}\mid H_{0}),P(\rho^{2}_{\rm L}\mid H_{0})$ were fit using

$\alpha_{i}$ and $\beta_{i}$ are fit to the background computed from time-slides in the region close to the $(\rho^{2}_{\rm H},\rho^{2}_{\rm L})$ combination of the event. We find this approximation robust in all cases where the event is close to the detection threshold and when the difference between $\rho^{2}_{\rm H}$ and $\rho^{2}_{\rm L}$ is not big.

using the table constructed in Section III.11. Here, \mathcal{R}_{>100}=\mathcal{R}\big{(}\rho^{2}_{\rm H}+\rho^{2}_{\rm L}>100\,\big{|}\,H_{1},n_{\rm H},n_{\rm L}\big{)} is the astrophysical rate of detecting gravitational wave mergers in the event’s bank, with the detector sensitivity at the time of the event. Because $\mathcal{R}_{>100}$ can be easily estimated and updated using a list of known astrophysical events, it is assumed to be known. We then provide the estimate for the event’s astrophysical origin to be:

For ease of future interpretation of the results, we report in Section IV both $W/\mathcal{R}{({\rm event}\mid H_{0})}$ and the computed $p_{\rm astro}$ using our best knowledge of $\mathcal{R}_{>100}$ at the time of writing.

IV Results of the BBH search

Here we report all the signals and sub-threshold candidates found in the search. We report the FAR in units of “O1” to reflect the fact that there was a volumetric correction factor in the coherent score. If we assume the sensitivity of the first observing run to be roughly constant, then the “O1” unit can be converted to roughly 46 days, the effective coincident time we used in the analysis (that has some variation across banks due to differences in the data flagging thresholds). There was no background trigger with a better coherent score than GW150914, GW151012 and GW151226 in their respective banks, so we obtain only an upper limit on the FAR of $1/($20\,000$\,\text{O1})$ for all of these events, with an effective $p_{\rm astro}=1$ for all of them. We report their recovered squared SNR for each detector. We further found an additional event, GW151216, with a FAR of $1/($52$\,\text{O1})$, reported in greater detail in a companion paper Zackay et al. (2019). These and two additional sub-threshold candidates with FAR of approximately 1/O1 are reported in Table 3.

V Conclusions and Discussion

In this paper we presented an overview of a new and independent pipeline to analyze the publicly available data from the first observing run of Advanced LIGO. We used this pipeline to identify a new gravitational merger event in the O1 data. In companion papers we will provide additional details of our techniques and implementation choices and further characterize our search by providing simple estimates of the space-time volume searched as a function of parameters.

There are several areas for future development and improvements in this pipeline, including precise determination of the merger rate/sensitive volume, analysis of single detector triggers, and triggers with subthreshold candidates in the other detector. For future runs, it also remains to incorporate more than two detectors into the ranking of coincident triggers in our pipeline.

Acknowledgment

We thank the participants of the JSI-GWPAW 2018 Workshop at the University of Maryland, and the Aspen GWPop conference (2019) for constructive discussions and comments.

This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.

TV acknowledges support by the Friends of the Institute for Advanced Study. BZ acknowledges the support of The Peter Svennilson Membership fund. LD acknowledges the support by the Raymond and Beverly Sackler Foundation Fund. MZ is supported by NSF grants AST-1409709, PHY-1521097 and PHY-1820775 the Canadian Institute for Advanced Research (CIFAR) program on Gravity and the Extreme Universe and the Simons Foundation Modern Inflationary Cosmology initiative.

References