###### Abstract

$\sim 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.

###### keywords:

neutron stars; pulsars; glitches

xx 1 5 Pulsar Glitches: A Review Shiqi Zhou 1,5,6, Erbil Gügercinoğlu2*,
Jianping Yuan 3,7, Mingyu Ge 4 and Cong Yu 1,5,6 Shiqi Zhou, Erbil
Gügercinoğlu, Jianping Yuan, Mingyu Ge and Cong Yu Email:
<EMAIL_ADDRESS>(E. G.)

##  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 Hewish (1975). 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 Hartnett and Luiten
(2011). 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 Jenet et al. (2009), has achieved
unprecedented development. Two best known applications of this project are the
detection of ultra-low frequency ($\sim 10^{-9}-10^{-8}\ \rm Hz$)
gravitational waves Hobbs et al. (2010) and establishment of a pulsar-based
timescale Hobbs et al. (2020). Besides, binary millisecond pulsars constitute
most extreme physical laboratories and give opportunity to test several
theoretical models including the general theory of relativity Bisnovatyi-Kogan
(2010).

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
Harding (2018); Melrose et al. (2021); Philippov and Kramer (2022).
Nevertheless, the simple magnetic dipole model allows a basic understanding of
pulsar magnetosphere and the coherent radio emission in pulsars Gunn and
Ostriker (1969); Narayan and Ostriker (1990); van den Heuvel (2006). 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 Tchekhovskoy et al. (2013). These particles
are accelerated to relativistic velocities, and emit $\gamma$-ray photons by
either curvature radiation or inverse Compton scattering on lower energy
photons Ruderman and Sutherland (1975); Arons (1983); Daugherty and Harding
(1986). And then, the photons interacting with the magnetic field forms
electron-positron pairs Sturrock (1971). 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_{\rm 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 Baym et al.
(1969); Goldreich and Reisenegger (1992); Tauris and Manchester (1998); Tong
(2016). 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.

Figure 1: The conventional magnetic dipole model: a schematic view of the key
features of a pulsar magnetosphere Lorimer and Kramer (2005).

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 Kramer et al. (1999), 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 Bansal et al.
(2020). 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 Edwards et al. (2006):

$\phi(t)=\phi_{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. $\phi_{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 D’Alessandro
(1996). Timing noise is a slow, long-term, discernible, stochastic wandering
in the periods of pulsars. In 1969, two independent groups, Radhakrishnan &
Manchester Radhakrishnan and Manchester (1969) and Reichley & Downs Reichley
and Downs (1969) 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 Basu et al.
(2022) and 64-m Parkes Lower et al. (2021). 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
ATNF111https://www.atnf.csiro.au/people/pulsar/psrcat/glitchTbl.html $-$
accessed 2022 September Manchester et al. (2005) and also by Jodrell Bank
(JBO)222http://www.jb.man.ac.uk/ pulsar/glitches/gTable.html $-$ accessed 2022
September Espinoza et al. (2011). After combining entries from both the
catalogues there are a total of 719 glitches identified in 239 pulsars
Arumugam and Desai (2022). 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
Edwards et al. (2006):

$\displaystyle\phi_{\rm g}(t)=$ $\displaystyle\Delta\phi+\Delta\nu_{\rm
p}(t-t_{\rm g})+\frac{1}{2!}\Delta\dot{\nu}_{\rm p}(t-t_{\rm g})^{2}$ (2)
$\displaystyle+\sum_{i}^{N}\Delta\nu_{\mathrm{d}}^{(i)}\tau_{\mathrm{d}}^{(i)}\left[1-\exp\left({-\frac{(t-t_{\rm
g})}{\tau_{\mathrm{d}}^{(i)}}}\right)\right],$

whose value is zero if $t\leq t_{\rm g}$. $\Delta\phi$ is an offset in pulse
phase. The positive or negative sign of the permanent jumps in the spin
frequency $\Delta\nu_{\rm p}$ at the time of a glitch $t_{\mathrm{g}}$
represents usual glitches or anti-glitches, respectively.
$\Delta\dot{\nu}_{{\rm p}}$ is the permanent change in the frequency
derivative (spin-down rate) relative to the pre-glitch value. The fouth term
models the $i$ 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_{\rm g}$ and $\Delta\dot{\nu}_{\rm g}$ in a
glitch are described as:

$\Delta\nu_{\rm g}=\Delta\nu_{\rm
p}+\sum_{i}^{N}\Delta\nu_{\mathrm{d}}^{(i)},$ (3) $\Delta\dot{\nu}_{\rm
g}=\Delta\dot{\nu}_{\rm
p}-\sum_{i}^{N}\frac{\Delta\nu_{\mathrm{d}}^{(i)}}{\tau_{\mathrm{d}}^{(i)}}.$
(4)

In addition, the relative glitch amplitudes in frequency and frequency
derivative are $\Delta\nu_{\rm g}/\nu$ and $\Delta\dot{\nu}_{\rm
g}/\dot{\nu}$, respectively. The degree of recoveries are quantified by the
parameter $Q={\sum_{i}^{N}\Delta\nu_{\mathrm{d}}^{(i)}}/\Delta\nu_{\rm g}$.
Noteworthy, the maximum change of the spin frequency $\Delta\nu_{\rm max}$ and
$\Delta\dot{\nu}_{\rm max}$ its first time-derivative 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 Hotan et al. (2004) provide functionality for
generating TOAs; and then, these TOAs are used to fitted with a model that
contains a set of known parameters to achieve a phase-coherent timing
residuals ($R=(\phi-N)/\nu$, where $N$ is the nearest integer to $\phi$) Hobbs
et al. (2006), using softwares such as tempo Taylor et al. (1979), psrtime
Hobbs et al. (2006), tempo2 Hobbs et al. (2006), and pint Luo et al. (2021).
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. Freire and Ridolfi (2018) and Phillips et al. Phillips
and Ransom (2022) 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 the good post-glitch timing solutions. Melatos et al.
Melatos et al. (2020) 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_{\rm p}$ and spin-down frequency
$\Delta\dot{\nu}_{\rm p}$ during glitches with a HMM. Dunn et al. Dunn et al.
(2022) employed HMM to analyze low-cadence timing data to constrain the
magnitude of missing glitches in pulsars. Singha et al. Singha et al. (2021)
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 automatically
glitches. Soon, the searching techniques of glitches with massive human
intervention will be consigned to the past.

Figure 2: The most commonly employed data reduction routine for radio pulsar
glitches.

In China, some scientific achievements have been made in the area of glitches
research. Table 1 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
Zou et al. (2004) and significant changes in pulse profile associated with PSR
B2035+36’s glitch activity Kou et al. (2018) were discovered by NSRT. Ge et
al. Ge et al. (2020) and Zhang et al. Zhang et al. (2018) 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
continuously 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) Cong et al. (2017) and JRT (Jingdong 120-m Radio Telescope)
Qian et al. (2020) will provide an opportunity to track more pulsars, more
frequently, allowing glitches to be investigated in more detail.

Table 1: Chinese telescopes to perform a search for glitch events. Telescopes | Location | Diameter | Frequency | Pulsars | Start | No. of | Ref.
---|---|---|---|---|---|---|---
|  |  |  |  |  | Glitches |
NSRT | Nanshan | 26 m | 1.54 GHz | $\sim 300$ | 2000 | 103 | Zou et al. (2004); Yuan et al. (2010); Dang et al. (2020); Yuan et al. (2010); Kou et al. (2018); Liu et al. (2018); Zhou et al. (2022); Yuan et al. (2013); Wang et al. (2013); Zou et al. (2008); Yuan et al. (2017); Wang et al. (2001); Yuan et al. (2016); Wang et al. (2012)
TMRT | Shanghai | 65 m | 2.25/4.82/8.60 GHz | $\sim 100$ | 2013 | 5 | Liu et al. (2019, 2021, 2022, 2020)
Insight-HXMT | Space | – | 1-250 keV | $\sim 10$ | 2017 | 3 | Zhang et al. (2022); Ge et al. (2020)
KMRT | Yunnan | 40 m | 2.25/4.85 GHz | $\sim 90$ | 2008 | 1 | Xu et al. (2019)
HRT | Luonan | 40 m | 1.40 GHz | $\sim 10$ | 2014 | 1 | Luo et al. (2020)
XPNAV-1 | Space | – | 0.5-10 keV | $\sim 26$ | 2016 | 1 | Zhang et al. (2018)
FAST | Pingtang | 500 m | 1.25 GHz | $\sim 350$ | 2016 | – | Han et al. (2021)
QTT | Qitai | 110 m | – | – | – | – | Cong et al. (2017)
JRT | Jingdong | 120 m | – | – | – | – | Qian et al. (2020)

##  2 Properties of Pulsar Glitches

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\sim 10^{-9}$ do not lose phase coherence over several
hundred days, while those with sizes $\Delta\nu/\nu\sim 10^{-6}$ may see a
change in residuals of a large fraction of the pulse period in just a few days
Yu et al. (2013). 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 Janssen and Stappers (2006)
found that the minimum glitch size could be detected with a threshold of
$\Delta\nu/\nu\sim 10^{-11}$ by performing Monte Carlo simulations to test the
glitch detection method on PSR J0358+5413 Millhouse et al. (2022). The
millisecond pulsar J0613–0200 has suffered a glitch with a size
$\Delta\nu/\nu\sim 2.5\times 10^{-12}$, which is the smallest that has been
recorded McKee et al. (2016). Several studies have indicated that the
previously published pulsar timing results did not include all glitches in
individual pulsar or in particular datasets Janssen and Stappers (2006); Yu
and Liu (2017); Jankowski et al. (2019). The following reasons may result in
incomplete observations of glitches:

  * 1.

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 Shaw et al. (2021).

  * 2.

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. Shaw et al. (2018); Dang et al. (2020).

  * 3.

Glitches are too small to be resolved due to being below the present limit of
detectability Espinoza et al. (2011).

  * 4.

Pulsars that exhibit high levels of timing noise may have small glitches that
go undetected in the data Yuan et al. (2010).

  * 5.

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 Wong et
al. (2001); Espinoza et al. (2014).

  * 6.

When a glitch is identified to occur in a significant observation gap, it is
impossible to distinguish between a single glitch and multiple, smaller
glitches Yu and Liu (2017); Janssen and Stappers (2006).

Figure 3: Period–period-derivative ($P-\dot{P}$) diagram showing the pulsars
for which glitch is detected based on the ATNF Pulsar Catalogue Manchester et
al. (2005). The various sub-classes of glitching pulsars are represented by
the markers in the legend.

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. The known glitching pulsars are shown on the
the period–period-derivative ($P-\dot{P}$) diagram in Figure 3. Clearly,
glitches occur predominantly in young pulsars (characteristic ages $\sim
10^{3}-10^{7}\,\rm yr$ Yu et al. (2013)) located in the top right of the
diagram, but the size of and time between glitches do not depend only upon the
pulsar’s position in the $P-\dot{P}$ diagram Hobbs et al. (2002). The curious
glitches phenomena have been discovered from ordinary pulsars, to super-
powerful magnetic field harbouring magnetars Younes et al. (2020), to
millisecond pulsars Cognard and Backer (2004), 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 Antonopoulou et al. (2018); Ferdman et al. (2018),
which has undergone 53 glitches by September 2022 Ho et al. (2020). 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\sim 10^{-6}$ for most

of them. The post-glitch decay starts out exponential relaxation, with
different timescales from 1 minute 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 Dodson et al. (2007); Yuan et al. (2010). In
contrast, the glitch sizes of the Crab pulsar ($\Delta\nu/\nu\sim
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}$ Lyne et al. (1993); Yuan et al. (2010). 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 McKee et al.
(2016) and B1821–24 Cognard and Backer (2004)), 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. McKee et al. (2016)
indicated that the rate of glitches for MSPs is significantly lower than that
of the general population, with a probability of $\sim 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 Dib and Kaspi
(2014). 31 magnetars have been dicovered, 7 of which have exhibited a total of
22 glitches333http://www.physics.mcgill.ca/ pulsar/magnetar/main.html $-$
accessed 2022 September Olausen and Kaspi (2014); Jawor and Tauris (2022).
Chukwude & Urama Chukwude and Urama (2010) found that glitches have amplitudes
$\Delta\nu/\nu\gtrsim 10^{-9}$ and associated decreases in $\dot{\nu}$, 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. Yu et al. (2013), after a glitch, part of the sudden change in both
$\nu$ and $\dot{\nu}$ often recovers 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}$ McKee et al. (2016) to $10^{-3}$ Serim et al.
(2017). Near-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(2)\rm\ d$ in Vela puslar (PSR J0835–4510)
occurred in July 2010. The rapid spin-up timescales of glitches can be assumed
to be less than $5\rm\ s$ Palfreyman et al. (2018), $30\rm\ s$ Dodson et al.
(2007), and $40\rm\ s$ Dodson et al. (2002), 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 Howitt et al.
(2016); Yu et al. (2013) . 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 impermanent 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 Urama (2002); Yuan et al. (2017); Xu et al. (2019);
Zubieta et al. (2022). More specifically, Flanagan Flanagan (1990) reported
that the December 1988 Vela glitch best explained with three exponential terms
and Dodson et al. Dodson et al. (2002) confirmed four exponential terms
occurred in the January 2000 Vela glitch. In the early days following the
first glitch observation in the Vela pulsar Baym et al. Baym et al. (1969)
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 first evidence for
the existence of a superfluid interior of the neutron star Buchner and
Flanagan (2011); Baym et al. (1969). Soon after Anderson & Itoh Anderson and
Itoh (1975) 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 pulsar magnetosphere, but also the
dynamical process operating inside a neutron star.

