A semicoherent glitch-robust continuous gravitational wave search

Continuous gravitational-wave signals from isolated non-axisymmetric rotating neutron stars may undergo episodic spin-up events known as glitches. If unmodelled by a search, these can result in missed or misidentified detections. We outline a semicoherent glitch-robust search method that allows identification of glitching signal candidates and inference about the model parameters.


I. INTRODUCTION
Continuous gravitational wave (CW) searches for rotating neutron stars typically assume an underlying signal model (a template) for the signal observed in the detector and then perform a matched-filter analysis (see, e.g., Abbott et al. [1,2]). These templates assume the phase evolution of the source is well modelled by a spin frequency, and several frequency derivatives. On the contrary, observations of pulsars demonstrate that neutron stars are subject to low frequency timing noise [3] and can also undergo sudden spontaneous increases in their rotation frequency and frequency derivatives known as glitches [4,5]. While the former effect is unlikely to have a substantial negative impact for searches of data lasting less than a year [6,7], typical glitches seen in the pulsar population may adversely affect current and ongoing CW searches. In Ashton et al. [8], we provide a statistical analysis of pulsar glitches and demonstrated that for a fully coherent matched-filter analysis, a glitch can cause a substantial relative loss of signal-to-noise ratio (SNR); semicoherent searches (in which the data is segmented, searched coherently and then recombined; for a review see Prix [9]) will suffer smaller relative losses of SNR by comparison, but, during the follow-up process 1 , a glitching candidate's SNR will not increase as expected potentially resulting in dismissal of the candidate.
In this work, we introduce a glitch-robust detection statistic in which the template also models the size and epoch of one or more glitches. Standard CW searches (by which we mean those using a non-glitch-robust detection statistic) are already computationally constrained, so adding additional parameters will increase the computational load and also reduce the significance of results due to the increased number of trials. Moreover, wideparameter space searches typically begin with a semicoherent stage, which, as previously mentioned, is more * E-mail: gregory.ashton@ligo.org 1 A follow-up refers to the process whereby a candidate from a semicoherent search is subjected to a series of searches, each of which increases the coherence time until the candidate is detected with a fully coherent search robust to glitches. Therefore, we do not advise modifying the initial semicoherent blind search strategies, but that during the follow-up and vetting of candidates, a glitch-robust statistic be used to guard against dismissal of glitching signals. We begin in Sec. II by defining a glitch-robust detection statistic. Then in Sec. III we give a discussion and comparison of how the statistic could be applied in a gridor Markov chain Monte Carlo (MCMC)-based search for CW candidates. In Sec. IV we discuss how to perform a model selection between glitching and non-glitching signals. Since most glitching candidates will initially be identified by a standard-CW search, in Sec. V we discuss how glitching signals might manifest in such searches and what simple steps can be taken to identify them. We conclude with overall discussion in Sec. VI.

