{'fraction_non_alphanumeric': 0.03608216642283462, 'fraction_numerical': 0.015854059057116413, 'mean_word_length': 4.699081163859112, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accurate modeling and mitigation of overlapping signals and glitches in gravitational-wave data

The increasing sensitivity of gravitational-wave detectors has brought about an increase in the rate of astrophysical signal detections as well as the rate of"glitches"; transient and non-Gaussian detector noise. Temporal overlap of signals and glitches in the detector presents a challenge for inference analyses that typically assume the presence of only Gaussian detector noise. In this study we perform an extensive exploration of the efficacy of a recently proposed method that models the glitch with sine-Gaussian wavelets while simultaneously modeling the signal with compact-binary waveform templates. We explore a wide range of glitch families and signal morphologies and demonstrate that the joint modeling of glitches and signals (with wavelets and templates respectively) can reliably separate the two. We find that the glitches that most affect parameter estimation are also the glitches that are well modeled by such wavelets due to their compact time-frequency signature. As a further test, we investigate the robustness of this analysis against waveform systematics like those arising from the exclusion of higher-order modes and spin-precession effects. Our analysis provides an estimate of the signal parameters; the glitch waveform to be subtracted from the data; and an assessment of whether some detected excess power consists of a glitch, signal, or both. We analyze the low-significance triggers (191225_215715 and 200114_020818) and find that they are both consistent with glitches overlapping high-mass signals.


I. INTRODUCTION
Gravitational-wave (GW) analyses require accurate models for both the astrophysical signals and the detector noise [1]. The majority of source properties inference for transient signals such as compact binary coalescences (CBCs) is based on three assumptions about the detector noise that inform the functional form of the likelihood function: (i) the detector noise is uncorrelated between the detectors, (ii) it follows a Gaussian distribution, and (iii) it is stationary, i.e. its mean and covariance do not change with time. Violation of these assumptions could impact detection and inference efforts. For example, Schumann resonances in the Earth's large-scale magnetic field could cause correlated detector noise and affect detection and interpretation of a stochastic GW background [2][3][4][5][6][7][8]. Additionally, the detector Gaussian noise is stationary over short timescales [9], but longer signals might be subject to noise nonstationarity which has motivated relevant studies [10][11][12][13].
Transient noise artifacts, i.e., glitches, in a detector violate the assumption of Gaussianity and could bias parameter inference when they overlap with signals [14][15][16]. In the recent third observing run (O3), LIGO [17] and Virgo [18,19] have detected an astrophysical event approximately once every five days [20,21]. Glitches, however, appear in the detectors far more frequently. The average rates for glitch transients with signal-to-noise ratio (SNR) > 6.5 in the first and second half of the third observing run were 0.3 min −1 in LIGO Hanford (LHO), 1.13 min −1 in LIGO Livingston (LLO), and 0.75 min −1 in the Virgo detector [20,22]. Overall, in O3 a total of 18 events required some form of glitch mitigation [20,22]. Glitches are most likely to intersect lower-mass, longduration events such as binary neutron star (BNS) mergers; indeed, both such detected events have overlapped with a glitch and required mitigation [23,24]. As detector sensitivity improves not only will the event rate increase, but also the glitch rate might increase as weaker glitches that are currently below the noise floor could emerge above it. Additional detectors such as KAGRA [25] in the next observing run (O4) also increase the likelihood that a glitch will appear in at least one detector. In order to have an accurate and unbiased catalog of GW events, effective and generic methods for separating signal and glitch power are necessary. Proposed approaches for glitch subtraction include removing the contaminated data [23,[26][27][28][29][30] or subtracting detector noise based on data from auxiliary channels [13,[31][32][33][34][35][36][37][38][39]. The glitches discussed in this paper are those that remain after the noise mitigation described in [35,39].
A complementary analysis was proposed in Ref. [40] based on the BayesWave algorithm [41][42][43]. This analysis expands glitch-mitigation techniques already applied to LIGO and Virgo data [20,22,23], where it was used to subtract glitches in 15 out of 18 O3 events that required glitch mitigation [20,22]. The analysis of [40] simultaneously models the signal and glitch using waveform templates and wavelets respectively. The waveform templates are models for CBC signals that are obtained arXiv:2205.13580v2 [gr-qc] 17 Aug 2022 as solutions to the 2-body problem in General Relativity. The glitch model is based on a sum of sine-Gaussian wavelets and it is flexible enough to be able to reliably describe a wide range of potential glitch morphologies. This first study presented a number of examples of data containing binary black hole (BBH) signals and common glitch types and demonstrated the ability of the analysis to reliably separate the two [40]. In this study we expand upon this analysis by considering a wider rage of glitch classes, more instances of each class, and CBC injections of varying masses.
Our analysis results in a posterior distribution for the parameters of the glitch and the CBC signal. Depending on the exact placement of the signal in relation to the glitch, correlations between the two might exist and the resulting CBC posteriors might not be identical to those from data with no glitches. This is unsurprisingly most prominent for signals and glitches with similar timefrequency morphology. Despite this, the CBC parameters are correctly estimated within the extent of the posterior. As a point of comparison for each case, we also examine the bias incurred on CBC parameters by performing a standard analysis that ignores the presence of the glitch in the data entirely.
Our process allows us to obtain both a model for the glitch that can be subtracted from the data and parameter estimation results for the CBC signal, though the latter are restricted by the assumption of alignedspins in the current algorithm implementation. The glitch-subtracted data are ready for downstream analyses with more sophisticated waveform or detector calibration models. We demonstrate that the lack of spinprecession and higher-order modes in our CBC waveform models does not hinder accurate CBC-glitch separation. We also test whether we can do glitch subtraction on single-detector data. The additional examples and checks presented here demonstrate that the analysis is ready to tackle incoming data in O4 [44].
Finally, we consider some low-significance triggers and attempt to distinguish between CBC signals, Gaussian noise, glitches, or a combination of all three present in the data. Standard detection efforts consider only the possibility that either a CBC signal or a glitch is present in some data. It is therefore possible that the significance of a CBC signal could be impacted if it overlaps with a glitch as this could make the data inconsistent with our CBC model. We find that trigger 191225 215715 (hereon S191225) [20] and trigger 200114 020818 (hereon S200114) [45] are consistent with the presence of both a glitch and a CBC signal. The above results are subject to the caveat that we use a CBC waveform model without spin-precession or higher-order modes.
The paper continues as follows. In Sec. II, we describe the BayesWave algorithm, focusing on the CBC (templated) and the glitch models. In Sec. III we present our injection and recovery scheme. In Sec. IV we present our results on an array of different glitch classes. We test the systematics and limitations of our CBC sampler in Sec. V by including injections containing higher order modes, spin-precession, and events in a single detector. Finally, in Sec. VI, we apply our analysis on triggers S191225 and S200114 in order to assess the presence of signals and/or glitches in the data. In Sec. VII we conclude.

II. GENERAL METHODOLOGY DESCRIPTION
In this section we describe the main features of our analysis that estimates parameters of GW events using templates while jointly modeling detector glitches using wavelets. We briefly introduce our inference scheme, then discuss the features of BayesWave relevant for this study including the models, priors, and the joint sampler.
A. Brief Introduction to Inference Scheme GW parameter estimation aims to compute p(θ h |d), the posterior probability distribution that a model h with parameters θ h describes the given data, d. The quantity h can contain any component of the data we attempt to model, for example a CBC signal, a glitch, or Gaussian noise. In the case of CBC signals, h is typically a waveform template and θ h are parameters that describe the binary. The posterior p(θ h |d) is proportional to the prior p(θ h ) times the likelihood of observing data d given the model h(θ h ), p(d|θ h ). The likelihood function encodes our assumption about the detector noise. For LIGO analyses and stationary, Gaussian noise with power spectral density (PSD) S n (f ) in each detector, the likelihood is given by Eq. 3 in [43]. Under this formulation, the multidimensional parameter posterior p(θ h |d) is typically explored with stochastic sampling methods such as Markov Chain Monte Carlo (MCMC) and nested sampling [46,47].

B. BayesWave Analysis
BayesWave [41][42][43]48] is a flexible data analysis algorithm that models various components in GW data including signals, glitches, and noise. BayesWave explores the joint posterior distribution of its models using a combination of MCMC and reversible jump (RJ) MCMC [49] samplers. The algorithm has a wide range of applications including an unmodeled search pipeline [50,51], on-source PSD generation [9,20,21], tests of General Relativity [52][53][54][55] and consistency tests [56], and various studies of poorly modeled or unmodeled sources including binary neutron star (BNS) postmergers [57][58][59], eccentric BBH mergers [60], and supernova [61]. BayesWave has also been used to perform glitch subtraction [15,20,21,23]. With this flexibility to model a wide range of signal and glitch morphologies in hand, in Ref. [40] we extended the algorithm capabilities to include a model for CBC signals based on CBC waveform Time since t trig [s] templates. The algorithm implementation is described in detail in [41][42][43]62] and we briefly summarize the most relevant features here. Whereas most GW parameter estimation analyses [46,47] use single, point estimates of the PSD and assume that the data contains no glitches, BayesWave can relax both assumptions by modeling the CBC signal, the noise PSD, and any potential transient noise all at once. 1 The full BayesWave model for this study is the union of a CBC (waveform template) model, a glitch model, and a noise PSD model. The different BayesWave models are described below.
• The noise model, or PSD model, expresses the noise PSD as a sum of splines and Lorentzians. The splines capture the smooth underlying broadband noise whereas the Lorentzians capture the sharp spectral peaks. Within this paper, we use the color grey to represent the noise model.
• The glitch model expresses excess detector transient noise as the sum of sine-Gaussian wavelets: accordingly, each detector has its own, independent, set of glitch wavelets. The set of all sine-Gaussian wavelets form an overcomplete basis over smooth function space, and are thus able to describe a glitch of any morphology provided that it is sufficiently loud. The number of wavelets and hence the model dimensionality is not fixed and wavelets can be added or subtracted as needed. Each wavelet is 1 The effect of PSD uncertainty in CBC analyses has also been explored in [63][64][65][66][67][68][69][70][71].
described by five parameters (θ glitch :) central time, frequency, quality factor, amplitude, and phase. The quality factor is related to the damping time of the sine-Gaussian, and together with the frequency determines the duration of each wavelet. The functional form of the wavelet is given by Eq. 4 in [43].
Within this paper, we use gold to represent the glitch model.
• The CBC model uses waveform templates to capture the CBC signal in a manner similar to traditional GW parameter estimation [46,47]. Details of the implementation of the CBC model in BayesWave are given in [62]. For this analysis we restrict to quasicircular CBC signals whose spins are aligned with the orbital angular momentum. Such signals are characterized by up to 13 parameters, namely four intrinsic CBC parameters (the two masses m 1 , m 2 and two dimensionless spins χ 1 , χ 2 projected onto the Newtonian orbital angular momentum) and seven extrinsic parameters (a time and phase, the right ascension and declination, the luminosity distance D L , the inclination ι, and the polarization angle). For binary neutron stars (BNSs), we also have two tidal parameters Λ 1 , Λ 2 .
In what follows, we express the masses through the chirp mass M c ≡ (m 1 m 2 ) 3/5 (m 1 + m 2 ) −1/5 which determines the GW phase evolution to leading order and the mass ratio q ≡ m 2 /m 2 < 1. We also express the spin through χ eff ≡ (m 1 χ 1 +m 2 χ 2 )/(m 1 + m 2 ), which is conserved approximately throughout the binary inspiral [72]. Within this paper, we use pink to represent the CBC model.
• Though not used here, for completeness we also mention the signal model that fits for coherent, excess power ("unmodeled" astrophysical signals) using again sine-Gaussian wavelets. Unlike the glitch model, the signal model enforces that the wavelets must be coherent across the detector network as a genuine astrophysical signal would be. Both the signal and the CBC model have the potential to capture a CBC source though the former is more flexible, and thus less sensitive, particularly to weak or long-duration signals.

C. Priors
The priors for the glitch and CBC model parameters remain mostly unchanged compared to [40,62]. However, for some combinations of glitch and CBC signals, we find that additional flexibility is required in the time placement of glitch wavelets and the CBC template. By construction, BayesWave analyzes data of duration T obs around some trigger time t trig . The prior on the central time of the glitch wavelets then has support within a "window" of length T w < T obs around t trig , while the CBC time prior is by default (t trig −0.5s, t trig +1.5s). Here we relax the requirement that glitches and CBC signals have a time prior around a common time t trig and allow for them to be placed in different time intervals, though still within T obs . The priors for the noise PSD model are all unchanged compared to [40,42]. Figure 1 shows the relation between the different time intervals.
Though the wavelets that model the glitch can have central times only within T w , the current BayesWave implementation employs a preprocessing "cleaning phase". During this phase, the algorithm is run with only the glitch model activated and T w = T obs , i.e., wavelets can be placed anywhere in the analyzed data segment. At the end of the cleaning phase and before proceeding to the main analysis, wavelets with time centered within a specified interval are permanently subtracted from the data. We typically use a 1s cleaning window at the beginning and end of the data segment. This procedure removes any glitch power that might be present in the analysis segment but not necessarily close to the CBC signal itself as well as data artifacts caused by the finite segment duration and that could bias the PSD estimation.

D. Sampler
To characterize the joint glitch, CBC, and noise posterior, BayesWave uses a combination of MCMC and RJM-CMC samplers stringed together within a blocked Gibbs sampler. The blocked samplers give us the flexibility to trivially turn models on and off during an analysis. Here we sketch the workflow for the "CBC+Glitch" analysis, but other BayesWave running modes with different model combinations vary only in which samplers are active.

FIG. 2.
Visual depiction of the BayesWave workflow described in Sec. II. Each colored box represents a component sampler and displays its name (bolded) and its input data (in LaTeX). Pink (gold, grey) boxes indicate the sampler that searches over the CBC (glitch, PSD) parameter space. Underneath each colored box, the sampled (fixed) parameters are boxed (unboxed). The Gibbs sampler (i.e., the cycle of component samplers) is boxed in red. Before entering the Gibbs sampler, we use the Fourier domain strain data,d(f ), to obtain initial points for the noise and CBC samplers through the procedure described in [73]. We also optionally (dotted box) obtain an initial point for the glitch model [43], otherwise the glitch model begins with 0 wavelets. Within the Gibbs sampler, each component (RJ)MCMC goes through O(10 2 ) iterations before passing its values to the next sampler. After each Gibbs sampler loop, a single sample for all parameters is returned. The Gibbs sampler loops O(10 4 ) times.
Instead of sampling all parameters concurrently, we separate them into groups of related parameters that are sampled together while other parameters are held fixed. The order of the corresponding samplers and a breakdown of the fixed or varying parameters within each sampler is displayed in Fig. 2. Each component sampler runs for a predetermined number of iterations, typically O(100) and returns its last iteration to be used as fixed parameters by the other samplers. We iterate through the Gibbs sampler loop O(10 4 ) times before completing.
Before the Gibbs sampler begins, each model (i.e. PSD, CBC, and glitch) needs to be initialized. An initial estimate for the PSD is generated by the methods described in [43]. To initialize the CBC parameters we follow [73]. An optional GlitchBuster step finds initial parameters for the glitch wavelets by iteratively estimating the PSD, wavelet-transforming the data, and removing excess power wavelets [74]. The procedure is described in more detail in [43]. Without GlitchBuster, the glitch model begins with no wavelets.
The sampling procedure, and specifically the integration of the CBC sampler is described in detail in [62]. Below we briefly describe each individual sampler.
• The "wavelet (glitch) RJMCMC" updates the parameters of one of the current wavelets or adds/subtracts a wavelet. Details of the RJMCMC implementation are presented in [41,43].
• The "extrinsic MCMC" updates the extrinsic parameters of the signal, namely the distance, skylocation, inclination and polarization angles, and time 2 . Details are provided in [41,43,62].
• The "CBC MCMC" updates the intrinsic CBC parameters (masses, spins, tides) as well as parameters that are correlated with them, namely the distance, time, and phase. We use waveforms implemented in the LALSimulation suite of models [75]. The sampler can operate with any nonprecessing model available in LALSimulation and in this study we rely on IMRPhenomD [76,77]. We currently do not account for nonaligned spin degrees of freedom but plan to extend our analysis to include the effect of spin-precession in the future. The sampling proposals as well as the heterodyne procedure [73,78] used to speed up the likelihood calculation are described in [62].
• The "noise RJMCMC" updates the number and parameters of the splines and Lorentzians that describe the noise PSD. Details are provided in [42].

III. DATA AND ANALYSIS SETUP
To test the efficacy of our analysis, we use LIGO data from O2 and O3, available through the GW Open Science Center [79]. Since we do not have exact first-principles models for glitches 3 , we identify data containing genuine 2 Though not used in the study, the BayesWave signal model that describes a GW signal with coherent sine-Gaussian wavelets also makes use of the extrinsic sampler. 3 Phenomonological models for some glitch types have been constructed, for example [80]. detector glitches, which were classified through Gravity Spy [81]. We then inject a known CBC signal on top of the glitch, and analyze the data by simultaneously modeling the CBC, glitch, and noise. We label such analyses as "CBC+Glitch". In each case we also analyze the same data using only the CBC and noise models, i.e., ignoring the presence of the glitch. We label such analyses as "CBCOnly". To create a point of reference for the glitch, we also consider the original data with no CBC injection with a "GlitchOnly" run using only the glitch and noise models.
Going beyond the study of [40], here we analyze a more extensive set of glitch types (as classified by GravitySpy [82,83]), glitch instances, and a set of both BBH and BNS injections. The injected signals include high-mass (HM) BBH (GW150914-like [84]), low-mass (LM) BBH (GW170608-like [85]), and BNS (GW170817-like [23]) injections; their parameters are provided in Table I. Though initially we targeted specific glitches, in many cases secondary glitches (occurring in either or both detectors) were present in the data, speaking to the high occurrence rate of glitches. In those cases we do not discard the data; we analyze them nonetheless and attempt to model all glitches. These analyses also show that BayesWave can differentiate between signals and glitches even when they occur in multiple detectors. A spectrogram of one glitch per glitch family is shown in Fig. 3. Although the "worst case scenario" for CBC analyses is a glitch with SNR similar to or greater than that of the signal that has a similar time-frequency morphology, in our validation studies we include a wide range of combinations of signals and glitches.
We fix the intrinsic parameters and the sky location within injections of the same CBC type. For each com-FIG. 3. Spectrograms of representative glitches for each glitch family we consider in our study. We display two spectrograms for the same fast and slow scattering glitches to demonstrate the long-and short-term behavior of such glitch types. See Tab. IV for glitch GPS times.
bination of glitch and CBC type, we inject the signal at different times with respect to the glitch in order to vary the amount of overlap between the two. Each signal has a network signal-to-noise ratio (SNR) of 15 and we use the IMRPhenomD [76,77] (IMRPhenomD NRTidalv2 [86][87][88]) model for the BBH (BNS) signals for injection and recovery unless otherwise indicated. For computational efficiency and since the BNSs overlap with all glitches at low frequencies (below 40Hz), we use a low sampling rate that does not allow us to extract tidal parameters. The possibility of separating BNSs and glitches when they overlap in the high frequencies relevant for tidal effects was shown in [89] using the same analysis as here.
The recovered CBC parameters can be compared with the injected ones directly to assess the recovery reliability. However, no such comparison is possible for glitches as they are obtained from real data. We therefore use various overlaps (O, defined in Eq. 8 of [56]) and mismatches (one minus the overlap, M = 1 − O) between the recovered CBC and glitch to quantify the quality of the CBC-glitch separation. Definitions are given in Table II. Specifically, for each glitch we also analyze the original data with no CBC injection using only the glitch and noise models. This gives us access to a baseline estimate for the glitch that can be compared to estimates of the same glitch recovered in data coincident with a signal. We then compute the mismatch between the median glitch recovery from analyses with (g CBC+G ) and without (g G ) injected CBC signals as a way to test whether the recovered glitch subsumes part of the CBC signal power. A high mismatch could suggest that either the glitch model captures part of the CBC (and leads to it being inadvertently subtracted from the data) or the glitch model misses part of the glitch and fails to subtract all glitch power. Additionally, we compute the posterior for M(h rec CBC |h inj ), the mismatch between the injected and the recovered CBC signal whose expected value is a function of SNR [90].
In each case we also carry out a "CBCOnly" analysis that is akin to traditional parameter estimation (with the difference that we also marginalize over the PSD). Though such an analysis is not well motivated as the model cannot exactly match the data (and we encounter increased convergence issues), we find it instructive to compare its results to the full "CBC+Glitch" analysis. We do so by plotting both posteriors when discussing each glitch, and also via the Jensen-Shannon divergence in Sec. IV G. We denote the recovered CBC signal from such an analysis as h rec CBC and compute M(g G |h rec CBC ) as an estimate of how much of the glitch the CBC model can recover. We also compare M(h rec CBC+G |h inj ) against M(h rec CBC |h inj ) to test the effect of assuming pure Gaussian noise on CBC recovery in the presence of a glitch. An example result for the various overlaps is presented CBC |gG) (teal/gold, magenta/gold, and blue/gold; right panel) meaning that the CBC model is absorbing part of the glitch power, but doing less so when the glitch is part of the model. However, such overlaps are quite low, so the effect is small.
in Fig. 4 and described in detail in Table III Another potential test that was explored over the course of this analysis is the Anderson-Darling statistic [91]. This test can be used to assess the degree of Gaussianity in the residual and has been proposed in the context of PSDs [9]. Specifically, we explored the option of subtracting some point estimate (such as the median, or a fair posterior draw) for the CBC and the glitch models from the data and then computing the Anderson-Darling statistic. However, we found that the test is very forgiving and even fails to identify large amounts of residual glitch power. We attribute this to the fact that the form of the Anderson-Darling test we employ is better suited for identifying large non-Gaussian tails in distributions (and thus it is well suited for PSDs [9,64]) than for coherent non-Gaussian residual that affects only a few data points.

