Rotational evolution of the Vela pulsar during the 2016 glitch

TL;DR

This paper analyzes the 2016 Vela pulsar glitch, providing new constraints on the glitch rise time, evidence for a rotational overshoot and pre-glitch slow-down, offering insights into neutron star interior physics.

Contribution

It presents the first pulse-to-pulse analysis of a pulsar glitch, constrains the glitch rise time, and provides evidence for a rotational overshoot and pre-glitch slow-down.

Findings

Glitch rise time constrained to less than 12.6 seconds.

Evidence for a rotational-frequency overshoot and rapid relaxation.

Detection of a slow-down of the pulsar before the glitch.

Abstract

The 2016 Vela glitch observed by the Mt Pleasant radio telescope provides the first opportunity to study pulse-to-pulse dynamics of a pulsar glitch, opening up new possibilities to study the neutron star's interior. We fit models of the star's rotation frequency to the pulsar data, and present three new results. First, we constrain the glitch rise time to less than 12.6s with 90% confidence, almost three times shorter than the previous best constraint. Second, we find definitive evidence for a rotational-frequency overshoot and fast relaxation following the glitch. Third, we find evidence for a slow-down of the star's rotation immediately prior to the glitch. The overshoot is predicted theoretically by some models; we discuss implications of the glitch rise and overshoot decay times on internal neutron-star physics. The slow down preceding the glitch is unexpected; we propose the slow-down may trigger the glitch by causing a critical lag between crustal superfluid and the crust.

Tables

Table 1. Table 1: Definitions and Bayes factors for the four primary models tested in this Article, see Eq. ( 1 ) and the consistency relation. The final column is the log-Bayes factor between each model and the step-glitch model H s subscript 𝐻 𝑠 H_{s} . Uncertainties on the log-Bayes factors are  ≲ 0.2 less-than-or-similar-to absent 0.2 \lesssim 0.2 .

model	$f\left(t\right)$ parameters	${\log }_{10}B$
$H_s$	$\Delta f_r=\Delta f_d=0$	—
$H_1$	$\Delta f_r\ne 0$, $\Delta f_d=0$	$0$
$H_2$	$\{\Delta f_r,\Delta f_d\}\ne 0$	$0$
$H_{2+p}$	see Eq. (3)	$0$

Table 2. Table 2: Table of priors used throughout this work. For the parameters not introduced in the text, ϕ 0 subscript italic-ϕ 0 \phi_{0} is the phase parameter (number of rotations), we provide a wider prior to allow the reference pulse to not be zero; and σ 𝜎 \sigma is the standard-deviation of the Gaussian likelihood.

	Prior distribution	Units	Models
$f_0$	$\text{Uniform}(11.1854,11.1874)$	Hz	all
${\varphi }_0$	$\text{Uniform}(-5,5)$	—	all
${\sigma }_0$	$\text{Uniform}(0,0.01)$	s	all
$\Delta f$	$\text{Log-Uniform}({10}^{-8},{10}^{-4})$	Hz	all
$t_g$	$\text{Uniform}(-100,100)$	s	all
${\tau }_r$	$\text{Uniform}(0,1000)$	s	$H_1,H_2,H_{2+p}$
$\Delta f_d$	$\text{Log-Uniform}({10}^{-8},{10}^{-4})$	Hz	$H_2,H_{2+p}$
${\tau }_d$	$\text{Uniform}(0,1000)$	s	$H_2,H_{2+p}$
$\Delta t$	$\text{Uniform}(0,500)$	s	$H_{2+p}$
$\Delta f_p$	$\text{Half-Normal}(0,{10}^{-5})$	Hz	$H_{2+p}$

Equations

$$f_2(t)=f_0+H(t-t_g)\left[\Delta f+\Delta f_r{e}^{\frac{-(t-t_g)}{{\tau }_r}}+\Delta f_d{e}^{\frac{-(t-t_g)}{{\tau }_d}}\right],$$

$$\mathcal{B}\gtrsim 5.7\times 10^{-6}{\left(\frac{{\tau }_r}{12.6\mathrm{s}}\right)}^{-1}{\left(\frac{f_{\mathrm{sf}}}{11\mathrm{Hz}}\right)}^{-1}\left(\frac{I_{\mathrm{sf}}/I_{\mathrm{tot}}}{0.01}\right),$$

$$f_p(t)=f_2(t)-\Delta f_p\Pi (t;t_g-\Delta t,t_g),$$

