Dynamic Normalization for Compact Binary Coalescence Searches in Non-Stationary Noise

The output of gravitational-wave interferometers, such as LIGO and Virgo, can be highly non-stationary. Broadband detector noise can affect the detector sensitivity on the order of tens of seconds. Gravitational-wave transient searches, such as those for colliding black holes, estimate this noise in order to identify gravitational-wave events. During times of non-stationarity we see a higher rate of false events being reported. To accurately separate signal from noise, it is imperative to incorporate the changing detector state into gravitational-wave searches. We develop a new statistic which estimates the variation of the interferometric detector noise. We use this statistic to re-rank candidate events identified during LIGO-Virgo's second observing run by the PyCBC search pipeline. This results in a 7% improvement in the sensitivity volume for low mass binaries, particularly binary neutron stars mergers.


Introduction
The Laser Interferometer Gravitational-wave Observatory (LIGO) [1] and Virgo [2] have so far observed 12 distinct gravitational-wave signals, two from the merger of binary neutron star systems and ten from the merger of binary black hole systems [3,4]. These three detectors are currently operating at 60-75% of their design sensitivity in the third observing run (O3) [5]. So far ∼50 candidate gravitational-wave events have been identified [6] in this run, which will conclude in April 2020. After upgrades these detectors, in addition to KAGRA [7,8], will return online in late 2021 at design sensitivity.
The sensitivity of these interferometers to astrophysical signals is deduced by estimating the power spectral density (PSD) of the noise contributions to the measured strain [9]. These noise contributions are assumed to be a stationary Gaussian process. In reality however, the output from these instruments is highly non-stationary; the detector sensitivity can change on the order of seconds due to broadband sources of noise, whether instrumental or environmental in origin.
Searches for coalescing binaries using matched-filtering have contributed to the detection of all transient gravitational-wave signals [3,10,11]. Matched-filtering correlates the detector data with a set of possible compact binary coalescence (CBC) waveforms and reweights this correlation with the detector's estimated PSD. This produces a signal-to-noise ratio (SNR) timeseries for each waveform. In the presence of stationary and Gaussian noise the matched filter is the optimal detection strategy [12]. Typical searches estimate a detector's PSD over several minutes [13,14,15]. However, the detector noise can be highly variable over much shorter periods of time. This means the matched filter will not accurately capture the variable nature of interferometric noise, leading to a reduction in search sensitivity.
In this paper we investigate how incorrectly estimating a detector's PSD affects the detectability of gravitational-wave signals from merging neutron stars and black holes. This is similar to the work presented in Ref [16]. We develop techniques to account for the changing detector noise and incorporate these methods in to the PyCBC library [12,13,17]. We investigate the impact on the PyCBC search background on data from LIGO-Virgo's second observing run (O2) [18]. We find the sensitivity of PyCBC to binary neutron stars is improved by using techniques to dynamically normalize PyCBC's ranking statistic. However we find no significant improvement to the search's sensitivity for higher mass systems. This paper is organized as follows: Section 2 provides an overview of the matched filter in the presence of non-stationary detector noise. Section 3 details the methods used to track non-stationarity in real data. In Section 4 we identify and discuss noise variations in LIGO-Hanford during O2. Section 5 describes how the PyCBC ranking statistic is dynamically re-normalized to account for PSD variation. Finally, we demonstrate the improvements to the PyCBC search background and therefore sensitivity in O2 as a result of these methods.

