Measuring glitch recoveries and braking indices with Bayesian model selection

For a selection of 35 pulsars with large spin-up glitches ($\Delta{\nu}/\nu\geq10^{-6}$), which are monitored by the Jodrell Bank Observatory, we analyse 157 glitches and their recoveries. All parameters are measured consistently and we choose the best model to describe the post-glitch recovery based on Bayesian evidence. We present updated glitch epochs, sizes, changes of spin down rate, exponentially recovering components (amplitude and corresponding timescale) when present, as well as pulsars' second frequency derivatives and their glitch associated changes if detected. We discuss the different observed styles of post-glitch recovery as well as some particularly interesting sources. Several correlations are revealed between glitch parameters and pulsar spin parameters, including a very strong correlation between a pulsar's interglitch $|\ddot{\nu}|$ and $\dot{\nu}$, as well as between the glitch-induced spin-down rate change $\Delta\dot{\nu}_{\rm p}$ that does not relax exponentially and $\dot{\nu}$. We find that the ratio $\left|\Delta \dot{\nu}_{\mathrm{p}}/\ddot{\nu}\right|$ can be used as an estimate of glitch recurrence times, especially for those pulsars for which there are indications of a characteristic glitch size and interglitch waiting time. We calculate the interglitch braking index $n$ and find that pulsars with large glitches typically have $n$ greater than $3$, suggesting that internal torques dominate the rotational evolution between glitches. The external torque, e.g. from electromagnetic dipole radiation, could dominate the observed $\ddot{\nu}$ for the youngest pulsars ($\lesssim10^{4}\;\mathrm{yr}$), which may be expected to display $n\sim3$.

distinct components of the pulsar, i.e., the star's crust and the interior neutron superfluid (Haskell & Melatos 2015;Antonopoulou et al. 2022).The rotation of this superfluid is supported by vortices of quantised circulation whose density defines the superfluid rotation rate.Whereas the neutron star's solid crust and charged components of the core are coupled and spin down together under the external electromagnetic torque, the superfluid can follow only if vortices migrate outwards and their density reduces.As first proposed by Anderson & Itoh (1975), vortices can become pinned in some stellar regions, i.e. immobilised by their interactions with other components such as the inner crust nuclei.Where this happens, the superfluid cannot lose angular momentum at the same rate as the rest of the star and hence acts as a hidden angular momentum reservoir.Glitches are the result of a rapid exchange of angular momentum between this reservoir and the observed crust, for example due to a catastrophic unpinning event which enables a fast decrease of vortex density.The glitch spin-up sizes and time intervals between consecutive glitches can be used to infer the fractional moment of inertia of the superfluid component that causes a glitch (Link et al. 1999;Andersson et al. 2012).
Following a glitch, the spin-down rate increases, driving the spin frequency towards the pre-glitch values.Many glitches show some transient recovery, often described by one or more exponential terms (Lyne 1987;Basu et al. 2020).The post-glitch relaxation to the pre-glitch rotational state is often not complete; most commonly, persisting changes, Δ  p , are observed (Antonopoulou et al. 2018).Additionally, interglitch values of the frequency second order derivative  often lead to much higher braking indices ( =  /  2 ) than those expected from known spin-down mechanisms.Nonetheless, the cumulative effect of the permanent glitch changes, Δ  p , counteracts the high interglitch  so that the long term evolution of glitching pulsars is typically closer to a canonical  b = 3 braking (Shapiro & Teukolsky 1983), often with  ≲ 2 (Lyne et al. 1996;Espinoza et al. 2017).
In order to undertake statistical studies of glitches, it is vital to have consistent glitch measurements.Large catalogues have published the parameters of hundreds of glitches, but often lack detailed measurements of the transient recoveries (Espinoza et al. 2011;Basu et al. 2022).Yu et al. (2013) presented measurements for the exponential transient recoveries, but their work does not account for the effect of timing noise.Recently Lower et al. (2021) applied Bayesian techniques to measure the properties of glitches, including exponential recoveries, whilst also modelling timing noise.The glitch parameters were determined in a systematic fashion, and Bayesian inference was used to determine if exponential terms should be included in the model of the post-glitch recovery.
The Jodrell Bank observatory (JBO) routinely monitors the rotation of more than 800 pulsars and records all detected glitches (Basu et al. 2022).In this paper we adopt a similar technique to Lower et al. (2021) and apply it to JBO-monitored pulsars.Since it is difficult to distinguish transient recoveries in small glitches, we implement a cut off on the glitch size when selecting our sample.We focus on 35 pulsars which displayed large glitches, Δ/ ≥ 10 −6 , during our observation span.These sources have a total of 157 detected glitches, which we analyse together with their respective recoveries.The paper is structured as follows: in section 2 we briefly introduce the observational program at the JBO.In section 3, we explain the method we adopted for measuring pulsar glitch parameters, whilst the fitting results are discussed in section 4. We finally present a further analysis for the statistics of glitch parameters and their correlations in section 5, followed by a summary of our findings in section 6.