II. SEMI-COHERENT GLITCH-ROBUST DETECTION
In this section, we introduce the glitch-robust detection statistic, an adaptation of the semicoherent F-statistic for glitching signals. We begin by defining the standard-CW F-statistic and then describe the glitch-robust modification.
For an isolated CW signal, the gravitational wave signal template, h(t), has two sets of parameters: the amplitude parameters A = {h 0 , cos ι, ψ, φ 0 }, consisting of the CW amplitude h 0 , inclination angle ι, polarization angle ψ and initial phase φ 0 , and the phase-evolution parameters λ = {Ω, f,ḟ , . . .}, consisting of the sky location Ω, frequency f and higher-order frequency derivatives f (k) (cf. Prix [9] for a general review). One key component of defining h(t) is the source-frame phase evolution, which for a standard-CW signal can be written as (e.g., see Eq. (18) of Jaranowski et al. [10]) In this work, we model the th glitch by δf (k) , the permanent increment in the kth frequency derivative at an epoch t g . For a source with N g glitches, the glitching source phase evolution is then where H(t) is the unit step function. This is analogous to the method used in pulsar timing [11], except that we do not model any exponentially decaying components. We refer to {f (k) } as the set of the frequency and its derivatives up to s max , { {δf (k) } } as the set of all glitch magnitudes for all glitches, and {t g } as the set of all glitch epochs.
The fully coherent F-statistic, used by many wideparameter space searches as a ranking statistic, is the loglikelihood ratio for signal vs. Gaussian noise, marginalized over the amplitude parameters [10,12,13]. Using only data spanning times [t s , t e ] from the full set of data x, we write the fully coherent statistic as 2F(x; λ, t s , t e ). Often, wide-parameter space searches use a semicoherent approach in which the total data span T of x is divided into N contiguous segments. Defining {t } as the set of start times for each segment, the semicoherent F-statistic is Ideally, a glitch-robust statistic would modify the standard-CW fully coherent F-statistic with the glitching source phase evolution, Eq. (2), resulting in a fully coherent glitch-robust detection statistic.
However, we propose instead the following pragmatic approach: let us use a semicoherent detection statistic with the glitch-epochs t g partitioning segments. Then defining 2F x; λ, {δf (k) }, t g , t g +1 as the fully coherent detection statistic calculated between glitches and assuming the source phase model of Eq. (2), we can define a glitch-robust semicoherent F-statistic: For convenience, we also define t g 0 and t g Ng+1 to coincide with the start and end time of the data used.
In this semicoherent detection statistic there are N g +1 contiguous segments which is implied by the size of {t g }, with the first glitch occurring at t g =1 . This simplistic method leverages readily available and tested code. However, this approach is potentially suboptimal compared to a fully coherent glitch-robust detection statistic. By using the semicoherent statistic over glitches, we allow for independent amplitude parameters A in each inter-glitch segment. This allows not only for a jump in phase, which would be physically plausible for a glitch, but also in amplitude h 0 and polarization angles cos ι, ψ, which is more likely to be unphysical. In principle one could build a coherent glitch-robust statistic by allowing only a phase jump in each glitch as an additional search parameter to the f (k) -jumps, but it is unclear if this would gain much extra sensitivity, and we postpone this to a future study.

III. GLITCH-ROBUST SEARCHES AND PARAMETER ESTIMATION
A search using the glitch-robust statistic (i.e., Eq.(4)) could be implemented in any number of ways. Indeed, it could be added to any standard-CW wide-parameter space search. However, these searches (see, e.g., the recent all-sky searches in LIGO O1 data Abbott et al. [1,2]) already demand massive computing efforts and adding (at least) two additional search parameters t g , δf would decrease the sensitivity to standard signals. In this section we investigate how a glitch-robust search can instead be applied in the follow-up of potential candidates identified in such searches. Sec. V provides some examples of how a glitching-CW signal may appear in a standard-CW search.
We assume a candidate has been identified by a standard-CW search with some uncertainty on its phaseevolution parameters λ. We first discuss the prior ranges for the phase-evolution and glitch parameters. We then introduce the necessary tools to quantify the size of a given prior parameter space before comparing grid-and MCMC-based glitch-robust search methods.

A. Glitch-parameter priors
For the standard-CW phase-evolution parameters λ, the prior (in the absence of other information) is chosen as uniform over the parameter space of interest; for a glitch-robust follow-up these will primarily be determined by the candidate uncertainty. In addition, the glitch-robust search requires priors on the glitch epochs t g , the magnitude of the frequency jumps δf and spindown jumps δf (k) , and the number of glitches N g .
For the number of glitches, a prior could be formed using the glitch rate observed in the pulsar population. However, dynamically searching over the number of glitches, which determines the total number of parameters, can be difficult. For MCMC-based searches, this would require a reversible-jump MCMC algorithm [14]. Instead, we suggest searching over the number of glitches by hand, namely, perform the search for different numbers of glitches and compare the results. We will discuss in Sec. IV how to quantify this comparison.
For the glitch epochs t g , a uniform prior over the data duration makes intuitive sense; we also assert that t g < t g +1 ∀ . In this work we pragmatically bound t g between 0.1 and 0.9 of the fractional data duration. This avoids boundary issues where there is insufficient data to calculate the F-statistic in the first or last segment and also reduces the parameter space to the region of primary interest, e.g., where a glitch will cause the maximum loss of detection statistic.
Choosing a prior for the jump sizes { {δf (k) } } is more difficult. Clearly it should be informed by the glitches seen in the pulsar population and one option is to use fits to the observed set of glitches in the pulsar population (e.g., see Fuentes et al. [5], Ashton et al. [8]). However, these may be affected by observational biases. A simple option is to use a uniform prior on {δf (k) } between a minimum and maximum value. For δf (0) , one approach is to set the minimum at zero (excluding anti-glitches where δf (0) < 0, cf. Archibald et al. [15]) and the maximum at twice the maximum observed glitches in the pulsar population (∼ 5 × 10 −5 Hz, see, e.g., Livingstone et al. [16]). Similar approaches can be devised for higher-order spin-down components.

