Pulsar Glitches: A Review

Abstract

∼$6\%$ of all known pulsars have been observed to exhibit sudden spin-up events, known as glitches. For more than fifty years, these phenomena have played an important role in helping to understand pulsar (astro)physics. Based on the review of pulsar glitches search method, the progress made in observations in recent years is summarized, including the achievements obtained by Chinese telescopes. Glitching pulsars demonstrate great diversity of behaviours, which can be broadly classified into four categories: normal glitches, slow glitches, glitches with delayed spin-ups, and anti-glitches. The main models of glitches that have been proposed are reviewed and their implications for neutron star structure are critically examined regarding our current understanding. Furthermore, the correlations between glitches and emission changes, which suggest that magnetospheric state-change is linked to the pulsar-intrinsic processes, are also described and discussed in some detail.

1. Searching for Pulsar Glitches

Pulsars emit radiation across the whole electromagnetic spectrum, and the most plausible geometric explanation for the observed emission is the lighthouse model [1]. The electromagnetic radiation from pulsars is originated and beamed from above their magnetic poles. The magnetic axis of a pulsar has a certain angle with respect to its rotation axis. If the earth is located in the range of emission beams of the pulsar, a series of periodic pulse signals are received by telescopes as pulsars rotate. Among them, millisecond pulsars (MSPs) have the greatest long-term stability of period. The MSP PSR J0437–4715 is considerably more stable than standard atomic clocks [2]. In recent years, the international pulsar timing array project (PTA), which combines observations of a set of pulsars to search for correlated signatures in the arrival times of pulses [3], has achieved unprecedented development. Two best known applications of this project are the detection of ultra-low frequency (∼${10}^{-9}$–${10}^{-8}\mathrm{Hz}$) gravitational waves [4] and establishment of a pulsar-based timescale [5]. Besides, binary millisecond pulsars constitute most extreme physical laboratories and give opportunity to test several theoretical models including the general theory of relativity [6].

As of now, there is no fully self-consistent pulsar emission model capable of explaining the structure of the magnetic field and the time-variable phenomena [7,8,9]. Nevertheless, the simple magnetic dipole model allows a basic understanding of pulsar magnetosphere and the coherent radio emission in pulsars [10,11,12]. As illustrated in Figure 1, an induced electric field exists throughout the magnetised neutron star as a result of its rotation. Under this circumstance, charged particles are lifted out of the solid crust of the star, filling the magnetosphere with a dense plasma [13]. These particles are accelerated to relativistic velocities, and emit $\gamma$-ray photons by either curvature radiation or inverse Compton scattering on lower energy photons [14,15,16]. And then, the photons interacting with the magnetic field forms electron-positron pairs [17]. Near the neutron star’s polar caps, the flow of particles creates electric currents that produce the observed radio beams. Light cylinder is an imaginary surface at which co-rotation of the closed dipolar magnetic field lines with neutron star breaks-down and field lines open up to space. At light cylinder radius $R_{\mathrm{LC}}=cP/2\pi$ with P being the rotational period of the underlying neutron star the plasma speed reaches the speed of light c. The currents flowing out from the magnetosphere along the open field lines that cross the light cylinder, as well as the pulsar wind of plasma, exert a torque on the magnetic field lines which then slows the rotation of the neutron star. There is still a lack of reliable observational evidence and theoretical models to accurately describe the coupled magnetic, thermal and spin evolution of pulsars [18,19,20,21]. Fortunately, there is a powerful approach known as “pulsar timing”, which helps to obtain the accurate information on the pulsar spin-down and refine existing radiation models.

Pulsar timing is a diagnostic technique in the rotational behaviour by tracking each pulse’s time of arrival (TOA), and involves the evaluation and interpretation of the mismatch between measured and predicted values. Although the individual pulses measured during the observation are highly variable in shape and intensity [23], several hundreds or thousands of pulses can be summed in time, frequency and polarisation to produce a highly stable integrated pulse profile over the timescales of years [24]. In order to obtain a series of TOAs, each of the integrated pulse profile is cross-correlated with a high signal-to-noise ratio “standard” profile. Usually, rapidly rotating pulsars with narrow pulse profiles have higher accuracy of TOAs. After correcting for TOAs to the Solar System Barycentre in order to remove the Earth’s orbital motion effects, the evolution of the phase of the pulse sequence can be expressed as a truncated Taylor series [25]:

$$\varphi \left(t\right)={\varphi }_0+{\nu }_0(t-t_0)+\frac{1}{2!}{\dot{\nu }}_0{(t-t_0)}^2+\frac{1}{3!}{\ddot{\nu }}_0{(t-t_0)}^3,$$

(1)

where t is the observation time. ${\varphi }_0$, ${\nu }_0$, ${\dot{\nu }}_0$ and ${\ddot{\nu }}_0$ represent the phase, spin-frequency and its first and second derivatives as measured at the reference time $t_0$. Regardless of the effects of matter outflow and accretion, the spin-down of an isolated radio pulsar is mainly dominated by the magnetic dipole radiation, so that the accurate predictions of the arrival time of every one of its pulses can be made from the timing model given by Equation (1).

However in an actual long term observation, it is not always the case. The slow-down law of pulsars predicted by Equation (1) is generally perturbed by two types of timing irregularities: timing noise and glitches [26]. Timing noise is a slow, long-term, discernible, stochastic wandering in the periods of pulsars. In 1969, two independent groups, Radhakrishnan & Manchester [27] and Reichley & Downs [28] presented the first ever observation of a sudden increase in PSR J0835–4510’s (Vela pulsar) spin frequency $\nu$, known as “glitch”. Obviously, for some pulsars pulse arrival times are not well predicted by the simple spin-down model, reflected in the observed irregularity of the pulsar’s rotation. The discovery of glitches have naturally raised the curiosity to search for more similar events using telescopes around the world. Now, high cadence monitoring of an exceedingly large number of pulsars are conducted by the northern and southern hemisphere radio telescopes: 76-m Lovell [29] and 64-m Parkes [30]. Most of glitches are detected with these two telescopes. A database of all known pulsar glitches is regularly updated in two different locations: one at ATNF (https://www.atnf.csiro.au/people/pulsar/psrcat/glitchTbl.html—accessed on 30 September 2022) [31] and also by Jodrell Bank (JBO) (http://www.jb.man.ac.uk/~pulsar/glitches/gTable.html—accessed on 30 September 2022) [32]. After combining entries from both the catalogues there are a total of 719 glitches identified in 239 pulsars [33]. As will be seen in Section 3 below, these glitches can be further classified into four subclasses of normal glitches, slow glitches, glitches with delayed spin-ups, and anti-glitches.

These observed glitches are modelled using an additive function as below [25]:

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{g}}\left(t\right)=& \Delta \varphi +\Delta {\nu }_{\mathrm{p}}(t-t_{\mathrm{g}})+\frac{1}{2!}\Delta {\dot{\nu }}_{\mathrm{p}}{(t-t_{\mathrm{g}})}^2\hfill \\ \hfill & +\sum _i^N\Delta {\nu }_{\mathrm{d}}^{\left(i\right)}{\tau }_{\mathrm{d}}^{\left(i\right)}\left[1-exp\left(-\frac{(t-t_{\mathrm{g}})}{{\tau }_{\mathrm{d}}^{\left(i\right)}}\right)\right],\hfill \end{array}$$

(2)

whose value is zero if $t\le t_{\mathrm{g}}$. $\Delta \varphi$ is an offset in pulse phase. The positive or negative sign of the permanent jumps in the spin frequency $\Delta {\nu }_{\mathrm{p}}$ at the time of a glitch $t_{\mathrm{g}}$ represents usual glitches or anti-glitches, respectively. $\Delta {\dot{\nu }}_{\mathrm{p}}$ is the permanent change in the frequency derivative (spin-down rate) relative to the pre-glitch value. The fourth term models the ith component of N exponential recoveries with a transient frequency component $\Delta {\nu }_{\mathrm{d}}$ and decay time constant ${\tau }_{\mathrm{d}}$. $\Delta {\nu }_{\mathrm{d}}>0$ corresponds to the amplitude of the exponential recovery. $\Delta {\nu }_{\mathrm{d}}<0$ indicates the existence of a delayed spin-up component after the initial unresolved step. As for the criterion of slow glitches, a continuous increase in frequency and sudden decrease followed by a gradual increase in spin-down rate $|\dot{\nu }|$ are recognised from plots of $\nu$ and $\dot{\nu }$ versus time. Hence, the observed absolute jumps $\Delta {\nu }_{\mathrm{g}}$ and $\Delta {\dot{\nu }}_{\mathrm{g}}$ in a glitch are described as:

$$\Delta {\nu }_{\mathrm{g}}=\Delta {\nu }_{\mathrm{p}}+\sum _i^N\Delta {\nu }_{\mathrm{d}}^{\left(i\right)},$$

(3)

$$\Delta {\dot{\nu }}_{\mathrm{g}}=\Delta {\dot{\nu }}_{\mathrm{p}}-\sum _i^N\frac{\Delta {\nu }_{\mathrm{d}}^{\left(i\right)}}{{\tau }_{\mathrm{d}}^{\left(i\right)}}.$$

(4)

In addition, the relative glitch amplitudes in frequency and frequency derivative are $\Delta {\nu }_{\mathrm{g}}/\nu$ and $\Delta {\dot{\nu }}_{\mathrm{g}}/\dot{\nu }$, respectively. The degree of recoveries are quantified by the parameter $Q={\sum }_i^N\Delta {\nu }_{\mathrm{d}}^{\left(i\right)}/\Delta {\nu }_{\mathrm{g}}$. Noteworthy, the maximum change of the spin frequency $\Delta {\nu }_{\max }$ and its first time-derivative $\Delta {\dot{\nu }}_{\max }$ are derived from comparing the timing solutions away from the spin-up event to describe the slow glitch effect. Uncertainties associated with these parameters are obtained using standard error propagation methods.

For decades, the usual procedures of data reduction for radio pulsar glitches are as shown in Figure 2. The two key steps are the measurements of the TOAs and implementation of phase-connected timing residuals. The popular pulsar data analysis packages psrchive [34] provide functionality for generating TOAs; and then, these TOAs are used to fit with a model that contains a set of known parameters to achieve a phase-coherent timing residuals ($R=(\varphi -N)/\nu$, where N is the nearest integer to $\varphi$) [35], using softwares such as tempo [36], psrtime [35], tempo2 [35], and pint [37]. tempo2 remains by far the most popular software packages in the bunch of pulsar timing tools. Using these tools, the substantial manual efforts are required to obtain a satisfactory phase-connected timing residuals. Until recently, Freire et al. [38] and Phillips et al. [39] have developed new automated algorithms, DRACULA and APT, to determine phase-connected timing solutions for pulsars. Therefore, pulsar glitches can be identified via visual inspection with the sudden phase discontinuity in the timing residuals relative to a solution based on earlier data. Glitches manifest as timing residuals becoming increasingly negative with time and new models are needed for post-glitch behaviours. In the case of frequent occurrence of glitches with particularly large sizes, it is almost impossible to obtain good post-glitch timing solutions. Melatos et al. [40] presented a complementary method to this problem in the form of the Hidden Markov Model (HMM) algorithm, which measures the permanent changes in spin frequency $\Delta {\nu }_{\mathrm{p}}$ and spin-down frequency $\Delta {\dot{\nu }}_{\mathrm{p}}$ during glitches with a HMM. Dunn et al. [41] employed HMM to analyze low-cadence timing data to constrain the magnitude of missing glitches in pulsars. Singha et al. [42] presented a real-time automated glitch detection pipeline (AGDP) that enables statistical tests to determine phase coherence between TOAs and has been implemented at the Ooty Radio Telescope (ORT) to recognise glitches automatically. Soon, the searching techniques of glitches with massive human intervention will be consigned to the past.

In China, some scientific achievements have been made in the area of pulsar glitch research.

Table 1. Chinese telescopes to perform a search for glitch events.

Telescopes	Location	Diameter	Frequency	Pulsars	Start	No. of Glitches	Ref.
NSRT	Nanshan	26 m	1.54 GHz	∼300	2000	103	[43,44,49,50,51,52,53,54,55,56,57,58,59,60]
TMRT	Shanghai	65 m	2.25/4.82/8.60 GHz	∼100	2013	5	[61,62,63,64]
Insight-HXMT	Space	–	1–250 keV	∼10	2017	3	[45,65]
KMRT	Yunnan	40 m	2.25/4.85 GHz	∼90	2008	1	[66]
HRT	Luonan	40 m	1.40 GHz	∼10	2014	1	[67]
XPNAV-1	Space	–	0.5–10 keV	∼26	2016	1	[46]
FAST	Pingtang	500 m	1.25 GHz	∼350	2016	–	[68]
QTT	Qitai	110 m	–	–	–	–	[47]
JRT	Jingdong	120 m	–	–	–	–	[48]