IV. RESULTS ON OVERLAPPING SIGNALS AND GLITCHES
In this section we present results from our injection study. Each subsection corresponds to a different glitch family and injections using three different CBC sources, see Table I. In many cases, especially for the longer signal analyses, the data contain additional glitches beyond what was intended. We point these cases out, but our inability to easily find glitch-free data of duration 1 − 2 minutes speaks to their prevalence in LIGO data. Details about the times analyzed and the glitches that we encountered either intentionally or accidentally are provided in Table IV, together with analysis settings.

A. Blip Glitches
Blip glitches are short-duration glitches that occur in the most sensitive frequency band of the detectors with a frequency range from ∼ 32Hz up through 1024Hz. Because of their short duration and prevalence, they challenge analyses of high-mass events [92], especially for sources with unequal masses and high spins [93]. Their origin is largely unknown; < 8% of LHO glitches and < 2% of LLO glitches were identified with auxiliary channels during the first and second observing runs [92]. Details about the blip glitches used for this study are provided in Table IV while a spectrogram is provided in Fig. 3. By chance, the data around GPS time 1165578732 contain a second glitch in LHO, 1s later. The presence of the additional glitch did not require any modification of the priors since it occurs entirely after the signal and has a low SNR. Figure 5 presents our results, with runs labeled according to Table IV. Each row corresponds to the same data with a given glitch and the same CBC signal injected at various times relative to the glitch. The merger time of the CBC relative to the glitch center is given on the y axis.
The first and second column show various mismatches following the format of Fig. 4. The set in the first column measure how well the models were reconstructed: the lower the mismatch, the more faithful each model recovery is. Specifically in all cases we find M(g G |g CBC+G ) ≤ 0.01, suggesting that the recovered glitch model does not consume any significant amount of the injected CBC signal, nor does it miss part of the glitch due to the presence of the signal. The first column also shows that typically M(h rec CBC+G |h inj ) ∼ M(h rec CBC |h inj ), though we also occasionally find M(h rec CBC+G |h inj ) < M(h rec CBC |h inj ); in these cases the CBC model captures part of the glitch if we ignore the presence of the latter by analyzing the data assuming pure Gaussian noise.
In the second column we present information quantifying how similar the glitch the CBC are and, by extension, how difficult it is to separate them. Again the format is similar to Fig. 4. Since O(g G |h inj ) is evaluated on the injected CBC parameters and not maximized over CBC parameters, it does not directly correspond to how well a CBC template can recover the glitch, but it is a conservative estimate of the similarity between the glitch and injected signal. All overlaps are small, however these is some clear variation. Specifically, we find O(g G |h inj ) ∼ O(g G |h rec CBC+G ) ≤ O(h rec CBC |g G ) which means that in the "CBCOnly" analysis the CBC model absorbs part of the glitch power. This is not the   Table II. An example result for these mismatches is plotted in Fig. 4. Throughout, pink and light blue are used for CBC reconstructions with and without the glitch model respectively, gold is used for glitch reconstructions, and teal is used for CBC injections. The split colors indicate the two models used for the overlap or mismatch.  Fig. 1), the type of injected signal, the glitches encountered and the affected detector(s), the segment length T obs , the duration of the time window where glitch wavelets can be placed Tw, the maximum quality factor of the glitch wavelets Qmax, the maximum number of wavelets allowed Dmax, and the low-frequency cutoff f low . All BBH runs have a sampling rate of 2048Hz whereas the BNS runs used a sampling rate of 1024 Hz (which precludes any possibility of recovering tidal parameters). Horizontal lines group analyses that target overlapping data, though the exact center of the glitch window (GPS time) might be shifted according to long-and short-term glitch behavior.   case with the full "CBC+Glitch" model; indeed in all cases O(g G |h rec CBC+G ) is closer to the original value of O(g G |h inj ).
The remaining columns show the marginalized posterior distributions for the (detector frame) chirp mass M c , luminosity distance D L , and effective spin χ eff respectively. The light blue posterior correspond to h rec CBC ("CBCOnly") whereas the magenta posteriors correspond to h rec CBC+G ("CBC+Glitch"). In all cases the parameters recovered with the full "CBC+Glitch" model are consistent with the injected value. We find differences in the posteriors for the same CBC signal injected at different times with respect to the glitch. This is expected for two reasons. Firstly, the glitch and CBC posteriors are not completely uncorrelated, and hence the marginal CBC posterior will not be exactly the same as if there was no glitch. Secondly, each CBC is injected at slightly different times, and hence is subject to a different realization of the detector Gaussian noise. This distinction becomes more important for shorter signals, and indeed we find that the posteriors become more similar as the CBC signal duration increases from top to bottom. Additionally, we find numerous instances where the "CBCOnly" posteriors are significantly shifted and even inconsistent with the injected value. These cases are typically accompanied by M(h rec CBC+G |h inj ) < M(h rec CBC |h inj ) (first column) and/or O(h rec CBC |g G ) < O(g G |h rec CBC+G ) (second column). Similar biases were reported in [94] for extrinsic signal parameters computed in low latency.
Blip glitches are one of the most common glitch types and are very similar to high-mass BBHs. However, they are also one of the most straightforward glitch types to deal with due to their compactness in time and similarity to wavelets. The analyses presented here typically used the default BayesWave glitch settings (apart from cases where there were additional glitches in the data beyond blips), and could be easily automated.