$$\mathcal{L}(t_i;\theta )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{(t_i-h(i;\theta ){)}^2}{2{\sigma }^2}\right],$$

Full text

{affiliations}

School of Physics and Astronomy, Monash University, VIC 3800, Australia

OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery, Clayton VIC 3800, Australia

Department of Physics and McGill Space Institute, McGill University, 3550 rue University, Montreal QC H3A 2T8, Canada

School of Natural Sciences, University of Tasmania, Australia

Pulsar glitches, rotational irregularities of otherwise stably rotating neutron stars, are believed to be caused by the complex interplay between micro- and macrophysical properties of the star’s internal components. One model posits that superfluid vortices in the inner crust suddenly unpin, transferring angular momentum to the star’s lattice crust[1]. This is seen as an increase in the frequency of pulsations. Such models are difficult to verify; the internal components of the star are shielded from view and their behavior has to be inferred indirectly. Until recently, radio observations of glitches were limited to observations before and after the glitch, but not during. The detailed morphology of glitch dynamics (e.g., the glitch rise time) is therefore not well constrained or understood. In 2016, the first pulse-to-pulse observations of a glitch were made using the University of Tasmania Mt Pleasant 26 m radio telescope[2]. Those observations showed variations in pulse shape of four pulses starting 20 rotations before the inferred time of the glitch, which are attributed to variations in the magnetospheric state. A preliminary estimate of $\sim 4.4\,\mathrm{s}$ for the glitch rise time was given.

In this Article, we provide a detailed pulse-to-pulse analysis of the glitch morphology. Our main results are three-fold. First, we constrain the glitch rise time to less than 12.6 s with 90% confidence. We connect this with internal neutron-star physics using a body-averaged description of the components participating in the glitch. Second, we show that a frequency overshoot —an increase in the rotation frequency above the post-glitch equilibrium value — and subsequent fast relaxation exist immediately following the glitch, in agreement with the only two previous high time-resolution observations of Vela glitches[3, 4], and as explained by several models of neutron-star glitches[5, 6, 7, 8]. Third, we show that the glitch may be preceded by an initial precursor slow-down, whereby the crust of the star slowed before rapidly speeding up; we speculate that this preceding slow-down of the pulsar’s rotation triggered the glitch.

In the superfluid vortex model of pulsar glitches, the star’s crust, and components tightly coupled to it, spin down due to external electromagnetic torques. The crustal superfluid is decoupled from the lattice as a result of vortex pinning, implying a lag develops between the superfluid’s angular velocity $\Omega_{\rm sf}$ and the crust’s angular velocity $\Omega_{\rm crust}$; $\Omega_{\rm lag}\equiv\Omega_{\rm sf}-\Omega_{\rm crust}>0$. When $\Omega_{\rm lag}$ reaches a critical value $\Omega_{\rm crit}$, some mechanism initiates the glitch, simultaneously releasing a large number of vortices. The superfluid’s excess angular momentum is transferred to the crust, producing the observed spin up. The specific trigger for the angular-momentum transfer is not well understood, although speculation abounds[9, 10, 11, 12, 13, 14, 15, 16].

We show there exist fluctuations of Vela’s spin period prior to the glitch, see Fig. 1, and speculate that $\Omega_{\rm lag}>\Omega_{\rm crit}$ is reached due to a single, relatively large stochastic fluctuation. This is independent of the glitch trigger mechanism, but provides a means for reaching the critical lag between superfluid and rigid crust; we discuss implications below.

1 Frequency evolution

For a model-agnostic view of the pulsar’s evolution, we fit[17, 18, 19] a constant-frequency model to 200 s-long data segments, sliding this window throughout the data (for details, see the Methods section). In Fig. 1 we show the median frequency and 90% credible interval for each window. The frequency evolution on the vertical axis is given as the difference between the inferred frequency and a nominal value of 11.186433 Hz. Times are relative to the solar-system barycentre time of the fitted glitch MJD 57734.4849906[2]. As we detail below, this analysis method allows for a powerful comparison with physical models of pulsars’ rotational evolution. The glitch can clearly be seen near time zero. We caution that one cannot use this plot to find the glitch rise time as the sliding window time averages the frequency evolution, thereby elongating such features. The time window also suppresses the amplitude of features whose timescale is shorter than the window length. Despite this, Fig. 1 shows stochastic fluctuations in the frequency evolution for the $\sim 2000\,\mathrm{s}$ before and after the glitch. The largest fluctuations are immediately before, and immediately after the glitch. This motivates the more thorough analysis detailed below.