presents the information on the number of pulsars currently timed by Chinese telescopes and the number of glitches discovered by each telescope. Notably, frequent observations spanning 22 years with the Nanshan 26-m radio telescope (NSRT) have detected 103 glitches, which account for about $15\%$ of all published glitches. The first slow glitch of PSR B1822–09 [43] and significant changes in pulse profile associated with PSR B2035+36’s glitch activity [44] were discovered by NSRT. Ge et al. [45] and Zhang et al. [46] reported that a delayed spin-up process occurred before the normal recovery process in the November 2017 Crab (PSR B0531+21) glitch with Insight-HXMT (Hard X-ray Modulation Telescope) and XPNAV-1 (X-ray Pulsar Navigation satellite) data, respectively. Currently, about 350 pulsars have long been become a continuous timing target with the largest dish FAST (Five-hundred-meter Aperture Spherical Radio Telescope), which has the capacity to follow more than 2000 pulsars. In the future, next-generation telescopes QTT (Qitai 110-m Radio Telescope) [47] and JRT (Jingdong 120-m Radio Telescope) [48] will provide an opportunity to track more pulsars, more frequently, allowing glitches to be investigated in more detail.

2. Properties of Pulsar Glitches

The known glitching pulsars are shown on the the period–period-derivative ($P-\dot{P}$) diagram in Figure 3. Almost all published glitches were identified by visual inspection of the phase residuals, relative to a predicted model. Glitches with fractional sizes of $\Delta \nu /\nu$∼${10}^{-9}$ do not lose phase coherence over several hundred days, while those with sizes $\Delta \nu /\nu$∼${10}^{-6}$ may see a change in residuals of a large fraction of the pulse period in just a few days [69]. Once the glitch has been determined, groups of TOAs that overlap are fitted to demonstrate the time evolution of $\nu$ and $\dot{\nu }$ around the glitch. To illustrate the permanent and transient jumps in one glitch, the time-dependence of the frequency residuals $\Delta \nu$ relative to the pre-glitch spin-down model are generally presented as in Panel (a) of Figure 4. Also, an expanded plot of frequency residuals $\Delta \nu$ where the mean post-glitch value has been subtracted from the post-glitch data is derived to clearly show the exponentially changing components, as seen in Panel (b). The variations of spin-down rate $\dot{\nu }$ is visualized in detail as shown in Panel (c). Janssen & Stappers [70] found that the minimum glitch size could be detected with a threshold of $\Delta \nu /\nu$∼${10}^{-11}$ by performing Monte Carlo simulations to test the glitch detection method on PSR J0358+5413 [71]. The millisecond pulsar J0613–0200 has suffered a glitch with a size $\Delta \nu /\nu$∼$2.5\times {10}^{-12}$, which is the smallest that has been recorded [72]. Several studies have indicated that the previously published pulsar timing results did not include all glitches in individual pulsar or in particular datasets [70,73,74]. The following reasons may result in incomplete observations of glitches:

• Given the limited resources available to observatories and the large numbers of pulsars, timing observations of some pulsars are not carried out around the glitch [75].
• Making low-cadence, short dwell time observations near a $\dot{\nu }$ transition epoch may result in misinterpreting the timing behaviour, leading one to omit a glitch [50,76].
• Glitches are too small to be resolved due to being below the present limit of detectability [32].
• Pulsars that exhibit high levels of timing noise may have small glitches that go undetected in the data [49].
• If small glitches occur during the initial stages of recovery following a large glitch, they could be missed because of the need of timing fit [77,78].
• When a glitch is identified to occur in a significant observation gap, it is impossible to distinguish between a single glitch or multiple, smaller glitches [70,73].