B. Slow-Scattering Glitches
Glitches from slow-scattered light appear in the detectors as long duration, O(4s), arches in the time-frequency domain, evenly stacked in frequency, usually in the range 8 − 64 Hz. Each set of arches often recurs multiple times in the detector as shown in Fig. 3. Unlike blip glitches, they are not morphologically similar to CBC signals, yet they create long periods of non-Gaussianity and nonstationarity in the data, thus posing a challenge for noise PSD estimation. Their rate of occurrence increased during O3, when slow scattered light glitches overlapped with nine events throughout the observing run [20,21]. Due to their morphology, slow scattered light glitches required longer analysis segments and wavelets of higher quality factors compared to short duration glitches such as blip glitches. A few of the glitches also extended to lower frequencies than other glitch types, so we used f low = 8 Hz for come cases. See Table IV for run set-tings. Figure 6 presents our results. All recovered CBC posteriors from the full "CBC+Glitch" analysis are consistent with the injected values. Unlike the blip glitch case discussed in Sec. IV A, the "CBCOnly" analysis that ignores the presence of the glitch returns largely unbiased posteriors as well, exhibiting mostly small shifts. This is likely due to the fact that fast scattering glitches are morphologically very different than the types of CBC signals we consider.
The left column shows some variation between the median recovered glitch reconstructions. Though M(g G |g CBC+G ) remains mostly low and around 0.01, it can reach ∼ 0.1 in some cases mostly for the second scatter glitch (runs whose label starts with S2). We explore this further in Fig. 7 where we plot the spectrum of the data, signal, and glitch as well as spectrogram of the data for the S2 HM BBH injection at −500ms compared to the glitch. In each row, we plot the spectrum (left panel) and subtract from the data (right panel) the median glitch reconstruction (top panel) or a fair draw from the posterior (middle panel). In the middle right panel, more of the glitch has been subtracted compared to the top right panel. This is due to the low SNR of the glitch (11.9 for S2 compared to 15.1 for S1; computed by the Omicron pipeline [95]), which results in some of the scattering arches residing in the threshold for reconstruction by the glitch model and thus not consistently included in the median. As a result, the median reconstruction has a large variation between different analyses, resulting in higher values for M(g G |g CBC+G ). In such cases, the glitch-subtracted data are sensitive to the choice of which glitch reconstruction to subtract (median or some fair draw) and additional case-by-case attention is needed.
Motivated by the low SNR of the S2 glitch, we also considered the effect of the default priors in BayesWave. The prior on the amplitude of the wavelets is broad, but peaks at SNR=5 per wavelet by default, see Fig. 5 of [43]. The bottom panel of Fig. 7 shows results with a wavelet amplitude prior that peaks at SNR=1 per wavelet. Clearly more of the low-SNR glitch is subtracted. We leave further tests of the prior tunings such as this to future work. We conclude that although analysis of high-SNR slow scattering glitches can be potentially automated, lower-SNR instances will require some user attention.