OBSERVATIONS
The 76-m Lovell telescope at the Jodrell Bank Observatory (JBO) regularly observes around 800 pulsars.Typically, observations have a duration of ∼ 6 min depending on the pulsar's flux density, and are divided into a series of subintegrations lasting from 1 to 3 minutes (Hobbs et al. 2004).Data before year 2009 are recorded with an analogue filterbank (AFB), but in 2009 this was updated to a digital filterbank (DFB).The bandwidth for the AFB observations is 32 MHz, centered on 1400 MHz.A small number of observations are carried out at 500, 600 or 925 MHz.A bandwidth of 400 MHz, centred on 1520 MHz was used for DFB observations (Perera et al. 2019).
In each observation, psrchive (van Straten et al. 2012) is used to average the remaining sub-integrations, frequency channels, and polarisations to form a single integrated profile.The pulse timeof-arrival (ToA) can be determined by cross-correlating the integrated profile with a noise-free template profile representing an idealized realization of the expected pulse profile in the time domain.
The obtained ToAs are subsequently adjusted to the approximately inertial Solar system barycenter, employing the DE436 planetary ephemeris1 .This correction results in Barycentric arrival times, utilized in the subsequent timing analysis discussed below.Details on the observing setup and the JBO glitch database can be found in Basu et al. (2022).
The collective dataset of glitches detected by JBO was presented in Basu et al. (2022).It was shown therein that the distribution of all glitch sizes Δ can be modelled as two Gaussian components: a relatively narrow one centered at large glitch sizes, and a wide component spanning from small to intermediate/large glitches.For this work we focus on pulsars that have at least one large glitch, as these tend to have clearly detectable transient components in their recoveries.We therefore set a threshold based on the Gaussian mixture model of Basu et al. (2022) and select a subset of 35 pulsars that have at least one glitch with Δ/ ≥ 10 −6 .Details for the observational period and the total number of glitches for these pulsars are given in Table 1.Of these glitches, a subset -consisting primarily of intermediate to large size glitches -was used to explore post-glitch recovery signatures.

METHODOLOGY
In pulsar timing, we use the truncated Taylor series expansion to model the rotational phase as: where  0 ,  0 ,  0 , and  0 are the phase, rotation frequency, and its time derivatives at some reference epoch  0 .This timing model approximately describes the rotational phase of a pulsar over the full data span.The presence of glitches and timing noise introduces deviations from Equation 1.For each pulsar in our sample, an initial timing solution and preliminary estimates of glitch parameters are taken from Basu et al. (2022).Phase residuals are calculated as the differences between measured ToAs and those predicted by that model.Residuals are minimised by fitting for the rotational parameters of pulsars as well as any astrometric, interstellar medium (ISM), and (where applicable) binary parameters.
In this section, we introduce the glitch model and noise models we used to describe pulsar rotation.To investigate the post-glitch recoveries, we used Bayesian model selection to compare evidence of various models, such as models with a varying number of exponential relaxation terms.

Bayesian model selection
Bayesian inference offers a consistent approach to the estimation of a set of parameters Θ describing a model or hypothesis , given the data .The core of Bayesian analysis is Bayes' theorem, which states that  Our aim is to work out the most plausible glitch and postglitch recovery model (see equation 4).We tested for post-glitch exponential relaxation and glitch-induced permanent step changes in frequency second order derivative, and measure the posteriors of parameters.The 7 different post-glitch recovery models explored are summarised in Table 2. To decide on the appropriate model we adopt a slightly modified version of Jeffreys' scale (Jeffreys 1939), and use the logarithm of Bayes factor (ratio of the Bayesian evidence between model 2 and model 1) to access the strength of evidence.A value of ln  21 < 2.5 is interpreted as model 2 having insufficient evidence over model 1, 2.5 ≤ ln  21 < 5.0 is regarded as moderate evidence favouring model 2 over model 1, 5.0 ≤ ln  21 < 10.0 as strong evidence preferring model 2, and ln  21 ≥ 10.0 as decisive evidence for model 2. Since the prior odds ratio is unity among models of a specific glitch, the posterior odds ratio is equivalent to the Bayes factor.For example, a logarithmic Bayes factor ln  21 ≥ 2.5 means the odds ratio of model 2 over model 1 is larger than  2.5 ≃ 12.Note that only the ratios of the evidence between models matter, so one model might have much larger evidence than the other (for example a double recovery model compared to a model without any recovery terms), but both models may be incorrect.
In order to generate a single result for each glitch/pulsar, and to make plots, we select a single 'best' model for each glitch.If there is a single model that has at least moderate evidence over all other models then we select this as the best model.However, in the case where several models have similar evidence, we inspect all models with logarithmic Bayes factor ln  21 > −2.5 and choose from amongst these.In principle the statistical modelling suggests all these models are plausible, and we usually select the simplest model (the model labelled with smallest number), but in two cases (PSR J0631+1036 and the 1st glitch of PSR J1907+0631) we selected a different model based on visual inspection of the residuals.These special cases will be discussed in Section 4.1.

Glitch Model
We can describe the additional pulse phase change induced by a glitch (only applied to post-glitch data) by: where Δ =  −  g ≥ 0, and  g is the glitch epoch.Besides permanent changes in rotational parameters (denoted with index p), we assume that any transient post-glitch components can be described as a sum of  e exponentials, with amplitudes Δ d  and timescales   .In principle, one can include the quantity Δ ≡  0 −  ′ 0 to account for any discontinuity in phase due to uncertainties in the glitch epoch.However, in this work, we stipulate Δ = 0 to avoid the strong covariance between glitch epoch   and Δ.From the glitch model in Equation 4, the change in rotation frequency at the time of a glitch is Δ = Δ p +  Δ d  and the change in its first derivative is Permanent step changes Δ p and Δ  p are fitted for in all our glitch models of Table 2, whilst the inclusion of other terms varies.No exponential relaxation term is considered in Model 1 and 2. We fit for one exponential relaxation term in Model 3 and 4, two exponential relaxation terms in Model 5 and 6, and three exponential relaxation terms in Model 7. A permanent change in the second order frequency derivative Δ  p is fitted for in Model 2, 4, and 6.We always consider models 1-4, for every glitch.When we find evidence for exponential relaxation, then we also test for models 5 and 6.In rare cases when there is evidence in favour of multiple exponential recoveries, Model 7 is tested for as well.