B. The metric and the size of parameter space
In setting up any search, it is useful to have a metric to understand distances in the parameter space. Given a detection statistic d(θ) measured at some set of parameters θ, we first define a mismatch the fractional loss of detection statistic between the exact signal parameters θ s and some other point in the parameter space θ s + ∆θ.
For small mismatches, one may expand and approximate the full mismatch by the metric-mismatch where g ij is referred to as the metric and ∆θ i are the components of ∆θ. As discussed in the next section, the metric is useful in bounding the maximum loss of detection statistic when setting up grid-based searches. However, one should note that the metric mismatch is only a good approximation up to µ 0.3−0.5 [17][18][19]. Another useful application of the metric is in calculating N * , the approximate number of unit-mismatch templates covering the given parameter space [20], which can be understood as a proxy for the size of that parameter space.
Calculation of N * requires the ability to calculate the metric. We do not yet have the metric for the glitch-robust detection statistic defined in Eq. (4) (future searches may require this metric in, for example, a grid-based glitch-robust directed search). Nevertheless, it is still useful to calculate N * using the fully coherent standard detection statistic over the standard signal parameters, i.e., {f (k) } and Ω. This can be used a lower bound on the full N * for the full parameter space including the glitch parameters.

C. Grid-based glitch-robust search
Grid-based (or template-bank) searches compute the detection statistic over a number of prespecified points in parameter space with the grid of points covering the prior range. The grid spacing is selected to minimize both the maximum loss of detection statistic, bounded at some level, and the computing cost (i.e., to avoid oversampling the space). This spacing is determined using the metric; for the fully coherent and semicoherent F-statistic see Wette and Prix [18] and Wette [19] respectively. However, as previously discussed, we do not have the metric for the glitch-robust detection statistic. So, while we can apply the usual relations to any standard phase-evolution parameters used in the search and they should approximately hold, there is no simple way to determine the spacings in {t g } and { {δf (k) } } that guarantee a bound on the maximum mismatch.
In the absence of the relevant parameter space metric, we will employ a naive method here, simply dividing the full range of each search parameter into M steps. As such, the total number of grid points is M to the power of the number of search dimensions. This choice is not optimal (as would be the case if one were to derive and use the relevant metric), but captures many of the salient features of a grid-based search.
As an example of the grid-based method, we simulate a glitching signal in Gaussian noise with the properties given in Table I. Note that the glitch occurs in frequency alone, i.e., δf (k) = 0 for k > 0. We then perform a grid-based search over {f,ḟ , t g , δf } with M = 20; in this search the sky-location, Ω, is fixed to that of the simulated value. The prior ranges are given in Table II. In Fig. 1, we plot the semicoherent glitch-robust F-statistic in a grid-corner plot. This, as with the corner plots used in MCMC parameter estimation, displays the marginalized detection statistic for all one-and two-dimensional combinations.  Table II. Solid lines indicate the simulated signal parameters.  Table I and t g is defined from the start of the observation span. For the uncertainty in f andḟ , the number of fully coherent unit-mismatch templates is N * = 1000.