Effects of non stationarity in gravitational-wave signal extraction
The mis-estimation of the noise power spectral density due to non-stationarity affects the search for gravitational-wave signals. The output of a gravitational-wave interferometer is a time series s(t) such that: if a signal is present, where n(t) is the detector noise and h(t) is a gravitational-wave signal. The noise is assumed to be well described as stationary colored Gaussian noise with zero mean and one-sided PSD, S A (|f |), defined by where the angle brackets denote averaging over different ensembles of the noise. We use the variables t and f to indicate whether a quantity is in the time or frequency domain. Gravitational-wave transient searches estimate the detector noise to construct a SNR time series of the data and identify signals. This is done through the matched filter, which consists of the inner product of the data with the template h: which is appropriate for real-valued data, for which s * (f ) = s(−f ). Note that the estimated power spectrum of the noise S E (|f |) is an average of the noise spectrum over a certain time window. To obtain the SNR time series of the data we normalize the matched filter. If the data do not contain any signals s(t) = n(t), the mean value of the matched filter over different noise realizations cancels, i.e. (n|h) = 0 . In this case, the variance of the matched filter is given by where the second equality follows from the definition of S A , which is the actual spectrum of the noise. If S E ≡ S A , (5) reduces to (h|h), which motivates the definition of the SNR as It follows that if nothing but noise is present in the data the variance of the SNR is Any non-stationarity in the data will lead to a divergence between the estimated spectrum (S E ) and the actual one (S A ). An estimate of the SNR variance can therefore be used to develop a new statistic which measures the variation of the noise PSD. We will refer to this as the PSD variation statistic.
Dynamic Normalization for Compact Binary Coalescence Searches in Non-Stationary Noise4

Optimality of the matched filter
The mis-estimation of the noise spectrum affects both the inner product of the signal with the template and its variance in noise. Both must be taken into account. For instance, we could use a very narrow bandwidth to define the inner product. This would lower the variance of the inner product in noise, but the matched filter would become less effective overall due to the decrease in the inner product of the signal with the template. The optimal configuration is obtained when the ratio between the SNR squared and its variance in noise is maximized. To prove that this is consistent with and By considering a signal with unknown amplitude s = Ah and ignoring the noise contribution, the SNR squared is Using (9) and (10) we can rewrite the variance of the SNR squared in noise as Treating w(f )df as a measure dw, the SNR normalized by its variance becomes where the Cauchy-Schwartz inequality is used to set an upper limit to the ratio. The limit is saturated when u = 1, i.e. when the estimated spectrum is equal to the actual one. Accordingly, if data are non-stationary the inner product is not an optimal receiver and therefore the matched filter SNR does not provide the optimal detection statistic.

PSD variation of gravitational-wave data
We estimate the PSD variation with a similar approach as Ref. [16]. We begin applying a filter to the data, where N is a normalization constant. Here we use an approximate expression for CBC templates which captures the dominant amplitude behavior, giving [19]. We compute S E (f ) every 512 seconds using the Welch method on overlapping segments of 8 seconds. This is identical to how the PyCBC search estimates noise. To remove undesired frequencies, we combine F(f ) with a bandpass filter. The low frequency threshold (20 Hz) is determined by the LIGO detector sensitivity, while frequencies above 480 Hz are filtered out to remove strong spectral lines [1] and to reduce the computational cost. The combination of filters is then smoothed and convolved with the data. The spectrum of the resulting time series, integrated over frequency, gives an estimation of the variance of the SNR. Using Parseval's theorem, we can then compute the PSD variation, v s at a given time ∆t is chosen to match the typical time scale of non-stationarity and F(t) * s(t) is the convolution between the filter and the data. The resulting time series corresponds to the estimated SNR variance of the previous ∆t seconds, sampled over 1 second. Our estimate can be strongly affected by loud instrumental noise transients, known as glitches [20] [21], which corrupt the data over short time scales. To reduce this contribution we compute (15) over a short time (0.25 seconds) and filter out transients with a moving average. Using a set of simulated merger signals we verified that real gravitational-wave events do not affect the PSD variation statistic. We simulated a population of compact binary mergers [22,23,24] uniform in chirp mass and distance which spans the astrophysical signal space observable by LIGO. We added these signals to 4 hours of simulated stationary Gaussian noise and computed the PSD variation. The PSD variation distribution does not change due to the presence of the signals; none of the mergers are distinguishable from the noise.