Noise Model
In addition to the deterministic timing model (equations 1 and 4) and the time-uncorrelated instrumental noise, stochastic processes also contribute additional time-correlated noise in the observed ToAs (Cordes & Shannon 2010).In this work, we adopt the noise model described in Lentati et al. (2013) and solve for both the pulsar timing model, white noise parameters, and a stochastic Gaussian-process noise model simultaneously.
For the time-correlated timing noise, we use a Fourier basis Gaussian process, where the power spectral density (PSD) at a frequency  is given by where A is a dimensionless amplitude (Lentati et al. 2013).White noise is modelled such that a ToA with error   is weighted as if it had an error, where EFAC and EQUAD represent a scaling and quadrature addition term (Lentati et al. 2014).EFAC and EQUAD values are estimated for the AFB and DFB data separately as the systematic noise processes may be different for the higher resolution DFB data.

Implementation
In order to fit for the glitch parameters we have added them to run_enterprise (Keith et al. 2022), a single pulsar Bayesian toolkit based on the enterprise framework (Ellis et al. 2019), which can also simultaneously fit the noise model parameters.We adopt multinest (Feroz et al. 2009) as the primary solver for sampling and calculating Bayesian evidence.
To do the fitting, prior ranges for the fitting parameters need to be defined.We find that a uniform prior in Δ p and Δ  p results in high covariance among Δ p , Δ  p , and Δ d .Further, the majority of the parameter space yields phase predictions that are extremely divergent from the observed residuals, leading to numerical instabilities in the fitting, particularly as the time-correlated noise component tries to model the resulting residuals.Therefore, in order that the priors include only phase models sufficiently close to the observed residuals, we adopt an alternative parameterization for the glitches, replacing Δ p and Δ  p with Δ, the instantaneous change in the spin frequency, and Δ Δ , the additional glitch-induced change in spin frequency between the moment after the glitch and a time Δ later.The choice of Δ is arbitrary, but should be chosen such that the spin frequency at Δ can be well measured from the residuals.The relationship between the new and old parameterisation is The posterior distribution for the traditional glitch parameterisation, Δ p and Δ  p , can be reconstructed after the sampling is complete.We choose a Δ in the range 50 days ≲ Δ ≲ 700 days, where we find that the density of ToAs will allow for a good estimate of Δ Δ .Hence, we can estimate Δ and Δ Δ directly from the gradient of the residuals before the glitch, immediately after the glitch, and at Δ, in order to define a reasonable prior range on Δ Δ and Δ.Based on the estimate for the mean value of Δ and Δ Δ , and the standard deviation  for the gradient of the residuals, we set a uniform prior range around their mean value with a size of ±100.This range is usually narrow enough to exclude any implausible models that cause numerical issues but also wide enough as to not overly constrain the fitting.
We perform a preliminary parameter fit using tempo2 (Hobbs et al. 2006), and set uniform prior ranges for other parameters (glitch epoch, Δ d  , and Δ  p if applicable) around their estimated values.We pick a default prior range of ±2 days for the glitch epoch, ±0.8Δ p for the amplitude of the exponential transients Δ d  , and ±10  for Δ  p .
For   , we take a uniform prior in log space between 10 and 1000 days.A lower limit of 10 days is chosen because we do not have enough ToAs to constrain an exponential relaxation shorter than 10 days.The upper limit of 1000 days is set as it is hard to constrain an extremely long exponential relaxation.We stipulate  1 <  2 <  3 , and split the prior ranges into two or three parts for different relaxation components tested in models 5-7, so that 10 <  1 <  1 <  2 <  2 <  3 < 1000, where  1 and  2 are constants used to split the priors, in order to avoid covariance among exponential relaxations.
In some cases we find that the default prior bounds can still be too large and again lead to numerical issues, especially when fitting for a large number of parameters.If this occurs we selectively reduce the prior bounds to ensure that the fit converges cleanly.It is important to ensure that the priors are the same for all models of that glitch (for evidence computation) and for all cases we verify that the posterior distribution is small (≲ 10%) compared to the prior.Furthermore, the glitch parameters are not expected to be multi-modal, so further widening of the priors should not influence the posteriors.
In the model selection procedure, we normally include only ToAs around the target glitch.However, when the previous or next glitch is too close to the target glitch, or when nearby glitches have considerable effects on the timing evolution around the target glitch, we extend the ToA span to include these glitches as well.Small glitches without modelled recoveries are treated as step changes in  and  only, and all the remaining timing model parameters are extracted from the maximum likelihood solution using tempo2.
After the model selection procedure is completed for all target glitches in a pulsar, we produce a single timing solution that includes the maximum likelihood solution for the best model for each glitch.We then re-fit the red and white noise model parameters with multinest, and the basic timing model parameters,  0 ,  0 ,  0 , astrometric parameters RAJ, DECJ, parallax, proper motion, dis-persion measure with tempo2.We also allow the parameters of the smaller glitches (for which no model selection was performed and are simply modelled as step changes in  and ) to vary during this process.The final results consist of the time series of  and  and the best solution for all fitted glitch parameters (glitch epoch, Δ, Δ Δ , Δ  p , Δ d  ,   ), the Gaussian process noise model parameters (, ), and an EFAC and EQUAD for each observing system in the dataset.The permanent step changes Δ p , Δ  p , and the recovery ratio  = Δ d Δ are derived afterwards.