days
The grid spacing in this instance is sufficiently fine to provide reasonably good parameter estimation. For detection purposes it may even suffice to have sparser template coverage in the glitch time (where the signal appears quite wide compared to the prior range). At a fixed computing cost, this would allow for denser coverage in other parameters where the signal is narrower compared to the prior range.

D. MCMC-based glitch-robust search
MCMC-based standard-CW searches have already been used with success [22][23][24][25]. Recently we demonstrated [20] that this success relies on the size of the pa- rameter space being sufficiently small, as quantified by N * . Namely it was found that typically N * 1000 is a good guideline, but this can depend on the exact MCMC setup. For too-large parameter spaces, the MCMC algorithm tends to fail to converge to the signal peak in a reasonable amount of time.
For the follow-up of candidates from wide-parameter space searches, the size of the phase-evolution parameter space (i.e., the candidate uncertainty) is well constrained (or, this can be ensured by performing a refinement step). It is not possible to calculate N * for a glitch-robust detection statistic without the metric. However, in practice, we find that for typical glitch size and rates seen in the pulsar population [5,20] and typical observing spans, an MCMC-based glitch-robust search is effective at converging on simulated signals. For longer observing spans (or if allowing for larger glitches than those observed in the pulsar population) further work will need to be carried out to ensure the method is robust.
The advantage of an MCMC-based, instead of a gridbased, approach is that there is no requirement to predetermine the grid points. In effect, the ensemble MCMC sampler adapts to the topology of the maxima during the burn-in phase (for a more detailed overview of MCMCbased CW search methods see Ashton and Prix [20]).
To illustrate the results of an MCMC search, we run it on the same data set used to produce Fig. 1 (simulation properties are given in Table I) with the same uniform priors, as given in Table II. In Fig. 2 we plot the resulting corner plot. MCMC searches produce samples from the posterior, which, if the signal is successfully identified, usually occupies only a small fraction of the prior range. As a result, an MCMC search does not produce a posterior over the whole prior range, but only over the region of interest. As a consequence of this the range shown in Fig. 1 is much larger than that of Fig. 2: the latter shows the range of the posterior peak only while the former shows the entire prior range. Moreover, we note that in Fig. 1 for the grid-based search we plot the F-statistic, corresponding to the (marginalized) log-likelihood ratio. On the other hand, in Fig. 2 for the MCMC-based search we plot the estimated posterior which in this instance, where we use uniform priors, is proportional to the likelihood, corresponds to the exponential of the F-statistic. This is why the peak looks much narrower compared to Fig. 1 while showing in principle the same likelihood-function.

E. Comparing grid-and MCMC-based searches
In order to provide a simple comparison between gridand MCMC-based searches, we run a Monte-Carlo study. We produce 500 data sets containing a simulated signal with a single glitch in Gaussian noise. Such a signal, perfectly matched, has a predicted 2F of approximately 330. The noise, amplitude and standard phase evolution parameters are given in Table I, except that we jitter the frequency and spin-down, picking their value uniformly from within the inner half of the prior region given in Table II. We also select the glitch epoch from the distribution given in this table. Meanwhile, for the glitch magnitude, we sample from the observed pulsar population distribution [8]; while the aim of the section is to compare search methods, this choice of simulation distribution allows us to also verify that the naive priors are robust to a more astrophysically motivated population distribution.
Varying the required computation time (for the gridbased search, by varying the number of grid points; for the MCMC-based search, by varying the number of steps taken), in Fig. 3 we plot the relative difference between the recovered maximum 2F (for each method) and 2F inject = 2F(λ s ), the statistic measured at the simulated signal parameters λ s . Due to the presence of noise, the actual maximum 2F will typically not occur at the signal parameters λ s but be slightly offset and we therefore generally expect the maximum recovered 2F > 2F inject , provided the search method manages to localize the maximum 2F well enough.
From this figure, it is evident that at the same runtime, the MCMC-based search outperforms the gridbased search, with the majority of points finding a larger FIG. 3: Comparison of the relative maximum 2F found by each method compared to 2F inject , the value calculated at the simulated signal parameters λ s . All timings performed on an Intel Core i7-7820HQ CPU @ 2.90GHz processor.
detection statistic than 2F inject . This is to be expected since the MCMC search is operating optimally (i.e., the size of parameter space is sufficiently small). As such, the MCMC quickly converges to the maximum while a grid-based search spends most of the computing time calculating the detection statistic for points not near to the signal peak. For added context, Figs 1 and 2 both have an approximate run-time of 90s; for the grid-based search, the peak is only sampled a handful of times yet almost all of the MCMC samples (by design) come from the peak.

