On the distributions of pulsar glitch sizes and the inter-glitch time intervals

Pulsars are believed to be spinning neutron stars (Gold 1968), core remnants of supernova events (Baade & Zwicky 1934b,a). They are very compact and their central densities exceed the standard saturation nuclear density. This is a feature that makes them the second densest objects observed in the universe. The internal temperatures of radio pulsars are thought to be below the neutron superfluid transition temperature, suggesting that their interiors contain superfluid neutrons (Migdal 1959; Baym et al. 1969b; Shternin et al. 2011). The superfluid neutrons spin at a higher velocity than the observed velocity of the solid crust, enabling it to act as a reservoir of angularmomentum as we will discuss later.

A pulsar glitch manifests as a sudden increase (spin-up) in spin frequency (Δν) of an otherwise steadily spinning down pulsar, which is sometimes accompanied by a change in the magnitude of the spin-down rate ($\Delta \dot{\nu }$) (Chukwude & Urama 2010; Espinoza et al. 2011; Yu et al. 2013). For many pulsars, one glitch is usually followed by a recovery phase during which the pulsar returns to a steady spin-down state on a wide range of time scale (Wong et al. 2001; Dodson et al. 2007). These behaviors could be attributed to an internal process of a pulsar, which makes the pulsar glitch a feasible tool for studying the dynamic behavior of super-dense matter (such as neutron superfluid) inside the star.

For nearly half a century since the first known pulsar glitch was discovered, the causes of pulsar glitches are still speculative despite numerous works put forward in explaining these events. There could be external and internal origins for pulsar glitches, but for the fact that early known glitches in rotation-powered pulsars were not accompanied by radiative changes, the causes of glitches were strongly attached to the internal mechanism of pulsars (for a recent review, see Haskell & Melatos (2015)). These causes can be broadly grouped into two depending on the mechanism involved, namely: the starquake model and the angular momentum transfer model. The starquake model relies on a sudden reduction in the pulsar's moment of inertia arising from a change in the oblateness of the star (e.g. Alpar et al. 1994, 1996; Franco et al. 2000). When the oblateness of a spherically rotating body is reduced instantaneously, such a body will spin up to compensate for the reduction in moment of inertia. Such a model is convenient in explaining small-size glitches (SSGs) in young pulsars (≤ 10³ yr), such as the Crab pulsar (Lyne et al. 1993; Alpar et al. 1994). However, for a middle-aged pulsar ( ≈ 10⁴ yr) such as the Vela pulsar, the starquake model is not sufficient (Alpar et al. 1996; Warszawski & Melatos 2013). A more promising model for explaining pulsar glitches relies on the exchange of angular momentum between the star's components that are weakly coupled together (Baym et al. 1969a; Anderson& Itoh 1975; Alpar et al. 1984). The components are the solid crust coupled to the core and the superfluid neutrons permeating into the inner-crust (crustal fluid) of the star. The crustal fluid, viewed as an angular momentum reservoir, rotates via a quantized vortex array. These quantized vortices carry the circulation of the superfluid that mimics the rotation of normal fluid. The number of vortices contained in the superfluid determines its angular velocity. These vortices have the ability to transfer their momentum to the solid crust in the right conditions. Furthermore, these superfluid vortices are pinned in the lattice of the inner-crust (Anderson & Itoh 1975). The pinned vortices enable the crustal fluid to partially decouple from the other components spinning at the same rate with the solid crust. For as long as the vortices remain pinned, the superfluid velocity is conserved. In this frame, as the crust spins down, the crustal fluid maintains its natural velocity. The velocity difference between the two components (i.e. rotation lag) increases with time and could reach a point when the pinned vortices can no longer withstand the load on them due to the rotational lag. An as yet unclear mechanism triggers a spontaneous unpinning of some (or all of the) vortices. The unpinned vorticesmigrate outward transferring their momentum to the crust, the superfluid spins down and the crust spins up (i.e. glitch) (Anderson & Itoh 1975; Alpar et al. 1984). Glitches occurring in this way should follow a canonical order in the sense that the accumulation and release of momentum are a function of pulsar rotational parameters (i.e. the spin frequency,ν, and the spin-down rate, $\dot{\nu }$). This implies that there should be a form of correlation between glitch sizes and inter-glitch time intervals. However, this logic is not observed in many glitch events.