Figure 4: The PSR J0835–4510’s (Vela pulsar) July 2010 glitch with an
exponential decay timescale of $13(2)\rm\ d$ (marked grey area) observed by
the Fermi-LAT Gügercinoğlu et al. (2022); Lower et al. (2021): (a) the
variations in spin frequency after subtracting the pre-glitch timing model;
(b) expanded plot of $\Delta\nu$ where the mean post-glitch value has been
subtracted from the post-glitch data; (c) the evolution of $\dot{\nu}$. The
vertical line indicates this glitch epoch at MJD 55408(1).

###  3.2 Slow Glitches

A total of 31 slow glitches have been observed in 7 pulsars Zhou et al.
(2019); Singha et al. (2022), 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) Yuan et al. (2010); Zou et al. (2004); Shabanova (1998).
This process repeats periodically, and the peak value of the glitch amplitude
is the same 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\rm\
nHz$ and $0.15\times 10^{-15}-3.15\times 10^{-15}\rm\ s^{-2}$ ; the
corresponding relative maximum changes are $\Delta\nu/\nu\sim 0.9\times
10^{-9}-31.2\times 10^{-9}$ and $|\Delta\dot{\nu}|/|\dot{\nu}|\sim 1.8\times
10^{-3}-23.6\times 10^{-3}$. The occurrence of slow glitches and normal
glitches in the PSR B1822–09 is of special interest Shabanova (2009). Zhou et
al. Zhou et al. (2019) 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. Theoretically, the angular momentum exchange model
Alpar et al. (1989); Eya et al. (2017), which is currently widely accepted to
explain the physical mechanism of normal glitches, cannot be applied to the
process of slow events. Link & Epstein Link and Epstein (1996) believed that
the slow glitches are the result of a sudden increase in the temperature of
the inner crust of the neutron star. However, Hobbs et al. Hobbs et al. (2004)
suggested that the slow glitches are just a manifestation of timing noise,
which is significantly dominated by sustained random wandering in either the
phase, spin, or spin-down rate Boynton et al. (1972). That means they 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.

Figure 5: A very large slow glitch in PSR J1062–5100 Zhou et al. (2019): (a)
timing residiuals with respect to the pre-glitch spin-down solution; (b)
variations of the frequency residual $\Delta\nu$ after subtracting the pre-
glitch spin-down model. (c) observed variations of $\dot{\nu}$. The two red
dashed vertical lines indicate that the slow glitch occurred between MJD 54733
and 55330.

###  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 rapid decay of $\dot{\nu}$. At
present, delayed spin-up processes have been detected in the glitches of three
pulsars, Crab (six times) (Lyne et al., 1992; Wong et al., 2001; Shaw et al.,
2018; Ge et al., 2020; Shaw et al., 2021), 1E 2259+586 (once) (Woods et al.,
2004) and SGR J1935+2154 (once) (Ge et al., 2022). 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 post-glitch transient variation are
prevalent in glitches, the low and irregular cadence observations are mainly
why they are not detected. Although the upper limit of the rising time scale
for the Vela pulsar glitches have been constrained to within 5 s Palfreyman et
al. (2018) by southern hemisphere observers monitoring near-continuously, no
delayed spin-up components have been observed. This may be another evidence
that the glitches of Crab and Vela pulsars have a different origin Manchester
(2018). 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 2020 SGR J1935+2154 glitch magnitude including the delayed increase is
the largest among the all glitches observed to date, and the transient decay
timescale $\sim 8\rm\ d$ is much longer than that of the Crab’s events. Ge et
al. Ge et al. (2022) found that the pulse profile of SGR J1935+2154 changed
rapidly and drastically after the unusual large glitch. Combined with the
detection of repeated FRB 200428 after the glitch with a delay of likely 3
days, Ge et al. Ge et al. (2022) argue that internal changes in young
magnetars generate large glitch to significantly alter the structure of the
magnetosphere and trigger X-ray bursts and FRB 200428. Just recently Younes et
al. Younes et al. (2022) reported a possible connection between a new glitch
occurred on October 5, 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 Ge et al.
(2022). The finding is also interpreted as generally consistent with the
Gügercinoǧlu & Alpar Gügercinoǧlu and Alpar (2019) theory. Moreover, it is
impossible to tell whether delayed spin-ups are only related to large jumps,
because if the delay time scale is too short or small jumps are too small to
be observed by existing observation equipment. 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 magnetar polar cap which
effectively extracts angular momentum from the star Younes et al. (2022).

Figure 6: The simulated July 2019 Crab delayed spin-up behaviour as presented
by Shaw et al. Shaw et al. (2021): (a) timing residiuals relative to a pre-
glitch spin-down model; (b) the time-dependence of the frequency residuals
$\Delta\nu$; (c) the evolution of spin-down rate $\dot{\nu}$. The vertical
line signifies the glitch around MJD 58687.565(4).

###  3.4 Anti-glitches

The so-called anti-glitches are first coined when Archibald et al. Archibald
et al. (2013) observed a sudden negative jump of spin frequency in the
magnetar IE 2259+586, accompanied by a X-ray outburst and a significant spin-
down rate change. Recently, Ray et al. Ray et al. (2019) reported three large
anti-glitches in the accreting ultraluminous X-ray source NGC 300 ULX1 during
its most recent outburst. Ducci et al. Ducci et al. (2015) 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 Archibald et al.
(2013); Şaşmaz Muş et al. (2014); Younes et al. (2020); Içdem et al. (2012);
Ray et al. (2019); Pintore et al. (2016); Vurgun et al. (2019). An & Archibald
An and Archibald (2019) 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
magnetospheric activity Manchester (2018). Tong Tong (2014) used the wind
braking model to successfully interpret the anti-glitch behaviours of magnetar
1E 2259+586. Huang & Geng Huang and Geng (2014) suggested that anti-glitches
occur as a result of the collision of a small mass body with pulsars. Garcia &
Ranea-Sandoval Garcia and Ranea-Sandoval (2015) interpreted that anti-glitches
is an obvious consequence of the cumulative decay of the internal toroidal
magnetic field component resulting in the rearrangement of the stellar
structure. Howitt & Melatos Howitt and Melatos (2022) 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.

Figure 7: The simulated two anti-glitch in magnetar IE 2259+586 around 2012
Archibald et al. (2013) : (a) timing residiuals relative to a spin-down model
before anti-glitch; (b) the frequency residuals $\Delta\nu$; (c) the spin-down
rate $\dot{\nu}$ of before and after anti-glitch. The red vertical dashed
lines mark the glitch epochs: 56035 Archibald et al. (2013), 56088.4 Hu et al.
(2014).

##  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 sheds light into the multicomponent
structure of neutron stars and leads 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 Pines et
al. (1974). 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 Bisnovatyi-Kogan (1970).
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 superfluid flow 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 Warszawski and Melatos (2008,
2011); Warszawski et al. (2012); Warszawski and Melatos (2013); Khomenko and
Haskell (2018); Lönnborn et al. (2019); Howitt et al. (2020). In the early
days of pulsar glitch research several experiments were devised with
superfluid Helium II filled containers accelerated and decelerated
occasionally by hand to mimic neutron star rotational evolution Tsakadze and
Tsakadze (1972, 1975, 1979, 1980). In recent years numerical simulations have
been conducted to visualize several aspects of pulsar glitch recovery Larson
and Link (2002); Sidery et al. (2010); van Eysden and Melatos (2011, 2012);
Xie and Zhang (2013); van Eysden and Melatos (2013, 2014); Howitt et al.
(2016); Graber et al. (2017).

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

