Parameterised population models of transient non-Gaussian noise in the LIGO gravitational-wave detectors

The two interferometric LIGO gravitational-wave observatories provide the most sensitive data to date to study the gravitational-wave Universe. As part of a global network, they have just completed their third observing run in which they observed many tens of signals from merging compact binary systems. It has long been known that a limiting factor in identifying transient gravitational-wave signals is the presence of transient non-Gaussian noise, which reduce the ability of astrophysical searches to detect signals confidently. Significant efforts are taken to identify and mitigate this noise at the source, but its presence persists, leading to the need for software solutions. Taking a set of transient noise artefacts categorised by the GravitySpy software during the O3a observing era, we produce parameterised population models of the noise projected into the space of astrophysical model parameters of merging binary systems. We compare the inferred population properties of transient noise artefacts with observed astrophysical systems from the GWTC2.1 catalogue. We find that while the population of astrophysical systems tend to have near equal masses and moderate spins, transient noise artefacts are typically characterised by extreme mass ratios and large spins. This work provides a new method to calculate the consistency of an observed candidate with a given class of noise artefacts. This approach could be used in assessing the consistency of candidates found by astrophysical searches (i.e. determining if they are consistent with a known glitch class). Furthermore, the approach could be incorporated into astrophysical searches directly, potentially improving the reach of the detectors, though only a detailed study would verify this.