C. Fast Scattering Glitches
Fast-scattering glitches (also referred to as "crowns" [82]) are long duration glitches composed of many short bliplike bursts in frequencies from 10 − 60 Hz. Fast-scattering glitches have been linked to light scattered off the LIGO optical systems, particularly during ground motion [82]. They were the most common glitch type in LLO in O3, comprising 27% of all glitches [82]. Two spectrograms of fast-scattering glitches are given in Fig. 3 and display the long-and  Table IV for run settings and labels. Spectra and residual plots for the S2 HM BBH injection at −500ms from the glitch, described in Sec. IV B. The first column shows the power spectra of h inj (teal), the raw data (grey) and a point estimate for the PSD (black), h rec CBC+G (pink), and gCBC+G (gold). The second column shows spectrograms of the original data without the injection, while the third column shows spectrograms of the data where an estimate of the glitch has also been subtracted. The first and second row show results with the median glitch reconstruction and a fair draw from the posterior respectively. The third row shows results where the glitch amplitude prior set to peak at SNR=1 that aids identification of lower-SNR glitches. Due to the low SNR of the glitch, we find that all these choices impact the quality of glitch reconstruction and the subtraction.

Frequency [Hz]
short-term glitch behavior; relevant run settings are presented in Table IV.
Since a single fast-scattering glitch contains multiple, time-separated bursts, some adjusted settings are required. Such glitches create long-term nonstationarity, particularly at low frequencies so we increase the duration of the analyzed segments from 16 to 32 seconds for the low-mass BBHs. Despite the overall long duration of the glitch, we found that an increase in Q max is not necessary, as each glitch consists of individual shorter burst that are each reconstructed by a few low-Q wavelets. Figure 8 presents our results. Overall, we find similar results to the slow-scattering case of Fig. 6, as the full "CBC+Glitch" analysis is able to separate the signals  Table IV for run settings and labels. and glitches while reliably estimating the parameters of the former. The "CBCOnly" analysis returns mostly unbiased results, with the exception of isolated cases. Additionally, we find M(g G |g CBC+G ) < 0.01 in the first column which indicates that the glitch reconstruction is reliable even in the presence of a CBC signal. Due to the adjusted settings required for this analysis, automating analyses of fast-scattering glitches will be challenging.

