Subtracting glitches from gravitational-wave detector data during the third LIGO-Virgo observing run

The BayesWave algorithm assumes that in each interferometer, the detector data h(t) can be expressed as

$$\begin{align} h(t) = n(t) + s(t) + g(t) \end{align} \tag{ 1 }$$

where n(t) is Gaussian noise, s(t) is an astrophysical signal, and g(t) is a non-Gaussian instrumental glitch. BayesWave models non-Gaussian features (both signal and glitch) as a sum of sine-Gaussian (also called Morlet-Gabor) wavelets. In the time domain, the wavelets are given by

$$\begin{align} \Psi(t;\vec{\theta}) = A e^{-(t-t_0)^2/\tau^2}\cos\left(2\pi f_0(t-t_0)+\phi_0\right) \end{align} \tag{ 2 }$$

where $\vec{\theta} = \{t_0, f_0, Q, A, \phi_0\}$ are the parameters of each wavelet, and $\tau = Q/(2\pi f_0)$. The wavelet parameters, as well as the total number of wavelets used, is marginalized over via a transdimensional Markov chain Monte Carlo (MCMC).

For the glitch model, the wavelet parameters in each detector are independent. For the signal model, there is one common set of wavelets which are projected onto the detectors using a set of extrinsic parameters that describe the sky location and polarisation content. The power spectral density (PSD) is modelled using [50], which also uses a transdimensional MCMC to model the noise as a combination of a cubic spline and Lorentzians. Full details can be found in [29, 30].

The transient features in the data can be reconstructed using three different model assumptions:

• (a)  
  The data contain Gaussian noise with an astrophysical signal.
• (b)  
  The data contain Gaussian noise with an instrumental glitch.
• (c)  
  The data contain Gaussian noise with both a signal and a glitch.

To identify and remove glitches from detector data, we use either the glitch-only model (case II above), or signal-plus-glitch model (case III). The glitch-only approach is used when the glitch is sufficiently separated in time and/or frequency from the part of the signal that BayesWave can reconstruct. When using the glitch-only approach, we only need to consider data from the single detector in which the glitch appears. If the signal and glitch overlap significantly, we use the signal-plus-glitch model. The signal-plus-glitch model requires data from multiple detectors in order to separate the coherent signal power and the incoherent glitch power. By simultaneously fitting the for signals and glitches, we ensure that we do not subtract any signal power during glitch mitigation.

As BayesWave uses an MCMC to marginalize over the parameters, the end result is a posterior distribution of time series of the glitch, g(t). To actually perform the glitch mitigation, we select a fair draw from the posterior [31] by randomly selecting a single sample from the posterior, and subtract it from the detector data. The timescale to model glitches and produce glitch-subtracted data using BayesWave is up to multiple days depending on the type of glitch that is being analysed [30]. The code for both BayesWave and the BayesWave glitch subtraction methods can be found at [51].

The codebase used to subtract glitches using an auxiliary witness, gwsubtract, is the same linear subtraction algorithm as was presented in [40]. In this section, we provide a short description of this method, and note how the use case considered in this work differs from the original use case presented in [40]. The code associated with gwsubtract is publicly available at [52].

In order to subtract noise, we first assume that the measured strain in the gravitational-wave detector, h(t), is a linear combination of time series from different sources, $n_i(t)$. We further assume that at least one of these noise sources can be modelled as the convolution of a witness time series, a(t), and an unknown transfer function $c_{ah}(t)$. In the time domain, the gravitational-wave strain can therefore be written as

$$\begin{align} h(t) & = n_1(t) + n_2(t) + \ldots + n_N(t) \nonumber\\ & = n_1(t) + n_2(t) + \ldots + a(t) * c_{ah}(t) \nonumber\\ & = h^{^{\prime}}(t) + a(t) * c_{ah}(t) \end{align} \tag{ 3 }$$

where all noise sources except one are included in $h^{^{\prime}}(t)$. This can be conveniently cast into the frequency domain as

$$\begin{align} \tilde{h}(f\,) = \tilde{h}^{^{\prime}}(f\,) + \tilde{a}(f\,) \cdot \tilde{c}_{ah}(f\,). \end{align} \tag{ 4 }$$

Following the derivation from [35, 40], we can estimate the unknown transfer function by using the discrete Fourier transform of the strain data, $\tilde{Y}_h(f\,)$, and the witness time series, $\tilde{Y}_h(f\,)$. A single value is calculated for the frequency band $[\,f_1, f_2)$ via