Around $6\%$ of known pulsars have been observed to glitch at least once since the first discovery of glitch in 1969. Out of the 239 glitching pulsars, 144 have only glitched once or twice. Clearly, glitches occur predominantly in young pulsars (with characteristic ages ${\tau }_{\mathrm{c}}=P/2\dot{P}\sim {10}^3$–${10}^7$ year [69]) located in the top right of the $P-\dot{P}$ diagram, but the size of and time between glitches do not depend only upon the pulsar’s position in the diagram [79] (Figure 3). The curious glitch phenomena have been discovered from ordinary pulsars, to super-powerful magnetic field harbouring magnetars [80], to millisecond pulsars [81], and are predominantly seen from young and high spin-down rate pulsars. PSRs J0534+2200 (Crab), J0537–6910, J0835–4510 (Vela), J1341–6220, J1740–3015 are remarkable examples of displaying frequent glitches. A significant linear correlation between glitch size $\Delta \nu$ and time to next glitch has been established for the big glitcher PSR J0537–6910 [82,83], which has undergone 53 glitches by September 2022 [84]. The Vela and Crab pulsars are the most widely studied glitching pulsars. The Vela pulsar, first observed in 1968, has since then had 24 glitches, with a frequency jump of $\Delta \nu /\nu$∼${10}^{-6}$ for most of them. The post-glitch decay starts out exponential relaxation, with different timescales from 1 min to several hundred days, but it eventually exhibits linear recovery of spin-down rate in time with a constant second time derivative of spin frequency [49,85]. In contrast, the glitch sizes of the Crab pulsar ($\Delta \nu /\nu$∼${10}^{-7}$–${10}^{-9}$) are about at least an order of magnitude smaller than those for the Vela pulsar. There is usually a sudden, exponential decay in the post-glitch behavior over several days, as well as a persistent increase in $\dot{\nu }$ [49,86]. Long-term timing observations of Crab- and Vela-like pulsars revealed their glitch behaviours were similar to the two pulsars. More interestingly, two micro-glitches have been observed from the millisecond pulsars (MSPs J0613–0200 [72] and B1821–24 [81]), which are thought to be the most stable type of pulsar population. Even so, if sufficient post-glitch timing observations are available to model the glitches, the influence on the timing precision for pulsar timing array applications, such as gravitational wave (GW) detection, can be eliminated. McKee et al. [72] indicated that the rate of glitches for MSPs is significantly lower than that of the general population, with a probability of ∼50% that another glitch will be seen in a timing array programme of pulsar within 10 years. Moreover, glitches occurred in magnetars are well known to have large sizes, and sometimes coincide with emission changes and outbursts [87]. 31 magnetars have been discovered, 7 of which have exhibited a total of 22 glitches (http://www.physics.mcgill.ca/~pulsar/magnetar/main.html—accessed on 30 September 2022) [88,89]. Chukwude & Urama [90] found that glitches have amplitudes $\Delta \nu /\nu \gtrsim {10}^{-9}$ and are associated increases in $\dot{\nu }$ (i.e., $\dot{\nu }$ is more negative), while micro-glitches are smaller in amplitude ($\Delta \nu /\nu \lesssim {10}^{-10}$) and can have either positive or negative changes in $\dot{\nu }$. According to Yu et al. [69], after a glitch, part of the sudden changes in both $\nu$ and $\dot{\nu }$ often recover exponentially, and $\dot{\nu }$ then continues to decay linearly until the next glitch.

3. Classification of Glitches

3.1. Normal Glitches

To-date, more than 200 pulsars have collectively shown hundreds of (> 650) normal glitches. The fractional amplitude $\Delta \nu /\nu$ of these events are in the range of ${10}^{-12}$ [72] to ${10}^{-3}$ [91]. Almost-instantaneous rise times and lengthy recoveries are common features in these glitches. Figure 4 shows a typical glitch with an exponential decay timescale of $13\left(2\right)\mathrm{d}$ in the Vela pulsar (PSR J0835–4510) occurred in July 2010. The rapid spin-up timescales of glitches can be assumed to be less than $5\mathrm{s}$ [92], $30\mathrm{s}$ [85], and $40\mathrm{s}$ [93], according to the results of glitch analysis for the Vela pulsar carried out by the southern hemisphere observers. Often, the instantaneous rise is followed by a temporary increase in the spin-down rate. For the post-glitch behaviours, the response of the spin frequency is diverse. Across the population of glitches, the evolution of the post-glitch rotation rates is usually modelled as a combination of exponential and linear recoveries to the pre-glitch values [69,94]. The former last from a few days to several months and the latter persist on a time-scale of years. Perhaps a more striking aspect of these glitches is that the transient part has been observed to require multiple decaying exponentials for the best timing fit. PSRs J0835–4510 (ten times), J1119–6127 (once), J1757–2421 (once), J1803–2137 (once) and J2337+6151 (once) have shown evidence of two exponential decay terms in same glitch [57,66,95,96]. More specifically, Flanagan [97] reported that the December 1988 Vela glitch can be best explained with three exponential terms and Dodson et al. [93] confirmed four exponential terms to interpret the January 2000 Vela glitch satisfactorily. In the early days following the first glitch observation in the Vela pulsar, Baym et al. [98] proposed the “two component model” that normal glitches are the outcome of the relaxation of the neutron-star superfluid after crust quakes to the current equilibrium oblateness. Glitch recoveries are thought to be the first evidence for the existence of a superfluid interior of the neutron star [98,99]. Soon after Anderson & Itoh [100] suggested that normal glitches are driven by the transfer of angular momentum from unpinned vortex lines within superfluid to the solid crust. Arguably, normal glitches contain invaluable information on the interior of neutron stars. Therefore, general glitch behaviours can be used to investigate not only the structure of the pulsar magnetosphere, but also the dynamical process operating inside a neutron star.

3.2. Slow Glitches

A total of 31 slow glitches have been observed in 7 pulsars [102,103], for which the spin frequency $\nu$ slowly increases whereas the spin-down rate $\dot{\nu }$ decreases over several months before $\dot{\nu }$ relaxes back to its original state over a slightly longer time scale (Figure 5) [43,49,104]. This process repeats periodically, and the peak value of the glitch amplitude is the same for each time. The maximum changes in the spin frequency $\nu$ and spin-down rate $|\dot{\nu }|$ for these slow glitches range from 2.3–46$\mathrm{nHz}$ and $0.15\times {10}^{-15}$–$3.15\times {10}^{-15}{\mathrm{s}}^{-2}$; the corresponding relative maximum changes are $\Delta \nu /\nu$∼$0.9\times {10}^{-9}$–$31.2\times {10}^{-9}$ and $|\Delta \dot{\nu }|/|\dot{\nu }|$∼$1.8\times {10}^{-3}$–$23.6\times {10}^{-3}$, respectively. The occurrence of slow glitches and normal glitches in PSR B1822–09 is of special interest [105]. Zhou et al. [102] linked the variations of the integrated pulse profile to the detection of a very large slow glitch in PSR J1602–5100. At the beginning and end of this slow glitch there was the presence and absence of a new trailing component, respectively. Theoretically, the angular momentum exchange model [106,107], which is currently widely accepted to explain the physical mechanism of normal glitches, cannot be applied to the process of slow events. Link & Epstein [108] suggested that the slow glitches are the result of a sudden increase in the temperature of the inner crust of the neutron star after a quake. However, Hobbs et al. [109] suggested that the slow glitches are just a different manifestation of timing noise, which is significantly dominated by sustained random wandering in either the phase, spin, or spin-down rate [110]. That means slow glitch events cannot be classified as a new subclass of pulsar glitches, and may be more widespread throughout the pulsar population. However, slow glitches have not been detected in any other glitching pulsars so far. While curious, the lengthy rise time, periodic recurrence and uniform amplitude of these slow glitches strongly indicate that they are a separate, though possibly related, category of event.

3.3. Delayed Spin-Ups

Figure 6 demonstrates the evolution of the rotation frequency $\nu$ and its first derivative $\dot{\nu }$ in the delayed spin-up of the 2019 Crab pulsar glitch. It is obvious that glitch with delayed spin-up component can be recognized by a short period of the $\nu$ increase exponentially immediately after the glitch occurrence and followed by a rapid decay of $\dot{\nu }$. At present, delayed spin-up processes have been detected in the glitches of three pulsars, Crab (six times) [45,75,77,111,112], 1E 2259+586 (once) [113] and SGR J1935+2154 (once) [114]. The timescale of these fast rising transient components ranges between 0.5 days and 14.1 days. The capture of delayed spin-up component in the multiple glitches from the Crab pulsar are partly attributed to the high and regular observing cadence at the Jodrell Bank Observatory. On the contrary, if the post-glitch transient variations are prevalent in glitches, the low and/or irregular cadence observations are the main reasons for their not being detected. Although the upper limit of the rise timescale for the Vela pulsar glitches have been constrained to within 5 s [92] by southern hemisphere observers through monitoring near-continuously, no delayed spin-up components have been observed. This may be another evidence that the glitches of the Crab and the Vela pulsars have a different origin [115]. It is worth noting that these observed delayed spin-up behaviours were only detected in young pulsars during their big glitches, and there is a slow exponential recovery process after the delayed spin-up. In particular, the total glitch magnitude of the 2020 SGR J1935+2154 event including the delayed increase is the largest among the all glitches observed to date, and the transient decay timescale ∼$8\mathrm{d}$ is much longer than that of the Crab’s glitches. Ge et al. [114] found that the pulse profile of SGR J1935+2154 changed rapidly and drastically after this unusual large glitch. Combined with the detection of repeated FRB 200428 after the glitch with a delay of likely 3 days, Ge et al. [114] argue that the internal changes in this young magnetar generate the large glitch and significantly alter the structure of the magnetosphere and trigger X-ray bursts and FRB 200428. Just recently Younes et al. [116] reported a possible connection between a new glitch occurred on 5 October 2020 and three FRB-like radio bursts in SGR J1935+2154. There is a power-law ($\alpha =2.0$) correlation between $\Delta \nu /\nu$ and the extended spin-up timescale $\tau$ of 6 glitches in 3 pulsars suggesting that the mechanisms behind these unusual glitches are similar [114]. The finding is also interpreted as generally consistent with the Gügercinoǧlu & Alpar [117] theory. Moreover, it is impossible to tell whether delayed spin-ups are only related to large jumps, if the delay time scale is too short or small jumps are too small to be observed by existing observation equipments. In short, the detection of delayed spin-ups offers a unique opportunity to study the micro-physical processes governing interactions between the neutron star interior and the crust. It was also proposed that SGR J1935+2154 glitch may arise from the ejection of a particle wind emanating from close to the magnetar polar cap which effectively extracts the angular momentum from the star [116].

3.4. Anti-Glitches

The so-called anti-glitches are first coined when Archibald et al. [118] observed a sudden negative jump of spin frequency in the magnetar IE 2259+586 red(Figure 7), accompanied by a X-ray outburst and a significant spin-down rate change. Recently, Ray et al. [119] reported three large anti-glitches in the accreting ultraluminous X-ray source NGC 300 ULX1 during its most recent outburst. Ducci et al. [120] evaluated the possibility that anti-glitches occur in accreting pulsars, and found that these spin-down glitches may take place more often in the superfluid vortex scenario. Until now, only a special class of pulsars, magnetars exhibiting contemporaneous radiative changes, have been observed to undergo anti-glitches. In total 7 anti-glitches are detected in 3 pulsars [80,118,119,121,122,123,124]. An & Archibald [125] speculated a possible anti-glitch occurred in magnetar CXOU J164710.2–455216 from the interrupted phase of residuals. Undoubtedly, the unexpected occurrence of anti-glitches are strongly challenging the standard glitch theories. The close correlation observed between spin-down glitches and outburst activity suggest that these events are related to the magnetospheric activity [115]. Tong [126] used the wind braking model to successfully interpret the anti-glitch behaviours of magnetar 1E 2259+586. Huang & Geng [127] suggested that anti-glitches occur as a result of the collision of a small mass body with pulsars. Garcia & Ranea-Sandoval [128] interpreted that anti-glitches are an inevitable consequence of the cumulative decay of the internal toroidal magnetic field component resulting in the rearrangement of the stellar structure. Howitt & Melatos [129] performed simulations in vortex avalanches scenario to find that anti-glitches are caused by interrupting the secular increase of the angular velocity. Whatever the cause, rethinking of the mechanism behind the glitches of all pulsars is needed.

4. Models of Pulsar Glitches

The properties of rotational glitches in a given pulsar, especially the relaxation of the spin frequency to a value slightly less than its extrapolated original pre-glitch level shed light into the multicomponent structure of neutron stars and lead to the view that glitches represent the exchange of angular momentum between these components. The literature on the mechanisms responsible for pulsar glitches is immense. Immediately following the first detection, pulsar glitches have been attributed to a wide variety of mechanisms by invoking processes occurring both in the magnetosphere and internal to the neutron star, such as accretion, magnetospheric instabilities, crustquakes, corequakes and catastrophic unpinning of vortex lines [131]. It was also conjectured that pulsar glitches may result from neutronization in the envelopes of neutron stars following loss of angular momentum and shallow increase in matter density [132]. Another possibility for large scale vortex unpinning is the vortex avalanches wherein an unpinned vortex segment encounters a pinned vortex and causing it to unpin which unpins other pinned vortices due to induced extra superfluid flow around them and collective unpinning develops in this way. Such a scenario is observed in simulations albeit for low number of vortices and point defects, i.e., lattice nuclei, compared to a real neutron star case [133,134,135,136,137,138,139]. In the early days of pulsar glitch research several experiments were devised with superfluid Helium II filled containers, which are accelerated and decelerated occasionally by hand to mimic the neutron star rotational evolution [140,141,142,143]. In recent years numerical simulations have been conducted to visualize several aspects of the pulsar glitch recovery [94,144,145,146,147,148,149,150,151].

4.1. Basics of Superfluid Vortex Dynamics for Neutron Star Rotational Evolution

Superfluidity is at the heart of the most glitch models [100,108,152,153,154,155,156,157]. Moreover, pulsar timing noise may arise from the erratic coupling of the superfluid with the normal matter inside a neutron star [158,159,160,161]. We begin this section by reviewing some key concepts of superfluid vortex lines for neutron star rotational dynamics and glitches. Those readers who want to learn more about the superfluidity of neutron stars can refer to the general review articles in the literature [162,163,164,165,166,167,168].

A superfluid can maintain rotation only if it is perforated with vortex lines. The number, distribution and interaction of these vortex lines fully determine the rotational behaviour of superfluid inside a container. In the case of neutron stars the container is its solid crust.

The equation of the observed crust’s rotation rate is given by the torque equilibrium on the neutron star due to decelerating external braking torque $N_{\mathrm{ext}}$ and internal superfluid torque $N_{\mathrm{int}}$ as

$$I_{\mathrm{n}}{\dot{\Omega }}_{\mathrm{n}}=N_{\mathrm{ext}}+N_{\mathrm{int}}=N_{\mathrm{ext}}-I_{\mathrm{cs}}{\dot{\Omega }}_{\mathrm{cs}}-I_{\mathrm{core}}{\dot{\Omega }}_{\mathrm{core}},$$

(5)

where the subscripts “n”, “cs” and “core” refer to normal matter, crustal superfluid, and core superfluid respectively, and I and $\dot{\Omega }$ are the moments of inertia and spin-down rates of the corresponding components. Recall that the rotation rate and spin frequency are related to each other via $\Omega =2\pi \nu$. In order to find the spin-down rate of the superfluid components one has to consider the superfluid’s circulation due to velocity field around each vortex line:

$$\oint {\overrightarrow{v}}_{\mathrm{s}}.d\overrightarrow{\ell }=2\pi r^2{\Omega }_{\mathrm{cs}}=N_{\mathrm{v}}\kappa =\int \int n_{\mathrm{v}}\kappa dS={\int }_0^{r}2\pi r^{\prime }{n}_{\mathrm{v}}\left(r^{\prime }\right)\kappa d{r}^{\prime }.$$

(6)

Here $\kappa =h/2m_{\mathrm{n}}$ with h and $m_{\mathrm{n}}$ being Planck constant and bare neutron mass, respectively is the vorticity quantum attached to each line and $N_{\mathrm{v}}$ is the total number of vortex lines inside a closed contour. If one takes spatial derivative of Equation (6) and considers rigid rotation for the superfluid, then

$$\frac{\partial \left(2\pi r^2{\Omega }_{\mathrm{cs}}\right)}{\partial r}=4\pi r{\Omega }_{\mathrm{cs}}=2\pi r{n}_{\mathrm{v}}\kappa \to n_{\mathrm{v}}=\frac{2{\Omega }_{\mathrm{cs}}}{\kappa }.$$

(7)

Equation (7) implies that the areal density of vortex lines $n_{\mathrm{v}}$ are fully determined by the rotation rate of the superfluid and in order for superfluid to spin down the vortices should be expelled from a given surface area. The time evolution of the vortex number density is given by the continuity equation:

$$\frac{\partial n_{\mathrm{v}}}{\partial t}+\overrightarrow{\nabla }.\left(n_{\mathrm{v}}{\overrightarrow{v}}_{\mathrm{L}}\right)=0.$$

(8)

If one takes the time derivative of Equation (6) and uses Equation (7) then

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\left(\oint {\overrightarrow{v}}_{\mathrm{s}}.d\overrightarrow{\ell }\right)& =2\pi r^2{\dot{\Omega }}_{\mathrm{cs}}=\frac{\partial }{\partial t}\left(\int \int n_{\mathrm{v}}\kappa dS\right)=-\int \int [\overrightarrow{\nabla }.\left(n_{\mathrm{v}}\kappa v_{\mathrm{L}}\right)]dS\hfill \end{array}$$

$$\begin{array}{c}\hfill =-\int \frac{1}{r}\frac{d}{dr}\left(r{n}_{\mathrm{v}}\kappa v_{\mathrm{L},\mathrm{r}}\right)2\pi rdr,\end{array}$$

(9)

where $v_{\mathrm{L},\mathrm{r}}$ is the vortex line’s velocity component in the radial direction. Equation (9) leads to the spin-down law for a superfluid:

$${\dot{\Omega }}_{\mathrm{cs}}=-\frac{n_{\mathrm{v}}\kappa v_{\mathrm{L},\mathrm{r}}}{r}=-\frac{2{\Omega }_{\mathrm{cs}}}{r}{v}_{\mathrm{L},\mathrm{r}},$$

(10)

where in the last step we have used Equation (7) for the vortex line number density. Equation (10) tells us that in order for superfluid to keep up with the rotation rate of the normal matter crust and thus spin down its vortices should migrate radially outward. Therefore, the rate of the superfluid spin change is fully dependent on the $v_{\mathrm{L},\mathrm{r}}$. In order to find the radial vortex line migration rate we consider the Magnus force equation,

$$\overrightarrow{F}={\rho }_{\mathrm{s}}\overrightarrow{\kappa }\times \left(\overrightarrow{v_{\mathrm{s}}}-\overrightarrow{v_{\mathrm{L}}}\right),$$

(11)

together with the force due to magnetised vortex electron scattering of the form [169,170,171]

$$\overrightarrow{F}=C\left(\overrightarrow{v_{\mathrm{c}}}-\overrightarrow{v_{\mathrm{L}}}\right),$$

(12)

The solution of system of Equations (11) and (12) gives

$$v_{\mathrm{L},\mathrm{r}}=\frac{\left(v_{\mathrm{s}}-v_{\mathrm{c}}\right)}{\left[\left(\frac{{\rho }_{\mathrm{s}}\kappa }{C}\right)\right]+\left[\left(\frac{C}{{\rho }_{\mathrm{s}}\kappa }\right)\right]}.$$

(13)

From Equations (10) and (13) the core superfluid obeys the spin-down law

$${\dot{\Omega }}_{\mathrm{core}}=-\frac{2{\Omega }_{\mathrm{core}}}{r}{v}_{\mathrm{L},\mathrm{r}}=-2{\Omega }_{\mathrm{core}}\frac{\left({\Omega }_{\mathrm{core}}-{\Omega }_{\mathrm{c}}\right)}{\left[\left(\frac{{\rho }_{\mathrm{s}}\kappa }{C}\right)\right]+\left[\left(\frac{C}{{\rho }_{\mathrm{s}}\kappa }\right)\right]}\equiv -\frac{\left({\Omega }_{\mathrm{core}}-{\Omega }_{\mathrm{c}}\right)}{{\tau }_{\mathrm{core}}},$$

(14)

where the core superfluid’s relaxation time is defined as

$${\tau }_{\mathrm{core}}=\frac{\left[\left(\frac{{\rho }_{\mathrm{s}}\kappa }{C}\right)\right]+\left[\left(\frac{C}{{\rho }_{\mathrm{s}}\kappa }\right)\right]}{2{\Omega }_{\mathrm{core}}}.$$

(15)

For the above considered process of magnetised vortex scattering of electrons ${\tau }_{\mathrm{core}}\simeq$ (10–200) P seconds [172]. The observed fast relaxation of the overshooting post-glitch spin rate after the 2016 Vela glitch within a minute [173] implies a time variable coupling of the core superfluid to the observed crustal rotation rate [174]. The form of radial motion of vortex lines in the crustal superfluid depends crucially on their microscopic interaction with the lattice nuclei [152,175]. In the vortex creep model [106,152], $v_{\mathrm{L},\mathrm{r}}$ is given by a microscopic trial rate of vortex motion around lattice nuclei times the transition probability due to finite temperature of the crust proportional to the Boltzmann factors, i.e.,

$$v_{\mathrm{L},\mathrm{r}}=v_0\left[exp\left(-\frac{E_{\mathrm{p},\mathrm{out}}}{kT}\right)-exp\left(-\frac{E_{\mathrm{p},\mathrm{in}}}{kT}\right)\right],$$

(16)

where $v_0$∼${10}^5$–${10}^7$ cm s⁻¹ is the microscopic trial velocity of a vortex line around nuclear clusters throughout the neutron star crust [176] and, $E_{\mathrm{p},\mathrm{out}}$ and $E_{\mathrm{p},\mathrm{in}}$ are defined as

$$E_{\mathrm{p},\mathrm{out}}=E_{\mathrm{p}}-\Delta E=E_{\mathrm{p}}\left(1-\frac{\omega }{{\omega }_{\mathrm{cr}}}\right),$$

(17)

and

$$E_{\mathrm{p},\mathrm{in}}=E_{\mathrm{p}}+\Delta E=E_{\mathrm{p}}\left(1+\frac{\omega }{{\omega }_{\mathrm{cr}}}\right).$$

(18)

Here $E_{\mathrm{p}}$ is the pinning energy for vortex-nucleus interaction, and ${\omega }_{\mathrm{cr}}$ is the maximum velocity difference between the superfluid and the normal matter that a vortex line can withstand before it unpins. The bias $\Delta E$ determines the height of potential barrier for a vortex line to overcome with subscripts “in” and “out” denote for vortex outward and inward motion, respectively. Note that when the crustal temperature drops below a few ${10}^5$ K the vortex line motion against potential barrier sustained by nuclear clusters proceeds through quantum tunnelling rather than classical thermal creep [175]. From Equations (16)–(18) we find radial component of the vortex line velocity in the crustal superfluid as

$$v_{\mathrm{L},\mathrm{r}}=2v_{0}exp\left(-\frac{E_{\mathrm{p}}}{kT}\right)sinh\left(\frac{E_{\mathrm{p}}}{kT}\frac{\omega }{{\omega }_{\mathrm{cr}}}\right).$$

(19)

From Equations (10) and (19) the crustal superfluid spin-down rate becomes

$${\dot{\Omega }}_{\mathrm{cs}}=-\frac{4{\Omega }_{\mathrm{cs}}}{r}{v}_{0}exp\left(-\frac{E_{\mathrm{p}}}{kT}\right)sinh\left(\frac{E_{\mathrm{p}}}{kT}\frac{\omega }{{\omega }_{\mathrm{cr}}}\right).$$

(20)

The angular velocity difference $\omega \equiv {\Omega }_{\mathrm{cs}}-{\Omega }_{\mathrm{c}}$ between the rotation rates of the crustal superfluid ${\Omega }_{\mathrm{cs}}$ and the observed crust ${\Omega }_{\mathrm{c}}$ evolves as [106,177]

$$\dot{\omega }=\frac{I}{I_{\mathrm{c}}}\left[-\frac{N_{\mathrm{ext}}\left(t\right)}{I}-\frac{4{\Omega }_{\mathrm{cs}}}{r}{v}_{0}exp\left(-\frac{E_{\mathrm{p}}}{kT}\right)sinh\left(\frac{E_{\mathrm{p}}}{kT}\frac{\omega }{{\omega }_{\mathrm{cr}}}\right)\right],$$

(21)

from Equations (5) and (20).

The differential rotation between the normal and the superfluid components dissipates some of the rotational energy of a neutron star and heats up its interior. The heating rate is determined by the difference between the work done by the external braking torque on the neutron star and the change in the rotational energies of the superfluid and normal components as [152,178,179,180]

$$\dot{E}\left(t\right)=N_{\mathrm{ext}}{\Omega }_{\mathrm{c}}\left(t\right)-\frac{d}{dt}\left[\frac{1}{2}{I}_{\mathrm{c}}{\Omega }_{\mathrm{c}}^2\left(t\right)+\frac{1}{2}{I}_{\mathrm{cs}}{\Omega }_{\mathrm{cs}}^2(r,t)\right]=I_{\mathrm{cs}}\left|{\dot{\Omega }}_{\infty }\right|\omega (r,t),$$

(22)

where in the last step we used the fact that in rotational equilibrium all the components are decelerating at the same rate, i.e., ${\dot{\Omega }}_{\mathrm{c}}={\dot{\Omega }}_{\mathrm{cs}}\cong {\dot{\Omega }}_{\infty }$.

4.2. Vortex Creep Model

Vortex pinning becomes favourable if the energy cost per particle inside vortex normal core is smaller when the vortex overlaps with the lattice nucleus than the line stays outside it. Early calculations indicate that in the density range $3\times {10}^{13}$ to $2\times {10}^{14}$ g cm⁻³ vortex pinning is energetically favourable [181]. The dragging of the superconducting protons by the superfluid neutrons in the neutron star core generates a proton supercurrent circulating around each vortex line which makes vortex lines strongly magnetized [169]. This leads effective scattering of electrons from vortices and thus maintains tight coupling of the neutron star liquid core to the solid crust. Hence, the only component inside the neutron star that is weakly coupled to the crust and may be responsible for the long-term relaxation after the glitches is the inner crust superfluid in which vortex lines can get pinned to the lattice nuclei.

From Equation (20) we can define a linearity parameter as

$$\eta \equiv \frac{r|{\dot{\Omega }}_{\infty }|}{4{\Omega }_{\mathrm{cs}}{v}_0}exp\left(\frac{E_{\mathrm{p}}}{kT}\right).$$

(23)

For $\eta \le 1$ crustal superfluid is said to be in the linear regime and since $sinhx\cong x$ we obtain

$${\dot{\Omega }}_{\mathrm{cs}}=-\frac{E_{\mathrm{p}}}{kT}\frac{4{\Omega }_{\mathrm{cs}}}{r{\omega }_{\mathrm{cr}}}{v}_{0}exp\left(-\frac{E_{\mathrm{p}}}{kT}\right)\omega \equiv \frac{\omega }{{\tau }_{\mathrm{l}}},$$

(24)

where the linear creep relaxation time-scale ${\tau }_{\mathrm{lin}}$ is defined as

$${\tau }_{\mathrm{lin}}=\frac{r\varpi }{4{\Omega }_{\mathrm{cs}}{v}_0}exp\left(\frac{E_{\mathrm{p}}}{kT}\right).$$

(25)

In the last step $\varpi$ is shorthand notation for

$$\varpi \equiv \frac{kT}{E_{\mathrm{p}}}{\omega }_{\mathrm{cr}}.$$

(26)

In the linear creep regime the angular velocity lag assumes the form

$${\omega }_{\mathrm{lin}}\left(t\right)=exp\left(-\frac{t}{{\tau }_{\mathrm{lin}}}\right)\left[\omega \left(0\right)-\frac{1}{I_{\mathrm{c}}}{\int }_0^{t}exp\left(\frac{t^{\prime }}{{\tau }_{\mathrm{lin}}}\right)N_{\mathrm{ext}}\left(t^{\prime }\right)d{t}^{\prime }\right],$$

(27)

with $\omega \left(0\right)$ being the initial value of the lag at some time $t=0$, while the response of the observed crustal rotation rate to $\omega \left(0\right)$ is

$${\dot{\Omega }}_{\mathrm{c}}\left(t\right)=\frac{N_{\mathrm{ext}}\left(t\right)}{I_{\mathrm{c}}}+\frac{I_{\mathrm{cs}}}{I{\tau }_{\mathrm{lin}}}exp\left(-\frac{t}{{\tau }_{\mathrm{lin}}}\right)\left[\omega \left(0\right)-\frac{1}{I_{\mathrm{c}}}{\int }_0^{t}exp\left(\frac{t^{\prime }}{{\tau }_{\mathrm{lin}}}\right)N_{\mathrm{ext}}\left(t^{\prime }\right)d{t}^{\prime }\right].$$

(28)

Due to the existence of coupling between the external pulsar braking and internal superfluid torques Equation (28) may give rise to over-relaxation, i.e., the spin-down rate after the exponential recoveries are over is greater than the pre-glitch original value $Q\gtrsim 1$ [177], which accounts for the observed case seen after many magnetar glitches displaying radiative changes [87].

In the case of constant external braking torque $N_{\mathrm{ext}}=I\left|{\dot{\Omega }}_{\infty }\right|$, ${\omega }_{\mathrm{lin}}$ simplifies to the steady state, time independent equilibrium value

$${\left[{\omega }_{\infty }\right]}_{\mathrm{lin}}=\left|{\dot{\Omega }}_{\infty }\right|{\tau }_{\mathrm{lin}},$$

(29)

and the post glitch spin-down rate given by Equation (28) now becomes

$${\dot{\Omega }}_{\mathrm{c}}={\dot{\Omega }}_{\infty }-\frac{I_{\mathrm{cs}}}{I}\frac{\delta \omega }{{\tau }_{\mathrm{lin}}}exp\left(-\frac{t}{{\tau }_{\mathrm{lin}}}\right),$$

(30)

where $\delta \omega$ is the change of the lag at the time of the glitch and is given by

$$\delta \omega =\omega \left(0\right)-{\omega }_{\infty }.$$

(31)

Equation (25) is related to the exponential decay time-scale of the transient post-glitch increase in the spin-down rate of the observed crust. Thus, post-glitch observations of pulsars can be used to constrain microphysical properties of the neutron star crust and superfluid traits [106]. Note that the response of vortex pinning and creep across flux tubes to a spin up glitch also leads to exponential decay [182,183]. The exponential decay timescales of solitary pulsars are found to be in qualitative agreement with the response of neutron star core whereas the case of magnetars implies that direct Urca process mediated fast cooling may operate inside this class of objects [184]. Exponential decay time-scales for magnetars estimated by the model fit observations very nicely if the internal core temperature of these class of neutron stars are below the temperatures predicted by the standard cooling calculations [185,186]. Accelerated cooling maintained by the direct Urca process gives right post-glitch recovery time-scales. Since the observed post-glitch exponential decay time-scales are at most of order ${\tau }_{\mathrm{lin}}\approx 100$ d. [30,69,187], the steady state lag before a glitch is ${\left[{\omega }_{\infty }\right]}_{\mathrm{lin}}\ll {\omega }_{\mathrm{cr}}$ for reasonable threshold values for vortex unpinning ${\omega }_{\mathrm{cr}}\cong {10}^{-2}$–${10}^{-1}$ rad s⁻¹ [152,188] in these linear creep regions. Therefore, linear creep regions cannot be responsible for the glitches by initiating collective vortex unpinning avalanche [106].

At the other extreme, if $\eta$ given by Equation (23) is greater than one, vortex creep is in the nonlinear regime. For the nonlinear creep stage since $sinhx\approx exp\left(x\right)/2$ Equation (20) becomes

$${\dot{\Omega }}_{\mathrm{cs}}=-\frac{2{\Omega }_{\mathrm{cs}}{v}_0}{r}exp\left(\frac{\omega }{\varpi }-1\right),$$

(32)

and the time evolution of the lag is now given by

$$\dot{\omega }=-\frac{I\varpi }{2I_{\mathrm{c}}{\tau }_{\mathrm{lin}}}exp\left(\frac{\omega }{\varpi }\right)-\frac{N_{\mathrm{ext}}\left(t\right)}{I_{\mathrm{c}}}.$$

(33)

The steady state lag in the nonlinear regime is given by

$${\left[{\omega }_{\infty }\right]}_{\mathrm{nl}}={\omega }_{\mathrm{cr}}\left[1-\left(\frac{kT}{E_{\mathrm{p}}}\right)ln\left(\frac{4v_0{\tau }_{\mathrm{sd}}}{r}\right)\right]\simeq {\omega }_{\mathrm{cr}}.$$

(34)

The differential rotation between the pinned superfluid in the nonlinear creep regime and the crustal normal matter results in energy dissipation at a rate

$${\dot{E}}_{\mathrm{diss}}\cong I_{\mathrm{s}}{\omega }_{\mathrm{cr}}\left|\dot{\Omega }\right|.$$

(35)

Equation (35) leads to a effective surface temperature

$$T_{\mathrm{s}}={\left(\frac{{\dot{E}}_{\mathrm{dis}}}{4\pi R^2\sigma }\right)}^{1/4},$$

(36)

where R is the neutron star radius and $\sigma$ is the Stefan-Boltzmann constant. Equation (36) predicts that an old pulsar which radiated away its original heat content would shine in ultraviolet region of the thermal spectrum. The observations of a handful of pulsars confirm this prediction [189,190,191,192,193,194,195,196]. Note also that other internal heating mechanisms, in particular rotochemical heating add constructively for high surface temperature detection in old pulsars [189,197]. In nonlinear creep regions the steady state lag is very close to the critical threshold for vortex unpinning as pulsar ages and temperature drops accordingly. Small statistical fluctuations due to local temperature increase or shift due to a quake in the crust can easily raise the lag above the critical level so that the excess superfluid angular momentum is tapped by the collective unpinning and radially outward motion of vortex lines contained within. Kelvon waves induced on unpinned superfluid vortices and their coupling with the lattice phonons entail angular momentum exchange with ambient normal matter and thereby give rise to a rapid speed up of the observed neutron star crust. The resulting angular momentum transfer process and therefore glitch spin-up time lasts for when these free vortex lines repin to new pinning sites. Therefore, these nonlinear creep regions are agent for pulsar glitches. At the time of a glitch superfluid rotation rate decreases by $\delta {\Omega }_{\mathrm{s}}$, i.e., ${\Omega }_{\mathrm{s}}\to {\Omega }_{\mathrm{s}}-\delta {\Omega }_{\mathrm{s}}$ and crustal rotation rate acquires an increment of amount $\Delta {\Omega }_{\mathrm{c}}$, i.e., ${\Omega }_{\mathrm{c}}\to {\Omega }_{\mathrm{c}}+\Delta {\Omega }_{\mathrm{c}}$. The glitch magnitude $\Delta {\Omega }_{\mathrm{c}}$ is then found from the angular momentum conservation at the time of the glitch event:

$$I_{\mathrm{c}}\Delta {\Omega }_{\mathrm{c}}=I_{\mathrm{s}}\delta {\Omega }_{\mathrm{s}}.$$

(37)

Thus, glitch magnitude is determined by two factors: The number of vortices unpinned at the time of glitch, $\delta N_{\mathrm{V}}=2\pi r^2\delta {\Omega }_{\mathrm{s}}/\kappa$ and the radial distance traversed by vortices before repinning. $I_{\mathrm{s}}$ consists of two parts: The nonlinear creep regions with moment of inertia $I_{\mathrm{A}}$ wherein vortices unpinned and regions with moment of inertia $I_{\mathrm{B}}$ through which vortices do not creep at all but contribute to the angular momentum balance via traverse of lines within them at the time of a glitch. Due to the movement of vortices across a glitch rotational energy of the pulsar changes by,

$$\Delta E_{\mathrm{rot}}=\Delta J{\omega }_{\infty }=I_{\mathrm{s}}{\omega }_{\infty }\Delta {\Omega }_{\mathrm{c}},$$

(38)

which is dissipated in the form of heat. Here $\Delta J$ is the angular momentum transfer across the glitch, ${\omega }_{\infty }$ is the pre-glitch lag between the pinned superfluid and the crust, $I_{\mathrm{s}}$ is the moment of inertia of the pinned superfluid and $\Delta {\Omega }_{\mathrm{c}}$ is the observed spin-up in the crustal rotation rate. The thermal luminosity expected from afterglow of glitches is given by [198,199,200,201]

$$\Delta L=\frac{\Delta E_{\mathrm{rot}}}{{\tau }_{\mathrm{th}}},$$

(39)

where ${\tau }_{\mathrm{th}}\approx {10}^6$ s is the thermal conduction time in the crust.

For constant external braking torque around the time of a glitch the post-glitch spin-down rate becomes

$${\dot{\Omega }}_{\mathrm{c}}\left(t\right)={\dot{\Omega }}_{\infty }-\frac{I_{\mathrm{cs}}}{I_{\mathrm{c}}}{\dot{\Omega }}_{\infty }\left[1-\frac{1}{1+\left[exp\left(\frac{t_0}{{\tau }_{\mathrm{nl}}}\right)-1\right]exp\left(-\frac{t}{{\tau }_{\mathrm{nl}}}\right)}\right],$$

(40)

where we defined nonlinear relaxation time

$${\tau }_{\mathrm{nl}}\equiv \frac{\varpi }{|{\dot{\Omega }}_{\infty }|},$$

(41)

and offset time

$$t_0\equiv \frac{\delta \omega }{|{\dot{\Omega }}_{\infty }|}.$$

(42)

Equation (40) can be integrated explicitly by summing contributions of all non-linear creep regions into account with the assumption of linearly decreasing superfluid angular velocity during a glitch through a superfluid layer with radial extent $\delta r_0$, i.e., $\delta {\Omega }_{\mathrm{s}}(r,0)=(1-r/\delta r_0)\delta {\Omega }_0$ [202]

$$\Delta {\dot{\Omega }}_{\mathrm{c}}\left(t\right)=\frac{I_{\mathrm{A}}}{I_{\mathrm{c}}}{\dot{\nu }}_0\left(1-\frac{1-({\tau }_{\mathrm{nl}}/t_0)ln\left[1+(e^{t_0/{\tau }_{\mathrm{nl}}}-1)e^{-t/{\tau }_{\mathrm{nl}}}\right]}{1-e^{-t/{\tau }_{\mathrm{nl}}}}\right).$$

(43)

Here total moment of inertia of the non-linear creep regions affected from vortex unpinning is indicated by $I_{\mathrm{A}}$. For times ${\tau }_{\mathrm{nl}}\lesssim t\lesssim t_0$ Equation (43) reduces to

$$\frac{\Delta {\dot{\Omega }}_{\mathrm{c}}\left(t\right)}{{\dot{\Omega }}_{\mathrm{c}}}=\frac{I_{\mathrm{A}}}{I_{\mathrm{c}}}\left(1-\frac{t}{t_0}\right).$$

(44)

After initial exponential recoveries are over the vortex creep model parameters are related to the post-glitch pulsar observables with three simple equations [203]:

$$\frac{\Delta {\Omega }_{\mathrm{c}}}{{\Omega }_{\mathrm{c}}}=\left(\frac{I_{\mathrm{A}}}{2I}+\frac{I_{\mathrm{B}}}{I}\right)\frac{\delta {\Omega }_{\mathrm{s}}}{{\Omega }_{\mathrm{c}}},$$

(45)

$$\frac{\Delta {\dot{\Omega }}_{\mathrm{c}}}{{\dot{\Omega }}_{\mathrm{c}}}=\frac{I_{\mathrm{A}}}{I},$$

(46)

$${\ddot{\Omega }}_{\mathrm{c}}=\frac{I_{\mathrm{A}}}{I}\frac{{\dot{\Omega }}_{\mathrm{c}}^2}{\delta {\Omega }_{\mathrm{s}}}.$$

(47)

The pre-factor $1/2$ in Equation (45) accounts for the assumption of linear decrease of the superfluid angular velocity through nonlinear creep regions with moment of inertia $I_{\mathrm{A}}$. Equations (45)–(47) predicts time to the next glitch from the magnitudes of the previous one as

$$t_0=\frac{\delta {\Omega }_{\mathrm{s}}}{|\dot{\Omega }|}=2\times {10}^{-3}\left[{\left(\Delta {\Omega }_{\mathrm{c}}/{\Omega }_{\mathrm{c}}\right)}_{-6}/\left(\beta +1/2\right){\left(\Delta {\dot{\Omega }}_{\mathrm{c}}/{\dot{\Omega }}_{\mathrm{c}}\right)}_{-3}\right]{\tau }_{\mathrm{sd}},$$

(48)

where $\beta \equiv I_{\mathrm{B}}/I_{\mathrm{A}}$. The same set of equations gives an estimate for the interglitch braking index

$$n_{\mathrm{ig}}=\left(\beta +1/2\right){\left(\Delta {\dot{\Omega }}_{\mathrm{c}}/{\dot{\Omega }}_{\mathrm{c}}\right)}_{-3}^2/{\left(\Delta {\Omega }_{\mathrm{c}}/{\Omega }_{\mathrm{c}}\right)}_{-6}.$$

(49)

Equation (49) estimates as large as $n_{\mathrm{ig}}$∼100 which compares well with the observed values [30,50,69].

Vortex creep model has been applied to the post-glitch timing data of many pulsars [101,117,174,202,204,205,206,207,208,209,210]. From the vortex creep model fits to the post-glitch observational data many invaluable information probing into the neutron star internal structure and dynamics can be obtained. Among them, total moment of inertia of the superfluid regions participated to glitches provides the fractional moment of inertia of the crustal superfluid and in turn places constraint on the equation of state of neutron star matter. The superfluid recoupling timescale given by Equation (41) with an information on the temperature for a given pulsar helps to restrict related microphysics at subnuclear densities including superfluid pairing gap. The theoretical estimate of Equation (42) for the time to the next glitch agrees well with the observed interglitch time-scales especially for middle aged pulsars exhibiting Vela-like glitches [101]. Therefore such estimates for other pulsars can be used to plan targeted observations of next glitches for the corresponding sources.

Anti-glitches with reverse sign in the glitch magnitude can also be explained in terms of vortex avalanches under radially inward bias [129,211]. A crustquake as an external agent may drive some of vortices radially inward at the time of glitch and acts as the required inward bias [210]. According to the vortex creep model delayed spin-ups observed in the young Crab pulsar results from the response of the internal superfluid torque on the neutron star crust when vortex lines transported inward after a crustquake [117]. Slow glitches, on the other hand, may be different manifestation of new vortex current distribution in old pulsars after inward movement due to crust breaking quake.

Spin-down rate oscillations observed in the long term rotational evolution of some pulsars [212,213,214,215] as well as post glitch relaxation of the Vela pulsar [216] and PSR B2234+61 [51,217] may result from the response of vortex line oscillations under the competition between the pinning potential and vortex tension [218].

If pulsar glitches are originated from crustal superfluid alone like in the vortex creep model, then the permanent post-glitch increase in the spin-down rate after exponential recoveries are over can be used to place limits on the crustal thickness and equation of state of neutron star matter and in turn pulsar mass [180,219,220,221,222,223].

The tightest limit on the glitch rise time to date is obtained for the observation of 2016 glitch of the Vela pulsar and implies a spin-up of the neutron star less than 12 s [173]. Such short spin-up time-scales result from involvement of vortex lines via the coupling of kelvon modes induced on vortices with lattice phonons when they freed [224,225,226]. The relativistic extension of the mutual friction mechanism between the superfluid and normal matter [227,228], and interpretation of the overshooting of the rotation rate immediately following the glitch as time variable coupling of neutron star core to the crust [174,229,230] provide insights for neutron star internal structure and dynamics.

Vortex pinning to lattice nuclei is at the heart of many glitch models including the vortex creep model which invokes pinned crustal superfluid as the angular momentum reservoir. A vortex will pin on to a nucleus if the interaction between them is attractive while hold on interstices of the lattice if the vortex-nucleus interaction is repulsive [181,231]. As the neutron star crust spins-down under the combined action of the magnetic dipole radiation and pulsar wind emission, the vortex pinning stores angular momentum in some part of the interior superfluid component by fixing angular velocity of the pinned superfluid [100]. The amount of angular momentum maintained in the pinned superfluid which drives a glitch is determined by the strength at which vortex lines pin to lattice nuclei. It is shown that pinning at various strengths may occur depending on the interaction potential between vortex and nucleus, the lattice orientation with respect to the vortex line and the vortex tension [188,232,233,234]. In particular an amorphous crustal structure or a lattice with large impurities would result in strong pinning wherein the unpinning process may start [235,236]. A vortex approaching to lattice falls into interaction potential well provided by neighbouring nuclei, and waves will be excited on the vortex during this process. As these waves damp as a result of the energy imparted to lattice phonons, the final shape of a vortex would be determined by the competition between the gain due to being located at the bottom of the potential well and the cost due to lengthening after bending [224,225,237,238]. There has been substantial controversy in the literature concerning the strength of vortex-nucleus pinning interaction as well as the interaction sign. Quantum calculations of Avogadro et al. [239,240] employ a mean field Hartree-Fock-Bogoliubov formalism approach to vortex-nucleus interaction which results in ∼3 MeV for the pinning energy per intersection. According to results of their calculations the interaction is in repulsive nature throughout the sizeable portion of the inner crust. Pizzochero and his collaborators [241,242,243,244] have used local density approximation in their semi-classical calculations and obtained a repulsive interaction of magnitude ∼1–2 MeV below densities ∼2 $\times {10}^{13}$ g cm⁻³ and an attractive interaction of strength ∼5 MeV at higher densities. An exact density functional theory approach to the same problem has yielded a repulsive interaction with ∼4 MeV per nucleus for the density range (3–7) $\times {10}^{13}$ g cm⁻³ [233]. Jones [232] argues that for an infinite, single component body centered cubic (bcc) lattice vortex line orientation with respect to the different lattice planes would cancel out the pinning force on a vortex largely and concludes that vortex lines are in the co-rotation with the crustal lattice without pinning. However, the neutron star crust likely involves not such symmetric but multicomponent bcc lattice which renders pinning to be effective for specific lattice orientations [245]. Also neutron star crustal lattice is not expected to be an ideal bcc as thermal fluctuations in the era of crustal solidification may lead to formation of monovacancies and lattice defects [235,236]. It is speculated that above the neutron drip point the attractive force between the interstitial unbound neutrons and the lattice makes bcc lattice unstable [246]. However, when the Coulomb interaction between two nuclei is ignored the total attractive interaction is entirely due to exchange of phonons which results in cancellation of destabilisation effects [247]. Recently three dimensional simulations of vortex line-nuclei interaction have been performed to understand the vortex line-nucleus pinning problem and its implications for pulsar glitches [188,248,249]. As the density increases lattice nuclei get closer so that they almost touch each other. The equilibrium shape of nuclei deviates considerably from being spherical as a result of destructive Coulomb and nuclear forces between them, and rod-like, plate-like nuclei appear. Such a sequence of nuclear shapes is collectively dubbed as pasta phase [250]. On the other hand, calculations of vortex pinning to non-spherical nuclei in the pasta phases of deeper layers of the neutron star crust has only premature nature [251,252,253].

In the neutron star crust dissipationless coupling of the dripped superfluid neutrons with lattice nuclei restricts mobility of unbound neutrons and gives rise to a significant decrement for the effectiveness of superfluid neutrons to impart their angular momentum to the crust in glitches [254,255,256,257]. Band theory calculations of [258] show that this entrainment effect may lead to large effective to bare mass ratio for neutrons with average enhancement factor $<m_{\mathrm{n}}^{\ast }/m_{\mathrm{n}}>=5$ throughout the neutron star crust. Therefore, the inferred moments of inertia found from post-glitch spin-down fits should be multiplied by the same enhancement factor. However, if the lattice in the neutron star crust is disordered, then the effects of the entrainment become less important [259]. Uncertainties regarding pairing potential may lead to smaller enhancement factors [247,260,261,262,263]. The pasta phases in the deepest regions of the inner crust with the series of different geometries of the atomic nuclei brings about reduction in the strength of entrainment compared to the spherical nucleus case [264]. Also some parts of the neutron star outer core may be involved in the glitches [182,222], which provided new way out for crustal superfluid based models.

So far we have investigated the effects of vortex pinning and creep for laminar superfluid flow. In the literature the effects of superfluid turbulence on pulsar glitches have been studied in some extent [265,266,267,268,269,270,271]. Instabilities associated with the superfluid hydrodynamical flow were proposed to play an important role in initiating collective vortex unpinning cascade [272,273,274,275,276,277]. It is shown that the creep and pinning of vortices across flux tubes in the outer core stabilises the superfluid flow and inhibits large scale turbulence [278].

4.3. Crustquake Model

Crustquakes are failures of the solid neutron star crust, resulting from spin-down assisted continual reduction in the centrifugal force which deforms and stresses the crust until it breaks. Crustquakes both spin up the neutron star and dissipate mechanical energy inside the crust, thereby producing heat.

The outermost layer of a neutron star solidifies in an early epoch of the star’s life when it was spinning comparatively faster and thus had a relatively oblate shape. As a neutron star spins down its fluid interior adjusts its shape to instantaneous new spin rate while its rigid crust resits such changes. Then, stresses build up in the crust which will be released in discrete quake events when the yield point is exceeded [279,280]. The shape change after the crack reduces the moment of inertia of the crust with reference oblateness accompanying a jump. The conservation of angular momentum and torque equilibrium acting on the neutron star at the time of the glitch implies

$$\frac{\Delta \Omega }{\Omega }=\frac{\Delta \dot{\Omega }}{\dot{\Omega }}=-\frac{\Delta I}{I}\approx 2\frac{\delta R}{R},$$

(50)

where $\delta R/R$ denotes the shrinkage in the solid component of the neutron star and is approximately given by [279]

$$\frac{\delta R}{R}\approx \frac{95}{7}\frac{\mu R}{GM\rho }\Delta \varphi .$$

(51)

Here, R is the radius and M is the mass of the neutron star, $\rho$ is the average stellar density, $\Delta \varphi$ is the strain relieved in the quake and $\mu \approx (\Delta /R){\left(Ze\right)}^2{n}^{4/3}$ is the shear modulus with $\Delta /R\sim 0.1$ being the crustal thickness to neutron star radius ratio, $Ze$ being the nuclear charge and n being the number density of nuclei. Recent theoretical and computational calculations show that the more accurate shear modulus value for rigid neutron star crust does not differ significantly from the above quoted order of magnitude estimate [281,282,283,284]. In Ref. [285] it has been argued that since the only dimensionless parameter related to the problem is the fine structure constant $\alpha =e^2/\hslash c$, the critical strain angle for crust cracking should be of the order ${\theta }_{\mathrm{cr}}=Ze\alpha$. The recent molecular dynamical simulations, which take realistic shear motion of the neutron star crust into account, yield ${\theta }_{\mathrm{cr}}\simeq 0.04-0.1$ [286,287,288], and is similar in order of magnitude to the early theoretical estimate. Starquake model has a definite prediction for the repetition of glitches since the time between successive quakes is that required to build up the right amount of stress from its current value to the critical strain.

Baym and Pines [280] established the formulation of starquake model for the solid neutron stars. In their model neutron star involves quadrupolar stellar distortion giving rise to mechanical stress. It is determined by a single time-dependent parameter, the crustal oblateness which can be only partially relieved during a quake. For quadrupolar deformation the rotational energy of a neutron star can be written as

$$E=\frac{L^2}{2I_0}+A{ϵ}^2+B{(ϵ-{ϵ}_0)}^2,$$

(52)

where the oblateness $ϵ$ is defined by the relation $I=I_0(1+ϵ)$ with $I_0$ being the moment of inertia of the nonrotating spherical star and ${ϵ}_0$ is the reference oblateness. In Equation (52) A and B quantifies the gravitational potential energy and the elastic energy stored in the neutron star crust, respectively. Minimizing Equation (52) with keeping stellar angular momentum L and ${ϵ}_0$ fixed gives

$$ϵ=\frac{I{\Omega }^2}{4(A+B)}+\frac{B}{A+B}{ϵ}_0.$$

(53)

A quake on the crust occurs whenever the mean stress in the crust

$$\sigma =\frac{\partial E}{\partial ϵ}\frac{1}{V}=\mu ({ϵ}_0-ϵ),$$

(54)

with V being the crustal volume exceeds the critical value $\mu {\theta }_{\mathrm{cr}}$. In a quake both $ϵ$ and ${ϵ}_0$ reduces according to

$$\Delta ϵ=\frac{B}{A+B}\Delta {ϵ}_0,$$

(55)

where $\Delta ϵ=\Delta I/I=-\Delta \Omega /\Omega$ is the observed glitch magnitude. After the quake stress will accumulate again once more and the next quake will take place after a time

$$t_{\mathrm{q}}\approx \frac{\Delta {ϵ}_0-\Delta ϵ}{\dot{ϵ}}\approx \frac{A}{B}\frac{\Delta ϵ}{2ϵ}\frac{\Omega }{\dot{\Omega }},$$

(56)

with the time evolution of the oblateness is given by

$$\dot{ϵ}\approx \frac{\partial I}{\partial ϵ}\frac{\Omega \dot{\Omega }}{4A}.$$

(57)

There are also other means for sinks of stresses in the neutron star crust, most prominent being the plastic flow [289,290,291]. If the response of the neutron star crust to the applied stress is not in the form of plastic flow, the time to the next glitch is proportional to the amount of the stress relieved in the previous quake and is predictable for a given neutron star model [290]. The time interval between successive quakes is proportional to $\Delta {ϵ}_0$ and set by critical strain angle at which solid crust exceeds the yield point.

The small Crab-like glitches can be accounted for stiff neutron star equation of state while larger Vela-like glitches cannot be explained within the framework of the starquake model as $(A/B)(\Delta ϵ/ϵ)$ becomes larger than the typical interval between glitches [131,292].

Since the rigidity parameter $b=B/A$ is so small, change of the ellipticity associated with a quake is tiny [293,294]. The actual value of the rigidity parameter b is rather uncertain mainly reflected in computation of B. Franco et al. [295] considered a uniform density solid crust having constant shear modulus afloat on a fluid core which led them to obtain 5 times smaller b compared to Baym & Pines model [280]. Cutler et al. [293] obtained even smaller (50 times) rigidity parameter by taking compressible neutron star matter and equation of state dependent shear modulus into account but still working in Newtonian gravity. Since the glitch induced increase in the spin-down rate is proportional to b the correct evaluation of structural parameters is vital in crustquake models which is far from being complete [296,297,298].

In Baym & Pines model the coefficients A and B were evaluated for a completely solid, incompressible constant density neutron star with uniform shear modulus in Newtonian gravity. However, the deformation of a realistic neutron star with different core and crustal densities is qualitatively dissimilar from the case of equal density incompressible star. In the following years the Baym & Pines model is extended to include the effects of general relativity [299], magnetic field [295,300], compressible neutron star model with density dependent shear modulus [293,301].

Molecular dynamical simulations suggest large values of ${\theta }_{\mathrm{cr}}$ implying that the secular spin-down may only produce crustquakes just a few times in the whole life of a neutron star [286,287,288].

Crustquakes can also play the role of trigger for vortex unpinning events by forcing vortices to leave their equilibrium configuration as their ends are anchored to the moving broken crustal platelet [117,210,302]. Crustquakes turn some metastable mechanical energy into heat by dissipation. Such heating may mobilize a very large number of vortices in the inner crust, thus serving as trigger for collective vortex unpinning events [108]. The energy released in a crustquake induced glitch is always much larger than the superfluid originated glitch of the same size in which angular momentum transfer from a faster rotating superfluid in the inner crust to the normal matter dissipates rotational energy inside the star [303].

In the presence of a global stellar magnetic field the crust fractures asymmetrically and fault planes moves in a way such that the angle between magnetic dipole and rotational axes and torque acting on the neutron star both increase [304,305,306]. This is proposed as a viable explanation of persistent step increases observed after the Crab glitches [45,75,77,307].

Progressive failures in the solid crust suffice to account for at least the small glitches observed in pulsars of all ages. Giliberti et al. [308,309] computed the deformation and the strain due to reduction in the centrifugal force in time by modelling the neutron star a solid crust floating above a constant density core. This is a very simplifying assumption since stable density stratification is expected to be present in realistic neutron stars [310]. Most recent works on stresses in the crust rely on quantifying the displacement field with centrifugal force changes as a part of the rotational evolution [291,296,308,309].

Elastically deformed crust under stresses motivates searches for detection of gravitational wave emission from neutron stars [311]. It was shown that stresses in the neutron star crust may sustain an ellipticity much larger than implied by observational upper limits for continuous gravitational wave [296].

The crustquake model is also envisaged in terms of solid quark stars [312,313,314,315]. Quark stars are bound by the strong nuclear force rather than the gravitational force relevant for neutron stars. For pure quark matter the shear modulus is larger than neutron star matter so that one can obtain Vela like large glitches in typical repetition timescale of a few years without exceeding observational bounds put by surface temperature measurements. There are two types of glitches within the quake model scenario for quark stars: Type I glitches are bulk invariable events resulted from solely a change in stellar ellipticity and stellar volume does not alter so that the amount of energy release is insignificant. While Type II glitches are bulk variable events involving stellar shape change induced probably by accretion and accompanied by large amount of energy. Type I glitches are invoked to explain standard glitches whereas Type II glitches may account for magnetar glitches with radiative changes. In the quark star model slow glitches are explained as conversion of a collapsed superficial layer with a melted layer beneath the surface. Cooling of this viscous fluid layer depletes the accumulated stress and releases energy to effectuate slow rise of pulsar spin.

4.4. Vortex Line-Flux Tube Interaction Model

About a year or so after the birth of a neutron star, its core temperature is expected to fall below the transition value and the interior proton plasma becomes a superconductor [220]. Since the flux expulsion timescale accompanying this transition is substantially long the magnetic field inside a neutron star core organises into tiny flux tubes of cross section $\Lambda$∼${10}^{-11}$ cm each carrying a field of $B_{\mathrm{c}}$∼${10}^{15}$ G [316]. The number density of flux tubes

$$n_{\Phi }=\frac{B}{{\Phi }_0}=5\times {10}^{18}\left(\frac{B}{{10}^{12}\mathrm{G}}\right){\mathrm{cm}}^{-2},$$

(58)

with ${\Phi }_0=hc/2e\cong 2\times {10}^{-7}$ G cm² being the unit of magnetic flux in each tube, far exceeds the number density of vortex lines given by

$$n_{\mathrm{V}}=\frac{2{\Omega }_{\mathrm{c}}}{\kappa }={10}^5\left(\frac{{\Omega }_{\mathrm{c}}}{100{\mathrm{rad}\mathrm{s}}^{-1}}\right){\mathrm{cm}}^{-2},$$

(59)

in a typical neutron star. Therefore, intimately strong interactions between these two types of topological defects in superfluid-superconducting interior of a neutron star are expected on qualitative grounds. Another important feature of proton superconductivity inside neutron stars is that following the transition with the nucleation of magnetic field into flux tubes tensile force and in turn the stress associated with the existence of magnetic field increases enormously, $B^2/8\pi \to B{B}_{\mathrm{c}}/4\pi$ while the repulsion between the flux tubes diminishes greatly since in equilibrium configuration they are separated well apart [317,318]. When a vortex line and flux tubes are within a distance of $\sim \Lambda$ these two types of structures would strongly interact with each other and the velocity difference between them is the main factor leading to distinct dynamical consequences for pulsar glitches. Radially expanding vortices in a spinning down neutron star either cut through or carry along with them flux tube array for each encounter.

Ruderman and his collaborators [155,318] considered the effect of vortex line-flux tube pinning in the core of a neuron star on the crustal failure and resulting evolution for pulsar glitches. Several years ago Sauls [164] and Srinivasan et al. [319] speculated that on microscopic-scale density fluctuation related interactions between neutron superfluid vortex lines and proton superconductor magnetic flux tubes within a neutron star’s core could give rise to a strong coupling between the star’s magnetic field and its spin history and thereby may play an important role for neutron star magnetic field evolution and spin glitches, respectively. This seminal idea was later strengthened by including the effect of a neutron vortex line’s own strong magnetic field generated due to proton entrainment mass currents around them [320,321]. Since the footpoints of flux tubes are anchored in the crust base, the associated stress on the rigid crust from moving flux tubes would grow as a result of the migration of the vortex line-flux tube conglomerate during secular spin-down. Then, either the crustal strain will ultimately exceed the yield point, or the core’s vortex lines will be forced to cut through the flux tubes. The stress associated with the flux tube squeezing the crust base leads to failure of the solid and break it when the following condition is met [322]:

$$\frac{B{B}_{\mathrm{c}}}{8\pi }{L}^2\gtrsim \left(\mu {\theta }_{\mathrm{cr}}\Delta \right)L,$$

(60)

where $\Delta$ is the crustal thickness and L is the lengthscale on which crust cracks. The main uncertainties in the solid neutron star crustal properties are the maximum sustainable shear strain and the amount of sudden strain-relaxation if yield point is exceeded. An initially hot neutron star rapidly cools down below the temperature via neutrino emission and crust-solidification begins just 10 s after the birth. A key ingredient of the model is the relative velocity of a flux tube with respect to the that of a vortex line. The comparative velocity $v_{\mathrm{c}}$ at which expanding vortex line array can push the flux tubes arrangement through core electron-proton plasma with the maximum pace would be [155]

$$v_{\mathrm{c}}=\beta \left(\frac{\Omega }{100{\mathrm{rad}\mathrm{s}}^{-1}}\right)\left(\frac{{10}^{12}G}{B}\right){10}^{-6}{\mathrm{cm}\mathrm{s}}^{-1}.$$

(61)

Here the parameter $\beta$ depends heavily on the various properties of the neutron star core:

$$\beta =0.4ln\left(\frac{\Lambda }{\xi }\right)\left(\frac{B_{\mathrm{V}}}{{10}^{15}G}\right)\left(\frac{B_{\mathrm{c}}}{{10}^{15}G}\right)\left(\frac{60\mathrm{MeV}}{E_{{\mathrm{F}}_{\mathrm{e}}}}\right)\left(\frac{{10}^{36}{\mathrm{cm}}^{-3}}{n_{\mathrm{e}}}\right),$$

(62)

where $\Lambda$ is the London penetration depth, i.e., the lengthscale over which magnetic flux decays for a flux tube, $\xi$ is the core radius of flux tube, $B_{\mathrm{V}}$ is the magnetic field attached to vortices, $B_{\mathrm{c}}$ is the field strength of each flux tube, $E_{{\mathrm{F}}_{\mathrm{e}}}$ is the Fermi energy of relativistic degenerate electrons and $n_{\mathrm{e}}$ is the number density of electrons. Note, however, that the form of the driving force used in derivation of Equation (61) which leads to magnetic field diffusion in the neutron star core has been questioned in the literature [323,324,325,326].

The vortex line-flux tube pinning model predicts two distinct families of glitches for pulsars. The diamagnetic screening currents set after superconducting transition in the neutron star core allows flux tubes to move independently through electron-proton plasma which leads bunching for flux tubes when they are pushed by vortex lines [327]. This is the underlying reason why the flux tubes move with the vortex lines rather than remaining immobile inside neutron star core and being cut through or carry along by the radially migrating vortices during spin-down of the neutron star. In young neutron stars since the radially outward vortex velocity is high crust cracking should also give rise to a permanent offset in the perpendicular dipolar component of global magnetic field configuration with respect to the spin axis. This would lead to spin-down rate change, part of which does not heal back. In younger and thus more spinning down neutron stars like the Crab pulsar the horizontal motion of broken platelet relieves strain of amount $\Delta \theta$∼$s/R\ll {\theta }_{\max }$∼${10}^{-1}$ with s being the horizontal displacement of broken platelet. This sudden reduction in the shear stress on the crust by the flux-tubes just beneath it has a magnitude $(B{B}_{\mathrm{c}}/8\pi )(s/R)$. Thus, the glitch magnitudes in spin and spin-down rates are related by [328]

$$\frac{\Delta \Omega }{\Omega }\sim \frac{B{B}_{\mathrm{c}}}{8\pi \rho R^2{\Omega }^2}(s/R)\sim {10}^{-4}\frac{\Delta \dot{\Omega }}{\dot{\Omega }}.$$

(63)

As a pulsar ages and becomes middle aged (like the Vela pulsar) a different kind of glitch family emerges. If there were not very dense array flux tubes around vortex lines, the outward moving vortices would smoothly shorten as they approach to the neutron star equator and then disappear. However, the very existence of flux tubes prevents this to occur and vortices pile up inside an annulus neighbouring to stellar equator and thereby excess angular momentum is deposited there. Highly conductive crust resists the entry of flux tubes towards the crust-core interface. When the Magnus force on these accumulated vortex lines overwhelms the repulsive force between closely packed flux tubes a sudden transfer of angular momentum from vortices to crust takes place making a spin up glitch. As a pulsar matures, the volume occupied by comoving vortex line-flux tube network increases and the glitch size attains much larger magnitudes given by [328]

$$\frac{\Delta \Omega }{\Omega }\sim \left(\frac{5V_{\mathrm{A}}}{2V_{\mathrm{NS}}}\right)\left(\frac{{\tau }_{\mathrm{g}}}{{\tau }_{\mathrm{sd}}}\right),$$

(64)

where $V_{\mathrm{A}}$ is the volume of the annulus of the vortex line-flux tube interaction region and $V_{\mathrm{NS}}$ is the volume of the neutron star. The repetition time between successive Vela like giant glitches is given by [328]

$${\tau }_{\mathrm{g}}\sim 3\times {10}^{-4}\left(\frac{\Omega }{\dot{\Omega }}\right){|n-3|}^{-1}.$$

(65)

Equation (65) predicts ${\tau }_{\mathrm{g}}=0.7$ year for ${\tau }_{\mathrm{sd}}=5$ kyr old “Big Glitcher” PSR J0537–6910 which has interglitch braking index $n_{\mathrm{ig}}\sim 7$ and observed glitch interval ${\tau }_{\mathrm{ig}}=0.4$ year [82,83,84,329,330,331].

Flux tube-vortex line pinning model predicts that after each glitch the amount of the spin-down rate that does not heal is given by [155]

$$\frac{\Delta {\dot{\Omega }}_{\mathrm{p}}}{\dot{\Omega }}=\frac{(3-n)}{2}\frac{{\tau }_{\mathrm{g}}}{{\tau }_{\mathrm{sd}}},$$

(66)

where ${\tau }_{\mathrm{sd}}=\Omega /\left(2\right|\dot{\Omega }\left|\right)$ is the pulsar spin-down (characteristic) time-scale.

Delayed spin-up observed in several glitches of the Crab pulsar can be interpreted in terms of the vortex line-flux tube pinning model as the process of setting of vortex lines on new equilibrium position responding to unbalanced force on vortices due to retarding drag by flux tubes immediately after crust breaking.

According to this model, the glitch variety observed in individual and different pulsars may be related to the complexities associated with breaking properties of the stratified neutron star crust [322]. The various crustal layers with specified nuclei would have different lattice orientations with respect to the crystal symmetry axis. The vortex pinning strengths would be also altered by the crustal depth and the temperature evolution through the crust which further complicates the place where a lattice crack starts. All these effects add constructively to the diversity of observed glitches across the pulsar population.

In all comparisons of glitch observations with the vortex line-flux tube interaction model-based predictions no apparent disagreement have been found. The application of the vortex line-flux tube pinning model to the post-glitch relaxation observations is discussed in Refs. [332,333].

Both the frictional drag of flux tubes within electron-proton plasma and the inelastic cut-through or pinning with vortex lines generate heat which can affect the thermal evolution of a neutron star [184,334]. The resulting heat production rate places stringent limits on the vigour of the vortex line-flux tube pinning model [155].

Much of the important model ingredients like how the vortex lines move together with flux tubes, under what circumstances they cut through them, the associated drag forces acting on moving flux tubes and how the crust breaks under flux tube stress, remain to be investigated at more quantitative level.

5. Radiative Changes Associated with Glitches

Pulsar glitches are generally accepted to be caused by the interior dynamics of a neutron star [98,100,335]. If radiative changes (including long-lived flux enhancements, short bursts, and pulse profile changes) are associated with glitches, a link between the interior and magnetospheric state change would be naturally established. Based on available observational studies, concurrent variations in emission and the rotation of pulsars are a rather rare phenomena.

Table 2. Detailed parameters for pulsars exhibited correlation between glitches and emission changes.

Pulsar Name (PSR)	P (s)	${\mathsf{\tau }}_{\mathit{c}}$ (kyr)	${\mathit{B}}_{\mathit{s}}$ (${10}^{12}\mathbf{G}$)	$\dot{\mathit{E}}$ (${10}^{32}\mathbf{erg}/\mathbf{s}$)	$\Delta \mathsf{\nu }/\mathsf{\nu }$ (${10}^{-9}$)	$\Delta \dot{\mathsf{\nu }}/\dot{\mathsf{\nu }}$ (${10}^{-3}$)	Profile	Ref.
J0738–4042	0.3749	4320	0.727	10	0.36(4)	3(1)	$W_{50}$(↓) ${}^a$	[53]
J0742–2822	0.1667	157	1.69	1400	102.73(11)	2.1(5)	${}^b$	[336,337]
J1119–6127	0.4079	1.61	41	23,000	9400(300)	580(14)	${}^c$	[338]
					5740(80)	79(25)	X(↑)	[339]
B1822–09	0.7690	233	6.42	45	4.08(2)	0.08(1)	${}^d$	[63]
					7.2(1)	1.65(7)	${}^d$	[63]
J2021+4026	0.2653	76.9	3.85	1200	<100(–)	56(9)	$\gamma$(↓)	[340,341]
B2021+51	0.5291	2740	1.29	8.2	0.373(5)	–0.24(3)	$W_{10}$(↓) ${}^a$	[62]
B2035+36	0.6187	2180	1.69	7.5	7.7(8)	67(8)	$W_{50}$(↓) ${}^a$	[44]

^a: W₁₀ and W₅₀ are the full widths of the pulse profile at 10% and 50% of the peak pulse amplitude, respectively. ^b: represents that the appearance of additional pulse components is closely correlated with an unusual glitch. ^c: indicates the correlation between the ratio of the two-components in the profile and $\Delta \dot{\mathsf{\nu }}$, which rapidly increases after the glitch. ^d: implies the variations of the integrated mean pulse profiles in both the radio-bright (B-mode) and the radio-quiet (Q-mode) modes.

lists the detailed parameters of each observed correlation between glitch and emission change. It is seen that just a handful of normal pulsars display radiative changes that are concurrent with the glitch events. The emission change-rotation correlation is found to occur in both young and relatively old pulsars. Glitches in young mode-switching pulsars have typically large $\Delta \nu /\nu$ sizes, which follows the statistical trend [29].

Recently, Zhou et al. [53] reported that a micro-glitch ($\Delta \nu /\nu$∼$0.36\left(4\right)\times {10}^{-9}$) occurred in the middle-aged pulsar PSR J0738–4042 is coincided with variations in the power of the average pulse profile components. As demonstrated in Figure 8, the average pulse profile after the glitch is slightly wider than the pre-glitch one, with a reduction in the leading component and an enhancement of the middle component. The post-glitch $\dot{\nu }$ shows a steady increase during this correlation, and similar post-glitch features have been seen in a very large glitch ($\Delta \nu /\nu$∼33,250 $\times {10}^{-9}$) in PSR J1718–3718 [342] and a delayed spin-up event in the Crab pulsar [112]. The very young pulsar, PSR J1119–6127, is another intriguing example of mode-changing radio pulsars. A large glitch activity in PSR J1119–6127 with two exponential recoveries is linked to the appearance of additional pulse components that have an intermittent and RRAT-like behaviour [338]. Shortly afterward, Archibald et al. [339] observed another large glitch to be accompanied by X-ray bursts and X-ray pulsations in this pulsar. The outburst timing anomalies seen in PSR J1119–6127 are similar to those in magnetars, providing evidence for a connection between the high magnetic field radio pulsar and magnetar populations. The potential future for high-magnetic-field, rotation-powered pulsars is to turn into magnetars [343]. It is important to make timing observations more frequently for high-B young glitching pulsars, such as PSRs J1718–3718 [344] and J1734–3333 [345], as they might display magnetar-like mode-switching behaviours in coming years. The 2016 Vela glitch was accompanied by polarization changes and emergence of a null pulse [92]. It is proposed that changes in the state of magnetospheric emission reflects Alfven mode conversion close to the light cylinder triggered by energy release aftermath of a crustquake [346,347].

The mechanism behind these state-switching correlations is unclear. Despite this, the scenario that the glitch may change the magnetic field structure and hence the inclination angle $\alpha$ [44,210,348], is gradually becoming more accepted. The outcome of this event is a change in the pulsar emission profiles and spin-down torque. In the cases of PSRs B1822–09 and B2021+51, Liu et al. [62,63] proposed that the flux tube, which moves independently in the emission zone, changes position during a glitch, and as a result the pulse profile varies. Zhou et al. [53] performed an in-depth analysis for the simultaneous variations in profile shape and spin-down rate $|\dot{\nu }|$ in PSR J0738–4042, based on a combined model of crustquake induced platelet movement and vortex creep response. They evaluated that the glitch-triggered emission variations may be a result of the crustquake-induced change in the inclination angle and the movement of vortex lines across the glitch [53].

6. Statistics of Glitches

After they have been observed in sufficient numbers of pulsars since their first discovery in 1969, there have been numerous statistical studies on the magnitudes, repetition times and relaxation behaviours of the pulsar glitches [29,32,69,101,349,350,351,352,353,354]. The total number of glitches occurred and glitching pulsars over time is shown in Figure 9. With increased detection of new pulsars and follow up observations of more pulsars by many radio telescopes, a steady big increase in the number of pulsar glitches has been witnessed over the past 50 years. A curve fit is applied to find the cumulative number of the glitch discovery between MJDs 50,000 and 58,000 in Figure 9, leading to a quadratic function $N\propto t^2$. Several statistical studies in particular datasets have been conducted to reveal the glitch populations and activities, even though the existing datasets suffer from incompleteness due to various reasons.

6.1. Glitches Sizes Distributions & Waiting Times

It was pointed out by Espinoza et al. [32] that both the distributions of absolute and fractional frequency jumps, $\Delta \nu$ and $\Delta \nu /\nu$, of all observed glitches show bimodal nature. This behaviour had also been indicated from recent studies using a larger sample [29,30,69,355,356]. Basu et al. [29] used the best-fit two-component Gaussian Mixture Model to label the $\Delta \nu$ distribution using 543 glitches in 178 pulsars. The first, wide Gaussian component with small magnitudes centres at $\Delta \nu$∼$0.032\mathsf{\mu }\mathrm{Hz}$ and the second, narrow Gaussian component with larger magnitudes peaks at $\Delta \nu$∼$18\mathsf{\mu }\mathrm{Hz}$. With a total sample of 700 glitches, Arumugam and Desai [33] applied Extreme Deconvolution (XDGMM) method, which is a generalization of Gaussian Mixture Model, to classify the relative glitch amplitudes $\Delta \nu /\nu$ into two classes. Two pulsar glitch components of the bimodal distribution are ascertained, with the mean values $5.9\times {10}^{-9}$ and $1.3\times {10}^{-6}$, respectively (Figure 10). Moreover, there is a pronounced dip at ${10}^{-7}$ presented in the $\Delta \nu /\nu$ distribution, suggesting that the bimodality of the glitch sizes distribution may be produced by two mechanisms by which the glitch events occur [69,355]. As previously described in Section 4, largest glitches are caused by the sudden transfer of angular momentum from the faster rotating interior superfluid to the solid crust, whereas small glitches may be result from starquakes due to relaxation of the neutron-star crustal oblateness to the current equilibrium shape. Celora et al. [357] proposed the non-linear mutual friction force as a cause of the exchange of angular momentum between the neutron superfluid and the observable normal component in glitches, and this effect is able to explain the observed bimodal distribution. An alternative explanation is that two different classes of pulsars may be the underlying reason of the bimodal behaviour [358]. As for spin-down rate changes $\Delta \dot{\nu }$, negative values are seen in the majority of glitches (that means after each glitch the spin-down rate becomes more negative since its absolute magnitude $|{\dot{\nu }}_0+\Delta \dot{\nu }|$ is greater), and there is usually a greater change in $\Delta \dot{\nu }$ for larger glitches [29].

For most pulsars glitches take place at irregular intervals with random distribution without a preferred distribution of magnitudes. The exceptions are the Vela pulsar and PSR J0537–6910 which have preferred large sizes and regular repetition times obeying Gaussian distribution [359]. At the other extreme, the Crab pulsar glitches display temporal clustering at times [307,360]. The random distribution of glitch sizes and waiting times until the next event can be understood within the framework of the vortex creep model. As can be seen from Equations (45) and (48) both of the observed parameters depend on the two ingredients: number of vortices unpinned at the time of the glitch and the distance covered by them through inner crust which in turn determined by vortex scattering from pinning traps [174,211]. Shallow heating of crust via a quake may sustain different initial conditions for such scattering events. Even though inter-glitch time-scales and glitch sizes in a given pulsar differ significantly, the process may be regarded as scale invariant in terms of the sandpile effect. Radio pulsar glitches thus can be used to probe the far-from-equilibrium processes involving stress accumulation and relaxation in neutron star interiors [361,362]. For productive pulsars exhibiting frequent large glitches the distribution of waiting times is found to be fitted by exponentials while their glitch size distribution obeys power law or log-normal functions [359,363,364,365,366,367,368,369].

6.2. Glitch Activity

The number of glitches per year $\dot{N}$ is a way to simply quantify how often glitches occur for known glitching pulsars [29,32,69]. According to Lyne et al. [350] and Espinoza et al. [32], the observed $\dot{N}$ is significantly correlated with characteristic age ${\tau }_c$ and spin-down rate $|\dot{\nu }|$. Basu et al. [29] found that the maximum and minimum glitch rates for the 134 pulsars observed by the Jodrell Bank Observatory for a long time, were ${\dot{N}}_{\max }=1.07{\mathrm{year}}^{-1}$ and ${\dot{N}}_{\min }=0.02{\mathrm{year}}^{-1}$, respectively. The glitch activity trends in previous work were also confirmed by Basu et al. with a larger sample of glitching pulsars. $\dot{N}$ decreasing for pulsars with larger ${\tau }_c$ implies that young pulsars exhibit glitches more often than old pulsars. A higher $\dot{N}$ occurring alongside a greater $|\dot{\nu }|$ reinforces the notion that glitches are driven by the spin-down. However, this analysis does not take into account the size of the glitches, only the number; and the incomplete observations of glitches could skew these results. Thus, the integrated glitch activity is defined to objectively determine the cumulative effect of the glitches on the pulsar’s rotation [29,350,355]:

$${\dot{\nu }}_{\mathrm{g}}=\frac{{\sum }_i\Delta {\nu }_i}{T},$$

(67)

where $\Delta {\nu }_i$ is the spin frequency increment at the i-th glitch, and T is the duration over which the pulsar has been searched for glitches. It is intriguing how the integrated glitch activity correlates with key pulsar quantities, such as spin-down rate $|\dot{\nu }|$, spin-down age ${\tau }_c$, energy loss rate $\dot{E}$, and magnetic field B. To carry out a robust estimation for these correlations, the average glitch activity ${\overline{\dot{\nu }}}_{\mathrm{g}}$ for groups of pulsars with similar properties is examined [29,355]. The specific results are: the linear relation between ${\overline{\dot{\nu }}}_{\mathrm{g}}$ and $|\dot{\nu }|$ in the range $-14<log|\dot{\nu }|<-10.5$ is $0.018\left(3\right)$, implying that $|\dot{\nu }|$ undergoes a reversal of about 1.8 percent through glitches; ${\overline{\dot{\nu }}}_{\mathrm{g}}$ decreases with ${\tau }_c$; $\dot{E}$ are correlated positively with ${\overline{\dot{\nu }}}_{\mathrm{g}}$; pulsars with lower B also have lower ${\overline{\dot{\nu }}}_{\mathrm{g}}$, and an increase in ${\overline{\dot{\nu }}}_{\mathrm{g}}$ is seen in neutron stars with the strongest B (magnetars) [29,32,350,355].

The accumulated effect of spin-up glitches on the rotational evolution for the long-term observations of PSR J0537–6910, the Vela pulsar (PSR J0835–4510), PSR J1341–6220 and PSR J1740-3015 are shown in Figure 11. The step-like increases correspond to glitches and imply that spin-up contributions are not totally reversed by the ongoing secular slow-down due to electromagnetic dipolar radiation and pulsar wind. In terms of the vortex creep model [152,174,370] this observation corresponds to the fact that at the time of a glitch vortex lines move through vortex traps inside which vortices do not creep at all. Since inside these traps continuous angular momentum impart between glitches does not occur via creep process, the associated part of spin-up increase as a result of discrete angular momentum transfer by glitches persists for interglitch intervals. The regularity of the glitches, i.e., almost constant magnitude and time interval between the events, in PSR J0537–6910 and the Vela entails that glitch activity described by Equation (67) is approximately given by the slopes of lines in the graphs of Figure 11 for these pulsars. On the other hand, large coverage of glitch sizes and unequal repetition timescales for PSRs J1341–6220 and J1740–3015 suggest that different kinds of glitch mechanisms come into play at different times or a combination of two mechanisms appear with variable efficiency. If one assumes that some stellar component with moment of inertia $I_{\mathrm{gl}}$ stores angular momentum by average spin-down rate $\langle \left|\dot{\nu }\right|\rangle$ during long-run of observations containing many glitches, then the rate of angular momentum transferred to the crustal component with moment of inertia $I_{\mathrm{c}}$ gives an lower limit on the fractional moment of inertia of the neutron star component participated in glitches as [254,255,371]

$$\frac{I_{\mathrm{gl}}}{I_{\mathrm{c}}}\gtrsim \frac{{\dot{\nu }}_{\mathrm{g}}}{\langle \left|\dot{\nu }\right|\rangle }.$$

(68)

Equation (68) is applied to the entire pulsar glitch population [220,221,353,372,373] and gives the result of a few percent implying only crustal superfluid and probably some part of the outer core region of a neutron star are involved in glitches when crustal entrainment effect has been properly taken into consideration [182,183,222].

The importance of continuously searching for new pulsar glitches through timing observations is self-evident. With new surveys being conducted by existing and under-construction radio telescopes around the world, the number of known pulsars is expected to increase steadily. It’s definitely worth mentioning that the world’s largest single-dish radio telescope, FAST, has discovered more than 600 new pulsars since it was built in 2016. In addition, monitoring programs to capture a glitch “live” and improved glitch-finding algorithms are being developed. Taking all these together, in coming years many production of increasing numbers of new glitches discoveries will be seen.