IV. GLITCHING VS. STANDARD-CW BAYES FACTOR
We now discuss how to quantify whether a signal is glitching and how many glitches best explain the data. We do this using a Bayes factor, the ratio of likelihoods for data x under two hypotheses. If H N implies that the data contains only Gaussian noise while H S implies that it contains an CW signal in addition to noise, then It can be shown [12,13,20] that the signal vs. noise Bayes factor at fixed phase-evolution parameters λ is whereF(x; λ, N ) is the N -segment semicoherent Fstatistic defined in Eq. (3) andρ max is an arbitrary upper cutoff on the prior range in signal strength [13].
Similarly, defining H gS as the glitching-signal hypothesis, we see that the targeted (in the sense that it depends on the model parameter) glitching-signal vs. noise Bayes factor is where the exponent is the glitch-robust semicoherent Fstatistic, defined in Eq. (4).
After marginalizing the targeted Bayes factor we get the signal vs. noise Bayes factor; i.e., for the standard search, while for the semicoherent glitch-robust The arbitrary prior cutoffρ max makes it difficult to interpret either of these Bayes factors by themselves: one can tune the Bayes factor by arbitrary changes in the prior. However, if we define the glitching-CW vs. standard-CW Bayes factor, then the arbitrary prior cutoff cancels and we are left with an interpretable Bayes factor for whether the signal is glitching or not. Calculation of the Bayes factor can be done by either a grid-based (using numerical integration of a dense sampling of the posterior) or MCMC-based method (using thermodynamic integration [26]). In future, we intend to extend the functionality to include nested sampling [27] which will improve the robustness of the evidence calculation (see, e.g., Ref. [28] for a comparison).
To understand the behaviour of B gS/S (x, N g ) as a function of the glitch magnitude, we run a Monte Carlo study, simulating 100 data sets (for each δf ) with a glitching signal in Gaussian noise. We use the parameters given in Table I, except δf which we vary systematically over a relevant domain. For each data set, we run a glitch-robust semicoherent MCMC search with N g = 1 along with a semicoherent MCMC search with N = 2 and calculate the resulting Bayes factor. The MCMC parameters are chosen such that the log Bayes factors are estimated to within a few percent. In Fig. 4 we plot the mean and standard deviation calculated over all data sets. We see that for small glitches, the Bayes factor prefers the standard signal hypothesis. But, once glitches are sufficiently large, the glitching-signal hypothesis is preferred.
To determine the preferred number of glitches, B gS/S (x, N g ) can be calculated for different N g and interpreted as a posterior over N g . Large numbers of glitches, 10 say, may be difficult to handle and require some tuning of the MCMC sampler. Fig. 4 illustrates that the glitch-robust search (as a function of the glitch size) plateaus above a certain minimum glitch size. An approximate way to characterize this size is to use the averaged (over glitch-time) single-glitch metric mismatch expressions derived in Ashton et al. [8]. Note that this is the metric for a standard CW search of a glitching signal and not the metric for the glitch-robust statistic introduced in Sec II. For example, for a fully coherent search at a fixed sky-location over frequency and spin-down, the minimum average metric mismatch is given byμ = (πT δf ) 2 /630. Setting the metric-mismatch to unity and inverting gives a rough order-of-magnitude estimate of the glitch size for which a standard-CW search would be sufficiently affected by the glitch that the glitch hypothesis will be preferred. Similar results can be derived for jumps in higher-order frequency spin-downs using the corresponding components of the glitch metric. In Fig. 4, we plot the value of Eq. (13), given the 50day duration of data used. Notably, this agrees with the point at which the Bayes factor begins to plateau.