RESULTS
In this section we present the parameters of 157 glitches in the 35 pulsars of our sample.The basic pulsar parameters are given in Table 3: the second to fifth column are the  0 ,  0 ,  0 , and  0 in the phase model Equation 1for each pulsar.Table 3 also shows the spin-down luminosity calculated by the inferred surface magnetic field strength and the braking index  =  /  2 .We also include the glitch activity parameter,  g = 1  Δ, which is the cumulative spin-up due to multiple glitches over a total observation span  (see Basu et al. (2022) and Antonopoulou et al. (2022) for discussions of this parameter).

Overview
Based on the work of Basu et al. (2022), pulsars with 0.03s ≤  ≤ 1s and 10 −15 ≤  ≤ 10 −12 (characteristic age around 10 1 ∼ 10 3 Kyr) are more likely to undergo large glitches that satisfy our selection criterion of Δ/ ≥ 10 −6 .Hence, the pulsars that form the sample of our work are located in a relatively small region in  − .
Most pulsars in our sample (22 out of 35) have multiple glitches during our observations, whilst 13 pulsars have only one glitch.For instance, PSR B1737−30 has 36 reported glitches but, in a similar time interval, PSR B1930+22 has just one.
From the total of 157 glitches, 61 have Δ/ ≥ 10 −6 and hence meet our criteria for a 'large' glitch, and 19 have an intermediate size (10 −7 ≤ Δ/ < 10 −6 ).All glitches with Δ ≥ 10 −7 are included in the model selection process.Small glitches, on the other hand, are typically treated just a steps in  and  (calculated by a linear least square fit) because they do not commonly show detectable exponential recoveries.However, there are 5 small glitches for which preliminary analysis suggested the presence of additional recovery features, so these are included in the model selection process as well.These are the 1st glitch of PSR B0355+54, the 6th, 12th, 17th glitch of PSR J0631+1036, and the 6th glitch of PSR J2229+6114.We mark these glitches with asterisk in Table 4.
There are two pulsars (PSR J0631+1036 and PSR J1907+0631) where we did not strictly follow the model selection process in Section 3.1, specifically when comparing models with equivalent Bayes factors we did not choose the simplest model.In both these cases there is some ambiguity between glitch behaviour and timing noise, and sometimes equivalent models with additional terms improved the final timing solution in the full dataset.Specifically PSR J0631+1036 has a large number of glitches close together and the overall timing solution is improved by including a single exponential term in the 4th and 17th glitches, and Δ  in the 6th and 15th glitches.For PSR J1907+0631 excluding the exponential recovery in the 1st glitch results in the same exponential recovery behaviour being absorbed into the timing noise model, which seems less likely.The logarithm of Bayes factor for all the models we tested are given in Table B1.
Among the 85 glitches for which we performed model selection, 37 include at least one exponential recovery term in the best model (see Equation 4).The preferred model for most of these includes only a single exponential, with a two exponential model favoured only for 15 glitches.In our dataset, only the 2nd glitch (MJD 46469.9) of PSR B0355+54 and the glitch at MJD 53642.3 of PSR B2334+61 present a Bayes factor indicating evidence for a third exponential.The timescales of the exponential recoveries range from around 10 days (our lower prior bound) to 879 days.We note that the lower prior bound is needed to prevent the timescale converging to extremely small values (which are impossible to resolve with current ToA cadence), and hence timescales consistent with 10 days should be taken to mean that the timescale might be smaller than what we can measure.
We list our solutions for the glitch parameters, including any exponential recovery components, in Table 4.In some cases, two large glitches were fitted simultaneously -accordingly, the range of ToAs includes the previous or next large glitch, as shown in Table 4.This was done when the timing solution for one glitch could affect the derived parameters for the other due to their temporal proximity.For the small glitches, we adjust the glitch epoch until there is no perceivable phase discontinuity at the glitch.