On the other hand, recently, a possible relic from the early universe could also be a trigger of a pulsar glitch (Lai & Xu 2016). Such a relic, like "strange nuggets," which are believed to mostly exist as dark matter, could collide with a pulsar to trigger a glitch and resemble that arising from a starquake. Lai & Xu (2016) have shown that the collision rate of pulsars with strange nuggets is consistent with the occurrence rate of SSGs, though this model would not be able to explain glitches in pulsars with mixed glitch sizes.

Pulsar glitches are a rare phenomenon. Since the first observed glitch in the Vela and Crab pulsars (PSR J0835–4510 and PSR J0534+2200 respectively) (Radhakrishnan & Manchester 1969; Boynton et al. 1969), only a few hundred such events have been recorded across the pulsar population with the majority of them having a few glitches per pulsar¹. This low glitch rate, which characterizes the few multiple glitching pulsars, has hindered comprehensive statistical study of pulsar glitches. The prototype glitching pulsars (i.e. the Vela and Crab pulsars) have by far been the best studied. Thus, they have become a reference point in discussing pulsar glitch events. This is partially because they were first to be identified with such events, and mainly due to the characteristic size and frequency of the glitches they exhibit.

Glitch sizes are in the range of 10⁻¹¹ ≤ Δν/ν ≤ 10⁻⁵ (Espinoza et al. 2011), while the inter-glitch time intervals are in the range of 20 d ≤ t_i ≤ 10⁴ d. Most of the Vela pulsar glitches occupy the right end of the size distributions (with Δν/ν ≈ 10⁻⁶), while most of the Crab pulsar glitches occupy the left end (with Δν/ν ≈ 10⁻⁹). Across the population of glitching pulsars, most pulsars with multiple glitches present a broad range of glitch sizes. This has made it extremely difficult to use a single mechanism to explain all features of pulsar glitches across the population. Predicting glitch sizes from inter-glitch time intervals has been quite challenging, as most events appear to be random in size. Cumulative glitch sizes in the Vela pulsar are quite linear with time (McCulloch et al. 1987; Link et al. 1999; Marshall et al. 2004; Eya & Urama 2014; Eya et al. 2017), but this trend is not exhibited by glitches in the Crab and many other pulsars. However, evidence of linear dependence for cumulative glitch size with time exists in two other pulsars, namely, PSR J0357-6910 (Middleditch et al. 2006; Eya et al. 2017; Ferdman et al. 2018; Antonopoulou et al. 2018) and PSR J1420-6048 (Eya et al. 2017).

Other statistical analyses of glitch parameters in ensembles of pulsars have also generated interesting results. These include correlation of glitch activity with pulsar spin-down rate (McKenna & Lyne 1990; Urama & Okeke 1999; Espinoza et al. 2011) and spin-down luminosity (Fuentes et al. 2017), glitch size following a power law with respect to energy released in the event (Morley & Garcia-Pelayo 1993) and young pulsars having higher glitch rate than old pulsars (Espinoza et al. 2011). Melatos et al. (2008) analyzed glitch sizes and waiting times in individual pulsars using the cumulative distribution function (CDF) to show that glitch events could be reproduced in a glitch mechanism that involves avalanches of vortices. In the model, the basic assumption is that the neutron superfluid vortices are in a self-organized criticality (SOC) state. As such, the glitch size distribution is seen to follow a power law, while the distribution of waiting times is exponential. Following Melatos et al. (2008), Onuchukwu & Chukwude (2016) presented a similar result using a dataset of microglitches recorded by the Hartebeesthoek Radio Astronomy Observatory. Also, using the CDF plot, Eya et al. (2017) demonstrated that the cumulative spinup sizes in 12 pulsars follow a normal distribution. These pulsars are mostly pulsars in which the dispersions in their glitch size distribution are narrow.