Superfluidity is at the heart of the most glitch models Anderson and Itoh
(1975); Alpar et al. (1984); Jones (1993); Sedrakian et al. (1995); Link and
Epstein (1996); Ruderman et al. (1998); Sedrakian and Cordes (1999);
Pizzochero (2011). Moreover, pulsar timing noise may arise from erratic
coupling of superfluid with the normal matter inside a neutron star Alpar et
al. (1986); Jones (1990); Haskell (2011); Melatos and Link (2014). 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 Baym and Pethick (1979); Shaham (1980); Sauls
(1989); Pines (1999); Haskell and Melatos (2015); Haskell and Sedrakian
(2018); Sedrakian and Clark (2019).

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_{\rm ext}$ and internal superfluid torque $N_{\rm int}$ as

$I_{\rm n}\dot{\Omega}_{\rm n}=N_{\rm ext}+N_{\rm int}=N_{\rm ext}-I_{\rm
cs}\dot{\Omega}_{\rm cs}-I_{\rm core}\dot{\Omega}_{\rm 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\vec{v}_{\rm s}.d\vec{\ell}=2\pi r^{2}\Omega_{\rm cs}=N_{\rm
v}\kappa=\int\int n_{\rm v}\kappa dS=\int_{0}^{r}2\pi r^{\prime}n_{\rm
v}(r^{\prime})\kappa dr^{\prime}.$ (6)

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

$\frac{\partial\left(2\pi r^{2}\Omega_{\rm cs}\right)}{\partial r}=4\pi
r\Omega_{\rm cs}=2\pi rn_{\rm v}\kappa\rightarrow n_{\rm v}=\frac{2\Omega_{\rm
cs}}{\kappa}.$ (7)

Equation (7) implies that areal density of vortex lines are fully determined
by the rotation rate of the superfluid and in order to 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_{\rm v}}{\partial t}+\vec{\nabla}.\left(n_{\rm
v}\vec{v}_{\rm L}\right)=0.$ (8)

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

$\displaystyle\frac{\partial}{\partial t}\left(\oint\vec{v}_{\rm
s}.d\vec{\ell}\right)$ $\displaystyle=2\pi r^{2}\dot{\Omega}_{\rm
cs}=\frac{\partial}{\partial t}\left(\int\int n_{\rm v}\kappa
dS\right)=-\int\int[\vec{\nabla}.(n_{\rm v}\kappa v_{\rm L})]dS$
$\displaystyle=-\int\frac{1}{r}\frac{d}{dr}\left(rn_{\rm v}\kappa v_{\rm
L,r}\right)2\pi rdr,$ (9)

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