D. Tomte Glitches
Tomte glitches are similar to blip glitches in that they can resemble CBCs with high, unequal masses and high spins [93]. We again initially considered various instances of tomte glitches in LIGO data. However, we found the various tomte glitches to be morphologically very similar to each other, and therefore here restrict to a single instance. The glitch spectrogram is again given in Fig. 3 and results are presented in Fig. 9. We find broadly similar results to the blip glitch case.
Similar to blip glitches, the "CBCOnly" analysis leads to large biases for all CBC parameters, both intrinsic and extrinsic. This suggests that of all glitches analyzed so far, tomtes are the ones most morphologically similar to CBCs. This is also evident in the second column of Fig. 9, where O(h rec CBC |g G ) ∼ 1 and O(g G |h inj ) ∼ 0, which means that in the "CBCOnly" analysis the CBC model ignored the signal entirely in favor of the glitch. However, even in this challenging case, the joint "CBC+Glitch" analysis is able to separate the signal and the glitch and result in reliable CBC parameter estimates and glitch reconstruction.

E. Low-Frequency Blip Glitches
Low-frequency blips, as the name suggests, are similar in morphology to blip glitches except that they infect lower frequency bands, see the spectrogram in Fig. 3. Given their similarity to blips, we consider a single instance of a low-frequency blip glitch and show results in Fig. 10. We obtain similar results to the blip glitch case, Fig. 5, with small mismatches M(g G |g CBC+G ) ≤ 0.01 and recovered posteriors that are consistent with the injected values. However, low-frequency blips do not cause as significant a bias on the recovered CBC parameters when the glitch is not included in the model. This might be due to their low-frequency nature, which means that they do not significantly overlap with the CBCs in the most sensitive detector frequency band.