$$\begin{align} \tilde{c}_{ah}(f^{\,^{\prime}}) = \frac{df}{f_2-f_1} \sum_{f = f_1}^{f_2} \tilde{Y}_h(f\,) \tilde{Y}_a^{\,*}(f\,), \end{align} \tag{ 5 }$$

where df is the frequency resolution of the data.

By averaging nearby frequencies when calculating the transfer function, contributions to the measurement from chance correlations are reduced. Specifically, the fractional measurement error from this averaging method is expected to be $\sqrt{\frac{df}{f_2-f_1}}$ if both h(t) and a(t) are stationary Gaussian processes.

Although the basic description of the process for glitch subtraction is quite similar to [35, 40], the process has a few practical differences. In [35, 40], multiple witness time series were used to subtract contributions from a single noise source, requiring that correlations between each of the witness time series were also calculated. In this work, only a single witness time series is considered. As the glitch subtraction can be completed by gwsubtract with a few straightforward calculations, the timescale to produce this data for one event is generally less than multiple hours.

An additional key difference is the stationary nature of the time series. Data that contains glitches is, by definition, not a stationary Gaussian time series. Regardless of the use case, if the transfer function is still stationary and linear, then theoretically this should not prevent these linear subtraction methods from being used. However, it is possible that the true transfer functions may have some non-linear features that are not captured by this method. It is also possible that some of the methods used in [40] to reduce the impact of false correlations could prevent the transfer function from being measured accurately. For these reasons, the ideal type of glitch to subtract with these methods is one that occurs at a high rate over a short time period. This increases the probability that the transfer function will be accurately measured.

Due to these practical differences between glitch subtraction and the use case that [40] was originally designed for, we completed additional validation tests to confirm that the methods of [40] were able to be used for glitch subtractions. These tests are summarized in appendix A.

One of the key steps in the glitch subtraction process, after deciding whether to subtract excess power in the data, is determining whether this excess power was successfully subtracted. For these purposes, we employ a residual test inspired by [53] to evaluate the statistical significance of any non-Gaussian features in the data. This test relies upon comparing the PSD of the data during the time of interest against the PSD measured over 512 s. Any positive fluctuations in the short-duration PSD are considered excess power.

Following the description in [53], we aim to measure the PSD variation, ρ, which is

$$\begin{align} \rho = N \int_{t_0-\Delta t}^{t_0} |F(t) * s(t) |^2~\text{d}t \end{align} \tag{ 6 }$$

for some data, s(t), at a given time, t₀. $\Delta t$ is the timescale on which the PSD variation is measured, F(t) is a filter function to weight different frequencies that are being analysed, and N is a normalizing constant. In the frequency domain with discrete samples, this is

$$\begin{align} \rho = N \sum_{f = f_\text{low}}^{f_\text{high}} |F(f\,) s(f\,)|^2~\Delta f. \end{align} \tag{ 7 }$$

We approximate the squared frequency spectrum of the data, $|s(f\,)|^2$, as the PSD of the data measured using GWpy [54] over a short duration that overlaps the glitch. We refer to this as $S_G(f\,)$. Rather than designing a filter function based on a gravitational-wave inspiral signal (as was done in [53]), we choose to use a filter function that is weighted by the sensitivity of the detector. Specifically, $F(f\,) = 1/S_L(f\,)$, where $S_L(f\,)$ is the PSD of the data over 512 s of data measured using the 'median method' [55] with windows the same size as the short duration PSD and $50\%$ overlap. With these substitutions, the variation in the PSD, ρ is

$$\begin{align} \rho = \sum_f S_L(f\,) \times \sum_f \frac{S_G(f\,)}{S_L^2(f\,)}, \end{align} \tag{ 8 }$$

where the normalizing constant, N, is

$$\begin{align} N = \sum_f \frac{S_L(f\,)}{\Delta f}. \end{align} \tag{ 9 }$$

However, the variance of this statistic is dependent on both the time window used to estimate $S_G(f\,)$, $\Delta t$, and the 'effective bandwidth' considered, $\Delta f_\textrm{eff}$ [53], making it unsuitable for comparisons where these two quantities differ. Specifically, the variance is

$$\begin{align} \mathrm{Var}(\rho) = (\Delta t \Delta f_\textrm{eff})^{-1}. \end{align} \tag{ 10 }$$