To quantitatively analyse the frequency evolution, we define three models. The most general of these is $H_{2}$, where rotational changes are characterised by a constant term plus two exponentials,

[TABLE]

where $H(t-t_{g})$ is the Heaviside step function, $t_{g}$ is the glitch time, and $\Delta f,\Delta f_{r},\Delta f_{d}$ are the glitch magnitude and amplitudes of each exponential term. Equation (1) is supplemented by the relation $\Delta f+\Delta f_{r}+\Delta f_{d}=0$ ensuring the frequency evolution is continuous at the glitch. This relation, along with positive log-uniform priors on $\Delta f$ and $\Delta f_{d}$ imply $\Delta f_{r}\leq 0$; in Eq. (1), the first exponential therefore describes the “rise” in the frequency on time-scale $\tau_{r}$, while the second exponential describes a “decay” in frequency on time-scale $\tau_{d}$.

Models $H_{s}$ and $H_{1}$, specified in Tab. 1, are limiting cases of $H_{2}$. These phenomenological models are motivated by analytic solutions to coupled rigid-body problems, discussed later where we also introduce a final model $H_{2+p}$.

For each model, we integrate the frequency evolution to obtain the phase evolution, which we invert to obtain the model-predicted arrival time for each pulse. A likelihood of the model given the data is calculated by modelling the pulse arrival time as a sum of the deterministic arrival time predicted by the timing model, and a zero-mean Gaussian process with unknown variance. Using this likelihood and a suitable set of priors for the model parameters, we infer their posterior distributions and the evidence for the model using PyMultiNest[20, 21, 22]. Complete descriptions of the likelihood and prior are given in the Methods.

The simplest model in Tab. 1, $H_{s}$—step glitch, ignores the complex morphology of the glitch and models the frequency evolution as a simple step function of amplitude $\Delta f$ at time $t_{g}$. We use this as a base-model against which we compare all other models in Tab. 1. For the $H_{s}$ model, we infer a glitch magnitude $\Delta f={16.11}_{-0.04}^{+0.04}\,\mathrm{\mu Hz}$ and time, $t_{g}={-0.31}_{-2.78}^{+2.74}\,\mathrm{s}$ (given relative to the reported value) consistent with the initial obsevation[2].

2 Glitch rise time

We analyse the rise time using a simple, physically-motivated reference model; a body-averaged model with two uncoupled spinning components that suddenly couple, has equations of motion that can be integrated to give a model with a single exponential rise time $\tau_{r}$[23], corresponding to model $H_{1}$ in Tab. 1. In the limit where $\tau_{r}$ is much smaller than the pulse period, the $H_{1}$ model is equivalent to the step-glitch model $H_{s}$.

In Fig. 2, we show the $\tau_{r}$ posterior for the $H_{1}$ model. The posterior peaks at zero; we cannot resolve the rise time of the glitch. This is consistent with the Bayes factor between $H_{1}$ and $H_{s}$ being in favour of the simpler step-glitch model; Tab. 1. Nevertheless, the $\tau_{r}$ posterior gives a 90% upper limit of $\tau_{r}\leq 12.6\,\mathrm{s}$. This improves upon the previous best upper limit on the rise time of the 2004 Vela glitch of $\tau_{r}\lesssim 30\,\mathrm{s}$[4].

We also fit a model with a logistic function such that $f(t)\propto 1/\left(1+e^{-t/\tau_{r}}\right)$. While the functional form differs, it captures a similar idea of a rise in frequency. Results are similar to the $H_{1}$ model; the Bayes factor favours the step glitch, with an upper limit $\tau_{r}\leq$ 7.64 s. The evidence for the logistic and $H_{1}$ models are comparable, implying neither is substantially favoured by the data. The physically-motivated reference model $H_{1}$ and the logistic model are phenomenological; we expect the true evolution of the glitch rise to be more complicated, although we know of no robust predictions in the literature. Nevertheless, the use of these two models shows our derived rise time is relatively insensitive to the details of the mathematical model. Throughout this work, we quote the more conservative upper-limit rise time of the $H_{1}$ model.

Within body-averaged models, the glitch rise is described by a dimensionless mutual-friction coefficient $\mathcal{B}$, controlled by the underlying vortex dynamics[24, 25, 26, 8]. Invoking a simplified two-component model, where the superfluid in the inner crust provides the angular-momentum reservoir for the glitch, the rise time is an indirect measurement of the coupling strength between the superfluid and the crust, with moments of inertia $I_{\rm sf}$ and $I_{\rm crust}$, respectively. We derive a lower limit