Spin Evolution
It can be informative to visualise time-series of () and ().The glitch contribution can be computed by analytic derivatives of Equation 4, however, we need to also consider the spin evolution attributed to timing noise in our model.A smooth model of the timing noise can be derived directly from the Gaussian process noise model (e.g.Keith & Niţu 2023), however, we can also directly measure  and  over a given small time window from the residuals.Therefore, after obtaining the best model, we use the striding boxcar method described by Shaw et al. (2018Shaw et al. ( , 2021) ) in order to visualise the evolution of the spin frequency  and the spin frequency derivative .We fit for , , and  in small windows whose widths span 6 times the average cadence.We shift the fit window by half its width generating a series of spin parameter values over the entire data span.For any window near a glitch, we adjust the start or finish epoch, so that they are immediately after or before the glitch respectively, in order to avoid fitting the rotational parameters over the glitch.
Figures 1 to 6 present the rotational evolution of () and () from the stride fitting method and as derived from the maximum likelihood timing solution, including the best glitch models and the timing noise model.All figures are organised in four panels as described in the caption of Figure 1.
The  evolution for these pulsars is dominated by the large glitches that lead to their selection for the sample.Once the large permanent changes are subtracted, the evolution of the frequency is The scaling of vertical axes in the plots gives the impression that the amplitude of the exponential terms varies strongly over the sample.In fact, the decaying amplitude Δ d only varies by about an order of magnitude (10 −8 − 10 −7 Hz) over the sample, but this is seen against the background of the overall scale of Δ  p and , which both vary by seven orders of magnitude across our sample.The evolution of the glitch parameters over the  −  diagram will be discussed in section 5.
Among our sample, we highlight PSRs J1740+1000 and B1757−24 as two pulsars which exhibit examples of some common glitch features.PSR J1740+1000 (middle right panel of Figure 2) has a single large glitch over 13 years in our data span.Typically such pulsars have a relatively small | | and , and their post-glitch behaviour is often dominated by exponential relaxation.PSR B1757−24 (bottom right panel of Figure 2), on the other hand, has several large glitches of similar amplitude and interglitch waiting times.Such a regular pattern of glitches resembles what is seen in the Vela pulsar, and hence such pulsars are often termed 'Vela-like'.
After the step change in  and  at the glitch, the post-glitch recovery process resets most of the change in  before the next large glitch occurs (see Subsection 5.2 for details on this).Whilst exponentially decaying components are detected in some of its glitches, the linear post-glitch recovery has the dominant effect on the long-term  evolution.
The youngest pulsars in the sample, PSRs J0205+6449 and J2229+6114, show an interglitch  evolution that does not indicate a 'linear' recovery nor can be described well by exponential terms.The post-glitch () curve appears more like an 'S-shaped' curve with the gradient of  increasing smoothly a few hundred days after the glitch before flattening off again.As a consequence, our models describe glitches in PSR J0205+6449 (top left panel of Figure 1) with permanent changes in  and  only, while for PSR J2229+6114 (middle right panel of Figure 6) our models favour exponential recoveries with negative changes in  for the 2nd, 4th, and 9th glitch.The difference between the preferred model for these pulsars seems to come down largely to the timing noise model: the Bayesian code may attribute the post-glitch behaviour to a sum of exponentials or to the timing noise model.Such 'S-shaped' swings can also be seen in PSR J1709−4429 (Lower et al. 2021) and are likely to be glitch related.Future improvements to glitch modelling may reduce the amount of this glitch activity currently attributed to timing noise.Lyne (1990); 9. Zou et al. (2008); 10.Shemar & Lyne (1996); 11.Krawczyk et al. (2003); 12. Urama (2002); 13.Janssen & Stappers (2006);14. Jankowski et al. (2016);15. Basu et al. (2020);16. Wang et al. (2000); 17.Hobbs et al. (2002);18. Hessels et al. (2004);19. Gügercinoğlu et al. (2022);20. Yuan et al. (2010b).
PSR B0355+54 is a well studied (see e.g., Janssen & Stappers (2006)) pulsar which showed a very large glitch, preceded by a more moderate one by ∼387 days.Our model of the first glitch of PSR B0355+54 (top right panel of Figure 1) finds an unusual positive value of Δ  p , which is not reported previously.However, there is a sparsity of data between the two glitches and the exponential recovery of the first glitch does not fully recover before the second glitch occurs.Therefore we suspect that the positive value of Δ  p should be treated cautiously as it is not clear that there is sufficient data to unambiguously determine the 'permanent' effect of the first glitch.
Another pulsar of interest is PSR B1737−30 (top right panel of Figure 2), for which we have 36 glitches in our dataset.We do not find any exponential terms in the glitch recoveries of this pulsar, but we also note that the very high density of small glitches means that the fitting may struggle to distinguish between spin evolution or The epochs of ToAs by vertical blue lines; residual , which is  after subtracting  f and  f ( −  f );  after further subtracting the permanent changes due to glitches  gp = Δ p + Δ  p Δ + 1 2 Δ  p Δ 2 and 1 2  f ( −  f ) 2 , as well as some constant offset so that the point before the first glitch is zero;  as a function of time.The reference epoch  f , reference frequency  f , and its derivatives (  f ,  f ) are some constants given in Appendix A in Table A1.Glitches best-modelled as steps in  and  (without exponential recoveries) are indicated with dotted lines at their glitch epochs.Dash-dotted lines mark the epoch of glitches for which one exponential term was included in the best model, dashed lines for those with two exponentials, and solid lines for those with three exponentials.An analytical model derived from the best-fit glitch parameters is shown as a black-solid curve.We use green-solid lines to mark gaps longer than 10 times the pulsar's average cadence in ToAs.The stride fitting results for  and  are shown as the red and purple dots with errorbars respectively.MNRAS 000, 1-22 ( 2024)     MNRAS 000, 1-22 (2024) glitch recoveries and the parameters of the next glitch.Urama (2002) reported an exponential recovery after the fifteenth (MJD 50936.8)glitch, with a timescale of 9 days.We have no ToAs within 10 days after this glitch and only 2 ToAs within 10-20 days, so it not possible to verify this exponential with our data.Zou et al. (2008) found a big ( = 0.103) exponential recovery ( d ∼ 50 days) at the twentyfirst (MJD 52347.7)glitch, and a smaller ( = 0.0302) exponential recovery of 100 days at the twenty-sixth (MJD 53036) glitch.Our one-exponential model for the twenty-first glitch ( d ∼ 40 days,  = 0.04(4)) did agree with their results, but the Bayes factor for this model is not large enough (ln  21 = 2) and the exponential timescale is poorly constrained so we did not accept it.In addition, Liu et al. (2019) reported an exponential of Δ d = 9.5(6) × 10 −9 Hz with a timescale of 71( 6) days for the thirty-sixth (MJD 58232.4)glitch, but our Bayes factor does not favour models with exponential recovery for this glitch.
PSR J0631+1036 (middle right panel of Figure 1) also has very frequent glitches, and has a  = −1.265× 10 −12 Hz/s similar to PSR B1737−30 (  = −1.266× 10 −12 Hz/s).The 17 glitches of PSR J0631+1036 in our data span are described by a range of different glitch models: four of them have exponential recoveries, two of them additionally have a noticeable change in .The frequent small glitches and noisy () variations in PSR J0631+1036 and PSR B1737−30 make it hard to be certain that all the glitch recovery is fully captured or discernible from other noise-like processes.Other pulsars with similar frequency derivatives are PSR J1841−0524 (middle left panel of Figure 4), PSR J1850−0026 (bottom left panel of Figure 4) and PSR B1830−08 (bottom right panel of Figure 3).They do not glitch as frequently as PSR J0631+1036 and PSR B1737−30, but they also show some modulations in their  evolution.
PSR J1921+0812 (bottom right panel of Figure 5) has an increased  for ∼ 1000 days following the glitch around MJD 55349 but then  returns to its pre-glitch level.The particular Bayes factors do not support exponential recovery models nor a permanent change in .Instead, our algorithm prefers to describe the change in  as timing noise.
All pulsars with a measurable  show a positive value, except PSR B1830−08 (bottom right panel of Figure 3), which has a clear systematic linear increase in | | since its large 2nd glitch.Unfortunately there is insufficient data to unambiguously determine if  before this large glitch is also negative.Unlike other pulsars, this negative  and the effect of glitches work together to increase | |.Understanding the origins of a negative  will be important for any theory that hopes to explain the long-term evolution of pulsars.
Whilst large glitches are typically accompanied by a negative Δ , the glitch of PSR B0919+06 (bottom right panel of Figure 1) at MJD 55151.7 is an unusual one.Our best model suggests a negative decaying component Δ d = −2.2×10−8 s −1 and a positive Δ  = 2.4×10 −14 Hz/s.The negative exponential recovers on a very short timescale of about 10 days as shown in the bottom right panel of Figure 1.We note that this model is favoured mostly due to the two ToAs right after the glitch.So our model for this glitch is highly dependent on the accuracy of these ToAs.Apart from this negative exponential and the three in PSR J2229+6114 mentioned earlier in this section, all other exponential recoveries we measure are associated with positive changes in .
Many large glitches have significant exponential recoveries, however not every large glitch is necessarily associated with an exponential recovery.The exponential timescale is sometimes shorter than we can reliably measure, and in these cases we caution that the exponential model may not be correct as typically only one or two ToAs close to the glitch are giving rise to the requirement for a transient recovery term.In these cases it is worth noting that the uncertainty on the instantaneous change in  strongly depends on the exponential model.
Finally, we note that several pulsars (PSRs J0729−1448, B1727−33, B1821−11, B1830−08, and B1930+22) exhibit quasiperiodic oscillations in , which can be seen by eye in the  time series.Such behaviour is not uncommon across the pulsar population, though not fully understood (Hobbs et al. 2010;Lyne et al. 2010;Niţu et al. 2022).The quasi-periodic oscillations are thought to originate from changes in the magnetosphere, though there have been hints of a connection to glitch activity (Keith et al. 2013).Further investigation of quasi-periodic oscillations are left for future work.