F. Whistle
Whistle glitches are fairly loud glitches with a characteristic morphology depicted in the spectrogram of Fig. 3. Our chosen instance of this glitch has an SNR of ∼ 275.
Given their strength, more than 200 wavelets are required to model them accurately, which poses a considerable challenge for sampler convergence. Given this fact, we only attempted injections on short duration segments.
To aid convergence, we use GlitchBuster to initialize the glitch model, see Sec. II. Despite the short duration, the high frequency of the glitch results in a lot of waveform cycles, we therefore also increase the maximum quality factor of the wavelets. Finally, we also increase the number of iterations within the wavelet RJMCMC (see Fig. 2) from 10 2 to 10 3 as this glitch reconstruction requires upwards of 230 wavelets at every posterior sample. By default, we retain one out of 100 samples, and only update (i.e., add/remove/change) a single glitch wavelet at each sampler step. These default settings would therefore not lead to independent samples as not all wavelets have a chance to be updated before a posterior sample is retained. Details about the run settings are provided in Table IV.
Results are presented in Fig. 11 where we find that we are able to subtract the glitch consistently as well as estimate the CBC parameters. Despite the strength of the glitch, the "CBCOnly" analysis returns mostly unbiased parameter posteriors, possibly due to the fact that whistle glitches are not morphologically similar to high-mass BBHs. The quality of glitch modeling and subtraction is further explored in Fig. 12 for the analysis of Fig. 11, specifically the injection at 100 ms relative to the glitch. Comparison between the middle and right panel shows that we can efficiently subtract the glitch power. In this case, we also find that the median glitch reconstruction (top panel) results in an oversubtraction of the glitch. The fair draw (bottom panel) leads to data that look more consistent with Gaussian noise. For this and the reasons discussed in Sec. IV B we generally prefer working with fair posterior draws rather than median glitch reconstructions when making glitch-subtracted data.

G. Jensen-Shannon Divergence
As a final test, we compute a simple summary statistic for the differences between the "CBCOnly" and the "CBC+Glitch" posteriors: the Jensen-Shannon (JS) divergence. The JS divergence describes the similarity of two distributions with JS = 0 for identical distributions and JS = 1 for disparate distributions. We plot the median and maximum/minimum JS for D L , χ eff , and M c across CBC injections on the same glitch at different times in Fig. 13. With the exception of the tomte glitch (discussed in Sec. IV D) where the posteriors are completely different, we recover the general trend that the JS divergence is smaller for low masses, which implies that the "CBCOnly" and "CBC+Glitch" posteriors are more similar for lower-mass events. This again supports the previous conclusion that though glitches are more likely to overlap long duration events, glitch subtraction is more important for high-mass than low-mass signals. Among  Table IV for run settings and labels. The recovered CBC posteriors without simultaneous glitch modeling (light-blue violin plots) are heavily biased.  Table IV for run settings and labels.
the different parameters, the chirp mass is the one with the lowest JS on average, which is expected given the fact that it is the best measured intrinsic source parameter.

V. FURTHER VALIDATION STUDIES
The results of Sec. IV show that the full "CBC+Glitch" analysis can separate signals and glitches. In this section we provide some further validation tests regarding robustness against waveform systematics in the case of injections including higherorder modes or spin precession. We also assess the performance for signals that are observed by a single detector. Run settings and injection parameters for this section are presented in Tables V and VI respectively.

A. Waveform Systematics: Higher-Order Modes
The results of Sec. IV are based on IMRPhenomD [76,77], a waveform model that does not include higher-order    FIG. 13. Jensen-Shannon (JS) divergence between the "CBC only" and the "CBC+Glitch" marginalized one-dimensional posteriors for DL, χ eff , and Mc. We then plot the median (marker) and minimum to maximum values (error bars) over CBC injections at different times with respect to the same glitch. The breaks in the y axis of the plot indicate different mass ranges, increasing upwards from the origin. The general trend is that JS increases with the signal mass, again suggesting that glitches affect high-mass systems more.
Signal Injected Waveform Varied param. min. max. HOM IMRPhenomHM cos(ι) -1 1 SPIN IMRPhenomPv2 χp 0.23 0.60 TABLE VI. Parameters of the injected signals for the tests of higher-order modes (HOM) and spin precession (SPIN). For high-order modes (spin precession) we vary cos(ι) (χp) between a minimum and a maximum value (third and fourth columns) to modify the strength of the deviation between the IMRPhenomD recovery waveform and the injected waveform.
modes, i.e., power from spherical harmonics beyond the dominant l = |m| = 2 mode. Such modes change the waveform morphology and thus neglecting them will lead to biases, especially for high SNR, unequal-mass systems, observed "edge-on" (cos ι = 0) [96][97][98][99][100][101] 4 . Our current CBC sampler can work with waveform models that include higher-order modes such as IMRPhenomHM [102], however, we do not perform such runs here because we lack an implementation of the heterodyne procedure [73,78] that speeds up the likelihood calculation for such waveforms. Such an extension was described in [103] and we plan to implement it in the future. Because a real signal will inevitably contain some amount of higher-order modes, recovery with a waveform that neglects them could induce a systematic error in parameter extraction. Perhaps even more worrisome would the possibility that the glitch model subsumes some of the higher-order mode power which is then inadvertently subtracted from the data together with the glitch. We check for both effects by injecting the high-mass BBH signal from Table I using IMRPhenomHM [102] with varying inclination cos ι ∈ [−1, 1] onto one of our blip glitches and recover them again with IMRPhenomD. Figure 14 shows recovered parameters for different system inclinations and Fig. 15 shows the recovered CBC and glitch reconstruction for the case with cos ι = 0.
Compared to the top row of Fig. 5, i.e., the same injected CBC and glitch but without higher-order modes and different inclinations, we find that M(g G |g CBC+G ) has increased, but still remains ≤ 0.01. The other mismatches are comparable between Figs. 14 and 5. Despite the increase in M(g G |g CBC+G ), its value remains small and comparable to results from other blip glitches without higher-order modes, for example the second row of Fig. 5. We attribute this to the fact that the higher-order mode power is still too low to overcome the parsimony requirement of the glitch model to be picked up. Figure 15 further reinforces this picture, by showing that the main effect of higher-order modes is additional high-frequency power that cannot be captured by the template (compare the magenta reconstruction to the teal injection). However, this residual power is still two orders of magnitude below the glitch power in the same frequency range. Figure 14 also confirms that the amount of bias expected on the CBC parameters (notably D L ) is a function of the binary inclination. Furthermore, we find that not including the glitch in the model now leads to more pronounced parameter biases (blue violin plots). Regardless, even in the edge-on case the glitch and CBC signal can be separated sufficiently well as demonstrated by the mismatches in the first column. Higher SNR signals or, in general, a signal with more than SNR∼ 6 − 7 in the higher-order modes could have more noticeable deviations from waveforms without higher-order modes where their power could then be picked up by the glitch model. However, current events with detectable higherorder mode content are below this SNR threshold [104][105][106][107][108].