The effective bandwidth is dependent on both the bandwidth considered, $\Delta f = f_\textrm{high}-f_\textrm{low}$, and the shape of the PSD. This effective bandwidth is calculated via

$$\begin{align} \Delta f_\textrm{eff} \approx \int_{f_\textrm{low}}^{f_\textrm{high}} \frac{\textrm{max} [S_L(f\,)]}{S_L(f\,)}df. \end{align} \tag{ 11 }$$

To account for this dependence on $\Delta t$ and $\Delta f$, we normalize the statistic by the standard deviation for the parameters we use in each instance and subtract the mean of the distribution. Hence our new statistic, $\tilde{\rho}$ is

$$\begin{align} \tilde{\rho} = (\rho - 1) \times \sqrt{\Delta t \Delta f_\textrm{eff}}. \end{align} \tag{ 12 }$$

This new statistic has a mean of 0 and a variance of 1 for all choices of $\Delta t$ and $\Delta f$.

Although this new statistic is less susceptible to analysis choices for the reasons listed above, the true variance of the statistic may be different due to practical limitations in the measurement of the PSD. For example, it is not possible to calculate this statistic when $\Delta f \lt \Delta t^{-1}$ as there will be no points in the measured PSD in the bandwidth of interest. When $\Delta f \approx \Delta t^{-1}$, we also see additional variance in the statistic that is not accounted for in the normalization due to the large frequency resolution of the PSD. If the frequency resolution of the PSD is sufficiently coarse, only a small number of data points are the bandwidth of interest, inflating the measured variance.

To address this issue, we calculate a p-value for the measured value of $\tilde{\rho}$ for a specific choice of $\Delta t$ and $\Delta f$ using simulated data. In each case, we generate 2048 s of stationary, Gaussian data with the measured long-duration PSD using Bilby [56, 57]. We then repeat the same procedure over each sequential short-duration time window for the whole time period simulated. These results are used to generate the distribution from which the p-value is estimated.

An example distribution of this statistic based on data around GW190425 [23] is shown in figure 1. Before subtraction, the measured statistic value ρ was ≈180, much larger than any trial in the simulated background. After subtraction, the statistic dropped to 1.4, which was higher than $18\%$ of the simulated noise, corresponding to a p-value of 0.18. This can be seen from the cumulative distribution of the simulated noise, shown on the left, and a histogram of the measured statistics in the noise, shown on the right. The plotted histogram also demonstrates that the peak of this distribution is indeed close to 0. The distribution of statistics is asymmetric, with a longer tail towards larger values. This asymmetry, even for cases that should be well-behaved, motivates the use of the p-value rather than a threshold on the measured statistic based on an assumed variance.

Zoom In Zoom Out Reset image size

In order to ensure that the glitch-subtracted data used in downstream analyses do not introduce biases, a series of glitch-subtraction review tests were completed as a part of the O3 LVK analyses [2–4]. These tests focus on ensuring that either the data was consistent with Gaussian noise or that any residual instrumental excess power does not bias estimates of an event's source properties. A flow chart providing an overview of the steps in this process can be seen in figure 2. Previous investigations into the potential biases on PE when glitches overlap or are nearby signals [14–17] have not yet covered the entire parameter space of glitch types and where in the signal the glitch overlaps. We therefore use conservative thresholds for each test to minimize potential biases.

Zoom In Zoom Out Reset image size

In order to identify if data surrounding gravitational-wave events required glitch subtraction, we first utilize the results of event validation procedures completed for each event. These procedures for O3 are detailed in [10, 58]. The outcome of event validation specifically indicates if any glitches were nearby a gravitational-wave signal. Events that do not have any glitches flagged by these procedures are deemed not to require glitch subtraction.

For flagged events, the residual test described in section 2.3 is used to determine if the data is inconsistent with Gaussian noise. The time–frequency bounds to use for the residual test are determined by visual inspection of the data around events using spectrograms produced with the Q-transform [59]. Cases where the residual p-value was greater than 0.1 are deemed consistent with Gaussian noise and no glitch subtraction is attempted, cases with p-values less than 0.001 are passed on for glitch subtraction, and cases with p-values between these two thresholds require further follow up before a decision was reached. In the cases when the result is inconclusive, the overlap of the excess power and signal is used to determine if glitch subtraction is needed; glitch subtraction is generally recommended in cases where excess power directed overlapped the signal. As previously described, glitch subtraction was completed with BayesWave in all but 1 case where gwsubtract was used. The choice of algorithm for each event is made based on the suitability of each algorithm to subtract the glitches in glitches. Technical considerations are not a part of this choice as the time to completion and computing usage are not a limiting factor in the generation of reviewed, glitch-subtracted data.