[TABLE]

where $f_{\rm sf}$ is the rotation frequency of the superfluid and $I_{\rm tot}=I_{\rm sf}+I_{\rm crust}$.

3 Glitch overshoot and relaxation

To investigate the overshoot and subsequent relaxation, we include a second exponential in the frequency evolution. This $H_{2}$ model (Tab. 1) is a simplified version of the three-component neutron-star model from Ref.[8], where the star is separated into crustal superfluid, core superfluid and the non-superfluid crust component. The $H_{2}$ model assumes the three constituents are rigidly rotating and coupled via constant mutual-friction coefficients. While the specific equation chosen to model the overshoot and decay is motivated by this three-component model, we treat it as phenomenological for understanding the glitch dynamics. This phenomenological model could also be interpreted in terms of alternative physical models that also predict frequency overshoots[5, 6, 7].

We fit model $H_{2}$ to the data, and show the maximum likelihood alongside the time-windowed data in Fig. 3. We show the raw model (dashed blue curve) and the time-averaged frequency evolution (solid blue curve); the latter can be directly compared to the time-windowed data (black curve).

Comparing the overshoot-decay model $H_{2}$ and step-glitch model $H_{s}$ yields a Bayes factor $\log_{10}B=$0, providing marginal support for the overshoot model. However, the $H_{1}$ model (which compared unfavourably against $H_{s}$) is a special case of model $H_{2}$. The more relevant Bayes factor to understand the importance of the overshoot and relaxation is between $H_{2}$ and $H_{1}$, for which $\log_{10}B=$0 showing substantial evidence in favour of the overshoot and relaxation. Alternatively, we can compare $H_{s}$ with a modified step-function evolution including a single decaying exponential. This model, often used in glitch-timing to model long-term, $\mathcal{O}(\gtrsim 1\,\mathrm{day})$, relaxation was also fit to the data: the Bayes factor, $\log_{10}B=$0$$, demonstrates strong support for a relaxation component with $\tau_{d}\sim 1\,\mathrm{min}$, further confirming the existence of the overshoot. We remind the reader that a Bayes factor of $\log_{10}B>2$ is considered “decisive” support for a model, while $1\leq\log_{10}B\leq 2$ is considered “strong” support[27].

The existence of the overshoot is clear both visually and through our quantitative analysis. This is not the first identification of an overshoot: Refs.[3, 4] found a similar feature in the 2000 and 2004 Vela glitches, the only other pulsar glitches with high-time resolution data. However, these were not as well resolved as the telescope was less sensitive, requiring 10 s-folding of the pulses to achieve sufficient signal-to-noise ratio to calculate times of arrival. Comparing with the work herein, the pulse folding also likely explains the less constrained $\lesssim$30 s glitch rise time.

The maximum-likelihood overshoot-decay model $H_{2}$ shown in Fig. 3 has a decay timescale of $\tau_{d}{=}{65.97}_{-24.44}^{+59.38}\,\mathrm{s}$, magnitude $\Delta f_{d}{=}{9.36}_{-6.38}^{+12.20}\,\mathrm{\mu Hz}$, and the rise time is similarly constrained as in the $H_{1}$ model. The large uncertainty on the decay time is due to a strong correlation with the size of the overshoot: larger overshoots with shorter decay times are just as probable as smaller overshoots with longer decay times. We discuss this in more detail below, including Fig. 4, which shows the covariance between the size of the overshoot and the decay timescale using model $H_{2+p}$.

4 Slow-down preceding the glitch

We investigate the frequency slow-down preceding the glitch by extending the $H_{2}$ model to include a step-function in frequency before the glitch. This phenomenologically models a spin-down event sometime before the glitch. The frequency evolution for $H_{2+p}$, the precursor slow-down model, is given by

[TABLE]

where $f_{2}(t)$ is given by Eq. (1), $\Pi$ is a rectangle function such that the frequency decreases by $\Delta f_{p}$ for the period $\Delta t$ prior to the glitch. Constructing the model in this way, $\Delta f$ remains the long-term frequency change at the glitch.

Model $H_{2+p}$ is phenomenological and motivated by the data; Fig. 1. Furthermore, it is one of many simple phenomenological models that could be used; e.g., exponential or linear drift models. Rather than perform a systematic study, we focus solely on $H_{2+p}$ with the aim to motivate further research in this area. To this end, we speculate below about the causes of the slow-down. Until we have a more physically-grounded model, the significance of the slow-down is difficult to establish.

