Statistical characterization of pulsar glitches and their potential impact on searches for continuous gravitational waves

Continuous gravitational waves from neutron stars could provide an invaluable resource to learn about their interior physics. A common search method involves matched filtering a modeled template against the noisy gravitational-wave data to find signals. This method suffers a mismatch (i.e., relative loss of the signal-to-noise ratio) if the signal deviates from the template. One possible instance in which this may occur is if the neutron star undergoes a glitch, a sudden rapid increase in the rotation frequency seen in the timing of many radio pulsars. In this work, we use a statistical characterization of the glitch rate and size in radio pulsars to estimate how often neutron star glitches would occur within the parameter space of continuous gravitational-wave searches and how much mismatch putative signals would suffer in the search due to these glitches. We find that for many previous and potential future searches continuous-wave signals have an elevated probability of undergoing one or more glitches and that these glitches will often lead to a substantial fraction of the signal-to-noise ratio being lost. This could lead to a failure to identify candidate gravitational-wave signals in the initial stages of a search and also to the false dismissal of candidates in subsequent follow-up stages.


I. INTRODUCTION
Electromagnetic (EM) observations of pulsar glitches have long been one of the most fruitful sources of insight into neutron star physics. They are characterized by a sudden increase 1 in the rotation frequency, often accompanied by a jump in the frequency derivative and an exponential recovery of some fraction of the initial frequency jump. The events happen rapidly and are sufficiently disruptive that pulsar timing models often lose phase coherence over the event.
Two leading models exist to explain glitches. In the superfluid pinning model, some portion of the interior superfluid is pinned and does not participate in the smooth torque-driven spin-down of the rest of the crust (where "crust" refers to the actual crust, plus whatever other parts of the star that are strongly coupled to it). After some period, the crust will therefore have developed a frequency lag compared to the pinned superfluid. A glitch occurs when the two components recouple, transferring angular momentum from the pinned superfluid to the crust and producing a spin-up of the crust [2,3]. Alternatively, glitches could be caused by crust cracking as the crust readjusts to a minimum energy configuration brought about by the gradual decay of the spindown rate [4]. It is also possible that glitches result from a combination of these two models. In either case, it seems reasonable to assume that both the crust and the core will be involved.
Rotating isolated neutron stars can produce continuous gravitational waves (CWs) from nonaxisymmetric distortions. These distortions can be stationary as viewed from the frame rotating with the star, so-called "mountains" (supported by either elastic stresses in the crust or by magnetic fields), or persistent oscillation modes of the neutron star. In the case of mountains, the star emits a CW at a frequency f s that is twice the rotation frequency ν, i.e., f s ¼ 2ν. In the case of oscillation modes, the frequency relation will be different and generally depends on the equation of state and the type of oscillation (e.g., f s ≈ 4ν=3 in the case of r modes; see the review by Prix [5] for an overview of different CW emission mechanisms). In the following, we assume the mountain model for simplicity.
It is possible (see, for example, the work by van Eysden and Melatos [6] and Keer and Jones [7], or Singh [8]) that a glitch could trigger a quadrupolar quasinormal mode resulting in a transient burst of gravitational waves; search methodologies for such a signal have been considered by Clark et al. [9] for signals lasting OðmsÞ and by Prix et al. [10] for signals lasting Oðhours-weeksÞ. However, in this work, we are not concerned with the gravitational radiation triggered by a glitch but rather the impact glitches may have on searches for CW signals. Specifically, assuming an isolated rotating neutron star is producing a CW signal, if a glitch occurs and causes a sudden increase in the CW frequency, what will the impact be on our ability to detect the signal?
Estimates for the intrinsic gravitational-wave strain amplitude h 0 for canonical models of CW emissions (see, for example, the work of Abbott et al. [11]) suggest they are extremely weak compared to the noise level of advanced detectors [12]. To detect a signal, significant effort has been put into data analysis methods, which may be capable of identifying the putative signals. Many of these methods rely on matched filtering in which a template is correlated with the data in the hope of detecting the presence of the unknown signal similar to the template. The power of these methods lies in the fact that the signal-tonoise ratio (SNR) grows as the square root of the observation time (e.g., see Ref. [5] for an overview). Because of the longevity of CW signals, this allows the weak signal to be discerned from the noise by using a sufficiently long stretch of data.
These methods are powerful but harbor a vulnerability in any instance in which the signal lies outside of the regular CW template manifold (as defined in Sec. IVA). For CW signals from nonaxisymmetric distortions of neutron stars, we can expect that discrepancies from such a template may manifest in one of two ways. First, we know from radio pulsar timing that the spin-down of a pulsar differs from a smooth spin-down due to timing noise. This is a continuous low-frequency structure in the residual between the best-fit Taylor-series timing model and the observed pulsations (for a review, see the work of Hobbs et al. [13]). The effect of timing noise on CW searches was studied by Jones [14] and Ashton et al. [15], where it was found that its presence limits the coherent observation span over which signals may be detected if searched for using a smooth signal. However, one can attempt to mitigate this effect when performing a targeted search for a known radio pulsar by including the timing noise seen in the EM channel in the template; a method to do this is described by Pitkin and Woan [16]. In this work, we will address the second potential discrepancy between the signal and template: glitches.
There are two distinct questions to answer in the case of glitches: (1) How probable is it that a glitch will occur during our CW observation? (2) If a glitch does occur, what effect will it have on our ability to discover the CW signal? To answer these questions, we use known radio pulsar glitch statistics to estimate the size and rate of radio pulsar glitches for the parameter spaces considered in typical CW searches. In other words, we assume that the glitches we observe in the radio pulsars population are representative of those that we may expect to see in the population of CW signals. This is purely pragmatic in that we do not know of any other populations on which to base our assumptions. We then quantify the effect such glitches will have on current CW detection methods by calculating the mismatch, i.e., the relative loss of the squared signal-to-noise ratio. We do this by modeling a glitch as a piecewise Taylor-series expansion with a discontinuity at the glitch; we do not model the exponential recovery observed in some glitches, but we will discuss the significance this may have in Sec. V C. Ultimately, the goal of this work is to estimate the risk faced by current and ongoing CW searches to glitches in their target population.
In Sec. II, we will briefly describe current CW searches and how glitches may effect them. Then, in Sec. III, we investigate the statistical properties of the observed radio pulsar glitches providing fitting formulas for the glitch magnitudes and rates. In Sec. IV, we calculate the mismatch (relative loss of the squared SNR) that a single glitch will cause. Finally, in Sec. IV, we translate the observed glitches into a prediction for mismatches during a few selected current and future continuous-wave searches and discuss the risk faced by CW searches from glitches.

II. CONTINUOUS GRAVITATIONAL-WAVE SEARCHES
Searches that target a known pulsar making use of the observed EM emission (for example, the targeted search for the Crab and Vela pulsars by Aasi et al. [17]) are able to handle the epoch of a glitch, either by avoiding searches over the glitch or allowing for a jump in the timing solution at that point [18,19]. By this merit, such searches have a very low risk of being disrupted by a glitch coupled to the EM channel, provided the CW channel closely follows the phase evolution of the EM channel in between glitches.
In contrast, wide parameter-space CW searches that, by definition, search for signals without an EM counterpart do not have any such prior knowledge. This category of searches includes both directed searches, in which a single sky point in which a neutron star is believed to exist is searched (see, for example, the work by Aasi et al. [20], Wette et al. [21], and Zhu et al. [22]), and all-sky searches; in both instances, a band of frequencies and frequency derivatives are usually searched since they are inherently unknown. These searches use a matched filter against smooth templates built from a Taylor expansion in the phase; as such, they do not include glitches. If a neutron star emitting detectable levels of CW emission undergoes a glitch in the CW channel, then the matched-filtering method will not behave as expected because the template is a poor match to the real, glitching signal.
CW searches ideally employ a fully coherent search method that consists of matched filtering the template against all the data. However, such a search is typically computationally infeasible, and so CW searches use a semicoherent method: the total observation time T is divided into N seg segments of duration T seg . Each of these segments is fully coherently analyzed and then recombined incoherently to give a semicoherent measurement, which is insensitive to phase jumps between segments. This method provides more sensitive searches at fixed computing cost [23,24]. Typically, a semicoherent search is performed first; then, interesting candidates are followed up with longer coherent integration times by reducing the number of segments, aiming to eventually confirm a signal with a final fully coherent search; see the work by Shaltev and Prix [25] for a discussion of a two-stage follow-up procedure and by Papa et al. [26] for a multistage application.