In addition, the histograms of glitch size distribution have consistently been bimodal, (e.g. Wang et al. 2000; Espinoza et al. 2011; Yu et al. 2013; Eya & Urama 2014; Eya et al. 2017) thereby signaling a dual glitch mechanism. Nonetheless, the relatively large size glitches (LSGs) are mostly from a few middle aged pulsars. In fact, it was suggested that this bimodal nature of the distribution could result from two different classes of pulsars or via a glitch mechanism, which advances as the pulsars grow old. However, this bimodal distribution of glitch sizes, which has been widely reported, is in contrast to that of the spinup sizes (Δν), which tend towards a multimodal distribution (Espinoza et al. 2011). With amixture of Gaussians fitted in the histogram of spin-up sizes, Fuentes et al. (2017) also reported on it. Recently, from the estimation of glitch sizes and waiting time distributions, there is no persuasive evidence for such a multimodal distribution (Howitt et al. 2018). Using the distribution of fraction of neutron star components participating in a glitch, Eya et al. (2017) showed that the missing glitches that caused the dip in the histogram of glitch size distribution is accommodated well within the range that is feasible with the angular momentum transfer models. As such, the bimodal distribution of glitch sizes may not be connected to the glitch mechanism. These contrasting views may linger longer than expected due to paucity of glitch events, though it is expected that a clearer picture will emerge any time glitch data improve significantly, especially for individual pulsars. In all these analyses, efforts have not been made to study the distribution of inter-glitch time intervals in ensembles of pulsars or in analyzing the glitch data corresponding to the peaks in the histograms of glitch size distribution. Such a study could present a clearer picture to ascertain whether inter-glitch time intervals are size dependent or whether the peaks in the histograms of glitch size distributions could be linked to glitch mechanisms or not. In this paper, we shifted from the previous approach by examining the distributions of time intervals between glitches and analyzing the inter-glitch time intervals corresponding to small (Δν/ν < 10⁻⁷) and relatively large (Δν/ν ≥ 10⁻⁷) size glitches using an ensemble of 483 glitches in 168 pulsars. Equally, we investigate the distribution pattern of $\frac{\Delta \nu /\nu }{{t}_{i}}$ in both small and LSGs, where t_i is the inter-glitch time interval preceding a given glitch.

The glitches for this analysis were taken from the JBO Glitch Catalogue². This catalog contains 483 glitches in 168 pulsars, which are known at the time of this analysis. The glitch sizes are in the range of 10⁻¹¹ ≤ Δν/ν < 10⁻⁴ while the corresponding change in the spin-rate is in the range of ${10}^{-6}\le |\Delta \dot{\nu }/\dot{\nu }|\le {10}^{0}$.

We assume that the observed sudden increases in pulsar spin frequency during glitches result from transfer of angular momentum stored in the angular momentum reservoir over a glitch-free period preceding the glitch (hereafter called the inter-glitch time interval, t_i). With this, we evaluate the inter-glitch time interval for the n^th glitch as ${t}_{i}^{{\rm{th}}}={t}_{n}-{t}_{n-1}$, where t_n is the Modified Julian Date for the n^th glitch and t_n–1 is the Modified Julian Date for the (n – 1)^th glitch. The histograms of the distributions of glitch sizes, Δν/ν, and inter-glitch time intervals, t_i, are shown in Figure 1. Both distributions are continuous: the distribution of Δν/ν (top panel) is bimodal as has been widely reported (e.g. Wang et al. 2000; Espinoza et al. 2011; Yu et al. 2013; Eya & Urama 2014; Eya et al. 2017) with peaks at Δν/ν ≈ 2.14 × 10⁻⁹ and Δν/ν ≈ 1.17 × 10⁻⁶. On the other hand, the distribution of t_i is unimodal with peak at t_i ≈ 2.96 × 10³ d. The histogram of t_i approximates that of a normal distribution as suggested by the fit for it. The mean and median t_i are coincident in the bin before the modal bin. From the histograms (Fig. 1) there is no sign of larger t_i corresponding to larger Δν/ν or vice versa.

Zoom In Zoom Out Reset image size