The Bayes factor shows the precursor slow-down model $H_{2+p}$ is the preferred of all models tested here: comparing with the overshoot-decay $H_{2}$ model (the next most preferred), $\log_{10}B=$0$$. This suggests the data supports a slow down of the rotation prior to the glitch, in addition to an overshoot and decay. In Fig. 3, we show the best-fit slow-down model in red; the dashed curve represents the raw frequency evolution, and the solid curve the time-averaged best-fit model.

In Fig. 4, we show the posterior for the size of the overshoot $\Delta f_{d}$ and the overshoot relaxation timescale $\tau_{d}$ using the $H_{2+p}$ model. As mentioned, these are inversely correlated, implying a wide range of equally-likely values for both parameters.

In Fig. 5, we show the posterior of the precursor slow-down $\Delta f_{p}$ and the time before the glitch at which this occurs $\Delta t$. Although the size of the slow-down is not well constrained with $\Delta f_{p}={5.40}_{-2.05}^{+3.39}\,\mathrm{\mu Hz}$, we note this is a significant fraction of the actual glitch size $\Delta f={16.01}_{-0.05}^{+0.05}\,\mathrm{\mu Hz}$ for the $H_{2+p}$ model.

We use a half-normal prior distribution on $\Delta f_{p}$ (see Table 2, mentary Material), which places the maximum prior probability at zero, gives reasonable support over values $\lesssim 10^{-5}\,\mathrm{Hz}$, and exponentially disfavours larger positive values. If, on the other hand, we use uniform priors, then another local maxima in the posterior distribution becomes present at significantly larger $\Delta f_{p}$ and at a shorter time preceding the glitch. We find it difficult to physically motivate a precursor slow-down many times larger than the glitch size; this explains our choice of the half-normal prior. However, as all our models are phenomenological, we leave open the possibility that this larger mode exists, and leave that for future exploration.

The Vela pulsar has been seen glitching three times with high-time resolution observations[3, 4, 2]. The first two observations were with a less-sensitive instrument, implying pulses must be folded to obtain sufficient signal-to-noise. The method detailed herein can be used on those data, and we encourage reanalysis of that data in search of precursor slowdowns and overshoot-decays; such features have been hinted at in Ref.[3].

We divide models of the physical nature of the precursor slowdown into two groups. In the first, the fluctuation seen prior to the glitch in Fig. 1 is a large incarnation of the stochastic fluctuations seen preceding and following the glitch. In the second group, the slow-down is unrelated to this stochastic noise.

In future work, we will probe the hypothesis that the stochastic fluctuations cause the glitch by developing a statistical analysis of the Vela noise using data away from the glitch. This requires development of a spin-down model in the absence of a glitch that takes into account the stochastic fluctuations, with subsequent analysis of a large number of off-glitch data segments. Comparing the amplitude of the frequency slow-down ($\sim 5\,\mathrm{\mu Hz}$) with the noise distribution, could indicate if the slow-down is a statistical outlier from the typical noise, potentially allowing falsification of the idea that the spin-down event caused the glitch.

The cause of the stochastic fluctuations in Fig. 1 are also of interest. While they could be a manifestation of jitter noise, they may be due to fluctuations in the rotation rate caused by instabilities, or to extrinsic effects such as fluctuations in dispersion measure or scattering[28]. These ideas can be explored by looking at other sets of high time-resolution data, using the same method used to produce Fig. 1.

We hypothesise that the slow-down may be due to some intrinsic mechanism in the star, although this is far from certain. Reference[2] reported short-timescale variations in pulse shape during the glitch, including a null pulse and unusual pulse shapes before and after the null. These could indicate the glitch or preceding slow-down are magnetospheric in origin, although it is difficult to disentangle cause and effect. Moreover, the large change in spin period on such short timescales is difficult to quantitatively explain without catastrophically changing the magnetic-field topology, which would likely be accompanied with a long-term change in pulse shape as observed following glitches in high magnetic-field pulsars and magnetars[29, 30, 31]. Such a long-term change in pulse shape has not been observed following the 2016 Vela glitch[2].

5 Did the spin-down event cause the glitch?

That the glitch is preceded by a spin-down event is intriguing. Many models[32, 33, 34, 35, 6, 8] posit that glitches are triggered when sufficient lag is built up between the superfluid component of the inner crust and the lattice crust. Some fraction of the fluctuations in Fig. 1 may be due to intrinsic, stochastic variations of the star’s rotational period. Those variations may take place on timescales faster than the coupling timescales of the crust and internal components. If the slow-down event is one such stochastic variation, albeit a large one, we hypothesise this may trigger the glitch by spinning down the crust and driving the lag above its critical value.