III. STATISTICAL PROPERTIES OF THE OBSERVED GLITCHES
In this section, we study the properties of glitches in the observed radio pulsar population using the glitch catalog maintained by Espinoza et al. [27] and available at www.jb .man.ac.uk/pulsar/glitches.html. Our goal is to make a statement about how often glitches occur and their magnitudes for the types of neutron star which may be emitting CWs in the parameter space of typical CW searches. This task is made difficult since many searches look for young, rapidly spinning-down stars for which we only have a small sample of observations or no observations at all. Therefore, we must extrapolate the glitch properties for the population of CW-emitting neutron stars from the observed radio pulsar population.
Radio pulsar timing methods detect glitches by fitting a piecewise Taylor-series expansion in the phase to either side of the event, with a modeled jump in between (see the work by Edwards et al. [28] for a detailed discussion). The glitch catalog [27] reports 472 events from 165 isolated pulsars (as of June 27, 2016); for each of these events, a value is reported for the frequency jump δν and frequency derivative δν, if it can be measured. We cross-reference the glitch catalog with the Australia Telescope National Facility (ATNF) [29] pulsar catalog available at www.atnf.csiro.au/people/pulsar/psrcat/ in order to obtain the glitching pulsar's timing properties.
Of the 472 listed glitches, we find 15 with no ATNF cross-reference, one with δν < 0 and four with no measuredν in the ATNF catalog; these pulsars are removed from our data set. Additionally, we find 54 glitches that have either no measured δν or a measured value consistent with zero; these will be included.

A. Glitch magnitudes
Espinoza et al. [27] argued that the glitch catalog contains glitches from two distinct subpopulations of pulsars. There is the main population with δν magnitudes ranging from 10 −9 to 10 −5 Hz (which we will refer to as the "normal" glitches) and a second, smaller population with larger magnitudes of δν, referred to as "Vela-like" because the pulsars undergoing these glitches have similar characteristic ages and magnetic field strengths to the Vela pulsar (PSR B0833-45). We reproduce the evidence for this finding in Fig. 1 where we plot the histogram of all observed δν values. This illustrates the bimodality found by Espinoza et al. [27].
To check that the bimodality is not an artifact of the histogram bin sizes, we estimate the probability density function using a Gaussian kernel density estimate (KDE). Specifically, we use the Jones et al. [30] implementation. This is also plotted in Fig. 1 and shows two distinct peaks, although the lower peak could also be interpreted as two peaks close together. By eye, it is clear that there are at least two modes to the histogram and possibly more. We investigate this in Appendix A by applying a Bayesian model comparison for Gaussian mixture models (see the work by Gelman et al. [31] for a review), varying the number of components, and also allowing for a skew as described by O'Hagan and Leonard [32]. We find that all models with two or more components fit the data decisively better than a single component. Marginal gains are found by allowing the models to be skewed and have four or more components, but no single model is outstanding among the others. For this reason, we choose to use a two-component model with skew; this provides a good empirical description of the data and is pragmatic in that we limit the number of components to two for interpretability. We note that this description is empirical, and we do not intend to make any substantive claim regarding the underlying physics of the two components.
The mixture components and individual distributions for the two-component skewed model are plotted in Fig. 1, and in Table I we provide the resulting mean, standard deviation, weights, and skewness of the two components in log space. This method identifies the two subpopulations in a manner FIG. 1. The distribution of observed glitch magnitudes δν (from the glitch catalog). This is given as both a binned histogram and a Gaussian KDE, as discussed in the text. The solid and longdashed lines mark the two components of the two-component skewed Gaussian mixture model (GMM). consistent with the observations by Espinoza et al. [27], and notably the Vela-like component suffers a significant skew.
We use the best-fit two-component skewed Gaussian mixture model to label each data point as originating from one of two skewed Gaussian distributions. Specifically, to each data point, we assign the label based on the maximum probabilities of each of the two components, given the maximum posterior model parameters derived in the fitting process. In Fig. 2, we plot histograms for δν and δν along with the raw data in a scatter plot. We have separated the data into the individual subpopulations, as labeled by the two-component skewed Gaussian mixture model. Several pulsars of interest are picked out using contrasting markers. It is interesting that not all of the Vela glitches are categorized by this method as Vela-like, which can be seen by looking at the distribution of Vela glitches in Fig. 2.

B. Overview of the population of glitches
To give an overview of all observed glitches in the context of the wider population of observed radio pulsars listed in the ATNF catalog, in Fig. 3, we plot two copies of the familiar ν-ν diagram. In panel (a), for each pulsar that has been observed to glitch, we add a filled circle with an area proportional to the number of glitches seen in that pulsar. In panel (b), we mark each pulsar that has been observed to glitch with a filled circle, but here the area of the filled circles marks the pulsar's average glitch magnitude. For both plots, different shapes have been used to partition the Vela-like and normal glitches (note that some pulsars display glitches from both populations). Finally, dashed lines mark isoclines of constant characteristic age, τ age ¼ jν=νj.
While the bulk of observed glitching pulsars are from the main pulsar population, the fraction of young pulsars (τ age < 10 5 yr) that glitch is proportionally higher than in the normal population. Vela-like glitches occur predominantly in the young pulsars with none seen in pulsars with τ age > 10 7 yr. It is also noticeable that younger pulsars display a greater number of glitches. Note that, since we have not observed all pulsars for the same duration, one cannot infer the relative glitch rate from the number of glitches alone.
For the normal-glitch population, Espinoza et al. [27] noted that "pulsars with τ age < 5 × 10 3 yr undergo small-or medium-sized glitches (δν < 10 −5 Hz)." It is postulated that the higher temperatures in younger pulsars prevent the glitch mechanism from working effectively. This effect is consistent with Fig. 3(b); the pulsars with the largest average glitch sizes have τ age ∼ 10 5 yr, while younger pulsars tend to exhibit smaller glitches on average.