In order to estimate a possible demarcation between relatively large and SSGs, we fitted twin Gaussians to the glitch size distribution. These Gaussians are fitted so that their modes are within the modal peak of each side of the glitch size distribution. Then the measure of central tendency for the Gaussians, the component of the glitch size distribution they demarcate and other shape parameters of the component are evaluated. This helps to determine the portion of the glitch size distribution, which has a similar characteristic with the Gaussian. In addition, our interest is in the fits in which the difference between the mode of the component it demarcates and the mode of the Gaussian is not more than $\pm \frac{1}{2}w$, where w is the binwidth of the histogram of the glitch size distribution³. The result is summarized in Table 1 and the distributions with the best fits are shown in Figure 2. From Figure 2, it is readily observed that the fit in the bottom panel gives a better demarcation line between LSGs and small ones than other panels. The mean of the Gaussians coincides with the mean of the components they demarcate. The standard deviation of the Gaussians and that of the component they demarcate are within the same range. Also the kurtosis⁴ of the components is closer to that of a normal distribution than others. The twin Gaussians intersect at approximately Δν/ν = 10⁻⁷. Consequently, glitches with Δν/ν < 10⁻⁷ are regarded as SSGs and LSGs are taken as those with Δν/ν ≥ 10⁻⁷.

Zoom In Zoom Out Reset image size

Table 1. A measure of the central tendency of the two peaks of the glitch size distribution with respect to the Gaussian fits in Fig. 2.

Statistic	Top panel		Middle panel		Bottom panel
	LP	RP	LP	RP	LP	RP
${\bar{\mu }}_{p}$	–9.40	–6.01	–9.06	–5.96	–8.54	–5.93
${\bar{\mu }}_{g}$	–8.79	–6.21	–8.54	–6.23	–8.54	–5.93
μp	–9.22	–5.73	–8.96	–5.95	–8.55	–5.93
μg	–8.79	–6.21	–8.54	–6.23	–8.54	–5.93
ρp	–8.67	–5.93	–8.67	–5.93	–8.67	–5.93
ρg	–8.79	–6.21	–8.54	–6.23	–8.54	–5.93
σp	1.08	1.10	1.27	1.12	1.32	0.55
σg	1.20	1.36	1.43	1.39	1.40	0.60
ku	2.10	2.30	2.60	2.25	2.62	2.87

Notes: LP (RP) denotes the component of the glitch size distribution for which the Left Peak (Right Peak) is assumed to be its mode. ${\bar{\mu }}_{p/g}$ represents the mean of the component/Gaussian, μ_p/g the median of the component/Gaussian, ρ_p/g the mode of the component/Gaussian, σ_p/g the standard deviation of the component/Gaussian and ku signifies the kurtosis of the component demarcated by the Gaussian.

In order to investigate whether the inter-glitch time intervals⁵ are actually size independent as envisaged in Figure 1, we constructed two classes of inter-glitch time intervals with respect to the two classes of glitches and investigated with the two dimensional Kolmogorov-Smirnov (K-S) test. We denote those corresponding to SSGs with _St_i, and those corresponding to LSGs with _Lt_i. The K-S test enables us to determine the similarity in the distributions of the classes of inter-glitch time intervals (i.e. _St_i and _Lt_i.)

In the K-S test, the null hypothesis (i.e. h = 0) is that in a given dataset, x₁ and x₂ described by CDFs F₁ (x) and F₂ (x) respectively are drawn from the same continuous distribution x. An alternative hypothesis (i.e. h = 1) is that the two datasets are drawn from different continuous distributions. The K-S statistic D is the maximum distance between F₁(x) and F₂(x). The probability that the result is by chance is ascertained by the magnitude of the P-value. The null hypothesis is rejected if the P-value is small, but a small P-value is strong evidence for accepting the alternative hypothesis⁶. In this paper, the K-S tests are performed at the 5% significance level. If the null hypothesis is true, it is an indication that the two distributions are drawn from a common continuous distribution, implying that inter-glitch time intervals are size independent. In addition, it could be an indication that a parameter involving time is scale invariant and plays the same role in the process that generated the two classes of glitch sizes.