Once the glitch subtraction was completed, the data is reevaluated to determine if the excess power is sufficiently subtracted. The residual metric is recomputed for the glitch-subtracted data using the same time–frequency regions. The same p-value thresholds are also used to decide if the data is now consistent with Gaussian noise.

In addition to the residual tests, comparisons of PE results using different mitigation methods are used to evaluate if further glitch subtraction is needed or if a different glitch mitigation method should be used ¹¹ . To test if excess power in the data could bias analyses, we estimate the source properties of events using either the nominal lower frequency limit (20 Hz) or a higher lower limit that excludes the data containing glitches. In cases where the glitch subtraction is able to improve the stationarity of the data, analyses using the original and glitch subtracted data are also compared. If the glitch-subtracted data does not pass the residual test, different frequency bounds that exclude the flagged time–frequency regions are recommended if the different frequency bound changed the measured chirp mass ¹² posterior by more than 1σ and that the signal-to-noise ratio (SNR) is reduced by less than $10\%$. Other parameters, such as χ_p and ${\chi_\text{eff}}^{12}$ are also compared between analyses. Although not the deciding factor in determining which mitigation method is recommended, any differences in these parameters between analyses are flagged. Additional details about these tests can be found in appendix B.

Once either the data around an event passed the residual test or the identified glitches are deemed to not bias PE, the data is approved for use in downstream analyses. These procedures generally are completed for each event within a few days (for events not requiring glitch subtraction) to weeks (for events that require glitch subtraction). The timeline for this process is limited by the need for human input at multiple steps; in future observing runs, additional automation will be needed to significantly speed up these procedures.

The process of glitch subtraction can be broken into two stages: glitch modelling and glitch subtraction. The subtraction process is done simply by finding the difference between the glitch model and the original data. However, as described in section 2, the methods used by the tools described in this work to estimate the glitch model differ significantly. These differences are highlighted in the two representative events we have considered in this section.

The first data visualizations we will consider for our comparisons are spectrograms generated using the Q-transform [59]. Spectrograms of the data around GW190424 and GW200129 from LIGO Livingston can be seen in figure 3. The panels of this figure show the data before subtraction, after subtraction, and the difference of the two spectrograms. The power subtracted near GW190424 is confined to two distinct burst of power: one corresponding to power from the camera shutter glitch at ≈70 Hz, and the other one from the arches in the nearby scattered light glitch. These bursts of power are not directly overlapping the signal, but are near enough that they might bias PE results. Conversely, the power subtracted near GW200129 is more broadband, with visible glitch power below 70 Hz subtracted for multiple seconds that directly overlaps the signal.

Zoom In Zoom Out Reset image size

A comparison of the original strain data, the deglitched data, and the glitch models can be seen in figure 4. Both the glitch models and the strain data have been whitened using GWpy [54]. The same types of features that were visible in the figure 3 can be seen in the whitened time series. The glitch model for GW190424 is confined to two short bursts, while the glitch model for GW200129 is non-zero over the entire time considered.

Zoom In Zoom Out Reset image size

These results also demonstrate limitations of these two methods. BayesWave was not able to successfully subtract all of the excess power from the scattering glitches, while the linear subtraction was not able to subtract any power above 70 Hz. Scattered light glitches are a known challenge for BayesWave. Recent investigations [31] have shown that in order to fully model scattered light glitches with BayesWave, low-SNR wavelets must be used. These investigations also find that wavelets at the required SNRs are strongly disfavoured when using priors similar to those used by the analyses discussed in this work. The linear subtraction method is also limited by the input data used to witness the glitch power. In the case of GW200129, the auxiliary sensor used in the subtraction was negligibly coherent with the gravitational-wave strain above 70 Hz. In other cases, such as the broadband noise subtraction performed in O2 [40], the sampling rate of the auxiliary sensors can limit the range of frequencies where subtraction is effective. Ultimately the effectiveness of the linear subtraction method relies upon how faithfully the relevant auxiliary sensors are able to witness the noise source of interest.