B. Waveform Systematics: Spin Precession
In physical scenarios where the component spins are misaligned with the orbital angular momentum, the binary system experiences spin-precession which modu- and spectrum (bottom) for a high-mass BBH injected with higherorder modes and edge-on. We show medians and 90% credible intervals for the CBC signal (magenta), the glitch (gold), and the noise PSD (grey/black) from the "CBC+glitch" analysis. The injected signal is given with a dashed teal line. The higher-order modes result in oscillations in the inspiral spectral amplitude as well as additional power at high frequencies.
The CBC waveform used for the reconstruction does not include higher-order modes, however the reconstructed glitch model does not recover the excess power from higher-order modes. The CBC model from a "CBCOnly" analysis (light blue) is similar to the one from the "CBC+glitch" analysis.
lates the observed waveform [109]. The current implementation of the CBC sampler used here only accounts for spins aligned with the orbital angular momentum, so we assumed that the injected waveforms in our main analysis were also nonprecessing, motivated also by the lack of strong precession effects in event catalogues [20,110]. However, signals with large in-plane spins, unequal masses, and/or observed edge-on could exhibit strong precessional effects [111][112][113][114][115][116][117][118][119][120][121]. The CBC sampler will be extended to include misaligned spin degrees of freedom in the future.
Similar to our analysis of the impact of higher-order modes, we study the impact of using nonprecessing templates by performing injections of high-mass BBHs with misaligned spins on the tomte glitch. The tomte glitch family was selected for this study as it is similar morphologically to highly-spinning, massive BBHs [93], and also because it consistently leads to the largest biases in CBC parameters when mismodeled, see Fig. 9. For the injections we use IMRPhenomPv2 [122] and we recover the signals with the same IMRPhenomD waveform as before. The degree of spin-precession in a signal is commonly characterized by the parameter χ p which is proportional to a mass-weighted maximum (over the two compact objects) of the in-plane spin magnitude [123]. Its range is χ p ∈ [0, 1], where χ p = 0 describes a system with aligned spins (no precession) and χ p = 1 is a maximally precessing system. Figure 16 shows results for different values of χ p and Fig. 17 shows the recovered CBC and glitch reconstructions for the case the largest χ p .
Compared to Fig. 9 that shows results with the same signal and glitch but with a nonprecessing injection, we find that parameter recovery is increasingly biased as the injected χ p increases. Again the "CBCOnly" analysis (blue violin plots) leads to overwhelming parameter biases, which seems to be a generic feature of tomte glitches. This is again reflected in the fact that O(h rec CBC |g G ) ∼ 1 whereas O(g G |h inj ) ∼ 0. Importantly, however, M(g G |g CBC+G ) ≤ 0.01 and is similar to the  Fig. 15 for a spin-precessing high-mass BBH signal with χp = 0.592. Spin precession induces oscillations in the inspiral spectral amplitude, most visible at around 40 Hz. The CBC waveform template used for the reconstruction assumes the spins are aligned with the orbital angular momentum, however the reconstructed glitch model appears unaffected. The CBC model from a "CBCOnly" analysis instead attempts to recover the glitch.
corresponding mismatches of Fig. 9; this means that the excess power due to spin-precession is not recovered by the glitch model, as its power is too low to be significant to the glitch model.
A similar conclusion is drawn from Fig. 17. The recovered signal is noticeably different from the injected signal, most prominently shown in the bottom panel where the precession-induced amplitude oscillations are absent from the posterior. However, again the glitch reconstruction appears to be unaffected, most likely because the relevant residual power is 2-3 orders of magnitude below the glitch power in the relevant frequency range. Higher binary inclinations and higher signal SNRs might make the difference between precessing and nonprecessing waveforms more stark, causing the glitch model to capture any residual power due to spin precession. However, with current signal strengths and inferred amounts of spin precession, we find this to be unlikely. Though the CBC parameters recovered are clearly biased, the glitch modeling appears to be robust. Glitch-subtracted data can therefore be constructed and further analyzed with more complete waveform models. Nonetheless, we plan to extend the analysis to include spin-precession effects in the CBC sampler.

C. Single-detector signals
The joint analysis of CBCs and glitches assumes that the astrophysical signal is coherent across the detector network, while the glitch is not. Hence, the results presented so far are based on data from both LIGO detectors. However, single-detector candidates have been reported [20], and we therefore test here if our analysis could separate them from glitches. The CBC waveform template we employ is a fairly restrictive model and we indeed find that this allows us in some cases to separate them from glitches even in single-detector data. Such a separation would be inherently impossible for the previous BayesWave analyses that distinguished between signal and glitch solely via coherence within a detector network [43].
We revisit runs B1 BBH HM and B1 BBH LM from Ta-  Table V for run settings and labels.
ble IV and analyze now only the LHO data that contain the glitch. We decrease the distance so that the singledetector SNR is 15 for consistency with all other analyses. We present parameter results in Fig. 18 and find that we are largely able to seperate the signal and the glitch. Even in this single-detector case the full "CBC+Glitch" model outperforms the simpler "CBCOnly" analysis and the glitch reconstruction is consistent with the one from data with no CBC injections. An example reconstruction plot is shown in Fig. 19 where again the CBC and glitch components of the full "CBC+Glitch" model recover their corresponding data component. The "CB-COnly" analysis, on the other hand, largely mistakes the glitch for a signal. Although these preliminary results are promising, we remain cautious of such cases. If presented with a similar scenario during an actual observing run, our analysis would require additional case-by-case attention and testing.

VI. CLASSIFYING TRIGGERS
The joint "CBC+Glitch" analysis simultaneously models CBC signals and glitches, however, the priors for both the CBC and the glitch model allow for the possibility of no CBC and/or no glitch in the data. In the glitch case, this is straightforward, as the model allows for 0 wavelets in all detectors, as was for example recovered in the case of GW150914 in [40]. While the CBC pri- ors ensure the presence of a CBC in the data, the luminosity distance prior extends to 10Gpc, which effectively corresponds to a signal with negligibly small SNR. This suggests that our analysis can be used to assess whether certain detected excess power consists of a CBC signal, a glitch, both, or neither. We revisit two low-significance candidates from [21,45] and analyze them with the joint "CBC+Glitch" and the "CBCOnly" analysis. Analysis settings are provided in Table VII.