Our hypothesis has natural corollaries for glitch statistics of the pulsar population: if the stochastic variations are large or comparable to the change in spin period from dipole radiation, glitches will occur probabilistically when the combination of the spin down and variations takes the crust-core lag above the trigger threshold. In such cases, the time period between glitches would neither be regular nor Poisson-distributed, but would depend on the relative size of the variations with respect to the spin down. When the variations are large with respect to the spin down, the glitch recurrence time should be Poisson-distributed. Finally, if the variations are small compared to changes from magnetic spin down, the glitch recurrence time should only depend on the spin-down timescale of the system. Additionally, we expect distinct behavior for those stars where glitches are driving the system far from the critical point and spin down is required to return to the threshold before variations can initiate a subsequent glitch, versus the objects where this is not the case and variations can always trigger a glitch.

It is worth noting that, in reality, the critical lag will depend on the neutron-star density and not have a single value, but there is likely also a statistical distribution associated with the macroscopically-averaged value of the critical lag. This should not significantly effect the arguments presented above as they are mainly qualitative. However it should be taken into account when quantitatively defending this model.

Although there are two populations of pulsars according to their glitch recurrence statistics[36, 37, 38], more work is required to establish whether the two populations correlate with the relative magnitude of the pulsar’s stochastic variations and their spin-down timescales. Such a statistically rigorous study would potentially be difficult; it is not clear which fluctuations are related to intrinsic pulsar spin noise and, e.g., pulse-jitter noise. This would also be difficult to generalize to other pulsars.

6 Conclusion

During the 2016 glitch, the Vela pulsar first spun down. A few seconds later it rapidly spun up, before finally spinning down with an exponential relaxation time of $\sim 60\,\mathrm{s}$. This model is substantially favoured over a simple step glitch, or one with only a single spin-up event (see Tab. 1).

Testing the rise time alone, we constrain $\tau_{r}\leq 12.6\,\mathrm{s}$ (90% confidence), consistent with the estimated value of Ref.[2], and reducing the previous-best constraint of $\tau_{r}\lesssim$30 s for the 2004 Vela glitch[4]. Invoking a two-component neutron-star model, our new constraint translates into a lower limit for the mutual-friction coupling of $\mathcal{B}\gtrsim 5.7\times 10^{-6}$; Eq. (2).

We find a frequency overshoot and exponential relaxation with amplitude $\Delta f_{d}={17.77}_{-7.99}^{+13.68}\,\mathrm{\mu Hz}$ and decay time scale $\tau_{d}={53.96}_{-14.82}^{+24.02}\,\mathrm{s}$ for the $H_{2+p}$ model, a feature that can be theoretically explained[5, 6, 7, 8]. For example, within the three-component model of Ref.[8], the overshoot only exists if the crust mutual-friction coefficient (coupling the crustal superfluid and crust) exceeds the core friction coefficient (coupling the core superfluid and crust). Providing a more qualitative analysis of the internal physics, e.g., constraining moments of inertia and coupling coefficients, is difficult at this point as the phenomenological models studied herein do not produce sufficient information.

Finally, we find evidence for a slow-down, or possibly a precursor antiglitch, immediately before the glitch. To the best of our knowledge this has not been predicted. We hypothesize that it may be a statistical fluctuation consistent with the overall noise fluctuations and speculate such fluctuations drive the differential lag between the superfluid and the crust above its critical value, thus triggering the glitch. This suggests a large number of glitches could be preceded by a slow-down, providing testable predictions.

Analyses like that presented herein only assess the relative evidence of models. We focus on phenomenological, albeit physically-motivated models, in a bid to remain model agnostic. Even the best fitting models tested here do not explain all the features in the data, e.g., Fig. 3. Future explorations may uncover new descriptions that explain the data better than the models used herein. For example, further theoretical modelling may provide a more nuanced view of how the slow-down preceding the glitch should manifest; the method we developed is easily extendable to compare more complex models.

While direct modelling is one avenue of further investigation, model-agnostic approaches may also yield considerable insight. Figure 3 is a first step in this direction, although it has the subtlety that the time window distorts temporal and amplitude features. Another method could be e.g., shapelet-based models for the frequency evolution, providing a means to study the underlying frequency evolution without modelling constraints.