Modeling compact binary signals and instrumental glitches in gravitational wave data
Katerina Chatziioannou, Neil Cornish, Marcella Wijngaarden, Tyson B. Littenberg
MModeling compact binary signals and instrumental glitches in gravitational wave data
Katerina Chatziioannou,
1, 2, 3
Neil J. Cornish, Marcella Wijngaarden,
3, 5 and Tyson B. Littenberg Department of Physics, California Institute of Technology, Pasadena, California 91125, USA LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA eXtreme Gravity Institute, Department of Physics, Montana State University, Bozeman, Montana 59717, USA Mathematical Sciences and STAG Research Centre, University of Southampton, SO17 1BJ, Southampton, UK NASA Marshall Space Flight Center, Huntsville, AL 35812, USA (Dated: February 9, 2021)Transient non-Gaussian noise in gravitational wave detectors, commonly referred to as glitches,pose challenges for detection and inference of the astrophysical properties of detected signals whenthe two are coincident in time. Current analyses aim toward modeling and subtracting the glitchesfrom the data using a flexible, morphology-independent model in terms of sine-Gaussian waveletsbefore the signal source properties are inferred using templates for the compact binary signal. Wepresent a new analysis of gravitational wave data that contain both a signal and glitches by simul-taneously modeling the compact binary signal in terms of templates and the instrumental glitchesusing sine-Gaussian wavelets. The model for the glitches is generic and can thus be applied to awide range of glitch morphologies without any special tuning. The simultaneous modeling of theastrophysical signal with templates allows us to efficiently separate the signal from the glitches, aswe demonstrate using simulated signals injected around real O2 glitches in the two LIGO detectors.We show that our new proposed analysis can separate overlapping glitches and signals, estimatethe compact binary parameters, and provide ready-to-use glitch-subtracted data for downstreaminference analyses.
I. INTRODUCTION
During the first half of their third observing run (O3a),the advanced ground based gravitational wave (GW) de-tectors LIGO [1] and Virgo [2] observed an astrophysi-cal transient signal about every 5 days of data [3]. Thelarge detection rate increases the chance of observing anevent while one of the detectors experiences transient nonGaussian noise, also known as instrumental glitches. In-deed, this scenario has come to pass for one event fromthe second observing run (O2) [4] and 8 events from thefirst half of the third observing run [3].Such coincidences are expected to become even morefrequent in the coming years. Planned improvements inthe detectors’ sensitivity will be directly reflected by aneven larger rate of astrophysical discoveries [5]. More-over, O3a was characterized by an increase in the rateof glitch occurrence in the two LIGO detectors, a trendthat might persist during the fourth observing run (O4)as the decreased average detector noise might help revealweaker sources of transient noise. For example the rate ofglitches in the LIGO Livingston detector increased from0.2 per minute in O2 to 0.8 per minute in O3a [3].The presence of a non Gaussian noise feature in thedata, a glitch, poses challenges for nearly all inferenceanalyses. GW inference is based on a model for the de-tector noise, expressed through the likelihood function.In the absence of glitches, detector noise is colored andGaussian to a very good approximation [6], with a spec-trum that is described through the noise power spectraldensity (PSD). The above considerations give rise to aGaussian likelihood function whose variance is the noisePSD, a choice that is almost ubiquitous [7, 8]. Different choices for estimating the PSD or treating its uncertaintycan result in different functional forms for the likelihood,but they are all based on the assumption of colored Gaus-sian noise [9, 10].Since instrumental glitches violate the basic assump-tions of GW inference, they need to be effectively miti-gated before the data are analyzed. One option is to re-move the offending data all together [4, 11–13], which canbe done quickly, allowing for low latency estimation ofsource parameters that enable followup observations [4].The downside of this approach is that part of the astro-physical signal is lost making it prohibitive for binaryblack hole (BBH) signals whose duration is comparableto the glitch duration. In order to avoid signal, and thusinformation loss, another option is to model the glitchand regress it from the data, leaving behind not only theastrophysical signal but also the Gaussian noise. Thisapproach is the topic of the current study .The wide variety of glitch morphologies, and even vari-ations within a certain glitch type, make constructing ex-act models for glitches challenging [21]. A more flexibleapproach is based on BayesWave [22, 23] which modelsvarious components of the GW data in a morphology-independent way. Non-Gaussian features in the dataare modeled in terms of sums of sine-Gaussian waveletswhose number and parameters are marginalized over with An independent effort to mitigate the effect of broadband and/ornonstationary detector noise is based on information from auxil-iary sensors [14–20]. This approach does not remove entire datasegments either and thus is not expected to lead to loss of infor-mation. a r X i v : . [ g r- q c ] F e b a suite of Markov Chain Monte Carlo (MCMC) and Re-versible Jump MCMC (RJMCMC) [24] samplers. Coher-ent features (i.e. features that appear in all detectors in amanner consistent with an astrophysical signal originat-ing from a specific sky location) are modeled by a singlesum of wavelets that is projected onto the detector net-work; these features are interpreted as having an astro-physical origin. Incoherent features are instead modeledby independent sums of wavelets in each GW detectorand are interpreted as instrumental glitches. The PSDof the Gaussian noise is also modeled in terms of splinesand Lorentzians using an algorithm sometimes known as BayesLine [6, 25].
BayesWave and
BayesLine are fullyintegrated and we will refer to the combined analysis withthe name
BayesWave in this paper.Modeling instrumental glitches with
BayesWave andsubtracting them from the data in order to make ready-to-use data for downstream inference has been a stan-dard step of LIGO/Virgo analyses since O2 [3, 4]. TheGW signal from the first binary neutron star (BNS)coalescence detection, GW170817, overlapped with aglitch in the LIGO Livingston detector approximately1 .
1s before coalescence [4]. The glitch was modeled with
BayesWave ’s glitch model in terms of a sum of waveletsand removed from the data, a procedure documented andreleased in [26]. Despite the glitch overlapping with theactual astrophysical signal, the subtraction process wasrobust against inadvertently removing the signal togetherwith the glitch. The reason is that the specific glitch wasshort in duration (less than a second) and extended infrequency, unlike the signal that lasted for about 2 min-utes in the detector sensitive band. As such, the sine-Gaussian wavelets that would fit the glitch and the signalare distinct in terms of their time-frequency features; thewavelets that model the glitch are short and hence donot model the long-lasting BNS signal. This procedurewas further shown to not introduce biases in the astro-physical parameter inference of the underlying signal byanalyzing simulated signals injected on instances of thesame glitch type in LIGO Livingston data [27].Motivated by the success of this first attempt at glitchmitigation and in preparation for the increased detectionrate of O3, BayesWave was extended to be able to si-multaneously model both the signal and the glitch [23].Both signals and glitches are modeled with a sum of sine-Gaussian wavelets, the only difference being that the sig-nal is coherent across the detectors in the network, whilethe glitch is not. The analysis effectively uses data fromall detectors available to determine which part of the nonGaussian data are coherent (and would thus correspondto an astrophysical signal), and which part is incoherent(and would thus correspond to an instrumental glitch).The combined signal+glitch analysis was applied to oneO3a detection [3], enabling glitch mitigation even for datathat contained short-duration BBH signals.The signal+glitch analysis models compact binary co-alescence (CBC) signals in terms of wavelets, and is thusagnostic to the signal morphology. However, accurate models exist for CBCs in terms of solutions to the Ein-stein field equations that are routinely used both for de-tection and parameter estimation. In this paper we takeanother step toward efficient separation of CBCs andglitches by constructing an analysis that simultaneouslymodels the CBC signal in terms of CBC templates andthe glitch in terms of sine-Gaussian wavelets. Similar tothe initial glitch-only analysis and the subsequent sig-nal+glitch analysis, we also model and marginalize overthe detector noise PSD. We test our analysis using publicO2 data that contain common glitch types and simulatedCBC signals. We demonstrate that we can efficiently sep-arate the glitch from the CBC, estimate the CBC param-eters, and provide ready-to-use glitch-subtracted data fordownstream inference analyses.The rest of the paper is organized as follows. In Sec. IIwe describe the updates to the standard
BayesWave algo-rithm in terms of the CBC analysis. In Sec. III we applyour analysis to simulated signal overlapping with knowndetector glitches from O2 data. In Sec. IV we analyzea selection of detected signals, namely GW170817 andGW150914. Finally, in Sec. V we conclude and point tofuture work.
II. GENERAL ALGORITHM DESCRIPTION
The combined
BayesWave algorithm is presented in de-tail in [23] and here we describe only the features relevantto our study.
BayesWave simultaneously models signals,glitches, and Gaussian noise in GW data by means ofdifferent models. The signal model describes astrophysi-cal signals through a sum of Morlet Gabor wavelets thatare coherent across the detector network. The numberof wavelets and the parameters of each are marginalizedover, as are the extrinsic parameters that determine howthe signal is projected in each detector. The glitch model describes instrumental glitches with an incoherent sumof Morlet Gabor wavelets whose number and parametersare again marginalized over. Glitch power in each detec-tor is described by an independent sum of such wavelets.The noise model describes the Gaussian noise PSD witha broadband spline model and sharp Lorentzians. Asabove, the number of spline points and Lorentzians aswell as their parameters are marginalized over.In order to sample the multidimensional posterior den-sity of all models,
BayesWave uses a blocked Gibbs sam-pler that takes turns between sampling each model withcompletely independent MCMC or RJMCMC samplers.This includes (i) an RJMCMC that samples the signaland glitch wavelet parameters, (ii) an MCMC that sam-ples the signal extrinsic parameters, and (iii) an RJM-CMC that samples the splines and Lorentzians for thenoise PSD. Each sampler in turn updates its parame-ters for a predetermined number of iterations, typically O (10 ), while all other parameters are kept fixed. For ex-ample, the extrinsic sampler updates the extrinsic signalparameters while the wavelet parameters and noise PSDare kept constant. Once the predetermined number ofupdates has been reached, the extrinsic sampler returnsits current parameters and the noise sampler begins up-dating the noise model while keeping the wavelet andextrinsic parameters fixed. This process of alternatingsampling between different blocks of model parametersis repeated for O (10 ) iterations.The construction of the algorithm in terms of a blockedGibbs sampler makes adding further models and sam-plers straightforward. In the current version describedin [23], the astrophysical signal is modeled with coherentsine-Gaussian wavelets that allow us to describe signalswith a large level of flexibility. We extend BayesWave ’sblocked Gibbs sampler by adding one more element,namely a model of the signal in terms of quasicircularCBC waveforms. In fashion with the existing implemen-tation, the MCMC that samples the posterior distribu-tion for the CBC parameters is completely independentfrom the remaining code samplers. The result is a flexi-ble algorithm that can be used with any combination ofCBC, signal , glitch, and noise models for the detectordata.The CBC model is integrated with LALSimulation [28]and can operate with any nonprecessing model availablethere . The eleven parameters of a spin-aligned quasicir-cular CBC signal, namely the four intrinsic parameters(the two masses and spin magnitudes) and seven extrin-sic parameters (the time of coalescence, the phase of co-alescence, two sky location angles, the polarization angleand the inclination angle the distance), are updated inoverlapping blocks. Common to both blocks is the phaseof coalescence since BayesWave ’s extrinsic sampler up-dates the overall phase of the signal as described in [23].The CBC MCMC sampler updates the four intrinsic pa-rameters, the time of coalescence, the phase of coales-cence, and the distance. The existing extrinsic samplerin
BayesWave updates the two sky angles, the polariza-tion angle, the inclination angle, and the phase of co-alescence while holding all other parameters fixed. Weuse standard priors for all parameters: uniform over thedetector-frame masses and spin magnitudes, uniform intime and phase, and uniform in luminosity volume.The CBC sampler is custom and not based on anyexisting samplers used in LIGO-Virgo parameter estima-tion. The CBC sampler is taken from the recently devel-oped
QuickCBC [29] analysis pipeline. A closely related We retain the original model names in
BayesWave , hence the signal model refers to the wavelet signal model, while the
CBCmodel refers to the model in terms of CBC templates. Bothmodels target astrophysical signals. Since we do not use the signal model in the remainder of the paper, we trust that thiswill not lead to confusion. Both the sampling and the jump proposals for the CBC param-eters are constructed to expect the signal amplitude and phasefrom the waveform generator. There is therefore no fundamentallimitation to non-precessing signals and we plan to extend ouranalysis to include the effect of spin-precession in the future. sampler [30] has been developed for analyzing data fromthe future Laser Interferometer Space Antenna. TheCBC sampler is a replica exchange (parallel tempered)Markov Chain Monte Carlo (PTMCMC) algorithm thatuses a mixture of proposal distributions. The defaultcollection of proposals are: Gaussian jumps along eigen-vectors of the Fisher information matrix, scaled by thereciprocal of the square root of the corresponding eigen-value; differential evolution using a rolling history arrayat each temperature, updated every 10 iterations andholding 1000 past samples; and small, Gaussian jumpsalong each parameter direction. Each chain carries itsown Fisher information matrix, which is updated peri-odically. The Fisher and differential evolution proposalsare effective at exploring parameter correlations, whilethe small jumps prevent the chains from getting stuck inregions where the Fisher matrix becomes ill-conditioned.The CBC sampler is not optimized for blindly find-ing signals, so it is best to initialize the sampler witha good starting solution for the source parameters suchas the output from a CBC search pipeline, or the in-jected parameters for a simulated signal. Alternativelythe sampler can be initialized using a custom built CBCsearch algorithm from the
QuickCBC [29] analysis pipelinethat has been incorporated into the
BayesWave prepro-cessing steps. The search is broken into two stages, arapid network-coherent search with analytic maximiza-tion over extrinsic parameters, followed by a fast MCMCover the extrinsic parameters using a likelihood functionthat precomputes the waveform inner products [31]. Thisprocedure returns the starting point for all 11 CBC pa-rameters. More details about the initial search step anddiscussion of its robustness against instrumental glitchesare presented in [29].
III. SIMULATED SIGNALS
We test the efficacy of separating CBCs from glitcheswith our CBC+glitch model by selecting 3 common glitchtypes from O2 data [32, 33] that are known to have anadverse effect on searches for CBCs [3]. We then add sim-ulated CBC signals consistent with a BBH with detector-frame masses of 36 M (cid:12) and 29 M (cid:12) and vanishing spin atdifferent times with respect to the glitch. All simulatedsignals have a signal-to-noise (SNR) ratio of 15. We usethe IMRPhenomD [34, 35] waveform model both for simula-tion and recovery as implemented in
LALSimulation [28].We then analyze the data from the two LIGO detectorswith our CBC+glitch+noise model, where the coherentsignal is modeled by the CBC template, the glitch is mod-eled by incoherent wavelets, and the noise PSD is mod-eled with splines and Lorenzians. Spectrograms for the 3glitches are shown in Fig. 1: blip glitch (left), scatteredlight (middle), and blue mountain (right). Further de-tails and run settings for each type of glitch are shown inTable I. -0.04 0 0.04 0.08 0.12 0.16 0.2 0.24 0.28 0.32 t (s) from GPS 1168989748 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1172917776 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1165069533 f ( H z ) N o r m a li ze d e n e r g y FIG. 1. Spectrograms for the three glitches of different types studied here: blip glitch (left), scattered light (middle), bluemountain (right). The three types of glitches are characterized by very different time-frequency properties.Glitch GPS time (s) Detector Segmentlength (s) Samplingrate (Hz) f low (Hz) Q max CBC SNRBlip 1168989748 Hanford 4 2048 16 40 15Scattered light 1172917779 Livingston 8 2048 8 160 15Blue mountain 1165069536 Hanford 16 2048 16 40 15TABLE I. Settings for the runs of Sec. III. From left to right, columns correspond to the type of glitch, the GPS time, theaffected detector, the segment length, the sampling rate, to low frequency cut off, the maximum quality factor of the glitchwavelets, and the SNR of the injected signals.
A. Glitch type 1: Blip
Blip glitches are one of the most common glitch typesfor the two LIGO detectors. They are characterized byshort duration, and hence pose a challenge for the de-tection of high mass BBH signals [36]. Their origin islargely unknown. Figures 2-5 show our results for sim-ulated signals injected at different times with respect toa blip glitch in the LIGO Hanford detector during O2.Details about the glitch, including its GPS time, and therun settings are presented in Table I. A spectrogram ofthe data containing the glitch is given on the left panelof Fig. 1, where the short duration and large frequencyextent are shown.The whitened data and reconstructions for the CBCsignal and the glitch are shown in Fig. 2 where we plot the90% credible intervals for each reconstruction in LIGOHanford (top) and LIGO Livingston (bottom). Theglitch is easily visible in LIGO Hanford as a short dura-tion ∼ σ noise excursion. No glitch power is identifiedin LIGO Livingston at that time, but the CBC signal isclearly identified. This allows us to separate the corre-sponding coherent CBC signal in LIGO Hanford from theinstrumental glitch, even when the latter overlaps withthe merger phase of the signal (left panel). The glitch re-construction is also consistent across the three simulatedsignals, suggesting that the glitch model is not fitting anypart of the CBC signal.Source parameters for the simulated CBC are pre-sented in Figs. 3 and 4 both for the CBC+glitch+noiseanalysis and a CBC+noise analysis for selected recov-ered parameters for the leftmost simulated CBC sig-nal together with the injected values with black crossesor vertical lines as appropriate. Figure 3 shows themass ratio q , the effective spin χ eff , and the detector frame chirp mass M , while Fig. 4 shows the luminos-ity distance and the cosine of the inclination angle. Inall cases the posterior distributions recovered under theCBC+glitch+noise model are consistent with the in-jected parameters, though the marginalized posteriors donot peak at the injected values, as expected from infer-ence of signals in Gaussian noise. For reference, we showposteriors under the CBC+noise model in orange thatassumes that the data are consisted of just a CBC signaland Gaussian noise, without any provision for a glitch.Since this assumption is violated by the presence of theblip glitch, the resulting posteriors are expected to bebiased compared to the true parameters and the orangecontours in Figs. 3 and 4 quantify this bias. We find thatthe extrinsic parameters that are primarily determinedby the signal amplitude are more biased than the intrin-sic ones that are measured through the GW phase, asalso discussed in [37].The separation of the CBC signal from the glitchdemonstrated in Fig. 2 can be used to produce ready-to-use deglitched data for downstream inference analyses, aswas done in [3]. An estimate of the glitch reconstruction(the median or a fair draw from the glitch model poste-rior) is subtracted from the data to produce strain datathat contain only the CBC signal and Gaussian noise.The result of the glitch subtraction is shown in the spec-trograms of Fig. 5 that show the LIGO Hanford databefore (left) and after (middle) the subtraction of a fairdraw glitch reconstruction for the leftmost injection ofFig. 2. The left panel includes both the chirping signaland the blip glitch, while only the former is visible in themiddle panel. The right panel shows the data after a fairdraw from both the CBC and the glitch models has beensubtracted, resulting in residual Gaussian noise only. − − − σ n o i s e LIGO Hanford − − −
505 DataCBCGlitch − − − − . − . − . − .
02 0 .
00 0 .
02 0 .
04 0 . t (s) from GPS 1168989748.125 − − σ n o i s e LIGO Livingston − . − . − . − .
02 0 .
00 0 .
02 0 .
04 0 . t (s) from GPS 1168989748.125 − − − . − . − . − .
02 0 .
00 0 .
02 0 .
04 0 . t (s) from GPS 1168989748.125 − − FIG. 2. Credible intervals for the glitch (orange) and the CBC (blue) signal reconstruction for data containing a blip glitchin LIGO Hanford and a simulated CBC signal at 3 different times with respect to the glitch (left to right). Shaded regionscorrespond to 90% credible intervals for the whitened reconstruction, while in grey dashed lines we plot the data whitened witha fair draw PSD from our noise model posterior. The top row corresponds to LIGO Hanford and the bottom row correspondsto LIGO Livingston. . . . . CBC+glitch+noiseCBC+noise . . . . M ( M (cid:12) )
24 26 28 300 . . . . . . . . q − . − . . . . χ e ff
24 26 28 30 M ( M (cid:12) ) − . − . . . . − . − . . . . χ eff FIG. 3. One- and two-dimensional posterior distributions forselected source parameters of the simulated signal from theleft panel of Fig. 2 injected on top of a LIGO Hanford blipglitch. We include the mass ratio q , the effective spin χ eff ,and the detector frame chirp mass M posteriors, while blackcrosses or black vertical lines denote the true parameters ofthe injection. Blue (orange) contours and lines correspond tothe CBC+glitch+noise (CBC+noise) run. B. Glitch type 2: Scattered light
Glitches caused by scattered light in the interferome-ter became particularly prominent during O3 [3]. Unlikethe blip glitches studied above, scattered light glitcheshave a longer temporal duration of a few seconds andare characterized by arches in a time-frequency spectro-gram [38, 39], as depicted in the middle panel of Fig. 1.We inject simulated signals on an instance of such a glitchin LIGO Livingston and analyze the data from bothLIGO detectors with our CBC+glitch+noise model. De-tails of the glitch and the run settings are given in Table I.Due to the duration of the glitch and its low frequency − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . cos( ι ) D ( M p c ) CBC+glitch+noiseCBC+noise
FIG. 4. Two-dimensional posterior distributions for the lu-minosity distance and the binary inclination of the simulatedsignal from the left panel of Fig. 2 injected on top of a LIGOHanford blip glitch. A black cross at (1,1200Mpc) denotesthe true parameters of the injection. Blue (orange) contourscorrespond to the CBC+glitch+noise (CBC+noise) run. power we extend our analysis duration and bandwidth.The longer duration helps the noise model determine thelow-frequency Gaussian noise PSD and thus separate thelow frequency part of the glitch from Gaussian noise. Wealso increase the maximum quality factor Q max of thewavelets due to the glitch’s long duration.Figure 6 shows the data and reconstructed CBC andglitch models. We zoom in around the CBC signals,though the glitch extends beyond the time range plotted.In all cases the CBC signal is separated from the glitch,aided by the presence of a coherent signal in LIGO Han-ford. The glitch reconstruction is also consistent for all 3simulated signals, as expected for runs on the same glitch.The reconstruction exhibits oscillations at around 32Hzand 16Hz, consistent with expectations from the glitchspectrogram. Figure 7 shows posterior distributions forselected source parameters for the left-most injection in t (s) from GPS 1168989747.125 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1168989747.125 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1168989747.125 f ( H z ) N o r m a li ze d e n e r g y FIG. 5. Spectrogram of the LIGO Hanford data around the time of the blip glitch for the leftmost injection from Fig. 2. Leftpanel: data containing the blip glitch and the simulated CBC signal. Middle panel: data after a fair draw from the glitchmodel has been subtracted leaving behind only the chirping CBC signal. Right panel: data after a fair draw from the glitchand CBC models has been subtracted, leaving behind only Gaussian detector noise. blue, as well as the injected parameters. In all cases therecovered parameters are consistent with their injectedvalues. In orange, we plot results from a CBC+noise runand find small biases in the source intrinsic parameters,most notably the mass ratio.Finally, Fig. 8 shows the spectrogram of the data be-fore and after various components of the model have beensubtracted from the data. The left panel corresponds todata that contain both a signal and the glitch and thusboth the signal chirp and the characteristic glitch archesare visible. In the middle panel we plot data after afair draw from the glitch model has been subtracted, re-sulting in both the high and the low frequency arches ofthe glitch having been regressed, leaving only the chirp-ing signal behind. The right panel corresponds to datawhere a fair draw from the CBC model has further beensubtracted and is consistent with Gaussian noise.
C. Glitch type 3: Blue mountain
The final type of glitch we consider is the blue moun-tain; the spectrogram of the LIGO Hanford instance ofa blue mountain glitch we consider is shown in the rightpanel of Fig. 1. The glitch has a duration of multiple sec-onds and is characterized by higher frequencies ∼ Q max , as the glitch is composed ofshort individual bursts of power, each of which is mod-eled by individual wavelets with a small quality factor.Figure 9 shows the whitened data and credible inter-vals for the whitened CBC and glitch reconstruction ineach detector for each of the injected signals. Due to thelarge glitch duration, the signals are injected sufficientlywide apart that the reconstruction plots show non over-lapping parts of the data and the glitch. The glitch re-constructions are therefore not expected to match. As expected from the glitch spectrogram, the glitch is char-acterized by a series of short high frequency bursts, eachof which is modeled by different wavelets within our glitchmodel. Figure 10 shows posterior distributions for se-lected source parameters for the left-most injection inblue, as well as the injected parameters. In all cases therecovered parameters are consistent with their injectedvalues, suggesting that the presence of the glitch doesnot incur biases on the inferred source properties if thetwo are modeled simultaneously. As before, we also plotresults from a CBC+noise run that neglects the glitchin the data in orange and again find small biases by thepresence of the glitch in the source intrinsic parameters.The glitch subtraction process is detailed in Fig. 11that again shows spectrograms of the original data con-taining both the glitch and the signal (left), data after afair draw glitch model has been subtracted (middle), anddata after both the glitch and the fair draw CBC modelhave been removed (right). As before, data from the mid-dle panel could be used for further data processing. Theright panel shows data where a model for both the glitchand the CBC have been subtracted. Even though the ma-jority of the glitch power is absent (compare the left andright panels), some small non Gaussian power might beleft behind. The reason for this is that the blue mountainglitch is manifested as individual short bursts of glitchpower, which our flexible analysis attempts to modelcompletely independently. Indeed, the glitch model forthis run uses O (70) wavelets. Each of these wavelets,needs to model sufficient non-Gaussian power in the datain order to overcome the parsimony penalty incurred byadding more parameters to the model. As such, we ex-pect that some of the weaker “bursts” of the glitch willnot be recovered. Possible ways to alleviate this are dis-cussed in Sec. V. IV. GRAVITATIONAL WAVE EVENTS
As a further demonstration of our CBC+glitch+noisemodel, we also analyze two astrophysical events,GW170817 [4] and GW150914 [40] whose data are avail-able from GWOSC [32, 33]. Though not the main focus − σ n o i s e LIGO Hanford − − Data CBC Glitch − . . . t (s) from GPS 1172917779.7 − σ n o i s e LIGO Livingston − . . . t (s) from GPS 1172917779.7 − − . . . t (s) from GPS 1172917779.7 − FIG. 6. Credible intervals for the glitch (orange) and the CBC (blue) signal reconstruction for data containing a scattered lightglitch in LIGO Livingston and a simulated CBC signal at 3 different times with respect to the glitch (left to right). Shadedregions correspond to 90% credible intervals, while in grey dashed lines we plot the data whitened with a fair draw PSD fromour noise model posterior. The top row corresponds to LIGO Hanford while the bottom row corresponds to LIGO Livingston. . . . . CBC+glitch+noiseCBC+noise . . . . M ( M (cid:12) )
24 26 28 300 . . . . . . . . q − . − . . . . χ e ff
24 26 28 30 M ( M (cid:12) ) − . − . . . . − . − . . . . χ eff FIG. 7. One- and two-dimensional posterior distribution forselected source parameters of the simulated signal from theleft panels of Fig. 6 injected on top of a LIGO Livingston scat-tered light glitch. We include the mass ratio q , the effectivespin χ eff , and the detector frame chirp mass M posteriors,while black crosses or black vertical lines denote the true pa-rameters of the injection. Blue (orange) contours and linescorrespond to the CBC+glitch+noise (CBC+noise) run. of this paper, the analysis presented below also providesan estimate of the effect of marginalizing over the noisePSD has on the inferred astrophysical parameters. Moredetails about this effect will be presented in a separatestudy. A. GW170817
Perhaps the most known instance of a GW signal over-lapping with an instrumental glitch is GW170817 [4].Inference on the GW170817 source properties is per-formed on data where the glitch in LIGO Livingstonhas been modeled with
BayesWave ’s glitch-only model and subtracted. Analysis of simulated signals suggeststhat this procedure leads to unbiased inference, whileany analysis on data that contain the glitch results inhighly biased source parameters [27]. Both versions ofthe data are publicly available, both with and withoutthe glitch [26], so we analyze them both with differentmodels. We use data from the LIGO Hanford and theLIGO Livingston detectors and analyze 64s of data from16Hz to 2048Hz using the
IMRPhenomD NRTides wave-form model that includes finite-size effects [41]. We em-ploy our CBC+glitch+noise model on the data with theglitch and the CBC+noise model on data where the glitchhas already been subtracted. For the CBC+glitch+noisecase we use
GlitchBuster [23] to provide a quick fit tothe glitch and use that as a starting point for our glitchmodel during sampling.Credible intervals for the signal and glitch reconstruc-tions are shown in Fig. 12 for each detector for ∼ ∼ σ relative to the background detector Gaus-sian noise. The signal is not visible in the LIGO Liv-ingston data given the plotting scale. Figure 13 showsselected source parameters obtained from data both withand without the glitch. We find consistent results, show-ing that our combined CBC+glitch+noise analysis canfaithfully fit the CBC signal and the glitch simultane-ously, without the need for the two step process of firstremoving the glitch and then reanalyzing the data. B. GW150914
The first GW signal directly detected by the LIGOdetectors, GW150914 [40], did not overlap with an in-strumental glitch [42]. However, since it is one of thebest studied and loudest signals, we select it as a demon-stration of our analysis on data without glitches. Our t (s) from GPS 1172917778.7 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1172917778.7 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1172917778.7 f ( H z ) N o r m a li ze d e n e r g y FIG. 8. Spectrogram of the LIGO Livingston data around the time of the scattered light glitch for the leftmost injection fromFig. 6. Left panel: data containing the scattered light glitch and the simulated CBC signal. Middle panel: data after a fairdraw from the glitch model has been subtracted leaving behind only the chirping CBC signal. Right panel: data after a fairdraw from the glitch and CBC models has been subtracted, leaving behind only Gaussian detector noise. − σ n o i s e LIGO Hanford − − t (s) from GPS 1165069535.9 − σ n o i s e LIGO Livingston -0.4 -0.2 0 0.2 0.4 t (s) from GPS 1165069536.9 − Data CBC Glitch -0.4 -0.2 0 0.2 0.4 t (s) from GPS 11650695367.9 − FIG. 9. Credible intervals for the glitch (orange) and the CBC (blue) signal reconstruction for data containing a blue mountainglitch in LIGO Hanford and a simulated CBC signal at 3 different times with respect to the glitch (left to right). Shadedregions correspond to 90% credible intervals, while in grey dashed lines we plot the data whitened with a fair draw PSD fromour noise model posterior. The top row corresponds to LIGO Hanford while the bottom row corresponds to LIGO Livingston. glitch model has the flexibility to use no glitch wavelets,we therefore expect many samples in the glitch modelposterior to contain exactly zero glitch power. We an-alyze 4s of data starting at 16Hz and with a samplingrate of 2048Hz. We perform two runs, one with theCBC+glitch+noise model and one with the CBC+glitchmodel using otherwise identical settings.Relevant results are shown in Figs. 14 and 15 where asbefore we plot the CBC and glitch reconstructions of theCBC+glitch+noise model in the two detectors and therecovered source parameters. The CBC reconstruction ofFig. 14 is consistent with previous results [43]. The glitchreconstruction is too small to identify in the scale of theplot, as we find that 86% and 14% of our posterior sam-ples had exactly zero glitch wavelets in LIGO Hanfordand LIGO Livingston respectively. Figure 15 shows theposterior distribution for selected source parameters ofGW150914 obtained under the CBC+glitch+noise andthe CBC+noise models. The two posteriors yield con-sistent results, showing that the glitch model does notaffect the CBC parameters when no glitch is present inthe data.
V. CONCLUSIONS
We construct and validate an analysis of GW datathat simultaneously models astrophysical CBC signalsand instrumental glitches. We test the analysis againstreal instances of glitches in the two LIGO detectors fromO2 data and simulated CBC signals injected at differenttimes with respect to the glitch. We find that our analy-sis can separate the two, and provide both estimates forthe CBC source parameters and glitch-subtracted datafor subsequent analyses. The glitch model we employ isa sum of sine-Gaussian wavelets that is not tuned to anyspecific glitch type and morphology; it can thus handleeven novel glitch types that might first appear duringO4. Even though this flexibility is desirable given theunpredictable and evolving nature of glitches, the effi-cacy of glitch subtraction can be improved by employingtargeted priors for different glitch types. One such exam-ple would be a prior that anticipates arches at frequencymultiples in the case of scattered light glitches. We leavesuch targeted priors to future work.Our analysis considered only simulated BBH sig-nals, though we also present an analysis of the BNSGW170817. We expect overlapping CBCs and glitchesof similar duration to be a worse-case-scenario due to . . . . CBC+glitch+noiseCBC+noise . . . . M ( M (cid:12) )
22 24 26 28 300 . . . . . . . . . q − . − . . . χ e ff
22 24 26 28 30 M ( M (cid:12) ) − . − . . . − . − . . . χ eff FIG. 10. One- and two-dimensional posterior distribution forselected source parameters of the simulated signal from theleft panels of Fig. 9 injected on top of a LIGO Hanford bluemountain glitch. We include the mass ratio q , the effectivespin χ eff , and the detector frame chirp mass M posteriors,while black crosses or black vertical lines denote the true pa-rameters of the injection. their similar morphology [44, 45]. Given that, we plan tocarry out a larger scale study of our CBC+glitch analysisthat includes more glitch types and CBC classes, such asBNSs and lower mass BBHs. Additionally, the analysispresented here did not make use of GlitchBuster [23]to provide initial fits to the glitch, apart from theGW170817 case. In the future we plan to investigateinterfacing
GlitchBuster and
BayesWave in more de-tail, in the hopes that an efficient starting point for the glitch model during sampling will decrease the sam-pler’s convergence time and result in ready-to-use glitch-subtracted data more quickly. We hope that our analysiswill contribute to robust and efficient glitch mitigationagainst the increased event rate anticipated in O4; ourgoal is to facilitate analysis of as much data as possibleand maximize the science output of the upcoming obser-vations.
ACKNOWLEDGMENTS [1] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. ,074001 (2015), arXiv:1411.4547 [gr-qc].[2] F. Acernese et al. (Virgo Collaboration), Class. Quant.Grav. , 024001 (2015), arXiv:1408.3978 [gr-qc].[3] R. Abbott et al. (LIGO Scientific, Virgo), (2020),arXiv:2010.14527 [gr-qc].[4] B. P. Abbott et al. (LIGO Scientific Collaboration, VirgoCollaboration), Phys. Rev. Lett. , 161101 (2017),arXiv:1710.05832 [gr-qc].[5] B. P. Abbott et al. (LIGO Scientific Collaboration,Virgo Collaboration), Living Rev. Rel. , 1 (2013),arXiv:1304.0670 [gr-qc].[6] K. Chatziioannou, C.-J. Haster, T. B. Littenberg, W. M.Farr, S. Ghonge, M. Millhouse, J. A. Clark, and N. Cor-nish, Phys. Rev. D , 104004 (2019).[7] J. Veitch et al. , Phys. Rev. D , 042003 (2015),arXiv:1409.7215 [gr-qc].[8] B. P. Abbott et al. (LIGO Scientific, Virgo), Class.Quant. Grav. , 055002 (2020), arXiv:1908.11170 [gr-qc].[9] C. Rover, R. Meyer, and N. Christensen, Class. Quant.Grav. , 015010 (2011), arXiv:0804.3853 [stat.ME]. [10] C. Talbot and E. Thrane, (2020), arXiv:2006.05292[astro-ph.IM].[11] S. A. Usman et al. , Class. Quant. Grav. , 215004(2016), arXiv:1508.02357 [gr-qc].[12] S. Sachdev et al. , (2019), arXiv:1901.08580 [gr-qc].[13] B. Zackay, T. Venumadhav, J. Roulet, L. Dai, andM. Zaldarriaga, (2019), arXiv:1908.05644 [astro-ph.IM].[14] R. DeRosa, J. C. Driggers, D. Atkinson, H. Miao,V. Frolov, M. Landry, J. A. Giaime, and R. X.Adhikari, Classical and Quantum Gravity , 215008(2012), arXiv:1204.5504 [physics.ins-det].[15] V. Tiwari et al. , Class. Quant. Grav. , 165014 (2015),arXiv:1503.07476 [gr-qc].[16] G. D. Meadors, K. Kawabe, and K. Riles, Class. Quant.Grav. , 105014 (2014), arXiv:1311.6835 [astro-ph.IM].[17] J. Driggers et al. (LIGO Scientific), Phys. Rev. D ,042001 (2019), arXiv:1806.00532 [astro-ph.IM].[18] D. Davis, T. Massinger, A. Lundgren, J. Driggers, A. Ur-ban, and L. Nuttall, Class. Quant. Grav. , 055011(2019), arXiv:1809.05348 [astro-ph.IM].[19] G. Vajente, Y. Huang, M. Isi, J. C. Driggers, J. S. Kissel,M. J. Szczepanczyk, and S. Vitale, Phys. Rev. D , t (s) from GPS 1165069534.9 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1165069534.9 f ( H z ) N o r m a li ze d e n e r g y t (s) from GPS 1165069534.9 f ( H z ) N o r m a li ze d e n e r g y FIG. 11. Spectrogram of the LIGO Livingston data around the time of the blue mountain glitch for the leftmost injection fromFig. 9. Left panel: data containing the blue mountain glitch and the simulated CBC signal. Middle panel: data after a fairdraw from the glitch model has been subtracted leaving behind only the chirping CBC signal. Right panel: data after a fairdraw from the glitch and CBC models has been subtracted, leaving behind only Gaussian detector noise. − σ n o i s e LIGO Hanford − . − . − . − . − .
03 0 . t (s) from GPS 1187008881.4457 − − − σ n o i s e LIGO Livingston
Data CBC Glitch
FIG. 12. Credible intervals for the glitch (orange) and theCBC (blue) signal reconstruction for GW170817. Shadedregions correspond to 90% credible intervals, while in greydashed lines we plot the data whitened with a fair draw PSDfrom our noise model posterior. The top row corresponds toLIGO Hanford while the bottom row corresponds to LIGOLivingston. The LIGO Hanford plot zooms in to show thesignal that is invisible in the LIGO Livingston plot due to thesize of the glitch; note the y-scale difference in the two plots042003 (2020), arXiv:1911.09083 [gr-qc].[20] R. Ormiston, T. Nguyen, M. Coughlin, R. X. Adhikari,and E. Katsavounidis, Phys. Rev. Res. , 033066 (2020),arXiv:2005.06534 [astro-ph.IM].[21] S. Coughlin et al. , Phys. Rev. D , 082002 (2019),arXiv:1903.04058 [astro-ph.IM].[22] N. J. Cornish and T. B. Littenberg, Class. Quant. Grav. , 135012 (2015), arXiv:1410.3835 [gr-qc].[23] N. J. Cornish, T. B. Littenberg, B. B´ecsy, K. Chatziioan-nou, J. A. Clark, S. Ghonge, and M. Millhouse, (2020),arXiv:2011.09494 [gr-qc].[24] P. J. Green, Biometrika , 711 (1995),http://oup.prod.sis.lan/biomet/article-pdf/82/4/711/699533/82-4-711.pdf.[25] T. B. Littenberg and N. J. Cornish, Phys. Rev. D91 ,084034 (2015), arXiv:1410.3852 [gr-qc].[26] BayesWave Glitch Subtraction for GW170817, https://dcc.ligo.org/LIGO-T1700406/public . . . . . CBC+glitch+noiseCBC+noise (cleaned) . . . . M ( M (cid:12) ) . . . . q . . . χ e ff M ( M (cid:12) ) . . . .
00 0 .
05 0 . χ eff FIG. 13. One- and two-dimensional posterior distribution forselected source parameters for GW170817. We include themass ratio q , the effective spin χ eff , and the detector framechirp mass M posteriors. Blue curves show posteriors underthe CBC+glitch+noise model on the full data, while orangecurves correspond to the CBC+noise model on data where theglitch has already been subtracted. The two sets of results areconsistent with each other.[27] C. Pankow et al. , Phys. Rev. D98 , 084016 (2018),arXiv:1808.03619 [gr-qc].[28] LIGO Scientific Collaboration, Virgo Collaboration,“LALSuite,” https://git.ligo.org/lscsoft/lalsuite (2018).[29] N. J. Cornish, (2021), arXiv:2101.01188 [gr-qc].[30] N. J. Cornish and K. Shuman, Phys. Rev. D , 124008(2020), arXiv:2005.03610 [gr-qc].[31] N. J. Cornish, (2016), arXiv:1606.00953 [gr-qc].[32] Gravitational Wave Open Science Center (GWOSC), .[33] R. Abbott et al. (LIGO Scientific, Virgo), (2019),arXiv:1912.11716 [gr-qc].[34] S. Husa, S. Khan, M. Hannam, M. P¨urrer, F. Ohme,X. Jim´enez Forteza, and A. Boh´e, Phys. Rev. D ,044006 (2016), arXiv:1508.07250 [gr-qc].[35] S. Khan, S. Husa, M. Hannam, F. Ohme, M. P¨urrer,X. Jim´enez Forteza, and A. Boh´e, Phys. Rev. D ,044007 (2016), arXiv:1508.07253 [gr-qc].[36] M. Cabero et al. , Class. Quant. Grav. , 15 (2019),arXiv:1901.05093 [physics.ins-det]. − − σ n o i s e LIGO Hanford
Data CBC Glitch − . − . − .
02 0 .
00 0 .
02 0 . t (s) from GPS 1126259462.391 − − σ n o i s e LIGO Livingston
FIG. 14. Credible intervals for the glitch (orange) and theCBC (blue) signal reconstruction for GW150914. Shadedregions correspond to 90% credible intervals, while in greydashed lines we plot the data whitened with a fair draw PSDfrom our noise model posterior. The top row corresponds toLIGO Hanford while the bottom row corresponds to LIGOLivingston. Our glitch model recovers essentially no incoher-ent power coincident with the astrophysical signal and there-fore the reconstruction is not visible in the scale of the plot. . . . CBC+glitch+noiseCBC+noise . . . M ( M (cid:12) )
28 30 32 340 . . . . . . . q − . . . χ e ff
28 30 32 34 M ( M (cid:12) ) − . . . − . . . χ eff FIG. 15. One- and two-dimensional posterior distribution forselected source parameters for GW150914. We include themass ratio q , the effective spin χ eff , and the detector framechirp mass M posteriors. Blue curves show posteriors un-der the CBC+glitch+noise model, while orange curves corre-spond to the CBC+noise model. The two sets of results areconsistent with each other. [37] J. Powell, Class. Quant. Grav. , 155017 (2018),arXiv:1803.11346 [astro-ph.IM].[38] T. Accadia et al. , Classical and Quantum Gravity ,194011 (2010).[39] S. Soni et al. , arXiv e-prints , arXiv:2007.14876 (2020),arXiv:2007.14876 [astro-ph.IM].[40] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 061102 (2016), arXiv:1602.03837 [gr-qc].[41] T. Dietrich et al. , Phys. Rev. D99 , 024029 (2019),arXiv:1804.02235 [gr-qc].[42] B. P. Abbott et al. (LIGO Scientific, Virgo), Class.Quant. Grav. , 134001 (2016), arXiv:1602.03844 [gr-qc].[43] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 241102 (2016), arXiv:1602.03840 [gr-qc].[44] B. P. Abbott et al. (LIGO Scientific, Virgo), Class.Quant. Grav. , 065010 (2018), arXiv:1710.02185 [gr-qc].[45] D. Davis, L. V. White, and P. R. Saulson, Class. Quant.Grav. , 145001 (2020), arXiv:2002.09429 [gr-qc].[46] D. Macleod, A. L. Urban, S. Coughlin, T. Massinger,M. Pitkin, paulaltin, J. Areeda, E. Quintero, T. G. Bad-ger, L. Singer, and K. Leinweber, “gwpy/gwpy: 1.0.1,”(2020).[47] J. D. Hunter, Computing In Science & Engineering9