I. INTRODUCTION
The LIGO, Virgo, and KAGRA gravitational-wave detectors [1][2][3] have opened the door to the Gravitational-Wave Universe and provided our first insights into the mergers of black holes and neutron stars. These kilometer-scale interferometers, sensitive to signals in the tens of Hz to kHz band, survey the sky nearly isotropically during observing runs [4]. Between each observing run, the detectors are improved, reducing the level of noise and yielding improvements in the sensitivity. Therefore, each new observing run probes deeper than before, yielding an increase in the rate of detections.
The LIGO detectors contain noise which is often assumed to be stationary and Gaussian for analysis, but often exhibits non-stationarity and frequent non-Gaussian transient noise artefacts termed glitches [5][6][7][8][9]. All glitches are believed to be instrumental or environmental in origin. In some cases, the cause of a glitch class can be identified and improvements made to the detector to reduce or eradicate the class. In other cases, so-called witness channels can identify when such glitches occurred and veto coincident candidate events (see, e.g. Davis et al. [9]). Finally, glitches in classes without a known origin or witness channel remain a feature of the sciencemode data from interferometric gravitational-wave observatories. * gregory.ashton@ligo.org Glitches are detrimental to achieving the scientific goals of gravitational-wave observatories.
Suppose glitches occur in coincidence with an actual astrophysical signal. In that case, this may result in the signal being missed by searches, reducing the number of detectable signals, or biasing our astrophysical inferences about the signal itself. However, glitches that do not occur in coincidence with real signals also reduce the sensitivity of the astrophysical searches. To identify signals, search pipelines use information about waveform morphology and coincidences between detectors (see, e.g. [10][11][12][13][14]) to calculate a detection statistic. The presence of glitches in the detector data produces values of the detection statistic larger than the expectation from Gaussian noise alone. To estimate the non-astrophysical background of the detection statistic due to glitches, most search pipelines use a bootstrap approach, for example, by time-shifting the data from independent detectors. These approaches enable a significance to be attached to a candidate (for example, by the false alarm rate or FAR). Therefore, a candidate's significance is affected both by how alike it is to predictions for astrophysical signals and the level of background noise created by glitches. Vetoing glitches or utilising detection statistics that better identify transient noise outliers (see, e.g. [15]) reduce the level of background noise from glitches and hence improve the number of sources that can be identified at a given confidence level. But, in either case, the presence of glitches makes it more challenging to identify astrophysical signals.
Significant effort has been made to classify glitches into arXiv:2110.02689v1 [gr-qc] 6 Oct 2021 separate classes based on their morphology [7,[16][17][18][19][20][21]. Glitch classification can potentially identify the cause of glitches enabling detector improvements to reduce their impact. Classification can also help mitigate the impact of glitches on a search, for example, by building a glitchrobust detection statistic [22] or help qualitatively interpret putative astrophysical signals by placing them in the context of known glitch classes.
Davis et al. [23] investigated the impact of four glitch classes, Blips, Koi Fish, Scattered Light, and Scratchy glitches on the PyCBC search pipeline [12]. Calculating a rate at which the glitches pass a given detection threshold over the physical parameter space of sources, they demonstrate a method to calculate the significance of a possible signal that coincides with a glitch. An alternative approach explored in Ashton et al. [24] builds information about the distribution of glitches into a Bayesian odds, eschewing a bootstrapped estimation of the background. However it is done, detection approaches which reduce the impact of glitches are critical to improving the search sensitivity and hence identifying more astrophysical candidates.
In this work, we take a modelled approach to glitch classification. We aim to deliver simple probabilistic models for the glitch classes most harmful to astrophysical analyses. To do this, we analyse single-detector glitch triggers pre-classified by the GravitySpy citizen science project [16] from the LIGO Hanford and Livingston detectors [1]. We analyse each glitch using a model of a precessing binary black hole merger (see Merritt et al. [25] for an alternative approach developing parametric glitch models). Using an astrophysical black hole merger model enables us to project the properties of the glitch into the physical parameter space of astrophysical signals. Taking these projections, we then build simple parameterised population models for each glitch class and study this population with respect to known astrophysical signals.
The rest of this paper is structured as follows. In Sec. II, we introduce the methodology based on the principles of population inference. In Sec. III, we describe our results fitting population models to 4 different glitch classes. In Sec. IV, we compare the glitch population with the population of binary black-hole mergers described in the GWTC2.1 transient gravitational-wave catalogue. Finally, we conclude in Sec. V.

II. METHODOLOGY
The methodology used herein follows the principles of hierarchical population inference for compact binary coalescence (CBC) signals (see Refs. [26][27][28] and Thrane and Talbot [29] for a review). We aim to identify the population properties of 4 glitch classes Blip, Tomte, Fast Scattering, and Scattered Light which are known to most adversely impact search sensitivity [9,23]. In Fig. 1, we provides time-frequency spectrograms demonstrating typical examples of each glitch class. These classes were the most numerous during the O1 to O3 observing runs and will likely persist for future observing runs. We perform the population inference separately on glitches in the LIGO Hanford (H1) and LIGO Livingston (L1) detectors. This produces two sets of population properties that we can use to understand the consistency of glitch classes between the detectors. To infer the properties of glitches in a given class, we first identify a list of glitch times pre-categorised by GravitySpy with a confidence threshold of 95%. We check that these glitches would typically contribute to the noise backgrounds estimated by search pipelines, i.e. they would not be vetoed by the single-detector signalto-noise ratio (SNR) cuts; the details of this are given in the subsections of Sec. III. Each glitch is individually vetted by eye using a time-frequency visualisation to ensure it is consistent with the corresponding glitch class. We then analyse the data d i surrounding the i th glitch in our list using a stochastic sampling algorithm (specifically, the nested sampling dynesty [30] algorithm through the Bilby [31] package) and a model M . In this work, we use the IMRPhenomPv2 [32,33] model which includes the inspiral, merger and ringdown of merging binary black holes and has been a standard benchmark waveform for several observing runs [34,35]. IMRPhenomPv2, however, only models the dominant ( , m) = (2, 2) modes in the gravitational-wave radiation. Recent advances in waveform modelling have produced improved models, including the effect of higher-order modes [36][37][38][39][40][41]. Nevertheless, we use IMRPhenomPv2 as it is well standardised and fast to evaluate, enabling us to analyse hundreds of glitches without a high computational cost. However, this implies that we infer the properties of glitches given the IMRPhenomPv2 waveform model. It is well known that systematic uncertainty exists between waveform models and that models including higher-order modes produce quite different posterior distributions, especially when the system has asymmetric masses [42][43][44]. We, therefore, caution when comparing the glitch population obtained herein with candidate triggers that are analysed using waveforms other than IMRPhenomPv2, especially if those waveforms include the effect of higher-order modes.
The IMRPhenomPv2 model has an associated set of model parameters θ, which describe the mass, spins, location, and orientation of the source. For each glitch class, we produce an estimate of the posterior distribution p(θ|d i , M ) (in the form of a set of posterior samples) and the Bayesian evidence Z(d i |M ). During the analysis of individual glitches, we use a Bayesian prior π(θ) appropriate for analysing astrophysical binary black holes. Specifically, the prior is uniform in the detector-frame chirp mass 1 M det [45] between 8 and 200 M and uniform in mass ratio q between 1/20 and 1 and assumes an isotropic distribution on the black holes spins. For all other parameters, we use the standard non-informative astrophysical prior distributions defined in Veitch et al. [46].
To infer the population properties of glitches in a given class, we take the set of posterior samples and Bayesian evidence from each individual analysis and perform a hierarchical population inference. Specifically, for a parameter θ i ∈ θ, we define a population hyper-model H with associated model hyperparameters Λ such that we can evaluate π(θ i |Λ, H). Once defined, we then utilise a variation on importance sampling to construct the likelihood L({d i |Λ, H). (We describe the population hypermodels and specifics of the likelihood used in this work below). Given this likelihood and a suitable hyperprior π(Λ|H), we then use stochastic sampling to infer p(Λ|{d}, H), the population hyperparameters conditioned on all glitches in the given class. For this analysis, we use the Bilby-MCMC [47] Markov-Chain Monte-Carlo sampling algorithm; we verify that both the dynesty and Bilby-MCMC samplers produce equivalent results, but that Bilby-MCMC is more efficient in this instance where only the posterior distribution and not the Bayesian evidence is of interest.
In principle, a hyper-model can predict the population behaviour for the full 15-dimensional parameter space θ, including correlations between parameters. However, for a sub-set of parameters, we do not observe any population trends (e.g., the sky localisation of a population of glitch events tends to be isotropically distributed). We find that three parameters, the chirp mass M det , mass ratio q, and the dimensionless spin-magnitude of the primary (heavier-mass) object a 1 are the primary indicators for the glitch classes considered in this work. We will build population models for each of these parameters independently. There is a weak correlation between these parameters, but we choose to infer their properties sep-arately. This discards information about the correlation but greatly simplifies the population models. In future work, we will investigate whether including such correlations can produce a more accurate probabilistic model.
For other parameters, we do not observe strong trends in the population, except for the detector-frame luminosity distance, for which we observe a piling up of events at the lower bound (100 Mpc). The glitches selected for this analysis are, by design, loud glitches which pass a confidence threshold for categorisation by the GravitySpy (see Zevin et al. [16] for details). Hence, our set of glitches has an intrinsic selection bias in the loudness and hence luminosity distance. It is no surprise then that we observe a trend that all glitches tend to have small luminosity distances, as this is the dominant term in the signal amplitude and hence determines how "loud" a glitch is. We choose not to model the luminosity distance as we would need to properly calibrate the model for the intrinsic selection bias of our set of glitches. Other parameters, such as the system mass, may also be affected by the selection bias, but we model them nonetheless.
In this work, we will use three hypermodels to capture the broad-scale population properties of each glitch class for each parameter. These are three hypermodels from an infinite set of possible models. We arrived at these by studying the performance of a larger superset of hypermodels using posterior-predictive checks (described below). Where more complicated models failed to improve the posterior-predictive checks, we used simpler models. These hyper-models are designed to capture the broad features in a simple parameterised form. The three hypermodels are: M1: A power law distribution, it has a single hyperparameter Λ = {α} and density with bounds given by the physical bounds on θ i . This distribution is useful for modelling populations which tend to rail against the physical bounds (for example, the primary spin). For α, we use a normally distributed prior with zero mean and a standard deviation of 5. This removes the requirement for an arbitrary upper bound on the magnitude while exponentially suppressing very large values.
M2: A truncated skew-normal distribution [48,49], it has three hyperparameters Λ = {µ, σ, κ} and density where φ(x) is the standard normal probability density function and Φ(x) is the error function. We truncate the distribution again to the physical bounds on θ i . This distribution is useful for modelling unimodal populations with a distinct peak.
We use a uniform prior for µ over the prior support on θ i , a normal prior for κ with zero mean and a standard deviation of 10 and a truncated normal prior for σ with standard deviation given by the maximum support on θ i .

M3
: A mixture-model of two normal distributions with five hyperparameters Λ = {ξ, µ 0 , σ 0 , µ 1 , σ 1 } and density This is effective in modelling multi-modal distributions; we do not include a skew as this was found to be poorly constrained in the cases where the M3 model was required. As is done for the M2 model, we use a uniform prior for µ 0 and µ 1 and a truncated normal distribution for σ 0 and σ 1 . For ξ, we apply a uniform prior on [0, 0.5].
For each model, the likelihood is constructed from Equation (32) of Thrane and Talbot [29] using "recycling". This enables a two-step approach (first analysing individual events, then the population properties), which is critical to reducing the wall-time of the analysis.
For each hyper-model, we report the median values of p(Λ|{d i }, H). These, together with the model descriptions above, constitute the core product of this work: simple parameterised models of glitches in the LIGO detectors.
To verify that the hyper-model captures the broadscale features of the population, we use predictive posterior checking. First, to plot the "Measured" population, we bin each individual glitch-posterior in a histogram, then plot the mean and the 90% credible interval (between 5% and 95%). This provides a summary of the population properties, which includes the variation in the individual posteriors 2 . An example of this population visualisation can be seen in the right-hand column of Table II: the "Measured" curve shows the mean (solid grey line) and 90% credible interval (filled grey region) for the population of Blip glitches across separate parameters. Second, to plot the "Predicted" population, we repeat the process above but simulate the individual posteriors from the population model. Specifically, we draw θ i from p(θ i |Λ ) (where Λ is the median) 3 then simulate a posterior and repeat the steps above. The simulated posterior is derived by randomly sampling a posterior from the "Measured" population and aligning its mean with a random draw from the posterior population distribution. This is a non-optimal method, but the alternative (simulating the signal predicted by the posterior population distribution in real detector noise unaffected by other transient noise and performing Bayesian inference) is too computationally demanding to be feasible. Examples of the resulting posterior predictive checks can be found in the right-hand column of Table II: the "Predicted" curve shows the mean (solid orange line) and 90% credible (filled orange region) predicted by the hyper model across separate parameters. Here we see that the models capture the broad features (e.g. the number and location of modes and their width) but do not always capture the details (e.g. the Measured and Predicted curves do not overlap perfectly).

III. RESULTS
We now describe the results for each glitch class in the sections below. For each glitch class, we provide a table with a summary (the glitch model and median inferred population parameters) and the results of the posterior predictive checks. We also provide this information in a machine-readable table A [50]. For each model, we report the results of a single hypermodel. However, we fitted various models during development and used the posterior predictive checks as a qualitative guide to select the best fitting model.

A. Tomte glitches
In Table I, we report the hypermodel, median, and posterior predictive checks for the population properties in M det , q, and a 1 for Tomte glitches in the H1 and L1 detectors. For each detector, we analyse a set of 1000 Tomte glitches identified with GravitySpy. 4 . The SNR of these glitches ranges from 7.7 to 74; the distribution is peaked with a median SNR of 21, but 90% of the Tomte glitches have an SNR less than 36. 5 We also verify that we find similar conclusions about the properties when we downsample to 100 glitches; we use this, later on, to justify analysing smaller sets of the Fast Scattering and Scattered Light glitches.
For the chirp mass, the measured glitch population is unimodal with a peak at ∼ 75 M and some positive skew. We fit the population using the M2 Skew-Normal distribution and find a qualitatively acceptable fit for both detector populations. However, the model underpredicts the rate of glitches with M det ∼ [125, 175] in the H1 detector. While the broad features are consistent between H1 and L1, the amount of skew is noticeably different, and Tomte glitches in the H1 detectors have chirp mass values up to ∼ 150, while Tomte glitches in the L1 detector, M det 125.
For the mass ratio, the population is unimodal with some skew and the bulk of support on q 1/2. We fit the mass ratio distribution using the M2 Skew Normal model and find a good qualitative fit in H1 and L1. As with the chirp mass, we note some qualitative differences between the H1 and L1 detectors with greater support in the H1 detectors for mass ratio's up to 0.5, which is not seen in the L1 glitch population.
The primary spin, a 1 piles up at the extreme spin case a 1 = 1; we fit this using the M1 power-law model and find a reasonable agreement. However, we note a lack of events at the a 1 = 1 bound. Our fitted model does not capture this feature; however, it is sufficiently narrow that overall the model remains a reasonable fit. The power-law index α ∼ 5 − 6 provides a rough guide to the extent to which the spin "piles up"; this value is less extreme than in the Blip glitch population discussed shortly. For the spin, we do not observe significant differences between the H1 and L1 detector populations.

B. Blip glitches
In Table II, we report the hypermodel, median, and posterior predictive checks for the population properties in M det , q, and a 1 for Blip glitches in the H1 and L1 detectors. We analyse a set of 1000 Blip glitches in both H1 and L1 identified with GravitySpy. 6 . The SNR of these glitches ranges from 7.5 to 47. The distribution peaks at the minimum bound, and 90% of the glitches have an SNR less than 20. We also verify that our results are robust when using only a subset of 100 randomly selected glitches.
The measured glitch population is unimodal for the chirp mass, with a peak at ∼ 25M and some positive skew. We find that the M1 Skew-Normal model provides a good qualitative fit. Comparing the measured and predicted posterior checks, the distribution means are well aligned, though the Skew-Normal under-predicts the density in the 40 − 60 M region.
For the mass ratio, the population is bimodal, with identical features found in both detectors. The bimodality originates from bimodal features observed in the posterior distributions of individual events. Though we see significant scatter in the shape of individual posteriors, averaging over multiple glitches, we observe similar distributions in both the H1 and L1 detectors indicating a common origin. The bulk of the support for both modes is found in the region q 1/2, with some support up to the equal-mass bound. We apply the M3 normal mixture model, which provides a reasonable fit to the data, though it tends to underpredict the number of nearly equal-mass glitches in the L1 detector.
Finally, for a 1 , the magnitude of the primary spin, we find that the measured glitch values pile up onto the extreme spin case a 1 ∼ 1. We also found this feature in the Tomte population, but here it is more extreme. We find the simplistic M1 power-law model provides a suitable fit to this data; the extreme power-law index α ∼ 37 demonstrates the strength of the pile-up.
Across all three parameters, we find remarkably consistent glitch behaviour between the H1 and L1 detectors. This, of course, may be simply due to the robust classification of events by GravitySpy. However, we note that the Tomte glitches did not show such consistency.

C. Fast Scattering glitches
In Table III, we report the hypermodel, median, and posterior predictive checks for the population properties in M det , q, and a 1 for Fast Scattering glitches in the H1 and L1 detectors. We use a set of 100 Fast Scattering glitches identified with GravitySpy. We limit ourselves to 100 glitches based on our analysis of the Blip and Tomte glitches, which indicated 100 glitches is sufficient to represent the population properties. This ten-fold reduction in events correspondingly reduces the computational burden of the analysis. The SNR distribution of these glitches is strongly peaked at the minimum bound of 7.5, with 90% of the glitches have an SNR less than 13. However, the distribution has a long tail with a maximum value of 121.
For the chirp mass, we find that the population has support from 75 M up to the artificial prior bound of 200 M . This indicates that the actual population likely have support at larger masses still. Such extreme heavy systems are beyond the space of likely astrophysical signals that ground-based detectors will observe in the advanced era. As such, extending our prior beyond this artificial bound is unlikely to help distinguish glitches from signals in the near term, so we choose not to do so. The observed distribution is "lumpy" with multiple components. We choose to fit the distribution with a Skew Normal distribution. This does not capture the "lumpy" behaviour but does broadly capture the typical support.
For the mass ratio, the Fast Scattering population have extreme mass ratios with a median of ∼ 1/10 with some support up to q ∼ 0.8. We fit this distribution using the M2 Skew Normal model and find a reasonable agreement between the predicted and measured distribution.
The primary spin, as with previously studied glitches, peaks at the extremal a 1 = 1 bound. However, the mea-sured population has significant support for systems with a 1 right down to zero. We fit this distribution using the M1 power-law model and find a good fit.

D. Scattered light glitches
In Table IV, we report the hypermodel, median, and posterior predictive checks for the population properties in M det , q, and a 1 for Scattered Light glitches in the H1 and L1 detectors. We use a set of 100 Scattered Light glitches identified with GravitySpy. The choice to use just 100 glitches mirrors the motivation in Sec. III C. The SNR of these glitches ranges from 7.5 to 370; as with the Fast Scattering glitch type, the distribution peaks at the lower bound with a long tail. 90% of the glitches have an SNR less than 47.
For the chirp mass, we find that the population has support from 75 M up to the artificial prior bound of 200 M with significant "lumpiness". The behaviour in the chirp mass distribution is similar in form between the Fast Scattering and Scattered Light glitches. For the L1 detector, there is a significant overabundance at the equal-mass bound. We choose to fit the distribution with a Skew Normal distribution. This does not capture the "lumpy" behaviour but does broadly capture the typical support.
For the mass ratio, the Scatter Light population have extreme mass ratios with a median of ∼ 1/10. But, unlike the Scattered Light glitches, q 0.2. We fit this distribution using the M2 Skew Normal model and find a reasonable agreement between the predicted and measured distribution.
As with previously studied glitches, the primary spin peaks at the extremal a 1 = 1 bound but with little support for smaller spins. We fit this distribution using the M1 power-law model and find a good fit. The powerlaw index α ∼ 28 is nearly as large as that of the Blip glitches.

IV. COMPARISON WITH THE GWTC2.1 CATALOGUE
So far, 57 gravitational-wave signals have been reported and analysed by the LIGO, Virgo and KAGRA collaborations [34,35,52,53]. Predominantly, these have been binary black hole systems (BBH) except for two binary neutron star (BNS) and two neutron star black hole (NSBH) binaries. Here we compare our glitch population with the 47 events with p astro > 0.5 reported in either the GWTC2 [35] or GWTC2.1 [52]  We provide two-dimensional plots comparing our glitch population with the observed signal in M det -q (Fig. 2), M det -a 1 (Fig. 3), and q-a 1 (Fig. 4). In each plot, we give the 90% credible interval of the four glitch classes (for both H1 and L1) along with the combined equal-weighted glitch population in black. Projections of the overall glitch population are given for each one-dimensional plane. Overlaid on this, we then plot the symmetric 90% credible interval for all events with p astro > 0.5 reported in either the GWTC2 or GWTC2.1 catalogues.
Comparing the glitch population with the astrophysical BBH distribution, three features stand out: i) glitches have large spin magnitudes while astrophysical signals have median spin magnitudes 0.9; ii) glitches have extreme mass ratios while astrophysical signals tend to have more equal-mass; and finally, iii), the combined glitch population spans the entire chirp mass range considered in GWTC2, but which had a re-calculated pastro < 0.5. This superset, therefore, includes all events during O3a reported by the collaboration with pastro > 0.5. in this work, though individual glitch classes tend to cluster.
Glitches have extreme mass ratios and spin because the CBC signals models are a poor fit to the glitch morphology (see, e.g. Merritt et al. [25]). In the equal-mass zero-spin limit, the CBC waveform is sinusoidal with a characteristic monotonically increasing frequency evolution. Non-equal mass ratios and spins modify this behaviour introducing "twisting-up" effects which produce irregular waveform morphology's. Glitches are caused by terrestrial disturbances, which are not expected to characteristically chirp up in frequency. So, when fitting the CBC signal model, it is unsurprising that the best fits happen in the non-equal mass and large spin limits.
Of the BBH systems reported so far in O3a, just one, GW190403 051519, falls firmly inside the 90% credible interval for the Tomte population in chirp mass, mass ratio, and primary spin magnitude. Additionally, GW190929 012149 sits at the edge of the Tomte distribution, but no signals are consistent with the other 3 glitch classes. Several other signals fall within the interval for the chirp mass and primary spin magnitude but are clearly distinguished when looking at the mass ratio and primary spin magnitude. 8 However, it should be noted that both GW190403 051519 and GW190929 012149 have lower signal to noise ratios than any of the glitches included in our population analysis since we select only high-confidence glitches. However, Merritt et al. [25] argue that both quiet and loud glitches (of the same class) exhibit similar morphology.
In Fig. 5, we plot the time-frequency spectrogram of GW190403 051519 and a characteristic signal track. This illustrates that the event is primarily identified in the L1 detector, with no clear signal found in the H1 detector.
GW190403 051519 was first identified in the GWTC2.1 catalogue by the GWTC2.1 PyCBC-BBH search alone [52] with a probability of astrophysical origin (p astro ) of 0.61.
No other pipelines identified GW190403 051519 as a candidate. Meanwhile, GW190929 012149 was identified in the GWTC2 catalogue by the GstLAL search alone [35] with p astro =1.
However, in the GWTC2.1 catalogue, GW190929 012149 was subsequently identified by all four pipelines used with p astro values ranging from 0.14 to 0.87. GW190403 051519 is notable because the mass of the primary (heavier) object has more than mode content not modelled by the IMRPhenomPv2 waveform used in this work can result in significant shifts in the posterior distribution. However, we re-analysed GW190403 051519 using the IMRPhenomPv2 waveform and found the event remained consistent with the Tomte population. Between the GWTC2.1 analysis and our re-analysis, the median detector-frame chirp mass shifted from 86 to 56 M , the mass ratio from 0.31 to 0.19, and the primary spin magnitude from 0.93 to 0.85. Time-frequency spectrograms of the GW190403 051519 event. We overlay the predicted signal using the SEOBNRv4 ROM waveform approximant [54] and the median-estimated signal properties from our IMRPhenomPv2 re-analysis of the event. (We note a mismatch between the waveform used for analysis and the waveform used to plot the signal track. However, the signal track here is intended only to guide the eye and not a quantitative test.) 50% support for a mass greater than 100 M and is either inside or above the "mass gap" predicted by pair-instability supernova theory [55].
That GW190403 051519 and GW190929 012149 are consistent with the Tomte glitch population does not, of course, mean that the events are Tomte glitches. However, they have p astro values, which allow for a terrestrial cause (i.e. 1 − p astro 0.1). Tomte glitches were numerous during O3a and are known to impact search sensitivity adversely. Therefore, it seems reasonable to conclude that if GW190403 051519 or GW190929 012149 are not of astrophysical origin, of the four glitch classes considered here, they are most consistent with Tomte glitches. However, we add that this assumes the signal morphology of the glitch classes studied herein extends to the low signal to noise ratio regime of GW190403 051519 and GW190929 012149.
This work does not change the estimated p astro values for GW190403 051519 or GW190929 012149. The search pipelines that identify signals are robust to all relevant glitch classes, use information about the consistency between separate detectors, and apply glitch discriminators [15]. Using bootstrap approaches (e.g. timeshifting the data from independent detectors), the search pipelines estimate the background noise from all glitch classes, and this is intrinsically part of the p astro calculation. Nevertheless, Fig. 4 does demonstrate that to identifying signals which have large spins and unequal mass ratios, astrophysical searches must contend with an elevated background of glitches compared to other areas of the parameter space. Search pipelines account for the distribution of noise triggers (including glitches) over the searched parameter space in their detection statistics and the resulting significance estimates (both the False Alarm Rate, FAR, and p astro ), though this is done in a variety of ways [13,14,56,57]. Therefore, the variation in glitch rate across the parameter space indicates that search pipelines are likely not uniformly sensitive across the parameter space [23,58].