C. Extrapolating: Glitch magnitudes
We would like to be able to predict the glitch magnitude for the unobserved neutron star population targeted by CW searches. In particular, we need to extrapolate up to large spin-down rates −ν ∼ ð10 −9 -10 −7 Þ Hz=s searched for in many recent CW searches, which are larger than what has been observed in radio pulsars.
It has previously been found [27,[33][34][35] that the glitch activity (defined in the first of these references) correlates well with jνj and the characteristic age τ age . We choose not to combine the rate and magnitude information together into the activity but estimate both separately as these are of most direct relevance to CW searches.
We investigate correlations of the glitch magnitudes δν and δν with the frequency ν, frequency derivativeν, and characteristic age τ age , as shown in Table II. This is done for three groups: all the data together and individually for the normal population and the Vela-like population. For the normal population, both glitch magnitudes most strongly correlate with the spin-down rateν, although we recognize that the correlation with τ age is only marginally weaker. In contrast, δν for the Vela-like population has a weak correlation with all predictor variables, but δν correlates well withν and most strongly with the characteristic age. For simplicity, we choose to useν as a predictor variable for both the normal and Vela-like populations, making it simpler to interpret later results as the same predictor is used for both populations. In practice, our conclusions will be robust to either choice of predictor variable. In Fig. 4(a) and Fig. 4(b), we show scatter plots of glitch magnitudes against the spin-down rate of the pulsar to demonstrate the correlation. For both plots, we have added contrasting markers to label several interesting pulsars. These help to show that there can be almost as much variation in the glitch magnitude of a single pulsar as from the entire population.
Fitting a linear function in log-log space (see Appendix B for details), our resulting fitting formulas for the frequency jump for each of the two populations is where the last factor provides an estimate of the variability about the linear fit. These fits do not provide a precise statement about the magnitude of glitches but are sufficient to estimate the orders of magnitude that we might expect. Our data set includes one glitch from a millisecond pulsar, PSR B1821-24 [36]; it is interesting to note that the glitch magnitude is typical of normal pulsars with similar spindown rates. This is also true for the second observation of a glitch in a millisecond pulsar, PSR J0613-0200 [37] (not included in the data set as it was added to the glitch database after the analysis was completed). This indicates that, despite the much greater spin frequency of these millisecond pulsars, the glitch magnitude depends on the spin-down rate in much the same way as seen for normal pulsars.

D. Extrapolating: Average glitch rate
To estimate the average rate of glitches, Espinoza et al. [27] grouped pulsars by their spin-down rateν, including pulsars that have not yet been observed to glitch. From this grouping, the authors used the measured number of glitches N g to calculate a mean glitch rate hṄ g i. In Fig. 10 of their work, they show that, to a good approximation, in log space, the mean glitch rate depends linearly on the spin-down rate; FIG. 3. ν-ν plot of all pulsars in the ATNF catalog [29], overlaid with squares for the normal glitching pulsars and circles for the Velalike glitching pulsars. (a): For each pulsar observed to glitch, the area of a surrounding filled circle is proportional to the number of observed glitches from that pulsar. (b): For each pulsar observed to glitch, the area of the filled circle is proportional to the average glitch magnitude from that pulsar. Note that in the plot showing glitch magnitudes the relative scaling of the Vela-like and normal populations are not equal; the area representing the normal glitch magnitudes are scaled to be three times larger than the Vela-like glitch magnitudes. we reproduce this in Fig. 5 using the data from Table 4 by Espinoza et al. [27]. To extrapolate, we fit a linear function to the glitch rate and find the fitting formula whereν is measured in Hz=s. The exponent agrees with that found by the original authors (they do not provide the prefactor).

IV. SIGNAL LOSS DUE TO GLITCHES
After the discussion in the previous section of the characteristics of glitches seen in radio pulsars, we will now answer the question of what the effect of a glitch in a CW signal on matched-filtered searches is. In the next section, this will be combined with the predictions of the previous section to try and quantify the risk posed by glitches to CW searches. We begin this section by introducing the mismatch, followed by analytic and numerical estimates for the mismatch of glitching signals.

A. General mismatch definitions
We assume the template family consists of a regular CW phase model (without glitches) with phase parameters θ, defined at some reference time t ref : The mismatch μ ð0Þ is generally defined 2 as the relative loss of the (squared) signal-to-noise ratio ρ 2 due to an offset between signal parameters θ s and template parameters θ, namely, where ρ 2 ðθ s ; θÞ denotes the expected squared SNR in template θ in the presence of a signal 3 with parameters θ s and ρ s ≡ ρðθ s ; θ s Þ is the maximal "perfect-match" SNR. When considering the loss of the SNR due to a glitch (or timing noise) in the emitting neutron star, we are in the situation of signals that lie outside of the regular CW template manifold T , as originally discussed by Apostolatos [39] and Owen [40]. This effect can be quantified by the "fitting factor" FF [39] or equivalently the minimal mismatch μ and per year (right-hand axis). Black dots are the original data points, while the solid line and shaded region are our best-fit line and a measure of the variability as given in Eq. (5). 2 Obviously, using any quantity that is merely proportional to the SNR would also yield the same mismatch definition.
Evidently, if the signal lies inside the template manifold T , i.e., if there is a θ ¼ θ s ∈ T , then the minimal mismatch would be zero, while for signals outside of the template manifold, it would be nonzero. Therefore, any additional mismatch from a finitely spaced template bank will further increase the total loss of the SNR, but in the following, we are only interested in the minimal match for glitching signals, which is independent of template-bank spacing.
The metric mismatch μ is a common and useful approximation [40][41][42] to the (full) mismatch μ ð0Þ of Eq. (7) for small offsets Δθ ≡ θ − θ s between template and signal; i.e., by Taylor expanding, we obtain where g ij is referred to as the metric.
The metric mismatch of Eq. (10) is generally a good approximation for small mismatches, and empirically, one finds that deviations between μ ð0Þ and μ start to become noticeable for mismatches above μ ≳ 0.3-0.5 (e.g., see Refs. [38,43]). For larger offsets, the (unbounded) metric mismatch will tend to increasingly overestimate the actual (bounded) mismatch μ ð0Þ , which only slowly asymptotes to 1, corresponding to a total loss of the SNR. See the work by Wette [44] for a more detailed discussion of this effect and an empirical fit to extend the applicability of the metric mismatch. This distinction will be relevant in the following discussion, as glitches can result in large parameter-space offsets Δθ compared to the range of validity of the metric approximation, and therefore one needs to be careful in interpreting results obtained via the metric mismatch.

B. Glitching CW signals
The model of a signal containing N g glitches 4 is defined as a piecewise CW phase function of the form of Eq. (6) over N g þ 1 continuous "domains" labeled by a; b; … ∈ ½0; N g ; i.e., a ¼ 1 is the domain following the first glitch, etc. For each continuous stretch a, the signal phase function is therefore Þ; for t ∈ ½t ðaÞ ; t ðaþ1Þ Þ; ð11Þ and the length of each continuous stretch a is T ðaÞ ≡ t ðaþ1Þ − t ðaÞ . We further denote the fractional length of each domain as A neutron-star glitch would produce a jump δ 0 f in signal frequency and a jump δ 0ḟ in the spin-down. We will also consider a phase jump δ 0 ϕ. Such a jump in the phase is not normally considered in radio pulsar astronomy, as the uncertainty in the time t ref at which a glitch occurs is much greater than the pulsar spin period. This is related to the fact that glitches often contain initial transient jump components in frequency and spin-down that decay away on various time scales (hours to days or longer) [34]. Given that these transient glitch components are not included in our simplified glitch model, we can effectively consider the asymptotic persistent jump to also contain a phase jump δ 0 ϕ. Such phase jumps are potentially relevant to gravitational-wave searches, however, and so need to be included in our analysis.
We therefore define glitch a as a discontinuity at time t ðaÞ in signal parameters with offset fδ 0 ϕ ðaÞ ; δ 0 f ðaÞ ; δ 0ḟ ðaÞ g, while all other signal parameters (e.g., sky position ⃗n and higherorder spin-down termsf; … and binary orbital parameters) remain constant. We denote the phase parameters explicitly as θ ¼ fϕ; f;ḟ;f; …g, so the jump induced by glitch a would be δ 0 θ ðaÞ ¼ fδ 0 ϕ ðaÞ ;δ 0 f ðaÞ ; δ 0ḟ ðaÞ ; 0;0;…g. Successive jumps would be cumulative, and so we define where the cumulative effect is a simple sum because all phase parameters θ refer to the same fixed reference time t ref and are therefore constant in time (in between glitches). The offset between signal and template parameters in each stretch a is denoted as C. Single-glitch metric mismatches The general expressions for the metric mismatch for N g glitches are derived in Appendix C. Here, we only present the explicit results obtained for the special case of a single glitch, occurring after time T ð0Þ ¼ RT, where T is the total observation time and where we average over the unknown fraction R ∈ ½0; 1. These metric results are limited to the case of "directed searches" in which the SNR is maximized over a search space ff;ḟ…g, but the sky position of the source is assumed to be known.

Coherent searches
Consider three different search scenarios: a four-dimensional (4D) search covering ff;ḟ;f; f … g, a 3D search including up tof, and a 2D search covering only ff;ḟg. For these three cases, we explicitly (see Appendix C 3) find the R-averaged minimal coherent glitch mismatches as These mismatches are minimal in the sense defined by Eq. (8), i.e., the mismatches that remain once one has searched over an infinitely finely spaced template bank. These expressions illustrate the effect of correlations in parameter space: including more search dimensions lowers the minimal achievable glitch mismatch over the template bank. In particular, includingf as a search dimension compared to onlyḟ reduces the mismatch of a pure spindown jump δḟ by a factor of ∼10, a pure frequency jump δf by a factor of ∼3, and a pure phase jump δϕ by roughly ∼35%. Including both f … andf reduces the mismatch even more: a spin-down jump by a factor of 39, a frequency jump by a factor of ∼6, and a phase jump by ∼70% compared to a search up toḟ. Similarly, including sky position ⃗n as a search dimension would likely further reduce the achievable minimal mismatch as well, we investigate this in Sec. IV E 2, but there is no simple analytic way to estimate this effect due to the intrinsic complications of the sky metric [43].
Note, however, that as a consequence of this "compensation" effect there will be biases in the estimated parameters of a glitching signal, which can lead to problems when following up or interpreting such candidates without taking into account the fact that the signal has potentially undergone a glitch.
As an interesting application of these R-averaged metric expressions, we can derive maximal coherent observation times T for a given glitch jump in either f orḟ (using canonical values corresponding to the largest observed jumps) and a maximal tolerated glitch mismatch (using a canonical value of 0.1). For pure frequency jumps δf, we find for a 4D, 3D, or 2D search, respectively. In the case of a pure spin-down jump δḟ, we find

Semicoherent searches
In a semicoherent search, the total observation time T is divided into N seg segments of length T seg ¼ T=N seg , each segment l is analyzed coherently, and the results are combined incoherently by summing the "power" from each segment (e.g., see the work by Prix and Shaltev [23] and Brady and Creighton [24]).
For simplicity, we assume the glitch to happen on a segment boundary, after N ð0Þ ¼ RN seg segments. If a glitch happens inside a segment, the additional mismatch in that segment is bounded within [0, 1] and at most contributes 1=N seg to the total mismatch.
From the detailed derivations in Sec. C 4, we can obtain explicit R-averaged glitch mismatch expressions in the limit of a large number of segments N seg in the three different search spaces (4D ¼ ff;ḟ;f;fg, 3D ¼ ff;ḟ;fg and 2D ¼ ff;ḟg), namely, Here, we see again how parameter-space correlations act to reduce the minimal-glitch mismatch when including more search dimensions. Note that the phase jump δϕ does not appear in these expressions, as in a semicoherent search the SNR in each segment is separately maximized over the phase of the signal within that segment, so that phase jumps in between segments have no effect. Turning these R-averaged expressions into bounds on the coherent segment length for a pure frequency jump, we find while for a pure spin-down jump, we obtain D. Upper bounds on glitch mismatches

Coherent searches
The metric mismatches discussed in the previous section can be a useful tool, for example, to estimate bounds on the maximal safe observation span in the presence of a glitch, as given in Eqs. (16)- (18) and Eqs. (21)- (23). However, the metric is of somewhat limited use for predicting mismatches from glitches in general because the parameter jumps δθ will often not be small enough for the metric approximation of Eq. (C12) to be applicable. In fact, it is the regions of "large" mismatches ∼Oð1Þ that would be most relevant for a realistic assessment of the impact of glitches on the detectability of signals. While numerical extrapolations to the metric mismatch do exist for large template mismatches [44], it is unclear to what extent they would be applicable in the present case of off-manifold signal mismatches. However, we can also obtain simple analytic upper bounds on the mismatches, which can be useful to estimate the relevant "scale" of the expected mismatches.
Namely, we can always find a template that perfectly matches the signal over the longest time stretch a 0 , i.e., R ða 0 Þ ¼ max a R ðaÞ , during which the difference between the template and signal is zero. While it is in principle possible for other time stretches a ≠ a 0 to reduce the total SNR by noticeable negative contributions in the (complex) coherent matched-filtering amplitude [cf. Eq. (C9)], this is generally quite unlikely: offsets in frequency Δf ðaÞ or spin-downs Δḟ ðaÞ would result in a rapidly oscillatory matched-filtering integrand e iΔΦðtÞ [where ΔΦðtÞ is the difference between the signal phase and the template phase defined in Eq. (C8)], and so the corresponding contributions would typically be small, while only pure phase jumps Δϕ ðaÞ ∼ π could result in canceling contributions. Therefore, generally, the matched-filtering amplitude of Eq. (C11) will be X ¼ max a R ðaÞ (assuming the template is perfectly matched in this domain but contributes neither positively nor negatively elsewhere). And so we obtain an upper bound on the coherent glitch mismatch as where the second inequality is obtained by observing that max a R ðaÞ ≥ 1 N g þ1 .

Semicoherent searches
By regrouping Eq. (C39) over the interglitch time stretches, we can write the semicoherent glitch mismatch aŝ ðaÞ : Using the same argument as before, there will always be a template that perfectly matches the signal over the longest interglitch time stretch max a R ðaÞ T, where the mismatch would be zero, while it would be at most 1 the rest of time (but quite possibly less). This yields the (strict) upper glitch mismatch bound of which, contrary to the coherent case of Eq. (26) is linear in max a R ðaÞ instead of quadratic.

E. Numerical estimates of the mismatch
The metric estimates discussed in Sec. IV C describe the behavior of small glitch mismatches, while the results of Sec. IV D provide useful analytic upper bounds. To verify these expressions as well as interpolate between the two regimes of small and maximal glitch mismatches, we also use direct numerical computation of the mismatch. In this section, we present two different approaches to numerically estimate the mismatch: (i) a relatively simpler numerical evaluation of the simplified matched-filter amplitude of Appendix C 1, with a minimization over the template parameters, and (ii) direct evaluation of the full CW F statistic [45] over the template space using an efficient Markov-chain Monte Carlo (MCMC) search.
Understanding the intermediate regime between metric estimates and upper mismatch bounds will be important when quantifying the risk posed by glitches to real searches, as glitches will often be sufficiently large to cause substantial mismatches (≳0.3), as was seen already from the simple estimates in Eqs. (19)- (20) and Eqs. (24)-(25).

Mismatch from simplified matched-filter amplitude
Given a particular glitching signal, we can estimate the mismatch by numerically evaluating the simplified matchedfiltering amplitude X defined in Eq. (C8) over the template manifold. This approach therefore only involves the differences between the template phase model and the signal phase model, while still neglecting antenna-pattern and polarization parameters.
Calculating the fully coherent mismatch involves generating two time series, namely, the signal-and template-phase functions. The phase difference ΔΦðtÞ ¼ Φ s ðtÞ − Φ t ðtÞ is then used to compute the simplified matched-filtering amplitude X, and the fully coherent mismatch is estimated as 1 − jXj 2 , as discussed in Appendix C 1, and numerically minimized over the template search parameters using a Nelder-Mead minimization (the Jones et al. [30] implementation), which was found to be effective for small mismatches. For larger mismatches, however, the mismatch topology contains multiple minima, and this method requires careful selection of the initial guess.
In Fig. 6, we show the results of a simple Monte Carlo study using this approach in the context of a semicoherent search. For each simulation, we generate a signal lasting 100 days with a single fixed glitch of magnitude δf ¼ 5 × 10 −7 Hz, and we choose the time at which the glitch occurs uniformly over the entire data span. In the upper plot, we show the average minimal mismatch of a 2D search over ff;ḟg as a function of the coherent segment length T seg of the search (and hence the number of segments), comparing the result to the metric prediction of Eq. (23). The fractional difference between the two estimates is shown in the lower plot. For short segment lengths, the mismatch is small hμi ≲ 0.1, and so the numerical estimates agree well with the metric mismatch approximation. However, as the segment length increases, the mismatch grows, and the approximation starts to overestimate the numerical results.
In addition to verifying the behavior of the metric mismatch approximation, Fig. 6 also demonstrates that for a fairly typical glitch size (compare with Fig. 2, for example) 10% of the squared SNR can be lost for longer segments. Moreover, this also provides insight into what will happen to a signal during a follow-up procedure when the segment length is increased to test the significance of candidates. For short segment lengths, only a small amount of the SNR is lost, but this mismatch will increase as the segment length is increased. Therefore, the SNR will not increase as it is expected to for a CW signal, and the candidate might potentially be classified as "not following the presumed signal model."

Full F -statistic mismatch using a MCMC search
The second approach consists in directly computing the mismatch of Eq. (7) by evaluating the F statistic over the template manifold using a MCMC search. For this, we generate data containing a glitching signal (using lalapps_Makefakedata_v5 [46]) and then search for it over the template manifold (not including glitches) using the LALSuite [46] implementation of the F statistic. The mismatch is then obtained as the fractional difference between the recovered maximal F -statistic value and a perfectly matched signal (including the glitches). This method allows one to include the effect of a sky search (that is, to answer the question of whether searching over the sky position can further reduce the minimal-glitch mismatch similar to what was found when searching over higher-order spin-downs).
The F statistic is the log-likelihood ratio of the signal vs Gaussian noise, analytically maximized over the four amplitude parameters h 0 , cos ι ∈ ½−1; 1, ϕ ∈ ½0; 2π, and ψ ∈ ½−π=4; π=4, where ι is the inclination angle of the source, ϕ is the initial phase, and ψ is the polarization angle. We compute the mismatch averaged over these amplitude parameters by randomly drawing these amplitude parameters from their prior ranges. Note that we fix the amplitude of the signal to h 0 ¼ 1 and generate data without noise in order to obtain the expected F statistic, which is directly related to the SNRρ via E½2F ¼ 4 þρ 2 [45], which enters the mismatch expression Eq. (7).
To minimize the glitch mismatch, we use a MCMC minimization step with priors for each of the search parameters chosen in order to allow all correlations to be fully explored. A gridded search would yield equivalent results, provided the grid points were sufficiently dense such that the template-bank mismatch was negligible; a MCMC approach (which must equivalently be run for a long enough period to ensure a good approximation of the global maximum) was found to be simpler and computationally less demanding than a gridded search.
As a simple demonstration of this method, in Fig. 7, we show the results of a Monte Carlo study of a semicoherent  7. Semicoherent mismatch as a function of glitch frequency jump δf for three different parameter spaces searched (2D=ff;ḟg, 3D=ff;ḟ;fg, and "2D þ sky"). In this example, the search setup is fixed (T ¼ 60 days and T seg ¼ 15 days), as is the time of the glitch, which occurs 1=4 of the way through the observation. Thick lines denote the mismatch computed using the F -statistic method of Sec. IV E 2. The horizontal line marks the mismatch upper bound at 0.25, since max a R ðaÞ ¼ 0.75. Black dashed lines indicate the metric mismatch estimates of Appendix C 3 for the 2D and 3D cases. search (with T ¼ 60 days and T seg ¼ 15 days) for a signal containing a single glitch which occurs exactly one-quarter of the way through the search; this setup is chosen such that in the large-glitch limit the maximum mismatch will be exactly 0.25 where the template matches three-fourths segments (see the discussion in Sec. IV D for more details). We plot the mismatch (averaged over amplitude parameters) vs a realistic range of the jump size δf in frequency (e.g., see Fig. 2).
This search is repeated for three different parameter spaces: a 2D search over ff;ḟg, a ð2 þ 2ÞD search over ff;ḟ; skyg, and a 3D search over ff;ḟ;fg. In the case of searching over the sky, we found that the minimum mismatches were found up to ∼0.1 rad away from the true sky position, which indicates that glitches can induce substantial biases in the recovered parameters of a search over a nonglitch template bank. The plot also shows the metric mismatch predictions for comparison, as derived in Appendix C 4, without taking the large-N seg limit (the exact expressions are lengthy and not generally useful and have therefore been omitted).
This confirms and extends the findings of Sec. IV C that searching over more parameters can yield a reduction in mismatch, which also applies to searching over the sky. The glitch mismatch reduction due to minimization over the sky is found to be somewhat smaller but comparable to searching over an extra spin-down parameterf.
While the explicit F -statistic search requires greater computational resources than the simplified matchedfiltering approach of the previous section, it models the effects of realistic searches most fully and will therefore be used as the method of choice to estimate the risk to real CW searches in the next section.

V. RISK POSED BY GLITCHES TO CW SEARCHES
In this paper so far, we have provided an empirical study of glitches in the radio pulsar population (Sec. III) and a discussion on how one or more glitches in the CW signal will produce a nonzero mismatch in a matched-filtered gravitational-wave search (Sec. IV). We will now combine these two analyses to assess the possible impact of glitches on some past and possible future CW searches.
Our crucial assumption is that the statistical properties of glitches in the known pulsar population are a good indicator of those of the target population of gravitational-wave emitters. To be more precise, the fitting formulas of Sec. III, for glitch size and rate which used the spin-down rateν as the indicator variable. For the results to be presented here, assessing the possible impact on CW searches, we will need to extrapolate up to spin-down rates of the order of 10 −7 Hz=s. This is about 2 orders of magnitude larger than the largest spin-down rates seen in the radio pulsars; see Figs. 4 and 5.
Given that the physical mechanism(s) underlying glitches are not well understood, it is difficult to assess the safeness of this assumption, although it seems reasonable to suppose that, whatever the mechanisms are, they will apply to all spinning-down neutron stars (but see the discussion of gravitars, below). Of course, the eventual observation and statistical characterization of glitches in a population of such rapidly spinning-down CW sources would provide a test of this assumption. Note that for spin-down rates with jνj > 10 −7 Hz, the fitting formulas for the glitch sizes δν for the normal and Vela-like populations cross, casting added doubt over their validity at such high spin-down rates; see Fig. 4. The same is true for the δν jumps. Fortunately, we will not need to extrapolate into this regime.
Gravitational-wave astronomers sometimes postulate the existence of a population of gravitars, electromagnetically unseen neutron stars that spin down mainly through gravitational-wave emission [47,48]. Such stars would necessarily have very low external magnetic fields, or else they would be visible as pulsars and would be acted upon by significant electromagnetic spin-down torques.
If the internal magnetic field is also small, one might wonder if the glitch mechanism might somehow be suppressed, as compared to the radio pulsar population. We can simply note that for glitches due to superfluid unpinning, if the pinning takes place in the inner crust, the large-scale stellar magnetic fields probably play no significant role [2], so such glitches are still to be expected in gravitars. On the other hand, if, as has been suggested relatively recently [49], pinning takes place in the core, on magnetic flux tubes, the lack of a significant internal magnetic field might indeed suppress this glitch mechanism. For crust quakes, which can be expected to occur alongside whatever superfluid glitch mechanism might be operative, the magnetic field is not expected to play any significant role [4], implying that crust quakes should occur also in gravitars.
To sum up, we can reasonably expect the glitch mechanisms that apply to the radio pulsar population to apply in whole or (at least) in part to the hypothetical gravitar population.
A. Glitch rate and associated probability CW searches typically use stretches of data from tens of days to a few months with a small number of searches spanning longer than a year or two. The first question that must be answered to understand the risk is as follows: given that a signal does exist in the data, how probable is it that one or more glitches will occur during the data span?
In Eq. (5), we reproduced the glitch-rate fitting formula for hṄ g i given by Espinoza et al. [27], which provides an estimate of the glitch rate per second as a function of the source pulsar's spin-down rate and includes the effect of pulsars that have not been observed to glitch. From this rate, we can estimate the expected number of glitches given a span of data as a function of the spin-down rate. For the glitch magnitude fitting formulas, we split the population of glitches into two subpopulations (normal and Vela-like) to avoid overestimating the glitch magnitude for large spin-down rates. To use these fitting formulas, we will present results in this section similarly split by subpopulation. However, the fitting formulas for the glitch rate were calculated from the whole population; therefore we define hN g i normal ðν; TÞ ¼ w normal hṄ g iT ð29Þ as the expected number of normal and Vela-like glitches in which w normal and w Vela-like are the weights of the two populations as given in Table I. Implicit in this definition is a prior specification that the proportion of normal and Vela-like pulsars in the target population is the same as in the observed population. There is some evidence that, in fact, the proportion of Vela-like pulsars increases withν; this could be modeled by aν-dependent weighting, but we will ignore this effect here. The average number of glitches in a given search is a useful quantity, but it is not easy to interpret; for a low average number of glitches, there remains a significant probability of having zero glitches and hence no loss of signal. To better understand this risk, we will therefore apply a simple substantive model, a Poisson process. Melatos et al. [50] demonstrated that glitch waiting times are consistent with an avalanche process transferring angular momentum from the core superfluid to the crust. Choosing nine pulsars that had glitched five times or more, they found that seven of these were consistent with a constant-rate Poisson process such that each glitch event was statistically independent. In the remaining two, PSR J0537-6910 and PSR B0833-45 (Vela), they found that a quasiperiodic component (i.e., the glitches occur quasiperiodically) coexists with the Poisson process and accounts for about 20% of the events; this is suggestive that these periodic glitches originate from a different mechanism.
Of all events in the glitch catalog, PSR J0537-6910 accounts for 23 and PSR B0833-45 (Vela) for 17 of the total 472 events. Assuming that 20% of these are due to the quasiperiodic component, this is ∼1.7% of the total number of observed glitches. It is possible that other pulsars also exhibit a quasiperiodic component, so the total fraction of glitches from a quasiperiodic component may be larger than 1.7%, but it seems likely that a Poisson-like process should provide a good description of the probability of glitches occurring in general.
Assuming the Poisson process is responsible for all the glitches in the catalog, we can calculate the probability of one or more glitches occurring given the expected number of glitches. To do this, we take the estimated number of glitches during a typical search λ and sum the Poisson probability mass function from 1 to infinity: Note that, in practice, we truncate the summation at a finite level where the mass function is negligible. We will apply this separately to the normal and Vela-like populations with λ ¼ hN g i normal and λ ¼ hN g i Vela-like . The total expected number of glitches is the sum of the individual expectations, and from this, the total probability could be calculated. In Fig. 8, we illustrate how this probability varies with the data span for four fixed spin-down rates. This shows that for the normal population, if the data span is shorter than a month or so, glitches are unlikely even for large (absolute) spin-down rates. However, for more typical data spans ≳100 days, the probability of seeing one glitch rises above one-half with the effect correlating with the magnitude of the spin-down. The same general pattern is found by the Vela-like population, with the effect being marginally weaker due to our prior weighting.

B. Past and future searches
As discussed, neutron stars with large spin-down rates will have larger and more frequent glitches than those with smaller spin-down rates; see Figs. 4, 5, and 8. It follows that gravitational-wave searches that search over larger spin-down ranges and over longer spans of data will suffer greater mismatches due to glitches than searches over smaller ranges in the spin-down rate and over shorter observation spans. In this section, we aim to quantify the potential mismatch for some recent and ongoing CW searches.
For each CW search that we consider, we estimate the probability of glitches and the mismatch if one or more glitches did occur. We evaluate these quantities for signals with the largest spin-down rate (i.e., largest value of −ν) considered in each search. In this respect, it represents a worst-case scenario, but we note that the results would not change much if we had instead averaged over the spindown rate of the search; averaging uniformly over spindown would still correspond to the scale set by the largest spin-down rate. These numbers are therefore still qualitatively representative for the "expected" effects over the whole search space. For each search listed in Table III, we summarize the relevant search parameters followed by our estimates for the effect of glitches in this search. Splitting the results up by their source population (such that to get the total expected number of glitches one adds the normal and Vela-like expectations), we list the expected number of glitches, the probability for one or more glitches, and the semicoherent and fully coherent mismatch expectations if one or more glitches did occur.
These mismatches have been computed using an F statistic-based MCMC search as described in Sec. IV E 2, while the glitch parameters have been sampled as follows: given the expected number of glitches computed at the largest spin-down rate, we draw the actual number of glitches from a Poisson distribution, cf. Eq. (31). If the number of glitches is zero, the mismatch would be zero (since we assume an infinitely fine template bank). In the case of nonzero glitches, we pick the glitch times uniformly throughout the search span. For each glitch, we draw glitch sizes based on a Gaussian distribution with a mean and standard deviation (in log space) obtained from the fitting formulas in Eqs. (1)- (4). For each such glitching signal, we then numerically compute the mismatch using the F -statistic MCMC search described in Sec. IV E 2 over the search dimensions of the original search. Subsequently, we average these mismatches (excluding cases where no glitches occur) to produce the expected mismatch numbers given in Table III. In brackets, next to the averaged mismatch, we also provide the average upper bound of Eqs. (26) and (28), respectively. This number is computed directly from hN g i as follows. First draw a large sample of random Poisson variables N i g ∼ PoissonðhN g iÞ. For each instance N i g , generate a realization of the fractional time at which the glitches occur (i.e., draw N i g uniform random variables from [0, 1]), and record the maximum interglitch fractional duration max a R i ðaÞ , which fully determines the mismatch upper bounds. We note that for many cases in this table the   TABLE III. Predicted effect of glitches on some recent and planned CW searches; note that the "S6 E@H all-sky HFU Ã " search refers to stage 4 of the hierarchical follow-up. For each search, we summarize the search parameters, followed by estimates for the number and probability of glitches and the resulting averaged semicoherent (hμ ð0Þ iÞ and fully coherent (hμ ð0Þ i) F -statistic mismatches (see the text for details). These quantities are calculated at the largest spin-down rate considered in the search. We also provide the average value of the corresponding mismatch upper bounds (given in brackets); see the text for details.
S6 E@H all sky average expected mismatch is close to and in some cases even larger than the average upper bound. In the cases in which it is larger than the upper bound, this indicates that our MCMC minimization failed to always find the global minimum. Upon inspection, it appears that this happens most when there is a large number of glitches (as is the case for the S6 E@H Cassiopeia A search and O1 E@H Multidirected Cassiopeia A). In the cases in which it is close to the upper bound, we conclude that the glitch sizes are so large that the mismatch often saturates at the upper bounds. This table shows that for these searches considered there is both a substantial probability of glitches occurring, and these glitches would result in a significant fraction of the SNR being lost. The table also reiterates one of the earlier findings shown in Fig. 6: the mismatches will increase in the follow-up. In this case, we have estimates of the mismatch during the initial semicoherent stage and the mismatch if followed up fully coherently. Taking as an example the S6 E@H all-sky search, during the initial stage, only 0.16 of the squared SNR would be lost, and hence the signal (if it where strong enough) might be classified as a candidate. However, if immediately followed up fully coherently, the mismatch may increase to 0.40 such that the signal would appear to be weaker than it should be if it had matched the presumed signal model.

C. Including the recovery from glitches
In this work, we have used the glitch catalog [27], which provides δν and δν, estimates of the change in rotation frequency and spin-down rate of the pulsar at the glitch. However, in addition to this instantaneous behavior, some pulsars also undergo a short-term exponential relaxation of some fraction Q of the total glitch magnitude over time scales τ d that are typically tens to hundreds of days [34]. This may have an important effect on our estimates since, if a large fraction of the glitch is recovered in a time scale that is short compared to the observation time, we will overestimate the mismatch. On the other hand, if not explicitly included in the search template, the relaxation itself could also cause a mismatch, which would tend to vanish when τ d ≪ T but be maximal when τ d ∼ T.
One can ask what the available pulsar data tell us about the likely values of Q. Unfortunately, the issue of extracting the recovery parameter Q from pulsar data is a subtle one. As noted by Ref. [27] in relation to the glitch catalog, "The results presented in this paper do not involve fitting of short-term recoveries because their parameters depend so critically upon the usually poorly known glitch epoch." Furthermore, if the recovery time scale is short compared to the baseline on which pulsar timings are constructed, the recovery will not be apparent and will effectively be absorbed into the measured δν values. If, instead, the recovery is on a time scale that is long compared to the baseline on which the pulsar timings are made (but not so long that one cannot see the recovery by connecting different timing solutions), then the recovery will be visible in the data, and recovery values Q will be estimated.
Such procedures have been carried out, with recovery values quoted for many glitches, e.g., in the work by Lyne et al. [34] and Yu et al. [53]. The results of Lyne et al. [34] pointed to a correlation between Q and the spin-down rate, with the largest recovery fractions being found in the pulsars with the largest spin-down rates (see their Fig. 6). However, the more recent data of Yu et al. [53] point to a more complex picture, with the largest glitches having values of Q from ∼10 −3 up to Q ∼ 1.
Furthermore, a recent detailed study by Shannon et al. [54] concluded that eight glitches in the Vela pulsar are best described by permanent and transient changes in the frequency alone, with a dominant decay time scale of 1300 days and Q ∼ 0.3-0.8. If such findings are replicated in many other pulsars, we might need to revisit the modeling of glitches used here, including the transient component but removing the jumps inν; however, in doing so, we would not expect drastic changes to our overall conclusions, since the decay time is long compared to typical observation spans. Given that the glitch catalog modeled glitches with permanent offsets in frequency and derivative only, the analysis performed here is appropriate and consistent with the current literature.
Currently, it is unclear what Q and decay time scale we might expect for the target population of all-sky and directed gravitational-wave searches. It is certainly possible that in the stars of interest to us the recovery fraction may be large (Q ∼ 1) so that the results of Table III overestimate the effect of glitches on CW searches. However, on the basis of the available evidence, it is not possible to say more than this. The only way of settling the issue definitively would be to allow for relaxation in CW templates and see if this better explains the gravitational-wave data.

VI. CONCLUSIONS
This work investigates the effects of glitches in the CW signal on searches that matched filter the data against nonglitching templates. We started with an initial study developing empirical fitting formulas for the size and rate of glitches in the pulsar population (as cataloged by Espinoza et al. [27]). Subsequently, we developed different ways to estimate and evaluate the loss of SNR in fully coherent and semicoherent searches due to the presence of glitches. Finally, we used our fitting formulas to predict the loss of SNR for typical CW searches and the probability of one or more glitches occurring.
The work is developed with three motivations: (i) to understand if glitches pose a risk to our ability to detect CWs from isolated neutron stars in wide-parameter space searches, (ii) to help guide decisions about planned future searches in order to minimize any risk, and (iii) to make a statement on how the risk manifests and what can be done to mitigate it.
With regard to the first motivation, it is clear from Sec. V that glitches, if they occur in CW signals with properties similar to those observed in pulsars, can both be frequent enough and of a sufficient magnitude to result in large losses of SNR for both semicoherent and fully coherent searches. This is particularly true for searches spanning long stretches of data and searching large spin-down rates, two attributes usually associated with increased detectability.
The main uncertainty in our analysis lies in extrapolating glitch sizes and rates from the known pulsar population to regions of large spin-down rates appropriate to CW searches. In the results presented here, we have had to extrapolate by about 2 orders of magnitude in the spindown rateν. We have also not attempted to account for postglitch relaxation in our analysis, as there is not a sufficiently clear pattern of how the amount of relaxation correlates with other quantities. We simply note that glitch recovery could serve to significantly reduce the impact of glitches on gravitational-wave searches, if the amount of relaxation is sufficiently large and if it occurs on time scales of relevance to CW searches.
For semicoherent searches, the impact of glitches is reduced by using shorter segment lengths. However, this loss of the SNR is not the primary risk; even if a glitching signal has been identified as a candidate in the initial wideparameter space search, the greater risk potentially lies in the follow-up procedure of these candidates, which is often considered the test for a real signal [26]. During the followup, for standard CW signals, it is expected that the squared SNR grows linearly with increasing segment length [e.g., see Eq. (C7)]. However, for a glitching signal (as seen in Fig. 6) the mismatch will also grow, and so the SNR will not increase in the expected way, and the candidate may potentially be dismissed. This is of concern to both future and past searches for CWs from neutron stars. If the effect of glitches is ignored, detectable signals could easily be missed due to the presence of glitches. We therefore recommend that the setup design of semicoherent searches and follow-ups take the possibility of glitching signals into account, especially when searching at large spin-down rates. Furthermore, more work is needed to investigate and implement modifications to current follow-up procedures in order to account for the possibility of glitches and ideally be able to fully localize them in CW signals.
ACKNOWLEDGMENTS D. I. J. acknowledges support from STFC and also travel support from NewCompStar (a COST-funded Research Networking Programme). All authors are grateful to Christobal Espinoza for help with using the glitch catalog and Danai Antonopoulou, Wynn Ho, Graham Woan, Matt Pitkin, and members of the Continuous Waves group of the LIGO Scientific Collaboration and the Virgo Scientific Collaboration for useful feedback and discussion. This article has been assigned the document number LIGO-P1600201.

APPENDIX A: BAYESIAN MODEL COMPARISON: TEST OF MIXTURE MODELS
It seems clear by eye that the histogrammed magnitudes of the frequency change in a glitch, log 10 jδνj, as shown in Fig. 1, exhibit at least two distinct modes. To model this empirically, we will use a Gaussian mixture model (GMM) [31] with N components. This model assumes that the measured data are taken from a population with N subpopulations, each having a Gaussian distribution with separate mean, variance, and weight (μ i , σ 2 i , ω i ), where i ∈ ½1; N; note that P N 1 ω i ¼ 1. Furthermore, we can also allow each of the components to be skewed with a dimensionless skew parameter α i , which can be either positive or negative determining the direction of the skew or 0, for which there is no skew. Following O'Hagan and Leonard [32], then, the probability density function of the ith skewed Gaussian component is where N denotes the Gaussian distribution.
Let y be the set of measured values of log 10 jδνj and ϑ ¼ fμ i ; σ i ; α i ; ω i g be the collection of all model parameters. Then, the probability density for a GMM with N components is To compare different choices of N, we will perform a Bayesian model comparison [55] between each of the mixture models and the simplest hypothesis, a mixture model with N ¼ 1.
For each model parameter, we must specify a prior. We list these in Eq. (A3), having defined hyi, jyj, and stdðyÞ as the average, range, and standard deviation of the data, For the mean μ i , we use a uniform prior over a range of values containing all data points. For the standard deviation σ i , we will use a half-Cauchy distribution with zero mean as suggested by Gelman et al. [56]. A large standard deviation, as compared to the standard deviation of the data themselves, provides a weakly informative prior. Instead, we use a standard deviation of stdðyÞ 2 to favor GMM components with small standard deviations as compared to the data. That is, our prior disfavors models in which any of the components are wide and flat. The prior for ω i is uniform on [0, 1] and for α i is normally distributed with zero mean and a wide, weakly informative standard deviation. The choice of a zero mean favors nonskewed components. Note that the nonskewed models do not include α i as a model parameter and the GMM with N ¼ 1 does not include ω i . We use this choice of prior for the model parameters of each component in the GMM with N components. In this way, models with larger values of N have a larger "prior volume," and hence there is a natural Occam factor favoring the simpler models with fewer components; this prevents overfitting.
We will present results for the Bayes factor between a GMM with N components and the simplest model, a GMM with N ¼ 1 components. This is computed by We use the EMCEE [57] MCMC algorithm to sample from the posterior and thermodynamic integration to estimate the evidence integrals [58]. In Table IV, we provide the log 10 of the Bayes factor for several possible models. The Bayes factor between any two of the models given in Table IV can be calculated from their difference.
This table clearly shows that the data are decisive; a Gaussian mixture model with N ≥ 2 fits the data a great deal better than the simple N ¼ 1 GMM. This is unsurprising, given the distinct multimodal nature of the data. However, the differences between the other models is more subtle. No single model distinguishes itself by a decisive odds ratio compared to its neighboring models. We have checked that these results are robust to small changes in the prior specification.
To help illustrate the differences between these models, in Fig. 9, we plot the probability density for the maximum posterior model parameters found in a few selected models. It is clear from these plots that the N ¼ 2 model does not explain the number of glitches found in between the two primary subpopulations around log 10 jδν=Hzj ¼ −5.5; by comparison, the N > 2 models and the N ¼ 2 model that allows for skewness can explain these points, and this is reflected in the Bayes factor.
From this analysis, it is difficult to decide which model best fits the data. However, what is clear is that simply modeling the data as a GMM with two components with a skew provides a reasonable empirical model. For this reason, in our analysis of the glitch population, we will use this model and not any of the models with a greater number of components.
It is important to realize that this comparison is purely empirical, in that the result was not conditioned on a substantive physical model. It would be interesting to include such modeling; this may provide some insight into the appropriateness of the mixture model and the number of components.

APPENDIX B: LINEAR REGRESSION IN LOG SPACE
In Sec. III, we perform several linear regressions in log space in order to calculate power-law fits. This assumes that the observed values logðy i Þ depend on the predictor values logðx i Þ as where the ϵ i are independent and identically central normally distributed variables with a standard deviation σ. In this way, m and c are the linear fit free variables, while σ is a measure of the variability in the observation about this linear fit. We use a Bayesian linear regression in which we estimate the posterior distributions of all three parameters using a Markov-chain Monte Carlo algorithm; for the prior distributions, we use noninformative priors and test that these do not induce any bias. In all cases, we find the  9. The distribution of glitch sizes in frequency jump along with the predictions for the components of several GMM; the Bayes factor for these models can be found in Table IV. resulting posteriors to be Gaussian and so can take their mean values to get best-fit parameters. The advantage of this method compared to a simple least-squares linear regression is that we also estimate hσi, the variation about the linear fit. The linear fit can therefore be written as logðyðxÞÞ ¼ hmi logðxÞ þ hci AE hσi: ðB2Þ We can then rearrange this equation to give the corresponding power-law fit in linear space, where the last term gives the variability about the mean. Hence, neglecting this term gives the mean. This is an inherently problematic approach since many functions besides a power law can appear linear in a log-log plot and the assumption of Gaussian error may not be a good description. Nevertheless, we will still apply this approach since we need only order-of-magnitude estimates and can always check our predictions; we must be clear that the power-law fit gives a good empirical fit but is not intended to signify any substantive underlying model.

Coherent matched-filtering approximation
To study glitch mismatches, it is useful to review the (simplified) matched-filtering amplitude introduced by Prix and Itoh [59], from which the metric can easily be derived. The coherent matched-filter scalar product for narrowband, long-lasting signals can be written [45] as in terms of the time average where t 0 denotes the start time and T denotes the duration of the data span analyzed. The log-likelihood ratio ln Λ of the signal model vs Gaussian noise is found as (e.g., see Ref. [ We note that the dominant term in the phase function is ΦðtÞ ∝ 2πft, where typical signal frequencies in ground-based detectors are f ≳ 10 Hz and typical observation times T are of order of a day or longer, which results in Assuming the data to be dominated by a signal (neglecting noise), i.e., xðtÞ≈sðtÞ¼A s cosΦðt−t ref ;θ s Þ, and maximizing the log-likelihood ratio ln Λ over the template amplitude A [but not yet over the initial phase ϕ ≡ Φð0; θÞ, which, contrary to the standard treatment, is left as an explicit template parameter for now] yields the (partial) maximumlikelihood expression with the perfect-match signal squared SNR defined as and the (coherent) "matched-filtering amplitude" X defined as Xðθ s ; θÞ ≡ he iΔΦðθ s ;θÞ i; ðC8Þ with jXj 2 ≤ 1 and where ΔΦ ≡ Φðt; θ s Þ − Φðt; θÞ. Note that, contrary to the discussion by Prix and Itoh [59], here X still depends on the initial-phase parameter ϕ, and one can verify that maximizing X over ϕ yields again the expression found by Prix and Itoh [59], namely, ln Λ ML A;ϕ ¼ 1 2 ρ 2 s jXj 2 . From Eq. (C6), it is natural to define the "mismatched (ϕ-coherent) SNR"ρðθ s ; θÞ as Specializing this to the glitching signal model of Eq. (11), we obtain This "ϕ-coherent" form of the phase metric differs from the standard "coherent" expression (e.g., see the work by Prix [38] and Brady et al. [41]) as it still covers the initial-phase parameter ϕ. Minimizing the ϕ-coherent mismatch over ϕ, one can easily recover the usual form of the standard coherent metric, namely, where indices i; j ¼ 1; … run over the phase parameters excluding ϕ, while indices α, β label all phase parameters including ϕ. Substituting Eq. (C12) into Eqs. (C10) and (C11) and keeping only terms up to second order in Δθ ðaÞ , we obtain the ϕ-coherent metric mismatch approximation as ðaÞ αβ Δθ α ðaÞ Δθ β ðaÞ : ðC15Þ

Minimal coherent metric mismatch
We can use the general ϕ-coherent multiglitch mismatch expression of Eq. (C15) to express the minimal mismatch over an infinitely fine template bank of Eq. (8). For this purpose, we parametrize the offsets Δθ ðaÞ of Eq. (15) between the signal and template in terms of the offset in the first continuous stretch a ¼ 0 before the first glitch; i.e., we write Δθ ð0Þ ≡ Δθ;

Coherent metric mismatch for a single glitch
In the special case of a single glitch N g ¼ 1, we write δθ ð1Þ ≡ δθ and R ≡ R ð0Þ ¼ T ð0Þ =T and therefore R ð1Þ ¼ 1 − R, and T ð1Þ ¼ ð1 − RÞT, and so Eqs. (C20) and (C21) yieldμ min ¼ δθ αG αβ δθ β ; with ðC22Þ and assuming the glitch occurs with uniform probability in R ∈ ½0; 1, the expected mismatch for a single glitch can be found as To obtain an explicit estimate for the minimal expected glitch mismatch of Eq. (C24), we consider the case of a "directed" search for a target with known sky position and binary-orbital parameters, where the search space only includes frequency f and spin-downs fḟ;f; …g. Note, however, that generally there will be metric correlations also between sky-position or binary-orbital phase parameters and the set of "glitch parameters" fϕ; f;ḟg, which can result in further reductions in the minimal mismatch of Eq. (C20) compared to the example given here. Assuming a directed search, the phase model is simply where Δt ¼ t − t ref is the offset from the reference time t ref at which the phase parameters fϕ; f;ḟ;f; …g are defined. It is easy to see that changing the reference time to t 0 with the reference-time shift operator T ðτÞ given by For evaluating the mismatch of Eq. (C24), however, it will be most convenient to choose the reference time at the time of the glitch, i.e., t ref ¼ t ð1Þ , such that δθ refers directly to the instantaneous changes in parameters at time t ð1Þ . We can compute the corresponding metric simply by appropriately shifting the reference time of Eq. (C30), namely, g ð0Þ ðt ð1Þ Þ ¼ T ⊺ ðτ 0 Þ ·g ð0Þ mid · T ðτ 0 Þ; g ð1Þ ðt ð1Þ Þ ¼ T ⊺ ðτ 1 Þ ·g ð1Þ mid · T ðτ 1 Þ ð C31Þ with Combining these expressions and substituting into Eq. (C23) yields the following R-dependent expressions for the mismatch: where we assumed a constant per-segment perfect-match SNRρ s . By Taylor expanding in small offsets λ ¼ λ s þ Δλ around the signal location λ s , using Eq. (12) and dropping terms of higher order, we obtain the incoherent metric mismatch approximation aŝ whereg l ij is the usual coherent metric (i.e., minimized mismatch over ϕ) for segment l, as given in Eq. (C14).
We consider again a signal that undergoes N g glitches at times t ðaÞ with a ¼ 1; …N g , and for simplicity, we assume these glitches to fall at or near segment boundaries, i.e., t ðaÞ ≈ l ðaÞ T seg . Note that if a glitch happens inside a segment l, the mismatchμ ð0Þ l is bounded within [0, 1] and at most contributes 1=N seg to the total mismatch. One can always derive an upper bound on the mismatch by effectively removing the affected segment, assuming the glitch to happen in between segments, and adding 1=N seg to