Uncertainty in the estimation
Tracking variations of the PSD over short time scales increases the uncertainty of the estimation of the SNR variance. If the data are stationary Gaussian noise we can gather the statistical properties of our estimator analytically. In this case, it is convenient to factor the filter F(f ) into a normalization, a weighting and a whitening factor. The filter applied to the data in the frequency domain is therefore where the final factor is white Gaussian noise. Each frequency bin is chi-distributed with two degrees of freedom and our statistic is a weighted sum of the square of these.
The moment generating function of the distribution is S E (f ) and N is the number of frequency bins. Accordingly, the variance of the PSD variation distribution relative to its mean is where the lower limit is defined by the Cauchy-Schwartz inequality and is reached when the weights w j are all equal. ∆f ef f is the number of weights which contributes to the summation per second and it represents the frequency range where the detector is effectively sensitive to gravitational-wave signals. For the current LIGO detectors we expect the dominant weights to be in a frequency range of approximately 100 Hz, due to the detector sensitivity. To estimate the variance we fix the integration time, ∆t, to 8 seconds. This time window is consistent with the typical time-scale of the noise variations and it is sufficient to encompass the main power of typical CBC waveform templates. Assuming all the weights are equal in the frequency range we obtain a variance of the PSD variation of 1.25 × 10 −3 . To verify our predictions, we computed the PSD variation over 4 hours of simulated LIGO-Hanford detector noise. We measure Var(v s ) ∼ 1.7×10 −3 . The bandwidth of the Hanford detector is therefore approximately 70 Hz. The SNR variance is not suitable to track noise variations in lower bandwidth detectors, such as Virgo. Indeed, its narrow effective bandwidth requires an integration time longer than the time scale of PSD variations to obtain a reasonable uncertainty. Accounting for noise variations will be increasingly important in the future, when the ground based gravitational-wave detectors reach their design sensitivity. Moreover, the next generation of ground based detectors, such as the Einstein Telescope [25], will have a higher bandwidth which allows tracking of non-stationarity over shorter time scales.

PSD Variation Statistic over example O2 Data
We use the SNR variance to track times of non-stationarity of LIGO data collected during the second observing run. We consider LIGO-Hanford data between January 22th 2017 08:00:00 UTC and February 3rd 2017 16:20:00 UTC, corresponding to ∼ 10 days of observing data. During this period the noise spectrum of the detector was extremely unstable due to environmental disturbances, such as high micro-seismic noise due to bad weather conditions. Consequently, we see wide variations in the detector sensitivity. Figure 1 shows the distribution of the PSD variation during this period compared to the theoretical prediction from Gaussian noise. The Hanford data clearly shows non-stationarity which appears as a tail of high PSD variation values. In Figure 1 we have restricted the x-axis for visualization purposes. In fact, the tail of the distribution extends to PSD variation values of order of 10 3 . These instances are due to a class of extended clusters of transient noise which are relatively rare and are not the target of this method. Data which are not consistent with Gaussian noise represent ∼ 3% of the whole period.
The right plot of Figure 1 shows a time-frequency representation of example data picked from the tail of the distribution. The power excesses, shown in yellow, are clearly not associated to any gravitational-wave signal and their origin is not trivial. Within the tail are other forms of noise which can look very visually different to Figure 1. Our method has the advantage of identifying many different forms of noise that can impact a matched-filter search. We therefore developed a PSD variation monitor to assess the stationarity of the LIGO detector data in low latency to inform possible retractions of false LIGO-Virgo events. This low latency monitor has been running for the entirety of O3.