The results obtained from application of the K-S test to the distributions of the two classes of inter-glitch time intervals are as follows: h = 0, D = 0.07 and P = 0.88. The CDF is shown in Figure 3. The outcome of this preliminary test indicates that there is no significant difference between the distributions of _St_i and _Lt_i. The maximum distance between the two curves is very small (Fig. 3) with a large P-value, enabling us to accept the null hypothesis. This is an indication that, on average, the inter-glitch time intervals are size independent with respect to glitch size. Furthermore, a parameter that is strongly tied to time is scale invariant and plays a similar role in the process that produces SSGs and LSGs.

Zoom In Zoom Out Reset image size

Next, following McKenna & Lyne (1990) who define the quantity, $\frac{\Sigma (\Delta \nu /\nu )}{{t}_{g}}$ as the 'glitch activity', where t_g is total time for the cumulative glitch size, Σ(Δν/ν), we define the quantity $\frac{\Delta \nu /\nu }{{t}_{i}}$ as the 'unit activity'. As the distribution of t_i in both classes of glitch size is not significantly different from each other, the distribution pattern of $\frac{\Delta \nu /\nu }{{t}_{i}}$ in both classes should be similar if both classes of glitch sizes originate from a single process, which is time invariant with respect to size. To investigate this, the distributions of unit activities in LSG and SSG were subjected to a normality test⁷.

In this analysis, we used the Lilliefors test at the 5% significance level for the normality test. Unlike the one dimensional K-S test, the Lilliefors test does not require predetermining the type of normal distribution, instead the type of normal distribution is determined from the dataset (Lilliefors 1967, 1969). A null hypothesis in the Lilliefors test (i.e. h = 0) is that the elements in the dataset are drawn from a normally distributed population or not.

Regarding application of the normality test to the unit activities, the results indicate that the distributions approximate a normal distribution in both classes. The P-values are 0.90 and 0.78 for LSG and SSG respectively. The CDF plots are shown in Figure 4. The fit of the CDF plots is based on a normality test result to indicate how the data fit the distribution. The goodnesses of the fits were quantified numerically using the two dimensional K-S test. The tests yielded P/D-values of 0.99/0.01 and 0.80/0.04, for LSG and SSG respectively, enabling us to accept the null hypothesis. The outcome of this test signifies that a parameter, which depends on time, plausibly the spin-down rate of a pulsar, plays a similar role in the processes that regulate both small and LSGs. This is because dividing the glitch sizes by their corresponding inter-glitch time intervals produces the same effect on the distribution pattern of both large and SSGs. Equally, it reaffirms that the inter-glitch time intervals are size independent.

Zoom In Zoom Out Reset image size

Pulsar glitch sizes and inter-glitch time intervals are key elements in discussing glitch events. From the glitch data, many pulsars exhibit a mixture of glitch sizes, though the majority of them has few glitches. SSGs are usually attributed to mechanisms involving a neutron starquake, while LSGs are attributed to a mechanism involving transfer of angular momentum. The dip in the distribution of glitch sizes (Fig. 1) has remained one of the strongest evidences in favor of dual glitch mechanisms (Yu et al. 2013). However, this idea has an implication as it could suggest that more than one glitch mechanism operates in a pulsar with a mixture of glitch sizes. If one regards the SSGs as events produced by a physical mechanism other than the angular momentum transfer, such as a starquake, it could be trivial to state that two different independent physical processes (i.e. starquake and angular momentum transfer processes) within a closed system have similar time intervals between them.