Evolution of 𝜈 and braking indices across the pulsar population
The spin-down of a pulsar under a braking torque can be assumed to have the form  ∝ −  , with  called the braking index.The  of a pulsar can then be predicted to be For the case where the braking is purely due to electromagnetic dipole radiation, we expect  b = 3.However, as seen in the results of Table 3 we typically find that  for the linear interglitch phase is of order 10 to 100.Note that the fitted parameter  in these glitching pulsars does not reflect the long-term evolution of the pulsar spin down, which -as can be seen in Figures 1 to 6 -is in general slower (smaller effective ).It has been postulated that the high interglitch  observed in several glitching pulsars is the result of internal, superfluid, torques (Antonopoulou et al. 2022).
Our results indicate that the interglitch  does not substantially change between glitches of the same pulsar.Conversely,  can vary considerably between different glitching neutron stars.To explore this, we investigate the evolution of  across the  −  diagram.We do this by fitting  as a function of powers of  and  using a Bayesian solver emcee with MCMC (Monte Carlo Markov Chain) (Foreman-Mackey et al. 2013).Specifically we model log where ,  and  are fit parameters and  is a Gaussian noise term with variance 2 .The result is shown in Figure 7 where  = −0.2± 0.4,  = 1.4 ± 0.1,  = −6.1 ± 1.9, and  = 0.53 ± 0.07.Comparing the exponents of  and  with Equations 8-10, we see that  does not strictly track any of ,  c or  S .We find that although  generally tends to decrease with  c , with Spearman rank correlation coefficient   = −0.90(see e.g.Spearman 2010), it most strongly correlates with  directly,   = −0.93(here we only consider measurements with  > 0 and at least 3 confidence.).
The evolution of  across the  −  diagram is shown in Figure 8, which includes our sample of pulsar as well as all other glitching pulsars found in the pulsar catalogue 2 (Manchester et al. 2005).If we assume that the observed  (described by Equation 12) is the result of internal and external torques acting on the star, and that the external torque has an  b ≃ 3, then it is clear that for most