SNR normalisation
Gravitational-wave searches for coalescing binaries based on matched filtering identify signals by correlating the detector data against a set of CBC template waveforms. The searches begin by dividing the data into blocks of several minutes and compute the average noise spectrum of the detectors for each period. These are used to create a matched-filter SNR time series for each template for each detector. Every local maximum in the SNR time series which exceed a fixed threshold defines a single detector trigger. These triggers can be generated either by gravitational waves or noise artefacts. The PyCBC search pipeline down-ranks loud noise transients with a chi-squared test which checks that the accumulation of signal power as a function of frequency is consistent with the matching template waveform. The re-weighted SNR,ρ, a function of ρ and χ 2 given by (1) of Ref. [26], defines the detection statistic for a single detector that ranks the likelihood for a trigger to be due to a real signal versus noise.
Short term fluctuations of the noise spectrum affect the optimality of the match filter and can bias the distribution of triggers, even in otherwise Gaussian noise. However, we can not account for these variations directly during the match filtering because estimating the PSD with sufficient accuracy to maintain search sensitivity requires durations of data of order a few hundred seconds [12]. Moreover, in order to construct the matched filter the same PSD must be used over the duration of the signal. Calculating the PSD over short time periods (i.e. order of seconds) is therefore problematic for long signals like binary neutron stars mergers (i.e. on the order of minutes). A better approach is to re-rank the SNR time series with the PSD variation statistic during times of non-stationarity.
We implemented this approach in the PyCBC search to account for noise variations of O(10s) timescales. In Ref. [10] a correlation was demonstrated between the PSD variation and the rate of noise triggers above a given threshold inρ. The rate was found to be a function ofρv S (t) −κ , where κ is a constant allowing for deviation from the expected behavior ρv −1/2 S . Estimates of κ 0.33 were obtained, thus a re-scaled statisticρv S (t) −0.33 was employed to obtain search results in Ref. [10]. The non-ideal estimated value of κ likely resulted from a combination of thresholding and clustering effects on the trigger distribution, and from not accounting for PSD variation in the chisquared test calculation. In this work we consider the effects of non-stationarity in both the SNR time series and the chi-squared test. Since the SNR scales as S −1/2 E , we correct the trigger detection statistic dynamically re-scaling the SNR by a multiplicative factor v s (t) −0.5 . The chi-squared time series depends on the power of the signal, we therefore re-scale it by v s (t) −1 .
In Figure 5 we present how this correction affects the single detector background distribution during the same period considered in Section 4. Each plot shows the cumulative trigger rate as a function of the detection statistic with and without the PSD variation correction. We remove triggers generated by real signals to only examine the noise distribution. The distribution of triggers associated to long duration templates (12.6 -36.6 seconds) changes drastically for both Hanford and Livingston detectors. Short templates however are only slightly affected by our correction. Short duration templates also appear with a very high detection statistic. In the case of Gaussian noise, we would expect to obtain similar distributions between short and long templates. This asymmetry suggests the presence of short noise transients which are not correctly identified by the pipeline affecting the sensitivity of the search [27]. The background distribution for short template will require further investigations.
To quantify the improvement in the search sensitivity due to the PSD variation correction, we estimate the sensitivity volume of the search. Considering just the two LIGO detectors, we can naively estimate the network coincident trigger rate, R N as: where R H1 and R L1 are the trigger rate for Hanford and Livingston detectors and ∆t is the time-distance between them. Noticing that the performance are similar for the two detectors, we can then approximate R H1 R L1 . Fixing the network false alarm rate (FAR) to 1 per year we then obtain R H1 = R L1 = 1.3 × 10 −3 , corresponding to a single detector FAR of 1 per 7 days. We finally estimate the volume increase as: whereρ is the initial detection statistics for a single detector. We obtain a volume improvement of ∼ 1% for short templates, while for long templates the sensitivity increases by ∼ 7%. We also estimate the search sensitivity by adding 20,000 simulated signals to the data during the PyCBC analysis [13]. The fraction of injections identified by the pipeline at a given FAR determines the sensitivity volume of the search. This method shows volume improvements due to the PSD variation correction which are compatible with our previous estimates within uncertainties.

Conclusion
This paper presents a new approach to account for noise variations in interferometric gravitational-wave detectors which can limit the detection of compact binary coalescences. We developed a new statistic to track noise variations which is based on estimating the variance of the SNR. We have incorporated this method in the PyCBC search pipeline. In particular, we used our statistic to dynamically re-rank the detection statistics of triggers found during times of non-stationarity. Analyzing 12 days of O2 LIGO data, we found that our correction significantly changes the single detector background distribution for long waveform triggers. Consequently, we measure a 7% increase in the sensitive volume of the PyCBC search for binary neutron stars, neutron star-black hole and low mass black hole systems. In contrast, the sensitivity to higher mass systems did not show any significant improvement.
During O3 both LIGO detectors have been strongly affected by non stationary noise, such as scattered light [28,29]. For this reason we expect our approach to provide much larger sensitivity improvements in the analysis of O3 data. Accounting for noise variations will be even more essential in the future, when the gravitational-wave detectors will reach their design sensitivity. Moreover, the next generation of ground based detectors will have a higher bandwidth which allows to track non-stationarity over shorter time scales.