V. IDENTIFYING GLITCHING SIGNAL IN STANDARD SEARCHES
In order to identify when a signal candidate from a standard-CW semicoherent search might best be followed up using a glitch-robust method, we now discuss the behaviour of glitching signals in a standard-CW search.

A. Multiple modes
One indicator of a glitching signal in a standard-CW search is the existence of multiple peaks in the detection statistic resulting from the template matching different parts of the signals. How exactly this behaviour manifests depends on the magnitude and size of the glitches, the data span, and the search setup.
Considering a signal which undergoes a single glitch with a jump {δf (k) }, we can identify two limiting cases depending on whether the glitch size is smaller or larger than a critical glitch size δf c , the effect of the glitch is negligible within the search setup; c , the signal can be thought of as two transient CWs, and we will find two distinct peaks in the detection statistic corresponding to the pre-and postglitch signal parameters. Between these two extremes, when δf (k) ∼ δf (k) c , the resulting structure in the detection statistic can be quite complicated. If required, the critical glitch size can be estimated from the single-glitch metric mismatches derived in Ashton et al. [8].
In order to illustrate this intermediate case, in Fig. 5 we show a standard-CW fully coherent grid search over frequency and spin-down for a simulated glitching signal. The simulation properties are given in Table I, except that we set δf = 2.4 × 10 −6 Hz and δḟ = −6.9 × 10 −13 Hz/s. As the reference time and glitch time FIG. 6: Illustration of the frequency sliding-window, a useful diagnostic tool for identifying glitching signals. In this example, a 10-day window is slid in increments of 1 day.
coincide, for the fully coherent search, the pre-glitch frequency and derivative is f s ,ḟ s while the post-glitch frequency and derivative are f s + δf ,ḟ s + δḟ . Two distinct ellipsoid patterns can be observed centered on the locations of the pre and post-glitch signals, but the maximum does not coincide with either.
Having multiple peaks in the frequency and its derivatives might be expected, but we typically also find multiple peaks in the sky position, even though the sky position of the source does not vary over a glitch. This is because by allowing the sky position to vary, the standard template fit to the glitching signal can be improved; this can happen in multiple ways, resulting in multiple peaks and will in general result in biases in the recovered sky position.

B. Sliding windows
A sliding window can be another simple, but powerful diagnostic test for a glitching signal. Fixing all other values to those of the maximum posterior estimate (or a set of parameters sufficiently close to the peak), the detection statistic is computed for a range of frequencies in an overlapping sliding time-window over the total data span. One could also do this for the frequency derivative (or any other parameter). Stacking the results together into a colour plot, if the signal is sufficiently strong, the glitch can easily be discerned from the change in frequency. We provide an example in Fig. 6 using the same data set used to produce Fig. 5.

VI. DISCUSSION
We have described a semicoherent glitch-robust detection statistic for use in evaluating if candidates found in wide-parameter space searches are glitching signals. This simple method adapts standard search routines, using a signal model that includes glitches as a set of instantaneous changes in the frequency and higher-order spindowns at a set of glitch epochs.
Comparing grid-and MCMC-based search methods, we find that the MCMC-based search is a superior method for performing glitch-robust searches of candidates from wide-parameter space searches. For the same computing cost it is able to better identify the maximum and perform parameter estimation vital to interpretation. MCMC-based glitch-robust searches are, for a suitable candidate uncertainty level, computationally cheap to run and provide parameter estimation and evidence estimates. Moreover, an MCMC-based method does not require a pre-specified grid template. We therefore recommend that such glitch-robust MCMC-based methods be used in the follow-up of candidates identified in wideparameter searches.
The methods introduced in this paper have been implemented in the package pyfstat [29]. Source code along with all examples in this work can be found at https:// gitlab.aei.uni-hannover.de/GregAshton/PyFstat.