A guide to LIGO–Virgo detector noise and extraction of transient gravitational-wave signals

The power spectral density of the noise S_n(f ) is not known a priori and must be estimated from the data. One can perform a complex FFT of the entire data stream around some time to be searched for signals, but that yields only two samples (real and imaginary parts) per frequency bin, hence the variance in the estimate of S_n(f ) in any single frequency bin is large. To overcome this, either some form of averaging is used [63], or a fit is made to a physical model for the spectrum [64]. For example, Welch averaging [65] can be used to reduce the variance in the estimated power spectrum, but at the cost of either reducing the frequency resolution or requiring longer stretches of data. The spectral estimate used to whiten the data in figure 3 was found by applying a Welch average to 1024 s of data centered on GPS time 1126259462 (the nearest integer GPS time to the peak of the GW150914 signal). The data were broken up into overlapping 4 s long chunks, each spaced by 2 s. The data in each chunk was Tukey filtered and Fourier transformed. The power spectrum from all the chunks was then averaged.

Figure 4 compares the power spectrum of the Hanford data shown in figure 3, before and after applying the Tukey window, to the power spectrum estimated using Welch averaging. The non-windowed spectrum is swamped by spectral leakage, and follows a 1/f ² scaling. This scaling results from the abrupt step function at the beginning and end of the data to be Fourier transformed. This non-windowed data chunk arises from multiplying a longer stretch of data by a boxcar (or top hat) window. Thus, when it is Fourier transformed, the result is the convolution of the desired spectrum of the original data with the Fourier transform of the 4s-long boxcar window, i.e. a cardinal sine (sinc) function whose amplitude decreases as 1/f ². Since the noise spectrum rises much more rapidly than 1/f ² towards low frequencies, the entire visible frequency range is then dominated by the leakage from this low-frequency component.

Zoom In Zoom Out Reset image size

When the noise spectrum varies significantly over time other spectral estimation methods have to be used [20, 66]. One approach used in LVC parameter estimation studies is to fit a parametrized spectral model to the data that has a smooth spline component and a collection of Lorentzian lines [64]. In section 5 we also discuss in detail the issue of stationarity and non-stationarity of the data, and the effects this has on the data analysis.

In addition to causing spectral leakage, improper windowing of the data can result in spurious phase correlations in the Fourier transform. Figure 5 shows a scatter plot of the Fourier phase as a function of frequency for the same stretch of data shown in figure 4, both with and without the application of a window function. The un-windowed data shows a strong phase correlation, while the windowed data does not.

Zoom In Zoom Out Reset image size

The degree to which a time series is consistent with being stationary and Gaussian noise can be diagnosed by looking at the distribution of its Fourier transformed frequency samples. If the noise is stationary and Gaussian the real and imaginary parts of the whitened noise in each frequency bin will be a collection of independent and identically distributed (i.i.d.) random variables with zero mean and unit variance: $x \sim {{\mathcal N}}(0,1)$. Departures from stationarity result in correlations between samples in different Fourier bins, while departures from Gaussianity can be identified by comparing the distribution of samples to a unit normal distribution. Loud instrumental noise transients and loud gravitational-wave bursts do contribute to non-stationary and non-Gaussian features, but away from these transient disturbances the LIGO–Virgo data can be approximated as stationary and Gaussian. Figure 6 shows the whitened Fourier amplitudes for a quiet stretch of data from the LIGO-Livingston observatory.

Zoom In Zoom Out Reset image size

The Advanced LIGO [9, 55, 73, 74] and Advanced Virgo [10, 47, 75] detectors use feedback loops to keep the optical cavities on resonance. The strain calibration must thus include models and measurements of all readout electronics, as well as of electronics and transfer functions of actuation hardware that act on the mirrors through multiple points in the suspension systems [76]. As shown in figure 11, there are three main components of the differential arm control loop for Advanced LIGO: the actuation function $A(\,f)$, the sensing function $C(\,f)$, and the digital filters applied, $D(\,f)$. All three are measured and modeled as functions of frequency.

Zoom In Zoom Out Reset image size

The digital filters are known to great precision so the calibration error and uncertainty come from the differences between the model and measurement (including measurement error) of the actuation and sensing functions, A and C. To independently measure the actuation and sensing functions, a pair of beams from auxiliary lasers are reflected off of each test mass mirror, with their intensities modulated at a known frequency and amplitude to actuate with radiation pressure. These auxiliary laser assemblies are referred to as photon calibrators [77]. Once A and C are known, the true differential arm length is extracted and translated to a strain by dividing by the common length of the arm L (4 km for LIGO, 3 km for Virgo), per the following equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{calib} h(t) = \frac{1}{L}\left[\mathcal{C}^{-1}*d_{\rm err}(t)+\mathcal{A}*d_{\rm ctrl}(t)\right] , \nonumber \end{align} \tag{ 5 }$$

where $\mathcal{C}$ and $\mathcal{A}$ are time-domain filters derived from frequency-domain measurements of the actuation function $A(\,f)$ and the sensing function $C(\,f)$.

Note that the gravitational-wave strain can also appear in the common-mode arm length changes, and in changes to the lengths of all degrees of freedom in the detector. However, only the sensing of the differential change in the interferometer's arm lengths is engineered to have low enough instrumental noise to be sensitive to the strain induced by gravitational waves. The other optical lengths are controlled in order to maintain an optimal and linear response to the gravitational-wave strain in the differential degree of freedom.

Calibration measurements are made periodically in each observing run. In addition, to monitor time dependent parameters such as optical gain, cavity pole frequency and actuation strength drifts, several calibration lines are continously injected, at specific frequencies, by applying sinusoidal forces on the test mass mirrors using the photon calibrators; these lines will be present in the raw strain data. The calibration line frequencies are different amongst the detectors.

For the second observing run O2 and onwards, the calibration lines are removed from the calibrated $h(t)$ strain data channel (within the calibration accuracy) [47, 78]; the calibration lines were not removed from the O1 data [74]. Even for O1, the presence of the calibration lines does not affect the search for compact binary coalescence gravitational-wave signals as the amplitude of data at the frequencies of the lines is suppressed via the whitening [79, 80] of the data when the calculation of the detection statistic is made for the data from each detector (see section 8). Similarly, for parameter estimation the presence of the noise spectral density in the likelihood (see section 9) minimizes the influence of spectral lines including calibration lines. Because the spectral lines are narrow, this frequency-domain weighting has a negligible effect on signal searches and parameter estimation and does not lead to any spurious effects such as generation of false candidate events or parameter biases.