Correlations between waiting times and glitch parameters
Previous works such as Melatos et al. (2018) have focused on the relation between glitch size and glitch waiting time, and many of them suggested some correlation between them (Carlin & Melatos 2019;Fuentes et al. 2019).Such a relation has implications on the physical mechanism behind glitches.Moreover, establishing a welldefined statistical procedure for predicting a glitch in a given pulsar will be greatly beneficial for long-term planning of observations.
The most active glitching pulsar that we know of is PSR J0537−6910, which has over 3 glitches per year on average.This is the only source for which a clear, strong, correlation between glitch amplitude Δ and the time to the next glitch Δ + has been established (Middleditch et al. 2006;Antonopoulou et al. 2018).The same studies also showed a weaker correlation between the change Δ  and the time preceding a glitch Δ − .
Another reported correlation is between the instantaneous change of frequency derivative at a glitch |Δ | and the subsequent  + Δ + , where  + is measured over the post-glitch interval Δ + .Lower et al. (2021) modelled this relation as a power-law and used collectively the glitch parameters for 16 pulsars, finding  = −4.3+2.5 −2.6 and  = 0.80 ± 0.12 and a Spearman's rank correlation coefficient of 0.74.The glitches of PSR J0537−6910 also show this correlation, and follow the general trend suggested by the other pulsars (Ho et al. 2022).
We assume a functional form as in Equation 13 for a possible Δ -interglitch time interval relation, and investigate the presence of a correlation with either forward (Δ +  ) or backward waiting time (Δ −  ), since indications of both have been previously found (Middleditch et al. 2006;Antonopoulou et al. 2018;Lower et al. 2021).Using our sample, we compute the Spearman's rank correlation coefficient   , as well as the exponents  and , using a bi-variate Gaussian likelihood (Hogg et al. 2010) to account for the Gaussian  2022) (H22), and our data (JBO).The third column is the cut-off in Δ/, the glitches below this threshold are not included, and the waiting time are therefore recalculated between the remaining large or intermediate size glitches.The forth column is the glitch parameters we are fitting in the right hand side of Equation 13.In addition, we test the correlation for 7 'Vela-like' pulsars (PSRs B1727−33, B1757−24, B1800−21, B1823−13, J0205+6449, J2021+3651, and J2229+6114) in our sample as well.Δ −  represents the backward waiting time, while Δ +  is the forward waiting time.uncertainties in both Δ   and  i .The results are presented in Table 5. Figure 9 shows the data and the median power-law model for both Δ +  and Δ −  .We note that in our timing solutions   is typically the same for each interglitch interval (   = ), except in the very few cases for which a change Δ  was included in the glitch model.On the other hand, Lower et al. (2021) used   measurements for each interval.Nonetheless our results largely agree with theirs and we confirm the correlation between   and |Δ   |/Δ +  .The difference between our result and Lower et al. ( 2021) is due to the different sample of pulsars; for the glitches that are in both samples we get consistent measurements.Our sample has a larger number of glitches in pulsars with  < 2 × 10 −23 , whereas Lower et al. (2021) have only one, and we once again note the inherent bias in our sample as pulsars were selected on the basis of having at least one large glitch.