$\dot{\Omega}_{\rm cs}=-\frac{n_{\rm v}\kappa v_{\rm
L,r}}{r}=-\frac{2\Omega_{\rm cs}}{r}v_{\rm L,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_{\rm L,r}$. In order to find the radial
vortex line migration rate we consider the Magnus force equation,

$\vec{F}=\rho_{\rm s}\vec{\kappa}\times\left(\vec{v_{\rm s}}-\vec{v_{\rm
L}}\right),$ (11)

together with the force due to magnetised vortex electron scattering of the
form Alpar et al. (1984); Alpar and Sauls (1988); Sedrakian and Sedrakian
(1995)

$\vec{F}=C\left(\vec{v_{\rm c}}-\vec{v_{\rm L}}\right),$ (12)

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

$v_{\rm L,r}=\frac{\left(v_{\rm s}-v_{\rm
c}\right)}{\left[\left(\frac{\rho_{\rm
s}\kappa}{C}\right)\right]+\left[\left(\frac{C}{\rho_{\rm
s}\kappa}\right)\right]}.$ (13)

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

$\dot{\Omega}_{\rm core}=-\frac{2\Omega_{\rm core}}{r}v_{\rm
L,r}=-2\Omega_{\rm core}\frac{\left(\Omega_{\rm core}-\Omega_{\rm
c}\right)}{\left[\left(\frac{\rho_{\rm
s}\kappa}{C}\right)\right]+\left[\left(\frac{C}{\rho_{\rm
s}\kappa}\right)\right]}\equiv-\frac{\left(\Omega_{\rm core}-\Omega_{\rm
c}\right)}{\tau_{\rm core}},$ (14)

where the core superfluid’s relaxation time is defined as

$\tau_{\rm core}=\frac{\left[\left(\frac{\rho_{\rm
s}\kappa}{C}\right)\right]+\left[\left(\frac{C}{\rho_{\rm
s}\kappa}\right)\right]}{2\Omega_{\rm core}}.$ (15)

For the above considered process of magnetised vortex scattering of electrons
$\tau_{\rm core}\simeq(10-200)P$ seconds Sidery and Alpar (2009). The observed
fast relaxation of the overshooting post-glitch spin rate after the 2016 Vela
glitch within a minute Ashton et al. (2019) implies time variable coupling of
the core superfluid to the observed crustal rotation rate Gügercinoğlu and
Alpar (2020). The form of radial motion of vortex lines in the crustal
superfluid depends crucially on their microscopic interaction with the lattice
nuclei Alpar et al. (1984); Link et al. (1993). In the vortex creep model
Alpar et al. (1984, 1989), $v_{\rm L,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_{\rm L,r}=v_{0}\left[\exp{\left(-\frac{E_{\rm
p,out}}{kT}\right)}-\exp{\left(-\frac{E_{\rm p,in}}{kT}\right)}\right],$ (16)

where $v_{0}\sim 10^{5}-10^{7}$ cm s-1 is the microscopic trial velocity of a
vortex line around nuclear clusters throughout the neutron star crust
Gügercinoğlu and Alpar (2016) and, $E_{\rm p,out}$ and $E_{\rm p,in}$ are
defined as

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

and

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

Here $E_{\rm p,in}$ is the pinning energy for vortex-nucleus interaction and
$\omega_{\rm cr}$ is the maximum velocity difference between the superfluid
and the normal matter that a vortex line can withstand before it unpins. 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 Link et al.
(1993). From equations (16), (17) and (18) we find radial component of the
vortex line in the crustal superfluid as

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

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

$\dot{\Omega}_{\rm cs}=-\frac{4\Omega_{\rm
cs}}{r}v_{0}\exp{\left(-\frac{E_{\rm p}}{kT}\right)}\sinh{\left(\frac{E_{\rm
p}}{kT}\frac{\omega}{\omega_{\rm cr}}\right)}.$ (20)

The angular velocity difference $\omega\equiv\Omega_{\rm cs}-\Omega_{\rm c}$
between the rotation rates of the crustal superfluid $\Omega_{\rm cs}$ and the
observed crust $\Omega_{\rm c}$ evolves as Alpar et al. (1989); Gügercinoǧlu
and Alpar (2017)

$\dot{\omega}=\frac{I}{I_{\rm c}}\left[-\frac{N_{\rm
ext}(t)}{I}-\frac{4\Omega_{\rm cs}}{r}v_{0}\exp\left(-\frac{E_{\rm
p}}{kT}\right)\sinh\left(\frac{E_{\rm p}}{kT}\frac{\omega}{\omega_{\rm
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 Alpar et al.
(1984); Shibazaki and Lamb (1989); Alpar (1998); Larson and Link (1999)

$\dot{E}(t)=N_{\rm ext}\Omega_{\rm c}(t)-\frac{d}{dt}\left[\frac{1}{2}I_{\rm
c}\Omega_{\rm c}^{2}(t)+\frac{1}{2}I_{\rm cs}\Omega_{\rm
cs}^{2}(r,t)\right]=I_{\rm cs}|\dot{\Omega}_{\infty}|\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}_{\rm
c}=\dot{\Omega}_{\rm 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-3 vortex pinning
is energetically favourable Alpar (1977). 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 Alpar et al. (1984). 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 lattice nuclei.

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

$\eta\equiv\frac{r|\dot{\Omega}_{\infty}|}{4\Omega_{\rm
cs}v_{0}}\exp{\left(\frac{E_{\rm p}}{kT}\right)}.$ (23)

For $\eta\leq 1$ crustal superfluid is said to be in the linear regime and
since $\sinh{x}\cong x$ we obtain

$\dot{\Omega}_{\rm cs}=-\frac{E_{\rm p}}{kT}\frac{4\Omega_{\rm
cs}}{r\omega_{\rm cr}}v_{0}\exp{\left(-\frac{E_{\rm
p}}{kT}\right)}\omega\equiv\frac{\omega}{\tau_{\rm l}},$ (24)

where the linear creep relaxation time-scale $\tau_{\rm lin}$ is defined as

$\tau_{\rm lin}=\frac{r\varpi}{4\Omega_{\rm cs}v_{0}}\exp{\left(\frac{E_{\rm
p}}{kT}\right)}.$ (25)

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

$\varpi\equiv\frac{kT}{E_{\rm p}}\omega_{\rm cr}.$ (26)

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

$\omega_{\rm lin}(t)=\exp\left(-\frac{t}{\tau_{\rm
lin}}\right)\left[\omega(0)-\frac{1}{I_{\rm
c}}\int_{0}^{t}\exp\left(\frac{t^{\prime}}{\tau_{\rm lin}}\right)N_{\rm
ext}(t^{\prime})dt^{\prime}\right],$ (27)

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

$\dot{\Omega}_{\rm c}(t)=\frac{N_{\rm ext}(t)}{I_{\rm c}}+\frac{I_{\rm
cs}}{I\tau_{\rm lin}}\exp\left(-\frac{t}{\tau_{\rm
lin}}\right)\left[\omega(0)-\frac{1}{I_{\rm
c}}\int_{0}^{t}\exp\left(\frac{t^{\prime}}{\tau_{\rm lin}}\right)N_{\rm
ext}(t^{\prime})dt^{\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$ Gügercinoǧlu and Alpar (2017),
which accounts for the observed case seen after many magnetar glitches
displaying radiative changes Dib and Kaspi (2014).

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

$[\omega_{\infty}]_{\rm lin}=|\dot{\Omega}_{\infty}|\tau_{\rm lin},$ (29)

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

$\dot{\Omega}_{\rm c}=\dot{\Omega}_{\infty}-\frac{I_{\rm
cs}}{I}\frac{\delta\omega}{\tau_{\rm lin}}\exp\left(-\frac{t}{\tau_{\rm
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(0)-\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 Alpar et al.
(1989). Note that the response of vortex pinning and creep across flux tubes
to a spin up glitch also leads to exponential decay Gügercinoğlu and Alpar
(2014); Sourie and Chamel (2020). 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 Gügercinoğlu
(2017). 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 Yakovlev and Pethick (2004); Page et al. (2006). 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_{\rm lin}\approx 100$ d. Yu et al. (2013); Lower et al.
(2021); Liu et al. (2021), the steady state lag before a glitch is
$[\omega_{\infty}]_{\rm lin}\ll\omega_{\rm cr}$ for reasonable threshold
values for vortex unpinning $\omega_{\rm cr}\cong 10^{-2}-10^{-1}$ rad s-1
Alpar et al. (1984); Seveso et al. (2016) in these linear creep regions.
Therefore, linear creep regions cannot be responsible for the glitches by
initiating collective vortex unpinning avalanche Alpar et al. (1989).

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
$\sinh{x}\approx\exp(x)/2$ equation (20) becomes

$\dot{\Omega}_{\rm cs}=-\frac{2\Omega_{\rm
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_{\rm c}\tau_{\rm
lin}}\exp\left(\frac{\omega}{\varpi}\right)-\frac{N_{\rm ext}(t)}{I_{\rm c}}.$
(33)

The steady state lag in the nonlinear regime is given by

$[\omega_{\infty}]_{\rm nl}=\omega_{\rm cr}\left[1-\left(\frac{kT}{E_{\rm
p}}\right)\ln\left(\frac{4v_{0}\tau_{\rm
sd}}{r}\right)\right]\simeq\omega_{\rm 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}_{\rm diss}\cong I_{\rm s}\omega_{\rm cr}|\dot{\Omega}|.$ (35)

Equation (35) leads to a effective surface temperature

$T_{\rm s}=\left(\frac{\dot{E}_{\rm 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
Gonzalez and Reisenegger (2010); Becker et al. (2006); Zharikov et al. (2008);
Durant et al. (2012); Rangelov et al. (2017); Guillot et al. (2019); Abramkin
et al. (2021, 2022). Note also that other internal heating mechanisms, in
particular rotochemical heating add constructively for high surface
temperature detection in old pulsars Gonzalez and Reisenegger (2010); Köpp et
al. (2022). 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 entails angular momentum exchange with
ambient normal matter and thereby gives 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_{\rm s}$, i.e. $\Omega_{\rm s}\rightarrow\Omega_{\rm
s}-\delta\Omega_{\rm s}$ and crustal rotation rate acquires an increment of
amount $\Delta\Omega_{\rm c}$, i.e. $\Omega_{\rm c}\rightarrow\Omega_{\rm
c}+\Delta\Omega_{\rm c}$. The glitch magnitude $\Delta\Omega_{\rm c}$ is then
found from the angular momentum conservation at the time of the glitch event:

$I_{\rm c}\Delta\Omega_{\rm c}=I_{\rm s}\delta\Omega_{\rm s}.$ (37)

Thus, glitch magnitude is determined by two factors: the number of vortices
unpinned at the time of glitch, $\delta N_{\rm V}=2\pi r^{2}\delta\Omega_{\rm
s}/\kappa$ and the radial distance traversed by vortices before repinning.
$I_{\rm s}$ consists of two parts: The nonlinear creep regions with moment of
inertia $I_{\rm A}$ wherein vortices unpinned and regions with moment of
inertia $I_{\rm B}$ through which vortices do not creep at all but contribute
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_{\rm rot}=\Delta J\omega_{\infty}=I_{\rm
s}\omega_{\infty}\Delta\Omega_{\rm 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_{\rm s}$ is the moment of
inertia of the pinned superfluid and $\Delta\Omega_{\rm c}$ is the observed
spin-up in the crustal rotation rate. The thermal luminosity expected from
afterglow of glitches is given by van Riper et al. (1991); Cheng et al.
(1998); Tang and Cheng (2001); Hui and Cheng (2004)

$\Delta L=\frac{\Delta E_{\rm rot}}{\tau_{\rm th}},$ (39)

where $\tau_{\rm 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}_{\rm c}(t)=\dot{\Omega}_{\infty}-\frac{I_{\rm cs}}{I_{\rm
c}}\dot{\Omega}_{\infty}\left[1-\frac{1}{1+\left[\exp\left(\frac{t_{0}}{\tau_{\rm
nl}}\right)-1\right]\exp\left(-\frac{t}{\tau_{\rm nl}}\right)}\right],$ (40)

where we defined nonlinear relaxation time

$\tau_{\rm nl}\equiv\frac{\varpi}{|\dot{\Omega}_{\infty}|},$ (41)

and offset time

$t_{0}\equiv\frac{\delta\omega}{|\dot{\Omega}_{\infty}|}.$ (42)

Eq.(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_{\rm s}(r,0)=(1-r/\delta
r_{0})\delta\Omega_{0}$ (Alpar et al., 1996)

$\Delta\dot{\Omega}_{\rm c}(t)=\frac{I_{\rm A}}{I_{\rm
c}}\dot{\nu}_{0}\left(1-\frac{1-(\tau_{\rm
nl}/t_{0})\ln\left[1+(e^{t_{0}/\tau_{\rm nl}}-1)e^{-t/\tau_{\rm
nl}}\right]}{1-e^{-t/\tau_{\rm nl}}}\right).$ (43)

Here total moment of inertia of the non-linear creep regions affected from
vortex unpinning is indicated by $I_{\rm A}$. For times $\tau_{\rm nl}\lesssim
t\lesssim t_{0}$ Eq.(43) reduces to

$\frac{\Delta\dot{\Omega}_{\rm c}(t)}{\dot{\Omega}_{\rm c}}=\frac{I_{\rm
A}}{I_{\rm 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 Alpar and Baykal (2006):

$\frac{\Delta\Omega_{\rm c}}{\Omega_{\rm c}}=\left(\frac{I_{\rm
A}}{2I}+\frac{I_{\rm B}}{I}\right)\frac{\delta\Omega_{\rm s}}{\Omega_{\rm
c}},$ (45) $\frac{\Delta\dot{\Omega}_{\rm c}}{\dot{\Omega}_{\rm
c}}=\frac{I_{\rm A}}{I},$ (46) $\ddot{\Omega}_{\rm c}=\frac{I_{\rm
A}}{I}\frac{\dot{\Omega}_{\rm c}^{2}}{\delta\Omega_{\rm s}}.$ (47)

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

$t_{0}=\frac{\delta\Omega_{\rm s}}{|\dot{\Omega}|}=2\times
10^{-3}\left[\left(\Delta\Omega_{\rm c}/\Omega_{\rm
c}\right)_{-6}/\left(\beta+1/2\right)\left(\Delta\dot{\Omega}_{\rm
c}/\dot{\Omega}_{\rm c}\right)_{-3}\right]\tau_{\rm sd},$ (48)

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

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

Equation (49) estimates as large as $n_{\rm ig}\sim 100$ which compares well
with the observed values Yu et al. (2013); Dang et al. (2020); Lower et al.
(2021).

Vortex creep model has been applied to post-glitch timing data of a handful of
pulsars Alpar et al. (1984, 1993); Chau et al. (1993); Akbal et al. (2017);
Gügercinoğlu and Alpar (2020); Alpar et al. (1985, 1996); Gügercinoǧlu and
Alpar (2019); Alpar et al. (1988); Akbal et al. (2015); Gügercinoğlu et al.
(2022). From the vortex creep model fits to the post-glitch timing
observational data many invaluable information probing into neutron star
internal structure and dynamics can be obtained. Among them, total moments of
the 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 Gügercinoğlu et al. (2022). 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 Pines et al. (1980);
Howitt and Melatos (2022). Crustquakes as an external agent may drive some of
vortices radially inward at the time of glitch and acts as the required inward
bias Akbal et al. (2015). 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 Gügercinoǧlu and Alpar (2019). Slow glitches, on the
other hand, may be different manifestation of new vortex current distribution
in old pulsars after moving inward due to crust breaking quake.

Spin-down rate oscillations observed in the long term rotational evolution of
some pulsars Lyne et al. (2010); Kerr et al. (2016); Parthasarathy et al.
(2019); Shaw et al. (2022) as well as post glitch relaxation of the Vela
pulsar McCulloch et al. (1990) and PSR B2234+61 Yuan et al. (2010); Espinoza
et al. (2017) may result from the response of vortex line oscillations under
the competition between the pinning potential and tension Gügercinoğlu et al.
(2022).

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 Datta and Alpar (1993); Larson and Link (1999); Ho et al. (2015);
Pizzochero et al. (2017); Montoli et al. (2020); Shang and Li (2021).

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 Ashton et al. (2019). 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 Epstein and
Baym (1992); Jones (1992); Graber et al. (2018). The relativistic extension of
the mutual friction mechanism between the superfluid and normal matter Sourie
et al. (2017); Gavassino et al. (2020), and interpretation of the overshooting
of the rotation rate immediately following the glitch as time variable
coupling of neutron star core to the crust Pizzochero et al. (2020);
Gügercinoğlu and Alpar (2020); Montoli et al. (2020) 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 invoking 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 Alpar (1977); Link and
Epstein (1991). 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 Anderson and Itoh (1975). 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 Jones (1997); Seveso et
al. (2016); Wlazłowski et al. (2016); Link and Levin (2022). In particular an
amorphous crustal structure or a lattice with large impurities would result in
strong pinning wherein the unpinning process may start Jones (1999, 2001). A
vortex is forced to move through lattice which will cause it to fall into
interaction potential well provided by neighbouring nuclei, exiting waves on
the vortex. 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 gain due to being located at the bottom of the potential well and cost
due to lengthening after bending Bildsten and Epstein (1989); Epstein and Baym
(1992); Jones (1992); Hirasawa and Shibazaki (2001). 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. Avogadro et al. (2007, 2008) employ a mean
field Hartree-Fock-Bogoliubov formalism approach to vortex-nucleus interaction
and results in $\sim 3$ MeV for the pinning energy per intersection. According
to results of their calculations the interaction is in repulsive nature
throughout sizeable portion of the inner crust. Pizzochero and his
collaborators Pizzochero et al. (1997); Donati and Pizzochero (2003, 2004,
2006) have used local density approximation in their semi-classical
calculations and obtained a repulsive interaction of magnitude $\sim 1-2$ MeV
below densities $\sim 2\times 10^{13}$ gr cm-3 and an attractive interaction
of strength $\sim 5$ MeV at higher densities. An exact density functional
theory approach to the same problem has yielded a repulsive interaction with
$\sim 4$ MeV per nucleus for the density range $(3-7)\times 10^{13}$ gr cm-3
Wlazłowski et al. (2016). Jones Jones (1997) argues that for an infinite,
single component body centered cubic (bcc) lattice vortex line orientation
with respect to 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 Grill and Pizzochero
(2012). 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 Jones (1999, 2001). It is
speculated that above the neutron drip point the attractive force between the
interstitial unbound neutrons and the lattice makes bcc lattice unstable
Kobyakov and Pethick (2014). 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 Kobyakov
and Pethick (2016). 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 Link (2009); Seveso
et al. (2016); Antonelli and Haskell (2020). 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 Chamel and Haensel (2008). 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 Lazzari and de Blasio (1995a,
b); de Blasio and Lazzari (1996).

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 Andersson et al.
(2012); Chamel (2013); Li (2015); Delsate et al. (2016). Band theory
calculations of Chamel (2017) show that this entrainment effect may lead to
large effective to bare mass ratio for neutrons with average enhancement
factor $<m_{\rm n}^{*}/m_{\rm 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 Sauls et al. (2020). Uncertainties regarding
pairing potential may lead to smaller enhancement factors Kobyakov and Pethick
(2016); Martin and Urban (2016); Watanabe and Pethick (2017); Durel and Urban
(2018); Minami and Watanabe (2022). 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 Sekizawa et al. (2022). Also some parts of the neutron
star outer core may be involved in the glitches Gügercinoğlu and Alpar (2014);
Montoli et al. (2020), 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 Greenstein
(1970); Peralta et al. (2005, 2006); Melatos and Peralta (2007); Andersson et
al. (2007); Mongiovì et al. (2017); Haskell et al. (2020). Instabilities
associated with the superfluid hydrodynamical flow were proposed to play an
important role in initiating collective vortex unpinning cascade Andersson et
al. (2003); Mastrano and Melatos (2005); Sidery et al. (2008); Glampedakis and
Andersson (2009); Andersson et al. (2013); Khomenko et al. (2019). 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 van Eysden
and Link (2018).

###  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 spins up the star and
dissipates mechanical energy inside it, 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 Ruderman (1969); Baym
and Pines (1971). 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 Ruderman (1969)

$\frac{\delta R}{R}\approx\frac{95}{7}\frac{\mu R}{GM\rho}\Delta\phi.$ (51)

Here, $R$ is the radius and $M$ is the mass of the neutron star, $\rho$ is the
average stellar density, $\Delta\phi$ is the strain relieved in the quake and
$\mu\approx(\Delta/R)(Ze)^{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 Ogata and Ichimaru (1990); Strohmayer et
al. (1991); Baiko (2011); Zemlyakov and Chugunov (2022). In Ref.Alpar and
Pines (1985) it has been argued that since the only dimensionless parameter
related to the problem is the fine structure constant $\alpha=e^{2}/\hbar c$,
the critical strain angle for crust cracking should be of the order
$\theta_{\rm cr}=Ze\alpha$. The recent molecular dynamical simulations, which
take realistic shear motion of the neutron star crust into account, yield
$\theta_{\rm cr}\simeq 0.04-0.1$ Horowitz and Kadau (2009); Hoffman and Heyl
(2012); Baiko and Chugunov (2018), 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 Baym and Pines (1971) 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\epsilon^{2}+B(\epsilon-\epsilon_{0})^{2},$ (52)

where the oblateness $\epsilon$ is defined by the relation
$I=I_{0}(1+\epsilon)$ with $I_{0}$ being the moment of inertia of the
nonrotating spherical star and $\epsilon_{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 $L$ and $\epsilon_{0}$ fixed gives

$\epsilon=\frac{I\Omega^{2}}{4(A+B)}+\frac{B}{A+B}\epsilon_{0}.$ (53)

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

$\sigma=\frac{\partial
E}{\partial\epsilon}\frac{1}{V}=\mu(\epsilon_{0}-\epsilon),$ (54)

with $V$ being the crustal volume exceeds the critical value $\mu\theta_{\rm
cr}$. In a quake both $\epsilon$ and $\epsilon_{0}$ reduces according to

$\Delta\epsilon=\frac{B}{A+B}\Delta\epsilon_{0},$ (55)

where $\Delta\epsilon=\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_{\rm
q}\approx\frac{\Delta\epsilon_{0}-\Delta\epsilon}{\dot{\epsilon}}\approx\frac{A}{B}\frac{\Delta\epsilon}{2\epsilon}\frac{\Omega}{\dot{\Omega}},$
(56)

with the time evolution of the oblateness is given by

$\dot{\epsilon}\approx\frac{\partial
I}{\partial\epsilon}\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 Smoluchowski (1970); Cheng et al.
(1992); Gourgouliatos and Lander (2021). 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 Cheng et al.
(1992). The time interval between successive quakes is proportional to
$\Delta\epsilon_{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\epsilon/\epsilon)$ becomes
larger than the typical interval between the glitches Pines and Shaham (1972);
Pines et al. (1974).

Since the rigidity parameter $b=B/A$ is so small, change of the ellipticity
associated with a quake is tiny Cutler et al. (2003); Zdunik et al. (2008).
The actual value of the rigidity parameter $b$ is rather uncertain mainly
reflected in computation of $B$. Franco et al. Franco et al. (2000) 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 Baym and Pines (1971). Cutler et al. Cutler et al. (2003) 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
Rencoret et al. (2021); Gittins and Andersson (2021); Morales and Horowitz
(2022).

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 extended to include the effects of
general relativity Carter and Quintana (1975), magnetic field Franco et al.
(2000); Kojima et al. (2021), compressible neutron star model with density
dependent shear modulus Ushomirsky et al. (2000); Cutler et al. (2003).

Molecular dynamical simulations suggest large values of $\theta_{\rm cr}$
implying that the secular spin-down may only produce crustquakes just a few
times in the whole life of a neutron star Horowitz and Kadau (2009); Hoffman
and Heyl (2012); Baiko and Chugunov (2018).

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 Akbal et al. (2015); Akbal and
Alpar (2018); Gügercinoǧlu and Alpar (2019). 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 Link and Epstein (1996). 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 Epstein (1988).

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
increases Link et al. (1992, 1998); Epstein and Link (2000). This is proposed
as a viable explanation of persistent step increases observed after the Crab
glitches Wong et al. (2001); Lyne et al. (2015); Ge et al. (2020); Shaw et al.
(2021).

Progressive failures in the solid crust suffice to account for at least the
small glitches observed in pulsars of all ages. Giliberti et al. Giliberti et
al. (2019, 2020) 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 Braithwaite (2009). 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 Giliberti et al. (2019, 2020); Rencoret et
al. (2021); Gourgouliatos and Lander (2021).

Elastically deformed crust under stresses motivates searches for detection
gravitational wave emission from neutron stars Johnson-McDaniel and Owen
(2013). 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 Rencoret et al. (2021).

The crustquake model is also envisaged in terms of solid quark stars Peng and
Xu (2008); Zhou et al. (2014); Lai et al. (2018); Wang et al. (2021). 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 Ho et al. (2015). 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\sim
10^{-11}$ cm each carrying a field of $B_{\rm c}~{}\sim 10^{15}$ G Baym et al.
(1969). The number density of flux tubes

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

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

$n_{\rm V}=\frac{2\Omega_{\rm c}}{\kappa}=10^{5}\left(\frac{\Omega_{\rm
c}}{100\,\mbox{rad s}^{-1}}\right)\,\mbox{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\rightarrow BB_{\rm c}/4\pi$ while the
repulsion between the flux tubes diminishes greatly since in equilibrium
configuration they are separated well apart Easson and Pethick (1977);
Ruderman (1991). 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 Ruderman (1991); Ruderman et al. (1998)
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 Sauls (1989) and Srinivasan et al. Srinivasan et al.
(1990) 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 Jones (1991); Mendell (1991).
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
core’s 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 Ruderman (2009):

$\frac{BB_{\rm c}}{8\pi}L^{2}\gtrsim\left(\mu\theta_{\rm 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_{\rm c}$ at which expanding vortex line
array can push the flux tubes arrangement through core electron-proton plasma
with the maximum pace would be Ruderman et al. (1998)

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

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

$\beta=0.4\ln\left(\frac{\Lambda}{\xi}\right)\left(\frac{B_{\rm
V}}{10^{15}\,\mbox{G}}\right)\left(\frac{B_{\rm
c}}{10^{15}\,\mbox{G}}\right)\left(\frac{60\,\mbox{MeV}}{E_{\rm F_{\rm
e}}}\right)\left(\frac{10^{36}\,\mbox{cm}^{-3}}{n_{\rm 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_{\rm V}$ is the magnetic field attached to vortices, $B_{\rm c}$ is
the field strength of each flux tube, $E_{\rm F_{\rm e}}$ is the Fermi energy
of relativistic degenerate electrons and $n_{\rm 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 Ding et al. (1993);
Jahan-Miri (2000); Jones (2006); Bransgrove et al. (2018).

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 pushed by vortex lines Ruderman (2003). 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\sim s/R\ll\theta_{\rm max}\sim 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
$(BB_{\rm c}/8\pi)(s/R)$. Thus, the glitch magnitudes in spin and spin-down
rates are related by Ruderman (2005)

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

As a pulsar ages (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 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 Ruderman (2005)

$\frac{\Delta\Omega}{\Omega}\sim\left(\frac{5V_{\rm A}}{2V_{\rm
NS}}\right)\left(\frac{\tau_{\rm g}}{\tau_{\rm sd}}\right),$ (64)

where $V_{\rm A}$ is the volume of the annulus of the vortex line-flux tube
interaction region and $V_{\rm NS}$ is the volume of the neutron star. The
repetition time between successive Vela like giant glitches is given by
Ruderman (2005)

$\tau_{\rm g}\sim 3\times
10^{-4}\left(\frac{\Omega}{\dot{\Omega}}\right)|n-3|^{-1}.$ (65)

Equation (65) predicts $\tau_{\rm g}=0.7$ yr for $\tau_{\rm sd}=5$ kyr old
"Big Glitcher" PSR J0537–6910 which has interglitch braking index $n_{\rm
ig}\sim 7$ and observed glitch interval $\tau_{\rm ig}=0.4$ yr Middleditch et
al. (2006); Ferdman et al. (2018); Antonopoulou et al. (2018); Ho et al.
(2020); Abbott et al. (2021); Ho et al. (2022).

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 Ruderman et al. (1998)

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

where $\tau_{\rm sd}=\Omega/(2|\dot{\Omega}|)$ 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 Ruderman (2009). 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
pulsar population.

In all comparisons of observations with 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.
Jones (2002); Wang et al. (2000).

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 Sedrakyan and Sedrakyan (1993);
Gügercinoğlu (2017). The resulting heat production rate puts stringent limits
on the vigour of the vortex line-flux tube pinning model Ruderman et al.
(1998).

Much of the important model ingredients like how 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 Baym et al. (1969); Anderson and Itoh (1975); Pines and
Alpar (1985). 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
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 Basu et al. (2022).

Recently, Zhou et al. Zhou et al. (2022) reported that a micro-glitch
($\Delta\nu/\nu\sim 0.36(4)\times 10^{-9}$) occurred in 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\sim 33250\times 10^{-9}$) in PSR J1718–3718 Manchester and
Hobbs (2011) and a delayed spin-up event in the Crab pulsar Shaw et al.
(2018). 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 Weltevrede et al.
(2011). Shortly afterward, Archibald et al. Archibald et al. (2016) observed
another large glitch to be accompanied by X-ray bursts and X-ray pulsations in
this pulsar. The outburst timing anomalies phenomena 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 Espinoza et al. (2022). It is important to make timing observations
more frequently for high-B young glitching pulsars, such as PSRs J1718–3718
Zhu et al. (2011) and J1734–3333 Olausen et al. (2013), 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
Palfreyman et al. (2018). 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 Bransgrove et
al. (2020); Yuan et al. (2021).

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

Pulsar Name | $P$ | $\tau_{c}$ | $B_{s}$ | $\dot{E}$ | $\Delta\nu/\nu$ | $\Delta\dot{\nu}/\dot{\nu}$ | Profile | Ref.
---|---|---|---|---|---|---|---|---
(PSR) | (s) | (kyr) | ($10^{12}\rm\ G$) | ($10^{32}\rm\ erg/s$) | ($10^{-9}$) | ($10^{-3}$) |  |
J0738–4042 | 0.3749 | 4320 | 0.727 | 10 | 0.36(4) | 3(1) |  $W_{50}$($\downarrow$)a | Zhou et al. (2022)
J0742–2822 | 0.1667 | 157 | 1.69 | 1400 | 102.73(11) | 2.1(5) | b | Keith et al. (2013); Dang et al. (2021)
J1119–6127 | 0.4079 | 1.61 | 41 | 23000 | 9400(300) | 580(14) | c | Weltevrede et al. (2011)
|  |  |  |  | 5740(80) | 79(25) |  $X$($\uparrow$) | Archibald et al. (2016)
B1822–09 | 0.7690 | 233 | 6.42 | 45 | 4.08(2) | 0.08(1) | d | Liu et al. (2022)
|  |  |  |  | 7.2(1) | 1.65(7) | d | Liu et al. (2022)
J2021+4026 | 0.2653 | 76.9 | 3.85 | 1200 | <100(–) | 56(9) |  $\gamma$($\downarrow$) | Allafort et al. (2013); Takata et al. (2020)
B2021+51 | 0.5291 | 2740 | 1.29 | 8.2 | 0.373(5) | –0.24(3) |  $W_{\rm 10}$($\downarrow$)a | Liu et al. (2021)
B2035+36 | 0.6187 | 2180 | 1.69 | 7.5 | 7.7(8) | 67(8) |  $W_{50}$($\downarrow$)a | Kou et al. (2018)

  * a: $W_{\rm 10}$ and $W_{\rm 50}$ 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{\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.

Figure 8: Emission-rotation correlation in PSR J0738–4042 reported by Zhou et
al. Zhou et al. (2022): (a) variations of the frequency resuduals $\nu$; (b)
the variation of $\dot{\nu}$; (c) the post-fit residuals after adding the
glitch terms to the timing model. (d) the integrated normalized pulse profiles
at pre- (blue line) and post-glitch (red line) modes. The vertical line
indicates the glitch epoch at MJD $\sim 57359(5)$.

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$ Akbal et al. (2015); Ng et al. (2016);
Kou et al. (2018), 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. Liu et al. (2021, 2022)
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. Zhou et al. (2022) 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 Zhou et
al. (2022).

##  6 Statistics of Glitches

Figure 9: Cumulative number of the detected glitching pulsars (blue) and
glitches (red) over time. A quadratic fit (black) to the glitch discovery
$N\propto t^{2}$ for $50000<\rm MJD<58000$ is carried out. The discovery
epochs of glitching pulsars and glitches are collected from the combined
glitch catalogues at ATNF 1 and Jodrell Bank (JBO) 2.

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
Alpar and Baykal (1994); Lyne et al. (2000); Wang and Yuan (2010); Espinoza et
al. (2011); Yu et al. (2013); Eya et al. (2019); Montoli et al. (2021); Basu
et al. (2022); Gügercinoğlu et al. (2022); Eya et al. (2022). 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
50000 and 58000 in Fig. 9, obtaining 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. Espinoza et al. (2011) that both the
distributions of absolute and fractional frequency jumps, $\Delta\nu$ and
$\Delta\nu/\nu$, of all observed glitches show the bimodal nature. This
behaviour had also been indicated from recent studies using a larger sample Yu
et al. (2013); Fuentes et al. (2017); Ashton et al. (2017); Basu et al.
(2022); Lower et al. (2021). Basu et al. Basu et al. (2022) 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\sim 0.032\ \mu\rm Hz$ and the second,
narrow Gaussian component with larger magnitudes at $\Delta\nu\sim 18\ \mu\rm
Hz$. With a total sample of 700 glitches, Arumugam and Desai Arumugam and
Desai (2022) 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}$ (Figure 10), respectively. Moreover, there is
a pronounced dip at $10^{-7}$ presented in the $\Delta\nu/\nu$ distribution,
suggesting that the bimodality of the glitches sizes distribution may be
produced by two mechanisms by which the glitch events occur Yu et al. (2013);
Fuentes et al. (2017). 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. Celora et al. (2020) 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 Haskell et al.
(2018). 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 Basu et al. (2022).

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 Fuentes et al. (2019).
At the other extreme, the Crab pulsar glitches display temporal clustering at
times Lyne et al. (2015); Carlin et al. (2019). 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 of
pinning traps Pines et al. (1980); Gügercinoğlu and Alpar (2020). 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 differs 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 Carlin and Melatos (2020, 2021). 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 Melatos et al. (2008);
Fulgenzi et al. (2017); Howitt et al. (2018); Melatos et al. (2018); Melatos
and Drummond (2019); Fuentes et al. (2019); Carlin and Melatos (2019a, b).

Figure 10: Normalized histogram of the pulsar glitches amplitudes
log($\Delta\nu/\nu$) Arumugam and Desai (2022). There are 390 glitches which
belong to the first component and 310 to the second component. The narrower
Gaussian component (navy) centred at $1.3\times 10^{-6}$ represents the large
glitches, whereas the wider Gaussian component (red) centred at $5.9\times
10^{-9}$ represents the smaller glitches.

###  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 Espinoza et al. (2011); Yu et
al. (2013); Basu et al. (2022). According to Lyne et al. Lyne et al. (2000)
and Espinoza et al. Espinoza et al. (2011), the observed $\dot{N}$ is
significantly correlated with characteristic age $\tau_{c}$ and spin-down rate
$|\dot{\nu}|$. Basu et al. Basu et al. (2022) 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}_{\rm max}=1.07\ \rm yr^{-1}$ and
$\dot{N}_{\rm min}=0.02\ \rm yr^{-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 Lyne et al. (2000);
Fuentes et al. (2017); Basu et al. (2022):

$\dot{\nu}_{\rm 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}}_{\rm g}$ for groups of pulsars with similar properties
is examined Fuentes et al. (2017); Basu et al. (2022). The specific results
are: the linear relation between $\overline{\dot{\nu}}_{\rm g}$ and
$|\dot{\nu}|$ in the range $-14<\log|\dot{\nu}|<-10.5$ is $0.018(3)$, implying
that $|\dot{\nu}|$ undergoes a reversal of about 1.8 percent through glitches;
$\overline{\dot{\nu}}_{\rm g}$ decreases with $\tau_{c}$; $\dot{E}$ are
correlated positively with $\overline{\dot{\nu}}_{\rm g}$; pulsars with lower
$B$ also have lower $\overline{\dot{\nu}}_{\rm g}$, and an increase in
$\overline{\dot{\nu}}_{\rm g}$ is seen in neutron stars with the strongest $B$
(magnetars) Lyne et al. (2000); Espinoza et al. (2011); Fuentes et al. (2017);
Basu et al. (2022).

Figure 11: The accumulated $\Delta\nu/\nu\ (10^{-9})$ for the most frequent
glitching pulsars PSRs J0537$-$6910, J0835$-$4510 (Vela), J1341$-$6220 and
J1740$-$3015.

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 Alpar et
al. (1984); Cheng et al. (1988); Gügercinoğlu and Alpar (2020) 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 continues 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_{\rm gl}$ stores angular momentum by
average spin-down rate $\langle|\dot{\nu}|\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_{\rm c}$ gives
an lower limit on the fractional moment of inertia of the neutron star
component participated in glitches as Link et al. (1999); Andersson et al.
(2012); Chamel (2013)

$\frac{I_{\rm gl}}{I_{\rm c}}\gtrsim\frac{\dot{\nu}_{\rm
g}}{\langle|\dot{\nu}|\rangle}.$ (68)

Equation (68) is applied to the entire pulsar glitch population Ho et al.
(2015); Pizzochero et al. (2017); Basu et al. (2018); Parmar et al. (2022);
Montoli et al. (2021) 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 Gügercinoğlu and Alpar (2014); Montoli et
al. (2020); Sourie and Chamel (2020).

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.

Conceptualization, S.Q.Z., E.G., J.P.Y. and M.Y.G.; methodology, S.Q.Z. and
E.G.; software, S.Q.Z. and J.P.Y.; data curation, S.Q.Z.; writing—original
draft preparation, S.Q.Z. and E.G.; writing—review and editing, S.Q.Z., E.G.,
J.P.Y., M.Y.G. and C.Y.; visualization, S.Q.Z.; supervision, S.Q.Z., E.G.,
J.P.Y., M.Y.G. and C.Y.; project administration, C.Y., J.P.Y. and M.Y.G. All
authors have read and agreed to the published version of the manuscript.

This research was funded by the National Natural Science Foundation of China
(NSFC) grant via NSFC-11373064, 11521303, 11733010, 11873103, U2031121,
11873080, U1838201, U1838202, U1838104, U1938103, and U1938109. EG is
partially supported by NSFC programme via 11988101.

The authors declare no conflict of interest.

References

yes

## References

  * Hewish (1975) Hewish, A.  Pulsars and High Density Physics.  Science 1975, 188, 1079–1083.  doi:black10.1126/science.188.4193.1079.
  * Hartnett and Luiten (2011) Hartnett, J.G.; Luiten, A.N.  Colloquium: Comparison of astrophysical and terrestrial frequency standards.  Reviews of Modern Physics 2011, 83, 1–9, [arXiv:astro-ph.IM/1004.0115].  doi:black10.1103/RevModPhys.83.1.
  * Jenet et al. (2009) Jenet, F.; Finn, L.S.; Lazio, J.; Lommen, A.; McLaughlin, M.; Stairs, I.; Stinebring, D.; Verbiest, J.; Archibald, A.; Arzoumanian, Z.; Backer, D.; Cordes, J.; Demorest, P.; Ferdman, R.; Freire, P.; Gonzalez, M.; Kaspi, V.; Kondratiev, V.; Lorimer, D.; Lynch, R.; Nice, D.; Ransom, S.; Shannon, R.; Siemens, X.  The North American Nanohertz Observatory for Gravitational Waves.  arXiv e-prints 2009, p. arXiv:0909.1058, [arXiv:astro-ph.IM/0909.1058].
  * Hobbs et al. (2010) Hobbs, G.; Archibald, A.; Arzoumanian, Z.; Backer, D.; Bailes, M.; Bhat, N.D.R.; Burgay, M.; Burke-Spolaor, S.; Champion, D.; Cognard, I.; Coles, W.; Cordes, J.; Demorest, P.; Desvignes, G.; Ferdman, R.D.; Finn, L.; Freire, P.; Gonzalez, M.; Hessels, J.; Hotan, A.; Janssen, G.; Jenet, F.; Jessner, A.; Jordan, C.; Kaspi, V.; Kramer, M.; Kondratiev, V.; Lazio, J.; Lazaridis, K.; Lee, K.J.; Levin, Y.; Lommen, A.; Lorimer, D.; Lynch, R.; Lyne, A.; Manchester, R.; McLaughlin, M.; Nice, D.; Oslowski, S.; Pilia, M.; Possenti, A.; Purver, M.; Ransom, S.; Reynolds, J.; Sanidas, S.; Sarkissian, J.; Sesana, A.; Shannon, R.; Siemens, X.; Stairs, I.; Stappers, B.; Stinebring, D.; Theureau, G.; van Haasteren, R.; van Straten, W.; Verbiest, J.P.W.; Yardley, D.R.B.; You, X.P.  The International Pulsar Timing Array project: using pulsars as a gravitational wave detector.  Classical and Quantum Gravity 2010, 27, 084013, [arXiv:astro-ph.SR/0911.5206].  doi:black10.1088/0264-9381/27/8/084013.
  * Hobbs et al. (2020) Hobbs, G.; Guo, L.; Caballero, R.N.; Coles, W.; Lee, K.J.; Manchester, R.N.; Reardon, D.J.; Matsakis, D.; Tong, M.L.; Arzoumanian, Z.; Bailes, M.; Bassa, C.G.; Bhat, N.D.R.; Brazier, A.; Burke-Spolaor, S.; Champion, D.J.; Chatterjee, S.; Cognard, I.; Dai, S.; Desvignes, G.; Dolch, T.; Ferdman, R.D.; Graikou, E.; Guillemot, L.; Janssen, G.H.; Keith, M.J.; Kerr, M.; Kramer, M.; Lam, M.T.; Liu, K.; Lyne, A.; Lazio, T.J.W.; Lynch, R.; McKee, J.W.; McLaughlin, M.A.; Mingarelli, C.M.F.; Nice, D.J.; Osłowski, S.; Pennucci, T.T.; Perera, B.B.P.; Perrodin, D.; Possenti, A.; Russell, C.J.; Sanidas, S.; Sesana, A.; Shaifullah, G.; Shannon, R.M.; Simon, J.; Spiewak, R.; Stairs, I.H.; Stappers, B.W.; Swiggum, J.K.; Taylor, S.R.; Theureau, G.; Toomey, L.; van Haasteren, R.; Wang, J.B.; Wang, Y.; Zhu, X.J.  A pulsar-based time-scale from the International Pulsar Timing Array.  Mon. Not. R. Astron. Soc. 2020, 491, 5951–5965, [arXiv:astro-ph.IM/1910.13628].  doi:black10.1093/mnras/stz3071.
  * Bisnovatyi-Kogan (2010) Bisnovatyi-Kogan, G.  Binary recycled pulsars: a powerful physical laboratory.  Mem. Soc. Astron. Italiana 2010, 81, 258.
  * Harding (2018) Harding, A.K.  Pulsar Emission Physics: The First Fifty Years.  Pulsar Astrophysics the Next Fifty Years; Weltevrede, P.; Perera, B.B.P.; Preston, L.L.; Sanidas, S., Eds., 2018, Vol. 337, pp. 52–57, [arXiv:astro-ph.HE/1712.02409].  doi:black10.1017/S1743921318000509.
  * Melrose et al. (2021) Melrose, D.B.; Rafat, M.Z.; Mastrano, A.  Pulsar radio emission mechanisms: a critique.  Mon. Not. R. Astron. Soc. 2021, 500, 4530–4548, [arXiv:astro-ph.HE/2006.15243].  doi:black10.1093/mnras/staa3324.
  * Philippov and Kramer (2022) Philippov, A.; Kramer, M.  Pulsar Magnetospheres and Their Radiation.  ARA&A 2022, 60, 495–558.  doi:black10.1146/annurev-astro-052920-112338.
  * Gunn and Ostriker (1969) Gunn, J.E.; Ostriker, J.P.  Magnetic Dipole Radiation from Pulsars.  Nature 1969, 221, 454–456.  doi:black10.1038/221454a0.
  * Narayan and Ostriker (1990) Narayan, R.; Ostriker, J.P.  Pulsar Populations and Their Evolution.  Astrophys. J. 1990, 352, 222.  doi:black10.1086/168529.
  * van den Heuvel (2006) van den Heuvel, E.P.J.  Pulsar magnetospheres and pulsar death.  Science (New York, N.Y.) 2006, 312, 539—540.  doi:black10.1126/science.1125934.
  * Tchekhovskoy et al. (2013) Tchekhovskoy, A.; Spitkovsky, A.; Li, J.G.  Time-dependent 3D magnetohydrodynamic pulsar magnetospheres: oblique rotators.  Mon. Not. R. Astron. Soc. 2013, 435, L1–L5, [arXiv:astro-ph.HE/1211.2803].  doi:black10.1093/mnrasl/slt076.
  * Ruderman and Sutherland (1975) Ruderman, M.A.; Sutherland, P.G.  Theory of pulsars: polar gaps, sparks, and coherent microwave radiation.  Astrophys. J. 1975, 196, 51–72.  doi:black10.1086/153393.
  * Arons (1983) Arons, J.  Pair creation above pulsar polar caps : geometrical structure and energetics of slot gaps.  Astrophys. J. 1983, 266, 215–241.  doi:black10.1086/160771.
  * Daugherty and Harding (1986) Daugherty, J.K.; Harding, A.K.  Compton Scattering in Strong Magnetic Fields.  Astrophys. J. 1986, 309, 362.  doi:black10.1086/164608.
  * Sturrock (1971) Sturrock, P.A.  A Model of Pulsars.  Astrophys. J. 1971, 164, 529.  doi:black10.1086/150865.
  * Baym et al. (1969) Baym, G.; Pethick, C.; Pikes, D.  Electrical Conductivity of Neutron Star Matter.  Nature 1969, 224, 674–675.  doi:black10.1038/224674a0.
  * Goldreich and Reisenegger (1992) Goldreich, P.; Reisenegger, A.  Magnetic Field Decay in Isolated Neutron Stars.  Astrophys. J. 1992, 395, 250.  doi:black10.1086/171646.
  * Tauris and Manchester (1998) Tauris, T.M.; Manchester, R.N.  On the Evolution of Pulsar Beams.  Mon. Not. R. Astron. Soc. 1998, 298, 625–636.  doi:black10.1046/j.1365-8711.1998.01369.x.
  * Tong (2016) Tong, H.  Pulsar braking: magnetodipole vs. wind.  Science China Physics, Mechanics, and Astronomy 2016, 59, 5752, [arXiv:astro-ph.HE/1506.04605].  doi:black10.1007/s11433-015-5752-x.
  * Lorimer and Kramer (2005) Lorimer, D.; Kramer, M.  Handbook of pulsar astronomy; Vol. 4, Cambridge University Press, 2005.
  * Kramer et al. (1999) Kramer, M.; Xilouris, K.M.; Camilo, F.; Nice, D.J.; Backer, D.C.; Lange, C.; Lorimer, D.R.; Doroshenko, O.; Sallmen, S.  Profile Instabilities of the Millisecond Pulsar PSR J1022+1001.  Astrophys. J. 1999, 520, 324–334, [arXiv:astro-ph/astro-ph/9903048].  doi:black10.1086/307449.
  * Bansal et al. (2020) Bansal, K.; Taylor, G.B.; Stovall, K.; Dowell, J.  Detection of Echoes in PSR B1508+55 at Frequencies below 100 MHz Using the LWA1.  Astrophys. J. 2020, 892, 26, [arXiv:astro-ph.SR/2003.03606].  doi:black10.3847/1538-4357/ab76bc.
  * Edwards et al. (2006) Edwards, R.T.; Hobbs, G.B.; Manchester, R.N.  TEMPO2, a new pulsar timing package - II. The timing model and precision estimates.  Mon. Not. R. Astron. Soc. 2006, 372, 1549–1574, [arXiv:astro-ph/astro-ph/0607664].  doi:black10.1111/j.1365-2966.2006.10870.x.
  * D’Alessandro (1996) D’Alessandro, F.  Rotational irregularities in pulsars — A review.  Ap&SS 1996, 246, 73–106.  doi:black10.1007/BF00637401.
  * Radhakrishnan and Manchester (1969) Radhakrishnan, V.; Manchester, R.N.  Detection of a Change of State in the Pulsar PSR 0833-45.  Nature 1969, 222, 228–229.  doi:black10.1038/222228a0.
  * Reichley and Downs (1969) Reichley, P.E.; Downs, G.S.  Observed Decrease in the Periods of Pulsar PSR 0833-45.  Nature 1969, 222, 229–230.  doi:black10.1038/222229a0.
  * Basu et al. (2022) Basu, A.; Shaw, B.; Antonopoulou, D.; Keith, M.J.; Lyne, A.G.; Mickaliger, M.B.; Stappers, B.W.; Weltevrede, P.; Jordan, C.A.  The Jodrell bank glitch catalogue: 106 new rotational glitches in 70 pulsars.  Mon. Not. R. Astron. Soc. 2022, 510, 4049–4062, [arXiv:astro-ph.HE/2111.06835].  doi:black10.1093/mnras/stab3336.
  * Lower et al. (2021) Lower, M.E.; Johnston, S.; Dunn, L.; Shannon, R.M.; Bailes, M.; Dai, S.; Kerr, M.; Manchester, R.N.; Melatos, A.; Oswald, L.S.; Parthasarathy, A.; Sobey, C.; Weltevrede, P.  The impact of glitches on young pulsar rotational evolution.  Mon. Not. R. Astron. Soc. 2021, 508, 3251–3274, [arXiv:astro-ph.HE/2109.07612].  doi:black10.1093/mnras/stab2678.
  * Manchester et al. (2005) Manchester, R.N.; Hobbs, G.B.; Teoh, A.; Hobbs, M.  The Australia Telescope National Facility Pulsar Catalogue.  AJ 2005, 129, 1993–2006, [arXiv:astro-ph/astro-ph/0412641].  doi:black10.1086/428488.
  * Espinoza et al. (2011) Espinoza, C.M.; Lyne, A.G.; Stappers, B.W.; Kramer, M.  A study of 315 glitches in the rotation of 102 pulsars.  Mon. Not. R. Astron. Soc. 2011, 414, 1679–1704, [arXiv:astro-ph.HE/1102.1743].  doi:black10.1111/j.1365-2966.2011.18503.x.
  * Arumugam and Desai (2022) Arumugam, S.; Desai, S.  Classification of Pulsar Glitch Amplitudes using Extreme Deconvolution.  arXiv e-prints 2022, p. arXiv:2206.02751, [arXiv:astro-ph.HE/2206.02751].
  * Edwards et al. (2006) Edwards, R.T.; Hobbs, G.B.; Manchester, R.N.  TEMPO2, a new pulsar timing package - II. The timing model and precision estimates.  Mon. Not. R. Astron. Soc. 2006, 372, 1549–1574, [arXiv:astro-ph/astro-ph/0607664].  doi:black10.1111/j.1365-2966.2006.10870.x.
  * Hotan et al. (2004) Hotan, A.W.; van Straten, W.; Manchester, R.N.  PSRCHIVE and PSRFITS: An Open Approach to Radio Pulsar Data Storage and Analysis.  Publ. Astron. Soc. Austr. 2004, 21, 302–309, [arXiv:astro-ph/astro-ph/0404549].  doi:black10.1071/AS04022.
  * Hobbs et al. (2006) Hobbs, G.B.; Edwards, R.T.; Manchester, R.N.  TEMPO2, a new pulsar-timing package - I. An overview.  Mon. Not. R. Astron. Soc. 2006, 369, 655–672, [arXiv:astro-ph/astro-ph/0603381].  doi:black10.1111/j.1365-2966.2006.10302.x.
  * Taylor et al. (1979) Taylor, J.H.; Fowler, L.A.; McCulloch, P.M.  Measurements of general relativistic effects in the binary pulsar PSR1913 + 16.  Nature 1979, 277, 437–440.  doi:black10.1038/277437a0.
  * Luo et al. (2021) Luo, J.; Ransom, S.; Demorest, P.; Ray, P.S.; Archibald, A.; Kerr, M.; Jennings, R.J.; Bachetti, M.; van Haasteren, R.; Champagne, C.A.; Colen, J.; Phillips, C.; Zimmerman, J.; Stovall, K.; Lam, M.T.; Jenet, F.A.  PINT: A Modern Software Package for Pulsar Timing.  Astrophys. J. 2021, 911, 45, [arXiv:astro-ph.IM/2012.00074].  doi:black10.3847/1538-4357/abe62f.
  * Freire and Ridolfi (2018) Freire, P.C.C.; Ridolfi, A.  An algorithm for determining the rotation count of pulsars.  Mon. Not. R. Astron. Soc. 2018, 476, 4794–4805, [arXiv:astro-ph.IM/1802.07211].  doi:black10.1093/mnras/sty524.
  * Phillips and Ransom (2022) Phillips, C.; Ransom, S.  Algorithmic Pulsar Timing.  AJ 2022, 163, 84.  doi:black10.3847/1538-3881/ac403e.
  * Melatos et al. (2020) Melatos, A.; Dunn, L.M.; Suvorova, S.; Moran, W.; Evans, R.J.  Pulsar Glitch Detection with a Hidden Markov Model.  Astrophys. J. 2020, 896, 78, [arXiv:astro-ph.HE/2005.09388].  doi:black10.3847/1538-4357/ab9178.
  * Dunn et al. (2022) Dunn, L.; Melatos, A.; Suvorova, S.; Moran, W.; Evans, R.J.; Osłowski, S.; Lower, M.E.; Bailes, M.; Flynn, C.; Gupta, V.  Systematic upper limits on the size of missing pulsar glitches in the first UTMOST open data release.  Mon. Not. R. Astron. Soc. 2022, 512, 1469–1482, [arXiv:astro-ph.HE/2202.12442].  doi:black10.1093/mnras/stac551.
  * Singha et al. (2021) Singha, J.; Basu, A.; Krishnakumar, M.A.; Joshi, B.C.; Arumugam, P.  A real-time automated glitch detection pipeline at Ooty Radio Telescope.  Mon. Not. R. Astron. Soc. 2021, 505, 5488–5496, [arXiv:astro-ph.IM/2106.01942].  doi:black10.1093/mnras/stab1640.
  * Zou et al. (2004) Zou, W.Z.; Wang, N.; Wang, H.X.; Manchester, R.N.; Wu, X.J.; Zhang, J.  Unusual glitch behaviours of two young pulsars.  Mon. Not. R. Astron. Soc. 2004, 354, 811–814, [arXiv:astro-ph/astro-ph/0407483].  doi:black10.1111/j.1365-2966.2004.08241.x.
  * Kou et al. (2018) Kou, F.F.; Yuan, J.P.; Wang, N.; Yan, W.M.; Dang, S.J.  The spin-down state change and mode change associated with glitch activity of PSR B2035+36.  Mon. Not. R. Astron. Soc. 2018, 478, L24–L28, [arXiv:astro-ph.HE/1801.01248].  doi:black10.1093/mnrasl/sly068.
  * Ge et al. (2020) Ge, M.Y.; Zhang, S.N.; Lu, F.J.; Li, T.P.; Yuan, J.P.; Zheng, X.P.; Huang, Y.; Zheng, S.J.; Chen, Y.P.; Chang, Z.; Tuo, Y.L.; Cheng, Q.; Güngör, C.; Song, L.M.; Xu, Y.P.; Cao, X.L.; Chen, Y.; Liu, C.Z.; Zhang, S.; Qu, J.L.; Bu, Q.C.; Cai, C.; Chen, G.; Chen, L.; Chen, M.Z.; Chen, T.X.; Chen, Y.B.; Cui, W.; Cui, W.W.; Deng, J.K.; Dong, Y.W.; Du, Y.Y.; Fu, M.X.; Gao, G.H.; Gao, H.; Gao, M.; Gu, Y.D.; Guan, J.; Guo, C.C.; Han, D.W.; Hao, L.F.; Huo, J.; Jia, S.M.; Jiang, L.H.; Jiang, W.C.; Jin, C.J.; Jin, J.; Jin, Y.J.; Kong, L.D.; Li, B.; Li, D.; Li, C.K.; Li, G.; Li, M.S.; Li, W.; Li, X.; Li, X.B.; Li, X.F.; Li, Y.G.; Li, Z.W.; Li, Z.X.; Liu, Z.Y.; Liang, X.H.; Liao, J.Y.; Liu, G.Q.; Liu, H.W.; Liu, X.J.; Liu, Y.N.; Lu, B.; Lu, X.F.; Luo, Q.; Luo, T.; Ma, X.; Meng, B.; Nang, Y.; Nie, J.Y.; Ou, G.; Sai, N.; Shang, R.C.; Song, X.Y.; Sun, L.; Tan, Y.; Tao, L.; Wang, C.; Wang, G.F.; Wang, J.; Wang, J.B.; Wang, M.; Wang, N.; Wang, W.S.; Wang, Y.D.; Wang, Y.S.; Wen, X.Y.; Wen, Z.G.; Wu, B.B.; Wu, B.Y.; Wu, M.; Xiao, G.C.; Xiao, S.; Xiong, S.L.; Xu, Y.H.; Yan, W.M.; Yang, J.W.; Yang, S.; Yang, Y.J.; Yang, Y.J.; Yi, Q.B.; Yin, Q.Q.; You, Y.; Yue, Y.L.; Zhang, A.M.; Zhang, C.M.; Zhang, D.P.; Zhang, F.; Zhang, H.M.; Zhang, J.; Zhang, T.; Zhang, W.C.; Zhang, W.; Zhang, W.Z.; Zhang, Y.; Zhang, Y.F.; Zhang, Y.J.; Zhang, Y.; Zhang, Z.; Zhang, Z.; Zhang, Z.L.; Zhao, H.S.; Zhao, X.F.; Zheng, W.; Zhou, D.K.; Zhou, J.F.; Zhou, X.; Zhuang, R.L.; Zhu, Y.X.; Zhu, Y.  Discovery of Delayed Spin-up Behavior Following Two Large Glitches in the Crab Pulsar, and the Statistics of Such Processes.  Astrophys. J. 2020, 896, 55, [arXiv:astro-ph.HE/2004.00791].  doi:black10.3847/1538-4357/ab8db6.
  * Zhang et al. (2018) Zhang, X.; Shuai, P.; Huang, L.; Chen, S.; Du, Y.  X-Ray Observation of the 2017 November Glitch in the Crab Pulsar.  Astrophys. J. 2018, 866, 82.  doi:black10.3847/1538-4357/aade46.
  * Cong et al. (2017) Cong, Z.; Qian, X.; Na, W.; Binbin, X.  The Prime Focus Receiver Positioner Design of the Xinjiang QiTai 110m Radio Telescope.  Astronomical Research and Technology 2017, 14, 172–178.
  * Qian et al. (2020) Qian, L.; Yao, R.; Sun, J.; Xu, J.; Pan, Z.; Jiang, P.  FAST: Its Scientific Achievements and Prospects.  The Innovation 2020, 1, 100053, [arXiv:astro-ph.IM/2011.13542].  doi:black10.1016/j.xinn.2020.100053.
  * Yuan et al. (2010) Yuan, J.P.; Wang, N.; Manchester, R.N.; Liu, Z.Y.  29 glitches detected at Urumqi Observatory.  Mon. Not. R. Astron. Soc. 2010, 404, 289–304, [arXiv:astro-ph.GA/1001.1471].  doi:black10.1111/j.1365-2966.2010.16272.x.
  * Dang et al. (2020) Dang, S.J.; Yuan, J.P.; Manchester, R.N.; Li, L.; Wang, N.; Wang, J.B.; Hobbs, G.; Liu, Z.Y.; Kou, F.F.  Results of 12 yr of Pulsar Timing at Nanshan. I.  Astrophys. J. 2020, 896, 140, [arXiv:astro-ph.HE/2005.02200].  doi:black10.3847/1538-4357/ab9082.
  * Yuan et al. (2010) Yuan, J.P.; Manchester, R.N.; Wang, N.; Zhou, X.; Liu, Z.Y.; Gao, Z.F.  A Very Large Glitch in PSR B2334+61.  Astrophys. J. Lett. 2010, 719, L111–L115, [arXiv:astro-ph.HE/1007.1143].  doi:black10.1088/2041-8205/719/2/L111.