For the Advanced Virgo calibration in O2 it was necessary to account for the transfer function of the optical response of the interferometer. This requires a calibration of the longitudinal actuators for the mirrors, still based on the laser wavelength as length reference using the so-called free swinging Michelson configuration described in [81]. In addition, the interferometer's output power as determined by the readout electronics also requires calibration. Calibration measurements are made weekly in each observing run and have shown stable actuation strengths. To monitor the time dependent optical gain and cavity pole frequency, several calibration lines are continously injected, as in LIGO detectors, at specific frequencies, by applying sinusoidal forces on the test mass mirrors using the electro-magnetic actuators. By construction, the calibration lines are removed from the calibrated $h(t)$ strain data channel (within the calibration accuracy). For Advanced Virgo in O2 the gravitational wave strain reconstruction removed the contributions from control signals to the test mass mirror motion. A photon calibrator, namely an auxiliary laser that is used to reflect photons off a mirror and induce momentum transfer, was used to verify the calibration and confirm the sign of the strain channel $h(t)$. See [47] for more explicit details of the Advanced Virgo calibration system for O2.

Calibrated strain data for the LIGO and Virgo detectors are created online for use in low-latency searches. After the completion of the observing runs, final time-dependent calibrations were generated for each detector. The results presented in GWTC-1 [3] use the full frequency-dependent calibration uncertainties described in [47, 82, 83]. It is important to note that the detector strain channel $h(t)$ is only calibrated between 10 Hz and 5 kHz for Advanced LIGO and 10 Hz and 8 kHz for Advanced Virgo [47]; the channel is not a faithful representation of strain at lower or higher frequencies.

As described in sections 3 and 5, calibrated LIGO and Virgo data can be both non-stationary and non-Gaussian at certain times and frequencies. Glitches may mimic true transient astrophysical signals in individual detectors [48], while spectral lines such as those seen in figure 2 can blind searches for long-duration signals at those specific frequencies [49]. In this section we outline how we identify and characterize these noise features so that we can either exclude the bad data or assess the impact of remaining artifacts on searches for gravitational-wave signals.

Figure 12 shows an example of a glitch. Glitches with power comparable to detectable signals have historically occurred on the order of once per minute, with larger glitches occurring less frequently. Even in their nominal state, the detectors' data contain glitches introduced by behavior of the instruments or complex interactions between the instruments and their environment. Many of these glitches (but not all) can be associated with transient signals in auxiliary channels from various sensors which serve as 'witnesses' to environmental disturbances coupling into the interferometer. These associations allow us to identify and catalog certain classes of glitches. See [48] for a detailed presentation on the characterization of transient noise in Advanced LIGO, especially pertaining to the observation of the gravitational-wave signal GW150914. Descriptions of various noise sources for Advanced Virgo in O2 can be found in [84–87].

Zoom In Zoom Out Reset image size

In searches for transient gravitational-wave signals, identified glitches and periods of poor data quality are flagged [50, 88, 89]. Periods of data are vetoed at various levels or categories depending on the severity of the problems; the GWOSC open data releases make this information available [24]. Sections of strongly non-stationary data that would corrupt the noise power spectral density estimates are removed entirely from the searches. Times when noise sources with known physical coupling to the gravitational-wave strain channel of the detector are active, and thus likely to cause glitches, are identified, and candidates at or around these times may be removed (vetoed) from search results. For a more detailed explanation of the strategy for mitigating noise sources, see section 4 of [48]. For searches for long-duration signals, frequency bands known to be dominated by instrument noise are omitted from the analysis [49]. The strategy of vetoing times of probable glitches is expected to increase confidence in detection candidates that survive the application of vetoes [50, 90, 91], and may thus increase the number of astrophysical signals that can be confidently detected. Detection methods are discussed in detail in section 8.

Although the great majority of transient noise sources are of local origin and thus uncorrelated between detectors, some noise sources exist that are potentially correlated between detectors, such as electromagnetic pulses from lightning coupling inductively into the detectors [92]. A key feature of the LIGO and Virgo experimental design is an array of physical environment monitors designed to detect environmental disturbances, and to have greater sensitivity to those disturbances than the detector's gravitational-wave strain channel does. The LIGO environmental sensor array includes seismometers, microphones, accelerometers, radio receivers and magnetometers to monitor ambient noise [93]¹⁸². Virgo has a similar array of sensors [10]¹⁸³. The environmental sensors' sensitivities are verified via a suite of noise injections performed at the beginning and end of each observing run; acoustic, magnetic, radio frequency, and vibrational tests are done to quantify the coupling from ambient noise to the gravitational-wave strain data $h(t)$ [48]. These injections are conducted at multiple locations around the detector such that sensor coupling functions to $h(t)$ via multiple potential coupling paths are verified and well understood.

As shown in figure 2 of [48], the external transient electromagnetic coupling of the ambient noise to the gravitational-wave data channel $h(t)$ is on the order of a factor of 100 below the current strain level, such that any electromagnetic source would have to register in one of the magnetometers surrounding the detector an SNR of 100 before registering in the gravitational-wave signal channel. This is easily confirmed with the study of lightning strikes during nearby storms [48]. The description of the coherence between the detectors' output strain signal $h(t)$ and magnetometers about the detector for the AC power frequencies (50/60 Hz) is nontrivial and described in [94]. In Virgo, a detailed study of the electromagnetic coupling to the gravitational-wave data channel was recently carried out [84, 85].

In addition, a potential correlated noise source in searches for a stochastic gravitational-wave background is Schumann resonances, or low-frequency magnetic field resonances between the Earth's surface and the ionosphere excited by lightning [95, 96]. These resonances are also monitored with sensitive magnetometers, with the future goal of subtracting their effect on gravitational-wave strain data [97]. The effect of Schumann resonances on the measured gravitational-wave strain is below the current Advanced LIGO–Advanced Virgo noise floor [98].

The LVC searches for gravitational waves compare the null hypothesis, ${{\mathcal H}}_0$, that a given stretch of data contains only noise, to the signal hypothesis, ${{\mathcal H}}_1$, that the stretch of data contains both noise and a gravitational-wave signal¹⁸⁴. Most searches assume that general relativity correctly describes the gravitational-wave signals. The likelihood of observing the data under the two hypotheses can be written in terms of equation (6) by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p({\bf d} \mid {{\mathcal H}}_0) = p_{0}({\bf d}) \quad {{\rm and}} \quad p({\bf d} \mid {{\mathcal H}}_1) = p_{1}({\bf d}) \nonumber \end{align} \tag{ 10 }$$

where ${{\mathcal H}}_0$ assumes noise alone with no signal in the data, while ${{\mathcal H}}_1$ assumes a signal parameterized by $\boldsymbol\theta$, ${\bf h}({\boldsymbol\theta})$ is present in addition to the noise. For the present, we will assume that each detector data stream is being analyzed independently, and we discuss below how data from multiple detectors is combined in LVC searches. The probability of the signal hypothesis given the observed data, known as the posterior probability, is given by Bayes' theorem as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p({{\mathcal H}}_1 \mid {\bf d}) = \frac{p({{\mathcal H}}_1)\,p_{1}({\bf d})} {p({{\mathcal H}}_0)\,p_{0}({\bf d}) + p({{\mathcal H}}_1)\,p_{1}({\bf d})} = \frac{p_{1}({\bf d})}{p_{0}({\bf d})} \left[\frac{p_{1}({\bf d})}{p_{0}({\bf d})} + \frac{p({{\mathcal H}}_0)}{p({{\mathcal H}}_1)} \right]^{-1} \nonumber \end{align} \tag{ 11 }$$