V. DISCUSSION AND CONCLUSION
We present population models for the behaviour of non-Gaussian transient noise (glitches) in the Hanford and Livingston LIGO detectors during the O3a observing run. We first identify lists of Blip, Tomte, Fast Scattering, and Scattered light glitches pre-categorised by the GravitySpy citizen-science project and then vetted for class consistency using time-frequency visualisations. These glitch classes are known to most adversely search sensitivity. Using Bayesian inference, we then analyse each glitch class using a precessing binary black hole model. In effect, this allows us to project the properties of the glitch onto the parameter space of astrophysical signals. Finally, we use hierarchical population inference methods to fit simple parameterised models to the glitch populations. We do this separately for the chirp mass, mass ratio, and primary spin magnitude parameters; all other parameters we find to be only weakly informative at the population level. We verify the performance of the models using posterior predictive checks and confirm that they reproduce the bulk features of the population.
Finally, we compare the population of glitches with events observed with p astro > 0.5 during the O3a observing run. We find two events, GW190403 051519 and GW190929 012149, which have a chirp mass, mass ratio, and primary spin magnitude consistent with the Tomte glitch class. GW190403 051519 is a fascinating case as it has a mass that may be inconsistent with the predictions of pair-instability supernova theory.
The parameterised models developed in this work could be used by search pipelines to potentially better distinguish astrophysical signals from terrestrial glitches 9 . This could be implemented either additively in the detection statistic used by matched-filtering pipelines or as part of the astrophysical odds [24] (see [59] for the use of a simplistic glitch model). Already, many search pipelines use a version of this in which the detection statistics and background distributions are calculated over the template banks. This enables a lowering of the background rate in regions of parameter space less contaminated by glitches. Including the probabilistic models here into the detection statistic itself could refine this approach, providing more fine-grained modelled information. Ultimately, this may lead to a greater number of astrophysical signals being detected. However, the impact of this will depend on how much additional information is provided by these glitch models relative to the approaches already taken. To determine if such an approach can improve the detection efficiency, a large scale simulation study is needed, which is beyond the scope of the current work.
We also identify multimodality in the mass ratio of Blip glitches which may be used to identify this glitch class. We do not yet understand why this multimodality is exhibited, but it seems likely the waveform fits two separate parts of the glitch morphology in different regions of parameter space. This information may help formulate an improved glitch model (e.g., breaking the astrophysical waveform in a non-physical way that fits the glitch). Such a model could then be used to identify glitches and improve astrophysical searches.
At the same time, In the future, we hope to improve these models by including correlations between parameters, using higher-order mode waveform models and extending the study to sub-dominant parameters. These improvements will provide a more refined glitch consistency test and may reveal new glitch population features.

VI. ACKNOWLEDGEMENTS
We a grateful to Derek Davis for help with the script used to produce Fig. 1 and Fig. 5, to Andrew Williamson for insightful comments about the nature of the bimodality seen in the Blip glitch mass ratio, to Tom Dent for helpful feedback on the manuscript, and to the CBC and detector characterisation groups of the LIGO and Virgo Collaborations for helpful feedback on this work. We also think Richard O'Shaughnessy, Jason Hathaway, Ben Henderson, Brian Cowburn for useful discussions related to this work. GA and LKN thank the UKRI Future Leaders Fellowship for support through the grant MR/T01881X/1. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gwopenscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. 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. We are also grateful to computing resources provided by the LIGO Laboratory computing clusters at California Institute of Technology and LIGO Hanford Observatory supported by National Science Foundation Grants PHY-0757058 and PHY-0823459. The majority of analysis performed for this research was done using resources provided by the Open Science Grid [60,61], which is supported by the National Science Foundation award #2030508. This work makes use of the scipy [62], numpy [63][64][65], gwpy [66], and PyCBC [12] software for data analysis and visualisation.