As previously stated, the uncertainty in the model of a glitch overlapping a gravitational-wave signal introduces another potential source of uncertainty in the PE results. Although it is possible to sample over the glitch model as an additional component in the PE (as is possible with recent updates to BayesWave [19, 31]), this was not done in the O3 analyses. The impact of accounting for the statistical uncertainties in the BayesWave glitch model while performing gravitational-wave PE is explored in [31]. Instead of these methods, a single glitch model was subtracted from the data as a preprocessing step. Although this means that any glitch modelling errors are not included in the analysis, this is not expected to significantly bias the results if the chosen glitch model is at least representative of the true glitch model. In this subsection, we discuss how we assess errors in glitch modelling for both BayesWave and gwsubtract, and demonstrate that statistical errors in the chosen glitch models are qualitatively low.

As discussed in section 2.1, the Bayesian nature of BayesWave results in a posterior distribution of the glitch time series. Previous studies of both simulated signals [69] and real glitches [31] have shown that the accuracy of BayesWave is related to the duration of the data analysed; BayesWave more accurately reconstructs short-duration ($\lesssim$1 s) events. To produce the glitch-subtracted data, we use a fair draw from the posterior distribution. An alternative choice is to use the median of the glitch time series posterior, however in some cases this can lead to oversubtraction [31]. An example of the 90% credible interval of the glitch time series posterior for GW190424, along with the whitened glitch model (a fair draw from the posterior), is shown in figure 5. We see that although the fair draw is a good representation of a possible glitch model, we could get slightly different results depending on which glitch time series we happen to draw.

Zoom In Zoom Out Reset image size

Statistical errors in the glitch model from gwsubtract come from both chance correlations between the two data streams considered and errors in the transfer function measurement itself. As described in section 2.2, the errors from chance correlations can be calculated directly from parameters used in the glitch modelling. In the case of GW200129 the 1σ statistical errors from chance correlations should be $\frac{1}{\sqrt{2048}} \approx 0.022$. The transfer function measurement errors, however, are not directly estimated as part of the glitch modelling process. However, if we assume that the transfer function is stationary over a timescale longer than used to measure the transfer function, we can take multiple measurements to estimate this source of error.

In order to qualitatively estimate the significance of measurement errors, we recalculate the model of the glitches overlapping GW200129 after changing the parameters used in gwsubtract. Specifically, we shift the time window used to estimate the transfer function by $\{-1024, -512, 512, 1024\}$ s with respect to original run. The glitch model estimated for these four additional attempts is shown in figure 5 compared against the original model. The statistical errors on the original glitch model due to chance correlations is represented by the orange shaded region. All five runs show qualitative agreement in the amplitude and phase of the glitch, within the expected errors from chance correlations. We therefore conclude that the statistical errors from the glitch modelling procedure are low and should not introduce biases in the analysis.

It is important to note that these estimates of the errors do account for systematic errors due to modelling assumptions. Although the wavelet model used by BayesWave aims to be fully generic, this model, along with specific choices on the priors used in the analysis, is not able to fully model all glitch features for all glitches. Furthermore, specific glitch classes are known to be more problematic for BayesWave than others [31]. Likewise, gwsubtract is limited by the assumption that the witness sensor used to subtract the excess noise is a perfect witness for the noise source and that the transfer function between the witness sensor and the gravitational-wave strain is linear over the amplitudes of interest. These types of systematic errors are mitigated by investigations into the limitations of each glitch subtraction algorithm and the choice of which algorithm to apply in each scenario. This type of choice was one of the reasons that gwsubtract was used for one of the events in O3. Other problematic glitches impacting events in O3 were investigated using gwsubtract, but only the glitches around GW200129 were successfully mitigated with this algorithm.

The potential impact of unaccounted systemic error sources on PE analyses is not addressed in this work. However, studies of BayesWave-based PE suggest that systematic errors are low relative to other sources of error in estimation of gravitational-wave event source properties [31]. Systematic uncertainties in the case of gwsubtract are more difficult to probe, as these will vary with each noise source considered. Previous investigations into the use of gwsubtract and other tools for broadband subtraction did not identify biases [40, 70], but it is not clear if these studies are fully applicable in this case due to this differences explained in section 2.2. As the sensitivity of the detectors increases, these systematic uncertainties may have a larger impact on the measured source properties. Recent investigations into GW200129 have already shown that the evidence of orbital precession from this event may be impacted by systematic errors in the glitch modelling [65].