where $p({{\mathcal H}}_0)$ and $p({{\mathcal H}}_1)$ are our prior beliefs of whether a signal is absent or present in the data. Regardless of these prior beliefs, the posterior probability is seen to be monotonic in the likelihood ratio

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda({\bf d}|\boldsymbol{\theta}) = \frac{p({\bf d} \mid {{\mathcal H}}_1)}{p({\bf d} \mid {{\mathcal H}}_0)} = \frac{p_{1}({\bf d})}{p_{0}({\bf d})} \nonumber \end{align} \tag{ 12 }$$

and so this quantity is the optimal test statistic [54]. For Gaussian noise, the log of the likelihood ratio can be written in terms of the inner product of equation (8) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{e:lnlambda} \log \Lambda({\bf d}|\boldsymbol{\theta}) = ({\bf d} \mid {\bf h}(\boldsymbol{\theta}))-\frac{1}{2}({\bf h}(\boldsymbol{\theta}) \mid {\bf h}(\boldsymbol{\theta})) \, . \nonumber \end{align} \tag{ 13 }$$

Only the first term of this expression involves the data; it is then observed that the posterior probability is a monotonic function of $({\bf d} \mid {\bf h}(\boldsymbol{\theta}))$, a quantity known as the matched filter, which, therefore, is also an optimal test statistic.

In a matched-filter search for gravitational-wave signals, the signal parameters $\boldsymbol\theta$ will not be known in advance. The optimal detection statistic would be obtained by marginalizing [107] the likelihood ratio $\Lambda({\bf d}|\boldsymbol{\theta})$ over all unknown parameters by integrating the likelihood ratio over these parameters¹⁸⁵.