S200114
We begin with 200114 020818 (referred to as S200114 from now on), which was identified by Coherent Wave Burst [124] with a false-alarm rate of 0.058 yr −1 [45]. Despite the low false alarm rate, the conclusion was that although an astrophysical origin could not be excluded,  Fig. 15 for the B1 HM BBH SING injection at 100ms from the glitch. When using the full "CBC+Glitch" model, h rec CBC+G and gCBC+G recover the CBC signal and glitch respectively even when data from a single detector only are available. The CBC model from a "CBCOnly" analysis instead attempts to recover the glitch.
the trigger is of likely glitch origin since the estimated CBC parameters depended heavily on the choice of waveform model. Figure 20 shows the time-domain reconstruction in each detector and Fig. 21 shows a few marginalized parameter posteriors. The joint "CBC+Glitch" analysis is consistent with the presence of both a CBC signal (at the 90% credible level) and a low-frequency glitch in LLO. The morphology of the LLO glitch is consistent with a low-frequency blip or a tomte. This suggests that while there is excess power that is morphologically similar to a CBC and is coherent between LHO and LLO, there is also additional incoherent power in LLO that is captured by the glitch model. The CBC reconstruction is consistent with zero in Virgo at the 90% credible level. Additionally, the "CBCOnly" analysis finds a broadly consistent CBC reconstruction, with small differences at low frequencies. Differences between the "CBCOnly" and "CBC+glitch" CBC parameter posteriors are much smaller than the waveform systematics reported in [45].
The presence of some amount of coherent power is consistent with the low false-alarm rate reported by a coherent detection pipeline in [45]. The additional incoherent power in LLO could explain the inconsistent parameter estimation results between different waveform models in each detector. We show medians and 90% credible intervals for h rec CBC+G (magenta), gCBC+G (gold), and h rec CBC (blue). The data are consistent with the presence of both a CBC signal (i.e., coherent power that is morphologically similar to a CBC), and a glitch (i.e., additional incoherent power in LLO). Posteriors for select parameters for S200114 from the "CBC+Glitch" analysis (pink) to the "CBCOnly" analysis (blue). Differences between these posteriors are smaller than those between different waveform models reported in [45].
presented in [45], especially since the NRSur7dq4 [125], SEOBNRv4PHM [121], and IMRPhenomXPHM [119] waveform models used in [45] include the effects of spin-precession and higher-order modes. Thus they can account for more complicated morphologies [126,127] in the CBC signal than the PhenomD model used here. The glitch in LLO is similar to a tomte which we have found to be recovered as a CBC in "CBCOnly" analyses. To check whether a single tomte glitch could trick the "CBC+glitch" analysis into concluding that both a CBC and a glitch are present in the data, we revisit the tomte glitch from Table IV, perform no CBC injection, and carry out a "CBC+Glitch" analysis. Reassuringly, the sampler indeed converges to the correct answer as shown in Fig. 22, namely that no CBC is present in the data, rather only a glitch. This suggests that the results of Fig. 20 cannot be the outcome of a single tomte glitch and the "CBC+glitch" analysis does not recover coherent power when there is none.
However, we do not attempt to obtain a background estimate and thus cannot assess the probability that such a combination of coherent/incoherent power has a terrestrial origin. Such a full background estimate could result in the calculation of a false alarm rate similarly to the BayesWave analysis in [50,128] based on the signal and glitch models only. In our case, we would use the "CBC+glitch" analysis to estimate Bayes Factors for the various models of interest and compare them to similar results obtained from data that have been shifted in time between LHO and LLO.

S191225
We then consider 191225 215715 (labeled S191225 from now on), a low-significance LLO-Virgo trigger found by the PyCBC Live [129] and the PyCBC-IMBH [130] searches with false-alarm rates of 0.4 yr −1 and 0.47 yr −1 respectively in O3 [20,45]. This candidate was ultimately deemed a glitch due to similar detector behavior surrounding the event. The reconstructed time-domain signal is shown in Fig. 23 for each detector, where we find that the data are consistent with a very high mass CBC and a glitch.
When contrasting parameters from the "CBCOnly" analysis to the "CBC+Glitch" analysis. Fig. 24, we find that the former displays the telltale signs of a glitch; negative χ eff and unequal masses [93]. The latter still recovers some coherent power with recovered parameters instead pointing to a much higher mass binary (total mass > 300M ). Since the new recovered mass is much larger than the original one, we might expect the false-alarm rate for this event to increase, although a full background estimate is outside the scope of this study.
Overall, we find that S191225 is consistent with a glitch in LLO (possibly a tomte) atop of some coherent power, which agrees and expands upon with the conclusions of [20,45]. Whitened time-domain reconstruction for an analysis of data that contain only a single tomte glitch. We show medians and 90% credible intervals for h rec CBC+G (magenta), gCBC+G (gold), and h rec CBC (blue). The data are consistent with the presence of solely a tomte glitch in LLO and do not recover any coherent power (the magenta CBC reconstruction is consistent with zero) when there is none.

VII. CONCLUSIONS
The various models that form the BayesWave algorithm allow us to analyze GW data that include multiple components, specifically noise PSD, glitches, and a CBC signal. We present multiple tests with injected signals that overlap with real LIGO glitches and show that we can reliably separate signals and glitches, estimate the CBC parameters, and provide estimates for the glitch to be subtracted from the data. Runtime estimates for the various analyses are presented in App. A.
Our analysis is able to identify all glitches analyzed. It is particularly reliable for short-duration glitches such as blips, tomtes, and low-frequency blips which could be tackled in an automated way with default analysis settings. These are also the glitch types that cause the largest biases for CBC parameters when left unaccounted for. Long-duration glitches such as fast and slow scattering glitches are more challenging and in some cases need specialized settings. Crucially, the necessary settings (such as maximum wavelet quality factor or analysis segment length) differ even between glitches of the same family, suggesting that automation is not yet feasible. This effect is particularly prominent for fast-scattering glitches that create long periods of nonstationarity that challenge PSD estimation, especially at low frequencies in LLO (top) and Virgo (bottom). We show the median and 90% credible intervals for h rec CBC+G (magenta), gCBC+G (gold), and h rec CBC (blue). Some low-frequency coherent power is recovered by h rec CBC+G whereas the high-frequency power is largely incoherent and is recovered by the glitch model. Posteriors for select parameters for S191225 from the "CBC+Glitch" analysis (pink) to the "CBCOnly" analysis (blue).
(below 40 Hz). Luckily, these glitches incur smaller biases in CBC parameter, most likely for the same reason they are difficult to model: they have a large time-frequency footprint that does not resemble a CBC chirp. In this study, (with one exception) changes on glitch wavelet priors concerned their ranges, while their shapes were unaltered compared to default settings. However, dedicated glitch priors that target particularly problematic glitch families would improve glitch modeling. For example, a prior that favors wavelets with lower am- plitude can lead to improved results for the low-SNR slow-scattering glitch (S2). Further examples of dedicated priors include a trained model based on a principal component analysis for tomte glitches [80] or prior information about the frequency spacing of the scattering arches [131,132]. Depending on the characteristics of the most prevalent and problematic glitch types in O4, we plan to explore such dedicated priors in the future.
vice of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes. This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by NSF Grants PHY-0757058 and PHY-0823459. S.H. and K.C. were supported by NSF Grant No. PHY-2110111. D.D. was supported by the National Science Foundation as part of the LIGO laboratory T.L. was supported by funding from the NASA LISA study office. N.J.C. was supported by NSF Grant No. PHY-1912053. Software: gwpy [133], matplotlib [134].

Appendix A: Runtimes
Including a variable-dimensional glitch and PSD model comes at an additional computational cost compared to standard CBC analyses. We display the runtimes for the different glitch types and analyses presented in this paper in Fig. 25. There is an essentially linear dependence between the number of data points (the segment length times the sampling rate) and the run time. The whistle glitch runs are outliers as additional settings were required to subtract such high-SNR glitches (see Sec. IV F). Since runtime is (approximately) linear with the number of parallel chains and the number of iterations, we rescale time by the default settings of 20 chains and 4 × 10 6 iterations, though lighter settings can be used for expedited results. For example, the BNS runs for Fast-Scattering and Blip glitches used 10 chains for speed, and the whistle glitch was run with 25 chains for convergence. The rest were run with 20 chains. These estimates are further based on single-core runs and will be sped up accordingly by ongoing work to parallelize the chains.