Meanwhile, in the frame of angular momentum transfer, a recent analysis has demonstrated that the fraction of neutron star components participating in a glitch could account for the current range of glitch sizes, including the missing glitches that caused the dip in the distribution of glitch sizes (Eya et al. 2017). As such, mechanisms involving transfer of angular momentum could readily account for all sizes of pulsar glitches. The glitch size has been shown to depend on three main factors: the number of vortices involved, the location of the vortices and the distance the vortices travel before re-pinning (Warszawski & Melatos 2011). Conventionally, these vortices are located in the inner-crust superfluid. The maximum distance the unpinned vortices could migrate before re-pinning is the thickness of the inner-crust (assuming a spherically symmetric inner-crust, and the vortices migrate radially outward after unpinning). In this regards, the distribution of the inter-glitch time intervals preceding LSG and SSG being the same could be feasible. This is because the electromagnetic braking torque on the stellar crust, which induces differential rotation that culminates in glitches, is not primarily connected to the vortices. If its magnitude is fairly stable at much longer times compared to the inter-glitch time intervals, it could lead to a similar inter-glitch time interval for both large and SSGs. The difference between LSG and SSG could be related to the number of vortices involved or the distance the vortices migrate before re-pinning. LSG could result from a situation in which almost all the vortices were unpinned; they migrate the maximum distance before re-pinning. On the other hand, SSG could result from a situation in which a few of the vortices were unpinned at each event. The vortices move less distance before re-pinning as their motion could be impeded by other pinned vortices. In such a scenario, a given pulsar could have a mixture of glitch sizes depending on how the vortex unpinning trigger mechanism operates in it.

Dividing the glitch sizes in each class by their corresponding inter-glitch time interval just normalizes the two distributions with respect to time (Fig. 4). This observation reaffirms the similarity of inter-glitch time intervals in two classes of glitches. If one is to present this picture in a model induced by a similar physical process, a mechanism that could produce varying glitch sizes with similar time intervals should be searched for. In addition, considering radiative changes associated with glitches in magnetars, and a few high magnetic field pulsars (Dib et al. 2008; Weltevrede et al. 2011; Dib & Kaspi 2014; Kaspi & Beloborodov 2017) and recently in a radio pulsar (Kou et al. 2018), such a mechanism should take into account magnetospheric activity. Magnetar glitches are part of the LSG distribution. Ordinarily with respect to a rotation-powered pulsar with a similar characteristic age, glitches in magnetars should be smaller compared to what is observed. Magnetars are mostly young pulsars powered by the decay of their conspicuous magnetic field rather than the loss of their rotational energy (Gao et al. 2013; Zhu et al. 2016; Liu et al. 2017); as such, there could be a contribution of magnetospheric activity to their glitch sizes (Gao et al. 2016).

Magnetospheric activity could come into play in glitches due to external torque on the pulsar coupled to it. The gross torque on a pulsar is from two sources, namely, the external electromagnetic braking torque and the internal braking torque coming from the neutron superfluid vortex creeping (Yuan et al. 2010). The external torque on the whole star causes small or micro-size glitch/timing noise and noticeable radiative changes, while the internal braking torque causes LSGs, SSGs and even timing noise (Yuan et al. 2010). As such, there could be external and internal origins for pulsar glitches depending on the contributions of torques involved, or even a combined action leading to LSGs accompanied with noticeable radiative changes.

Moreover, the narrowing of pulse profile and abnormal evolution of spin-down rate ($|\dot{\nu }|$) associated with glitch activity in the rotation-powered pulsar PSR J2037+3621 (B2035+36) could not be explained by standard glitch models (Kou et al. 2018). This is also an indicator that in some glitches, factors due to an external mechanism could be significant. In order to explain this anomaly that was observed in PSR J2037+3621, Kou et al. (2018) assume that there could be a change in external braking torque accompanied by the glitch activity. As the dipole magnetic field strength increases, the pulsar's braking index decreases due to an increasing braking torque and vice versa (Gao et al. 2017b,a). A change in braking torque affects the magnetic field structure, which in turn affects the inclination angle⁸. If this effect results from fluctuations in particle density outflow in the magnetosphere (Kou & Tong 2015), fluctuation in the inclination angle should lead to changes in magnetic field structure, further leading to radiative changes accompanying a glitch.

In conclusion, the distribution of inter-glitch time intervals having a single mode and that of _Lt_i and _St_i being the same are clear indications that a global parameter involving time, plausibly the pulsar spin-down rate $(\dot{\nu }$), plays the same role in both classes of glitches. Parameters that strongly regulate the glitch size are not strongly tied to the time parameter, rather, they have much to do with the glitch mechanism. Therefore, in searching for an explanation for pulsar glitch events, a hybrid mechanism that is related with a neutron star's interior, taking account of radiative change as well as having the ability to produce different sizes of glitches during similar time intervals, should be searched for.