Since the log likelihood ratio is a linear function of the signal model, its exponential—the likelihood ratio itself—will generally be sharply peaked about its maximum, thus the maximum value of $\Lambda({\bf d}|\boldsymbol{\theta})$ over unknown parameters $\boldsymbol{\theta}$ is expected to be a good approximation to the marginalized likelihood ratio (up to a possible constant rescaling). This maximization procedure is equivalent to minimizing the residuals seen in the detector, which can be seen as follows. The log likelihood can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \log \Lambda({\bf d}|\boldsymbol{\theta}) = -\frac12({\bf d} - {\bf h}(\boldsymbol{\theta}) \mid ({\bf d} - {\bf h}(\boldsymbol{\theta})) +\frac12 ({\bf d} \mid {\bf d}). \nonumber \end{align} \tag{ 14 }$$

Now it is clear that the parameters $\hat{\boldsymbol\theta}$ that maximize the log likelihood ratio are those that minimize the residuals ${\bf d} - {\bf h}(\boldsymbol{\theta})$ in terms of the noise weighted inner product.

The parameters ${\boldsymbol\theta}$ describing the strain observed in a detector include the signal amplitude A observed in the detector (which is inversely proportional to the distance to a gravitational-wave source), the phase $\phi$ of the sinusoidally-varying signal observed in the detector, the arrival time t of the signal (usually defined by the moment when it reaches peak gravitational-wave amplitude at the detector), and other parameters $\boldsymbol\mu$ describing the physical parameters of the source such as the masses and spins of the components. We write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{e:hdecomp} {\bf h}({\boldsymbol\theta}) = A {\bf p}(t,{\boldsymbol\mu}) \cos\phi + A {\bf q}(t,{\boldsymbol\mu}) \sin\phi \nonumber \end{align} \tag{ 15 }$$

where ${\bf p}(t,{\boldsymbol\mu})$ and ${\bf q}(t,{\boldsymbol\mu})$ are in-phase (cosine) and quadrature-phase (sine) waveforms, normalized so that $({\bf p}\mid{\bf p})=({\bf q}\mid{\bf q})=1$, and which are orthogonal, $({\bf p}\mid{\bf q})=0$.

Maximization over the amplitude and phase can be done algebraically as follows: equation (13) can be rewritten using equation (15) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \log \Lambda({\bf d}|\boldsymbol{\theta}) = A \rho(t,{\boldsymbol\mu}) \cos(\phi - \varphi) - \frac12 A^2 \nonumber \end{align} \tag{ 16 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varphi \equiv \arctan\frac{({\bf d} \mid {\bf q}(t,{\boldsymbol\mu}))}{({\bf d} \mid {\bf p}(t,{\boldsymbol\mu}))} \nonumber \end{align} \tag{ 17 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho(t,{\boldsymbol\mu}) \equiv \sqrt{({\bf d} \mid {\bf p}(t,{\boldsymbol\mu}))^2 + ({\bf d} \mid {\bf q}(t,{\boldsymbol\mu}))^2} \nonumber \end{align} \tag{ 18 }$$

is the SNR time series for waveform templates with parameters $\boldsymbol\mu$. The log-likelihood $\log\Lambda$ is maximal for amplitude $\hat{A} = \rho$ and phase $\hat{\phi} = \varphi$ with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \max_{A,\phi} \log \Lambda({\boldsymbol\theta}) \equiv \log \Lambda(t,\hat{A},\hat{\phi},{\boldsymbol\mu}) = \frac12 \rho^2(t,{\boldsymbol\mu}) ~ . \nonumber \end{align} \tag{ 19 }$$

Peaks in this time series correspond to times at which a signal is most likely to be present. Under the signal (noise) hypothesis, and in the presence of stationary and Gaussian noise with a known power spectrum, $\rho^2(t_{{{\rm peak}}},{\boldsymbol\mu})$ follows a non-central (central) chi-squared distribution with two degrees of freedom.

The SNR time series can be conveniently expressed in terms of a complex time series as equation (8) of [66]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{complexsnr} z(t,{\boldsymbol\mu}) = 4 \int_0^\infty \frac{\tilde{d}(\,f)\tilde{p}^\ast(\,f,{\boldsymbol\mu})}{S_n(\,f)}e^{2\pi ift} {\rm d}f \nonumber \end{align} \tag{ 20 }$$

as $\rho = |z|$ and the phase $\varphi = \arg(z)$; $\tilde{p}(\,f,{\boldsymbol\mu})$ is the Fourier transform of the in-phase waveform (see equation 15). Equation (20) is the inverse Fourier transform of

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{z}(\,f,{\boldsymbol\mu}) = 4 \frac{\tilde{d}(\,f)\tilde{p}^\ast(\,f,{\boldsymbol\mu})}{S_n(\,f)} \Theta(\,f) \nonumber \end{align} \tag{ 21 }$$

where $\Theta(\,f)$ is the Heaviside step function.

Parameters such as the masses and spins of the binary components (represented by $\boldsymbol\mu$) change the morphology of the gravitational wave. In order to detect signals with a wide range of possible masses and spins, a bank consisting of large numbers of signal templates spanning the parameter space is produced, and each template in the bank is used as a matched filter. The template bank used to initially find signals in the data is constructed with a density in parameter space sufficiently high that the loss in SNR between a true signal and the best-fit template is less than 3%. See [2, 3, 110] for more details on the template banks used for the O1 and O2 searches.

An important property of the SNR should be noted: in equation (20) it can be seen that the integrand is proportional to $\tilde{d}(\,f)/S_{n}(\,f)$, so the data are not simply being whitened (which would have been the case if the denominator were S^1/2(f )), but in fact noisier parts of the frequency spectrum (including narrow lines) are suppressed in the matched filter. Equivalently, the SNR integral can be seen as correlating a whitened data time series with a whitened template. The SNR therefore provides a natural way to down-weight the frequency bands where the noise is large, and effectively notches out the various lines.

While the SNR is the optimal detection statistic in the case of stationary Gaussian noise, transient instrumental artifacts make it a non-optimal statistic with real detector noise. Although the matched filter naturally suppresses stationary noisy features in the data, glitches can cause certain templates to produce high SNR values [50, 111–113]. We address this in several different ways:

• (i)  
  As explained in section 6.2, we use witness sensors to identify times when the environment or the instrument introduces frequent glitches and we veto a subset of these times found to impact search performance from our analysis [50]. These sensors include those that monitor the physical environment about the gravitational-wave detector, as well as those that record signals from within the internal control systems of the interferometer.
• (ii)  
  We implement waveform consistency tests which characterize the deviation of the data ${\bf d}$ from the model ${\bf n}+{\bf h}$ [20, 21, 66, 114]. For signals from compact binary mergers, these tests are extremely powerful and allow us to reject many glitches which have not been identified and vetoed, though for short signals the discriminatory power of these tests is diminished [3, 20, 114, 115].The exact implementation of these signal consistency tests vary among search pipelines, but all are based on the following principle: if the gravitational-wave model waveform is subtracted from the data to produce residuals ${\bf d}-{\bf h}$, the residuals should be consistent with Gaussian noise if the signal hypothesis is true. These residuals are re-filtered with the matched filter over different time or frequency intervals to determine if non-noise-like features persist; evidence of such features suggest that the model waveform ${\bf h}$ is not a good match to the non-Gaussian feature in the data, and the detection ranking statistic is down-weighted accordingly. For example, the consistency test described in [66, 114] constructs a chi-squared test statistic by dividing the matched filter into n frequency bands as

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi^2 = \sum_{i=1}^{n} \frac{|({\bf d}-{\bf h} \mid {\bf p})_i - ({\bf d}-{\bf h} \mid {\bf p})/n|^2 + |({\bf d}-{\bf h} \mid {\bf q})_i - ({\bf d}-{\bf h} \mid {\bf q})/n|^2}{1/n} \nonumber \end{align} \tag{ 22 }$$

  where $({\bf a}\mid{\bf b})_i$ is the same as the inner product in equation (8) but with the integrand restricted to the frequency interval $f_{i-1}<f<f_{i}$ with f ₀  =  0 and $f_n=\infty$. Here the bands are chosen so that $({\bf p}\mid{\bf p})_i=({\bf q}\mid{\bf q})_i=1/n$. If the residual ${\bf d}-{\bf h}$ is Gaussian noise, $\chi^2$ is chi-squared distributed with $\nu=2n-2$ degrees of freedom; values of $\chi^2\gg\nu$ are indicative of residual non-Gaussian features in the data after the model has been subtracted. A re-weighted ranking statistic proposed in [20]

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{\rho} = \rho \times \left\{\begin{array}{@{}ll@{}} 1 & \chi^2 \leqslant \nu \nonumber \\ \left[\frac12 + \frac12(\chi^2/\nu)^3\right]^{-1/6} & \chi^2 > \nu \end{array} \right. \nonumber \end{align} \tag{ 23 }$$

  down-weights the SNR for large values of $\chi^2$. A similar time-domain based signal consistency test is described in [21] and is incorporated into a likelihood ranking statistic.
• (iii)  
  For all detections published to date we have required that gravitational-wave signals be identified via matched filtering in at least two independent detectors with consistent parameters. For example, the arrival times of the gravitational waves at each detector must differ by no more than the the maximum time-of-flight between the detectors, e.g. 10 ms for the LIGO Hanford–LIGO Livingston pair, with an extra 5 ms added in order to account for uncertainty in the inferred coalescence time at each detector. This 5 ms addition to the coincidence window is also used when searching for simultaneous events for the LIGO Hanford–Virgo pair (27 ms light travel time), and the LIGO Livingston–Virgo pair (26 ms light travel time) [3]. However, having now established the existence and frequency of gravitational-wave signals, it may now also be possible to make detections when only one detector is operating, and thus this time coincidence test is not available [116].

The matched-filter based searches employed by the LVC construct ranking statistics from the SNR and the waveform consistency test statistics [3]. In addition an astrophysical signal received in several detectors will have a common set of parameters $\boldsymbol\mu$ (within limits imposed by limited SNR) in all detectors, and, furthermore, the amplitude, phase, and time-of-arrival of the signals observed in each detector will be determined by the direction of propagation of the wave (i.e. from where on the sky the signal originates) and the polarization state of the signal. Since gravitational waves have two polarizations (in general relativity), referred to as the plus-polarization ${\bf h}_+$ and the cross-polarization ${\bf h}_\times$, the strain on detector I is determined by the detector's antenna response patterns F_+,I and $F_{\times,I}$ by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf h}_I({\boldsymbol\theta}) = F_{+,I}(\alpha,\delta,\psi,t) {\bf h}_+(t-\tau_I,D,\iota,{\boldsymbol\mu}) + F_{\times,I}(\alpha,\delta,\psi,t) {\bf h}_\times(t-\tau_I,D,\iota,{\boldsymbol\mu}) \nonumber \end{align} \tag{ 24 }$$

where $\alpha$ and $\delta$ are the right ascension and declination of the source of the gravitational waves, D is the distance to the source of the waves, $\iota$ is the inclination of the orbital plane of the binary system (which, for circular orbits and leading order quadrupole emission, determines the ellipticity angle), and $\tau_I=\tau_I(\alpha,\delta,t)$ is the travel time of the signal from the geocenter to the detector. Although a fully-coherent search for gravitational waves across a network of detectors is possible, we opt instead to perform searches independently in each detector and then demand that triggers seen in different detectors have consistent times of arrival and the same parameters ${\boldsymbol\mu}$ since this provides a powerful glitch rejection consistency test as described above. However, further signal consistency requirements are also possible. For the leading-order quadrupole emission from a circular binary, the ratios of the amplitudes seen in two detectors is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{A_I}{A_J} = \sqrt{\frac{F_{+,I}^2\left(\frac{1+\cos^2\iota}{2}\right)^2+F_{\times,I}^2\cos^2\iota}{F_{+,J}^2\left(\frac{1+\cos^2\iota}{2}\right)^2+F_{\times,J}^2\cos^2\iota}} \nonumber \end{align} \tag{ 25 }$$

while the difference in arrival time is $t_I - t_J = \tau_I - \tau_J$ and the difference in phase is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi_I - \phi_J = \arctan\left(\frac{F_{\times,J}}{F_{+,J}}\frac{2\cos\iota}{1+\cos^2\iota}\right) - \arctan\left(\frac{F_{\times,I}}{F_{+,I}}\frac{2\cos\iota}{1+\cos^2\iota}\right) . \nonumber \end{align} \tag{ 26 }$$

It is therefore possible to include an amplitude-phase-time consistency measure in likelihood-based ranking statistics [20, 50, 117]; this was done for the most recent searches for gravitational-wave signals in the LIGO–Virgo O1 and O2 data [3].

After the steps described above which mitigate the effects of noise transients, the probability of remaining transients occurring simultaneously (within a time window that takes into account the maximal travel time of a signal, for example 10 ms for the two LIGO detectors) in two detectors and producing a large joint ranking statistic value becomes extremely small. Different searches adopt different approaches for measuring this probability as a function of the ranking statistic [50]. The basic method is to examine the statistical properties of the non-simultaneous transients observed in each detector and to artificially treat them as if they did occur simultaneously.

The statistical significance of any candidate event observed in two or more detectors is quantified by its false alarm rate, which is the expected rate of events per time due to noise which would be assigned an equal or larger ranking statistic than the candidate. One approach to estimating false alarm rates is to shift one detector's data stream in time (by a time interval larger than the maximum time-of-flight between detectors) and repeat the search. The resulting 'time-shifted' coincidences are then treated as a background noise sample. This is done numerous times with different time shifts in order to obtain a probability distribution for the joint detector ranking statistics. Each coincident trigger is assigned a false alarm rate given by the number of background triggers with an equal or larger ranking statistic, divided by the total time searched for time-shifted coincidences. For example, in [1] it is found that the frequency of transients producing more significant events than GW150914 is less than once every 200 000 years in both of the matched-filter searches employed by the LVC.

Another similar approach is to accumulate single-detector triggers not having simultaneous (within the time-of-flight of gravitational waves) triggers in another detector and therefore likely not associated with gravitational-wave signals. The distribution of the ranking statistic under the background (noise only) hypothesis is then estimated by randomly drawing single detector triggers from the inferred single detector distributions and artificially treating them as if they were simultaneous when constructing the ranking statistic. The significance of an observed ranking statistic value can then be evaluated using this background distribution [21, 117]. These two independent methods of determining the significance of observed gravitational-wave candidates have both yielded high significance for the gravitational waves that have been identified.

Terrestrial noise sources that are potentially correlated between detectors are not taken into account by these background estimation methods. Thus, as discussed in section 6.2 above, a detailed examination of physical and environmental sensors which monitor such noise sources and an assessment of their coupling to the measured gravitational-wave strain channel is carried out in order to check the validity of a detection candidate.

The LVC also performs searches of LIGO and Virgo data without specific waveform models [22, 118]. These searches first identify periods of excess power in each detector's data stream and then builds a detection statistic based on the cross-correlation between the data streams. The significance of a particular detection statistic value is again assessed using time shift analyses. In [1] it was shown that the frequency of noise transients producing more significant events than GW150914 in such a generic transient search is once every 8400 years.

The fact that multiple searches, employing different methods, all found GW150914 to be a highly significant candidate bolsters our confidence that this event is not the product of coincident transient noise. Furthermore, the signal is well matched by the waveform predicted by general relativity for the coalescence of a binary black hole system. Various tests performed using the first ten binary black hole mergers detected by the LVC with high confidence have shown no significant deviation from general relativity models [119].

A final component of a search pipeline is the determination of the sensitivity of the search to astrophysical populations of signals. The purpose of doing so is two-fold: first, it provides a metric by which a search can be tuned to optimize its detection efficiency for particular classes of signals; second, it provides a means for interpreting the rate of signals detected by the pipeline to the rate at which signals are generated by the population.

The sensitivity of the search pipeline is normally measured via a Monte Carlo procedure (see, e.g. [120]) in which simulated signals drawn from a hypothetical source population (e.g. some distribution of binary component masses and spins, orientation angles, arrival times, and distance) are added to real detector noise, and the search is rerun to determine the fraction of signals from this population that are detected by the search pipeline. A simulated signal is considered to be detected if the search pipeline produces a trigger above some chosen ranking statistic threshold. The result is represented as a time- and population-averaged spacetime sensitivity $\langle VT\rangle$ for a fixed ranking statistic threshold which corresponds to a fixed false alarm rate threshold (e.g. a threshold could be chosen to be one false alarm per century of observation). Alternatively, threshold-independent methods of astrophysical rate estimates can also be employed [121–123].

Let us focus now on the choice of signal model M. The signal model M determines the functional form of $h(t;\boldsymbol{\theta})$ which is key to calculating the likelihood function. For definiteness, we will concentrate on parametric forms of $h(t;\boldsymbol{\theta})$, as obtained by solving Einstein's equations. For a discussion of non-parametric signal models see [128]. Exact analytic solutions of Einstein's equations are notoriously difficult to obtain; therefore data-analysis-ready models are either based on perturbative solutions, e.g. the Taylor family of waveforms [129] or the effective-one-body waveforms [130–134], or on hybrid/phenomenological approaches such as the Phenom family of waveforms [135–138]. We will not discuss further the details of the waveform models here, but we will restrict ourselves to the two main types of waveforms employed in the original analysis of GW150914. For GW150914, the two models used in [127] were SEOBNRv2 and IMRPhenomPv2 [127, 139]. Both waveform models are full inspiral-merger-ringdown models that succeed in reproducing numerical waveforms, especially in the region of approximately equal mass and moderate spins magnitudes. The main difference between the two models lays in the treatment of the spin dynamics. SEOBNRv2 models the dynamics of the component of the spin vectors along the direction of the orbital angular momentum while IMRPhenomPv2 includes also an effective treatment of the dynamics of the in-plane components of the spins, and thus includes an approximate precessing dynamics¹⁸⁶. Because GW150914 was nearly a face-off system (orbital angular momentum vector pointing away from the Earth), the LIGO instruments were not sensitive to the in-plane spin components, hence the two waveform models are essentially equivalent. In [143], the LVC has empirically shown that the inferred properties of GW150914 depend relatively weakly on a change in the waveform model. This finding was confirmed by the analysis presented in [134] using an independent effective-one-body implementation and by [144], in which numerical relativity solutions were directly compared with GW150914 data.

The final functions necessary for the application of Bayes' theorem, equation (27), are the prior probability distributions for the parameters of interest. These are all the parameters necessary to completely characterise the gravitational-wave signal emitted during a coalescence event. For quasi-circular orbits, these are:

• the component masses m₁ and m₂;
• the spin vectors $\vec {S}_1$ and $\vec {S}_2$;
• the polarisation angle $\psi$ and the angle $\theta_{jn}$ between the total angular momentum $\vec {J}$ and the propagation direction of the gravitational wave $\hat{n}$;
• the source luminosity distance D_L;
• the source right ascension $\alpha$ and declination $\delta$;
• a reference phase $\varphi_0$ and a reference time, typically the gravitational-wave strain peak time, t₀.

The functional form of the prior distribution must be specified for all parameters. In some cases the prior distribution is determined via invariance (symmetry) properties of the parameter space [145]; for instance, the prior for the source position $D_L, \alpha, \delta$ in the Universe is chosen from the requirement that the number density of sources is uniform in the cosmological co-moving volume in accordance with a Friedmann–Lemaître–Robertson–Walker cosmological model. Thus, the probability $p(D_L, \alpha, \delta|M\, I) \propto dV$ and for redshift $z \ll 1$ reduces to $p(D_L, \alpha, \delta|M\, I) \propto D_L^2 \vert \cos(\delta)\vert$. In cases where invariance arguments do not apply, we choose simple forms of prior distribution so that the resulting posteriors are easily interpretable. Similar arguments determine the prior for the spin vectors $\vec {S}_1, \vec {S}_2$ and orientation angles $\psi,\theta_{jn}$ to be uniform over the azimuthal angles ranging between 0 and $2\pi$ as well as uniform in the cosine of the polar angles ranging between  −1 and 1. Regarding the spin vector magnitudes, several possible priors are possible, e.g. $p(|\vec {S}_i||M,I) \propto |\vec {S}_i|^2$ or $p(|\vec {S}_i||M,I) \propto 1$. The main analysis of the events in the GWTC-1 catalog employed a uniform distribution over the norm of the spin vectors. For the component masses m₁ and m₂, the chosen prior distribution is uniform, thus $p(m_1, m_2 |M\, I) \propto 1$, but limited from below so that $m_1,m_2 > 1M_{\odot}$.

In addition to uncertainties induced by detector noise, the accuracy and precision of our source parameter estimates are also affected by uncertainties in the amplitude and phase response of the detectors. For this reason, this source of uncertainty on the data is modelled and included in the analysis. The calibration uncertainty model employed is based on empirical estimates of the error magnitudes in both amplitude and phase in specified frequency bands [47, 74, 82]. In particular, the model assumes the value of the errors to be distributed as a Gaussian distribution with zero mean and a variance given by the empirically determined error magnitudes. Calibration uncertainty curves are then constructed using third order spline interpolation over the data frequency space. Typically, this introduces a total of $O(10)$ additional parameters per detector (half for the phase uncertainty and the other half for the amplitude), which are sampled in concert with the physical parameters of the system. Technical details for the LVC calibration model can be found in [146].

The total number of parameters to be inferred is thus 15 for quasi-circular orbits and generic spin vectors, and 11 for models where spins are forced to be aligned with the orbital angular momentum. To the set of physical parameters, we must add the 10 parameters per detector necessary to specify the calibration uncertainty model. Hence, for a typical three-detector analysis, we sample a grand total of 45 parameters for a quasi-circular system with generic spin orientations. Parameter spaces of such high dimensionality cannot be efficiently explored with grid-based methods. Therefore, over many years members of the LVC developed the LALInference stochastic sampler library [105] which implements two algorithms, a parallel tempering Markov chain Monte Carlo [147] and a nested sampling [102]. The parallel tempering Markov chain Monte Carlo is designed to generate samples from the multidimensional posterior distribution (27), while the nested sampling instead is designed to calculate the evidence, equation (28) and generates samples from the posterior distribution as a by-product. More details are given in [105] and references therein. Other parameter estimation pipelines are routinely used by the LVC, such as rapidPE [148] and BILBY [149], but for the rest of the discussion we will focus on LALinference. However, the same considerations will apply to other Bayesian analysis methods.

The end products of the LALInference analyses are posterior samples for all parameters that characterise the gravitational-wave waveform. Of particular interest are posteriors for the intrinsic parameters, masses and spin vectors, which help ascertain the nature of the coalescing objects. For GW150914, the detector-frame masses (i.e. redshifted due to cosmological expansion) measured using the IMRPhenomPv2 waveform model were $38.5_{-3.6}^{+5.6} M_\odot$ and $32.2_{-4.8}^{+3.6} M_\odot$; see table I in [127]. The quoted numbers are an extremely concise way of summarising the full posterior distribution. The output of LALInference analyses are samples from the full posterior distribution. In particular, a number such as $38.5_{-3.6}^{+5.6} M_\odot$ comes from the marginalisation over 14 of the 15 physical source parameters of the full posterior, as well as the calibration parameters, to obtain a one-dimensional posterior from which the 90% credible region is then calculated. Naturally, correlations between different parameters are invisible in a one-dimensional representation. For a clearer picture, multidimensional posterior distributions help to display the information extracted from the analysis. Figure 13 shows the joint two-dimensional posterior distribution for the component masses m₁ and m₂ as an example. In particular, the bottom left panel shows the non-negligible correlation between the component masses.

Zoom In Zoom Out Reset image size

The full details of the multi-dimensional posterior distribution can be visualised in a compact way by computing the posterior distribution over the waveform itself in the time domain. This is done simply by computing the predicted waveform over each of the posterior samples. Let $\boldsymbol{\theta}_i$ be the ith posterior sample, the corresponding waveform will be $h(t;\boldsymbol{\theta}_i)$. The waveform samples are the set $\newcommand{\e}{{\rm e}} \{h(t;\boldsymbol{\theta}_i)\}_{i=1,\dots,N}\equiv \{h_i\}$. Each of the waveform samples h_i can be whitened, see section 3, and then used to compute credible intervals at every time t_j at which the original data were sampled. The result of this procedure is summarised in the presentation of figure 6 of [127]. Figure 1 in [1] is representing a different procedure; the second row in this figure shows a comparison between the reconstructed 90% credible region obtained by the procedure described above, and a numerical relativity solution that, while not corresponding to any of the computed posterior samples, is consistent with the reconstructed 90% credible region.

The results from Bayesian inference are only as good as the models used in the analysis. If the waveforms used in the signal model or the underlying assumptions of the noise model are inaccurate, the results will suffer from systematic bias. A multitude of tests are used to check for possible mis-modeling error and to quantify the impact on the analyses. As discussed earlier in section 9.1 the waveform models are compared to highly accurate numerical relativity simulations, and multiple waveform approximants are used in the analyses and cross-compared. The difference between the results found using different waveform models provides an estimate of the systematic error due to the signal model. The noise model can also be checked. Other checks include adding simulated signals with similar parameters to the astrophysical events into nearby stretches of data and checking that the parameters are properly recovered by the parameter estimation algorithms.

Over long stretches of LIGO–Virgo data, the noise is known to be non-stationary and non-Gaussian. The overall noise levels fluctuate, and there are frequent low-SNR glitches, and less frequently high-SNR glitches, see section 5. On the other hand, the gravitational-wave signals spend very little time in the LIGO–Virgo sensitive band—seconds or less for black hole binaries and minutes for neutron star binaries, and over these shorter stretches of time the noise is generally (but not always) well approximated as stationary and Gaussian. When a significant trigger is found by the search pipelines, the first thing the analysts look at are multi-resolution time-frequency scalograms of the data surrounding the trigger (known as Q-scans). Q-scans are qualitative checks which require visual inspection [151, 152]. These scans reveal whether there are any loud glitches in the data, as was the case with the binary neutron star GW170817 [8]. Once the parameter estimation analyses have been run, Q-scans of the residuals are closely examined to see if any any unmodelled noise features might have affected the analyses. Figure 14 shows Q-scans of the whitened data and residuals surrounding GPS time 1126259462. The scans of the data reveal the signal from GW150914, while the residuals after subtracting the maximum likelihood waveform from the parameter estimation studies [127] show no visible evidence of glitches or correlated signal power.

Zoom In Zoom Out Reset image size

In addition to these qualitative checks, more rigorous quantitative checks can be applied. One test that is routinely applied is to reanalyze the residuals using the wavelet-based BayesWave algorithm [128] which is able to identify any glitches and remaining coherent power. Coherent power in the residuals could be evidence of departures from general relativity, or evidence of shortcomings in the template models or the noise model used for parameter estimation. No significant coherent power was found in the residuals for any of the detected events. In the case of GW150914 the lack of a coherent residual was used to place interesting bounds on possible departures from general relativity [153]. In the case of the binary neutron star merger GW170817 [8], a loud incoherent glitch was seen to overlap the signal in the Livingston detector. The glitch was reconstructed and removed using the BayesWave algorithm. The glitch removal procedure has been shown to be safe in a study that injected simulated neutron star merger signals into data with similar loud glitches, followed by removing the glitches with BayesWave and accurately recovering the true signal parameters with the LVC parameter estimation algorithms [154].

Figure 15 shows histograms of the whitened Fourier amplitudes of the residuals in the LIGO-Hanford and LIGO-Livingston detectors following the removal of the maximum likelihood waveform for GW150914. The residuals are taken from the parameter estimation analysis published at the time of the discovery [127]. These residuals were used to test for residual coherent power, which can be framed as a test of general relativity if we assume the template to be a sufficiently accurate solution of the theory. Applying the Anderson–Darling test of normality to the residuals yields p-values of 0.15 for LIGO-Hanford and 0.11 for LIGO-Livingston, indicating that the residuals are consistent with the Gaussian noise model used to define the likelihood.

Zoom In Zoom Out Reset image size

Even if the residuals had failed the formal tests of stationarity and Gaussianity discussed here, it would not necessarily imply that the parameter estimation would be strongly biased. When the noise deviates from the model the analysis will suffer systematic bias. But for this bias to be significant it has to be large compared to the statistical spread in the posterior distributions. Extensive studies using simulated signals added to real LIGO–Virgo data have shown that systematic errors due to deviations from the noise models are generally negligible compared to the statistical uncertainties [104, 143, 154–157]. One exception is when the simulated signals cover or overlap the times of glitches, in which case the biases can be large [158]. When glitches are present, tools such as BayesWave need to be used to model and remove the glitches, ideally in concert with the parameter estimation.

Gravitational-wave templates exhibit a variety of parameter degeneracies whereby templates with different parameters can have very similar amplitude and phase evolution, and yield very similar likelihoods. One example of such a degeneracy is evident in the posterior distribution for the component masses of GW150914 shown in figure 13. The degree of similarity between templates with parameters $\boldsymbol{\lambda}, \boldsymbol{\theta}$ is measured by the match

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm M}(\boldsymbol{\lambda}, \boldsymbol{\theta}) = \frac{({\bf h}(\boldsymbol{\lambda}) | {\bf h}(\boldsymbol{\theta}))}{\sqrt{({\bf h}(\boldsymbol{\lambda}) | {\bf h}(\boldsymbol{\lambda})) ({\bf h}(\boldsymbol{\theta}) | {\bf h}(\boldsymbol{\theta}))}}\, . \nonumber \end{align} \tag{ 29 }$$

If the true signal is described by ${\bf h}(\boldsymbol{\lambda})$, then the expectation value of the log likelihood for template ${\bf h}(\boldsymbol{\theta})$, maximized over amplitude is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm E}[\ln\Lambda (\boldsymbol{\lambda}|\boldsymbol{\theta})] = \frac{1}{2} {\rm M}^2(\boldsymbol{\lambda}, \boldsymbol{\theta})\, {\rm SNR}^2 \, , \nonumber \end{align} \tag{ 30 }$$

where $\rm SNR$ is the optimal signal-to-noise ratio [159]. We see that signals that have similar morphology, as measured by the match, yield similar likelihoods. Now suppose we hold one parameter, $\theta^k$ fixed, then maximize the likelihood with respect to all the other parameters. Up to an overall constant, we have [159, 160]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{logLFF} \ln\Lambda ({\bf d}| \bar\theta^k) \equiv {\rm max}_{j\neq k}\ln \Lambda({\bf d}| \theta^{\,j}) \simeq \frac{{\rm SNR}^2}{2}\, {\rm FF}^2(\bar\theta^k)\, , \nonumber \end{align} \tag{ 31 }$$

where ${\rm FF}(\bar\theta^k)$ is the fitting factor, or maximized match, between waveforms with $\theta^k=\bar\theta^k$ and the maximum likelihood waveform. Figure 16 compares the maximum likelihood as a function of the primary detector-frame mass for GW150914 to the fitting factor. The fitting factor as a function of m₁ was computed by maximizing the match between the overall maximum likelihood waveform and waveforms with fixed m₁. Note that the posterior distribution for this event had a 90% credible interval of $m_1 = 38.5_{-3.6}^{+5.6} M_\odot$, but templates with primary masses outside this interval continue to yield large fitting factors because other parameters can be adjusted to partly compensate the effects of the change in primary mass on the waveform. For example, we find a fitting factor between the maximum likelihood template for GW150914 and a template with a primary mass of $m_1=70 \, M_\odot$ of ${\rm FF}= 0.95$.

Zoom In Zoom Out Reset image size

The possibility of achieving high matches, or correlations, between templates with large primary masses and the maximum likelihood template have been cited as evidence that the LIGO–Virgo parameter estimation analysis for GW150914 and other systems may be flawed [43]. It was further hypothesized that the credible intervals were underestimated due to the instrument noise not conforming to the likelihood model [43]. However, our confidence in the reconstructed credible regions comes from extensive simulations designed to compare the cumulative distributions of simulated populations against the cumulative distribution of reconstructed credible regions. The agreement between the two, see for instance figure 10 in [105], demonstrates that our algorithms are properly computing the credible intervals. Moreover, as we have shown, the noise properties for GW150914 are compatible with the likelihood model used in the parameter estimation studies, further reinforcing our confidence in the method used to compute the credible intervals.

There is no contradiction in having templates with parameters outside the quoted credible regions producing large fitting factors with the best fit model, since even small template mismatches come with a large penalty for high signal-to-noise ratio systems such as GW150914. For example, the difference in log likelihood between the signal with $m_1=70 \, M_\odot$ and the global maximum is $\Delta \ln \Lambda = -32$, which is what we expect to see for a ${\rm SNR} \simeq 25$ signal and a template with a fitting factor of ${\rm FF}= 0.95$. But the relative likelihood for the higher mass solution is $e^{\Delta \ln \Lambda}=10^{-13.9}$, thus while templates with large primary masses can produce relatively good matches to the data, the probability that the primary mass is this high is vanishingly small.

As introduced above, the physical parameters of the signal $\boldsymbol{\theta}$ determine the shape and amplitude of the gravitational-wave signal $h(t;\boldsymbol{\theta})$. Numerical relativity simulations can be used to generate reference templates [24] using intrinsic parameters taken from the Bayesian parameter estimation studies. However, the templates still need to be projected onto the detectors using an appropriate set of extrinsic parameters. In a single detector the projection is equivalent to time shifting, phase shifting and rescaling the reference template: $\tilde h(\alpha, \delta t, \delta \phi)(\,f) = \alpha \, \tilde h_{\rm ref}(\,f) e^{2 \pi i f \delta t + i\delta \phi}$.

Figure 17 shows the reference numerical relativity template from figure 2 of the GW150914 discovery paper [1] along with maximum likelihood projections onto each detector. A smooth taper has been applied to the start of the template to avoid spectral leakage when transforming to the Fourier domain. The data file for the template was taken from the original posting at the GWOSC [162] and originates from the simulation SXS:BBH:0305, calculated for a system with a mass ratio of q  =  0.819, spins aligned with the orbital angular momentum with dimensionless magnitudes $\chi_1=0.330$ and $\chi_2=-0.440$, and detector-frame total mass scaled to $M=74.6 M_\odot$. These waveform parameters are consistent with those eventually determined for GW150914, within uncertainties, but do not exactly maximize the likelihood globally. Using the maximization procedure described in section 8 one finds that the signal arrived at the LIGO-Livingston detector 7.08 ms before the LIGO-Hanford detector, had a larger amplitude projected onto the antenna response pattern in LIGO-Hanford by a factor of 1.24, and had a phase difference of  −2.9 radians. These are, however, based on finding maximum-likelihood matches to the detector data individually with a fixed waveform without constraining them to be consistent (for example, the relative time shift could in principle be greater than the maximum light travel time between the detectors), a simplified procedure compared to the simultaneous multi-detector likelihood maximization described in section 9. When a loud signal is present in the data the individual and joint maximization techniques yield consistent results.

Zoom In Zoom Out Reset image size

In figure 1 of the GW150914 discovery paper [1] the LIGO-Hanford data were inverted (corresponding to a phase shift of $\pm$$\pi$) and overlaid on the LIGO-Livingston data to illustrate the similarity of the signals in the two detectors with minimal processing of the raw data. In addition, the reference numerical relativity template described above was approximately matched to the LIGO-Hanford and LIGO-Livingston data by adjusting the relative phase, amplitude and time offset. These adjusted templates for each detector were passed through the same bandpass and notch (band-reject) filters as the data and were then subtracted to produce the residuals plotted in the third row of figure 1 in that paper. Because those 'Fig 1 PRL' residuals were not globally optimized and were calculated from filtered data, they produce a somewhat different result than minimizing the residuals in the whitened and bandpassed data, as we will see below.

Figure 18 compares the whitened data to the whitened numerical relativity templates, maximized over arrival time, amplitude and phase, in the LIGO-Hanford and LIGO-Livingston detectors. Also shown are the residuals produced by subtracting the templates from the data. Prior to the publication of the GW150914 discovery paper [1], multiple tests were applied to the residuals to verify they were consistent with noise. The whitened residuals in each detector were found to be consistent with a Gaussian distribution: the Fourier amplitudes pass the Anderson–Darling test (see figure 15 in section 9), and the Fourier phases were found to be randomly distributed. The residuals from Bayesian parameter estimation studies [127] were analyzed using a wavelet reconstruction algorithm [128] that is able to detect coherent signals of general morphology. The degree of coherence in the GW150914 residuals was found to be entirely consistent with noise [153].

Zoom In Zoom Out Reset image size

A simpler, though less sensitive, test for coherence is to cross-correlate the data in the two detectors. The cross-correlation can be computed either in the time domain or the frequency domain using the whitened residuals. The correlation in the time domain is defined as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C(\tau) =\frac{\int H(t-\tau) L(t) {\rm d}t}{\sqrt{\int H^2(t) {\rm d}t \int L^2(t) {\rm d}t}} , \nonumber \end{align} \tag{ 32 }$$

where $H(t)$ and $L(t)$ represent the data streams from LIGO-Hanford and LIGO-Livingston respectively. When working with a finite data segment of duration T the data may be taken to be periodic: $H(t) = H(t+T)$. The correlation measure is very sensitive to the positioning and duration of the time window and the bandpass filtering that is applied to the data. To make meaningful statements about the significance of the correlation we need to know the distribution of the correlation measure for uncorrelated white noise, and these distributions change depending on the duration and bandpass. When applied to uncorrelated, unit variance Gaussian noise, the correlation coefficients follow a zero mean Gaussian distribution with a variance that depends on the duration and bandpass. Following [42] we apply the correlation analysis to four different time windows. The standard deviations for white Gaussian noise are $\sigma=0.0870$ for the 0.2 s segment, $\sigma = 0.121$ for either 0.1 s segment and $\sigma=0.193$ for the 40 ms segment.

Figure 19 shows the correlations using the whitened data shown in figure 18 (bottom panel), and in addition, the bandpass/notch-filtered data used to produce the panels in figure 1 of the GW150914 discovery paper [1] (top panel). There is a clear anti-correlation peak in the LIGO-Hanford—LIGO-Livingston data at a time lag of  ∼7.3 ms, which is consistent with the time delay inferred for the gravitational-wave signal.

Zoom In Zoom Out Reset image size

In contrast, figure 20 shows the correlations in the residuals produced using the procedure described above. The residuals from figure 18 show no notable anti-correlation at  ∼7 ms (bottom panel), while those from figure 1 of the GW150914 discovery paper [1] have a slight dip at this time lag (top panel), reflecting the fact that the reference waveform used for illustration in that paper was not the maximum likelihood waveform. For the shortest integration interval, the residuals from figure 1 of the GW150914 discovery paper have a  ∼3 sigma anti-correlation at a time lag of  ∼7.45 ms, which while marginally consistent with noise, is evidence that the signal subtraction was imperfect. In contrast, the residuals produced using the amplitude/time/phase maximized NR waveforms and whitened data show no significant excursions, and are fully consistent with noise. This is also the case for the residuals from the Bayesian parameter estimation described in section 9. Independent analyses of the GW150914 data have also found no significant correlations between the residuals in the Hanford and Livingston detectors [32, 33].

Zoom In Zoom Out Reset image size