Dataset
The Δ   values used in the above investigation stand for the total change in  at a glitch.Of this, part recovers exponentially whilst the rest (Δ  p, ) is treated as 'permanent' in our timing solutions and recovers linearly.The  values represent the linear evolution of  in between glitches.In most cases, the exponential decay timescales are shorter than the interglitch waiting times, hence the sum of any Δ  p, terms has decayed to negligible levels, compared to Δ  p, at the time of the next glitch.Therefore, |Δ   |/  can be viewed as the time required for  to return to the value it had right before glitch .The situation where a following glitch happens when |Δ | of the previous glitch fully recovers (  = |Δ   |/Δ +  ) is represented by the black line in the top panel of Figure 9.As we can see, this scenario is not inconsistent with our best-fit parameters  and , given their uncertainties, and in fact the median  of our power-law fit is very close to one.
From the bottom panel of Figure 9 a trend is also apparent for the case where the backwards interglitch time intervals have been used, though this correlation is weaker (as reflected in the smaller correlation coefficient in Table 5).Such a relation would not be unexpected within the superfluid glitch model: a longer waiting time means that more of the superfluid that decoupled at glitch ( − 1) had time to recouple before glitch , and can thus decouple anew leading to a larger |Δ   |.
Since  describes the approximately linear evolution of , we find it informative to examine how the correlations change if we consider the permanent change, Δ  p , rather than the instantaneous change, Δ  (which is the sum of the permanent change and any exponentially-relaxing changes).If the superfluid components which relax exponentially have time to recouple almost fully before a subsequent glitch, then Δ −  will mostly determine Δ  p, instead of the total Δ   .We therefore include in Table 5 the results for the correlations -|Δ  p, |/Δ +  and -|Δ  p, |/Δ −  .In accordance to our earlier argument, including or not the Δ  d components has minimal impact on the correlation with the forward waiting time, whilst the correlation with the backwards waiting time gets significantly stronger when only the permanent changes |Δ  p, | are considered.
Following Antonopoulou et al. (2018) and Fuentes et al. (2019) we also repeat the analysis focusing only on intermediate and large glitches, as mostly these affect  (whilst, for example, the small glitches in PSRs B1800−21 and B1823−13 have no noticeable effect on ).We therefore ignore glitches with Δ/ < 10 −7 and find that both correlations improve considerably and their correlation coefficients are very close to unity.
Ideally, when exploring the above correlations,  +  (post-glitch ) and  −  (pre-glitch ) should be used with Δ +  and Δ −  respectively.Nonetheless, our use of the average  is justified because  is often not well defined over a single interglitch interval, and we do not find strong evidence for a Δ  in most glitches.Such small changes in  of a given pulsar can introduce some (vertical) scatter in Figure 9 but the overall correlation still emerges as  for the collection of pulsars varies by ∼ 4 orders of magnitude.Using the overall  means, however, that we cannot check whether the correlations hold for each pulsar of our sample individually, as is the case for PSR J0537−6910 (Ho et al. 2022).
It is possible, given the large scatter in Figure 9, that one correlation is the artefact of the other.The reason is that -for many of the pulsars included -the waiting times distribution is rather narrow (see for example PSR B1757−24 with a characteristic exponential terms were preferred, and for 2 glitches three exponential terms were preferred.Pulsars with higher spin-down rate | | typically present larger decaying amplitudes Δ d at their glitches.Even though our sample is small, the emerging distribution of exponential timescales appears powerlaw-like, therefore the common monitoring programs that observe pulsars biweekly or monthly are not adept at capturing the full post-glitch recovery signature.
We confirmed the power law relation between  i and |Δ   |/Δ +  , where Δ +  is the 'forward' waiting time after glitch  until glitch  + 1 takes place, proposed by Lower et al. (2021).However we also find that the correlation holds nearly as well when the 'backward' waiting time Δ −  is used instead.According to our results, a major contribution to the instantaneous changes |Δ | comes from the persisting component |Δ  p | which correlates with its rate of recovery .An even stronger correlation exists between |Δ p | and | |.
To cancel out the effects caused by the | |-| | and Δ p -| | correlations we compare the ratio Δ  p /  to waiting times Δ.The data for the entire pulsar sample loosely follows the Δ  p, /  ±  = Δ ±  line, however the scatter is considerable and most individual pulsars have too few glitches to confirm if the correlation holds for each one of them.Therefore, currently, the power of this relation to predict glitch epochs with accuracy is limited; it can however serve as a general guide for observations.Glitch characterisation and parameters definition with the use of Bayesian model selection continues, as new glitches are being detected by the JBO timing programme.Our results revealed a range of post-glitch recovery timescales, from few days up to several hundred days.Regular pulsar monitoring is necessary for detecting new glitches and to resolve their subsequent long-term recoveries.It is also, however, critical to design observing programs able to detect glitches promptly and adapt the monitoring strategy in real-time in order to resolve the early phase of post-glitch relaxation.

Figure 1 .
Figure1.The evolution of  and  in PSRs J0205+6449, B0355+54, J0611+1436, J0631+1036, J0729−1448 and B0919+06.The top mini panel and other three panels show respectively: The epochs of ToAs by vertical blue lines; residual , which is  after subtracting  f and  f ( −  f );  after further subtracting the permanent changes due to glitches  gp = Δ p + Δ  p Δ + 1 2 Δ  p Δ 2 and 1 2  f ( −  f ) 2 , as well as some constant offset so that the point before the first glitch is zero;  as a function of time.The reference epoch  f , reference frequency  f , and its derivatives (  f ,  f ) are some constants given in Appendix A in TableA1.Glitches best-modelled as steps in  and  (without exponential recoveries) are indicated with dotted lines at their glitch epochs.Dash-dotted lines mark the epoch of glitches for which one exponential term was included in the best model, dashed lines for those with two exponentials, and solid lines for those with three exponentials.An analytical model derived from the best-fit glitch parameters is shown as a black-solid curve.We use green-solid lines to mark gaps longer than 10 times the pulsar's average cadence in ToAs.The stride fitting results for  and  are shown as the red and purple dots with errorbars respectively.MNRAS 000,1-22 (2024)

Figure 8 .
Figure 8.  −  diagram showing the  of pulsars in our sample (stars) and glitching pulsars from the pulsar catalogue (squares).Two young pulsars PSRs J0534+2200 (the Crab pulsar, in red pentagon) and B1509−58 (in green diamond) with  ∼ 3 are highlighted.The Vela pulsar is also highlighted in cyan hexagon for comparison with the 'Vela-like' pulsars.The color-bar represents the long-term .We draw 100 samples from the posterior of , ,  and shade the region in blue where the modelled  is less than that expected from  b = 3 braking.The magenta line shows this /  2 = 3 threshold for the mean solution.The green dot-dashed lines correspond to constant  values according to the best fit of Equation.12.

Figure 14 .
Figure 14.The correlation between exponential decaying amplitude and other glitch parameters.Black points are the positive exponential terms, and red points represent exponentials with negative amplitudes.Top panel: The timescale of exponential recoveries  d in pulsars from our sample and their corresponding amplitude |Δ d |.Middle panel: The spin frequency  0 and the amplitude of exponential recoveries |Δ d |.Bottom panel: The absolute value of frequency derivative |  0 | and the amplitude of exponential recoveries |Δ d |.

Table 1 .
The pulsars used in this work, the MJD range over which they were observed at JBO, number of ToAs in this range, and the total number of glitches that occurred in this period.In the last column we show the number of glitches for which we analysed the post-glitch recovery.

Table 2 .
The details of post-glitch models we fitted for in model selection.

Table 3 .
The spin frequency, its derivatives, and fit epoch, obtained from fitting the polynomial timing model (Equation1), and the derived physical properties for the pulsars in our sample derived from the data introduced in Table1.The braking indices of pulsars with unconstrained  are omitted.
typically dominated by either the transient, exponential, recovery, or by timing noise.The  time series are typically dominated by glitches and their recoveries, though depending on the pulsar, this may be the 'linear' recovery of  that is measured as , or the exponential recoveries.For a handful of low | | pulsars, the evolution of  mostly shows noise.

Table 4 .
The glitch epochs, sizes, recovery amplitudes, and corresponding time scales of 85 glitches in 35 pulsars from our sample for which model selection was performed.The data span used in the model selection for each glitch is given in the ToA span column, with the corresponding number of ToAs in the next column.Errors are given in parentheses in units of the last quoted digit.For glitches with exponential recoveries, we list the amplitude, timescale, and the recovery ratio  of the recovery components in the last three columns.The values quoted are the mean of posterior with standard deviation.The five small glitches with asterisks are also included in the model selection process.