Statistical Estimates of the Pulsar Glitch Activity

Abstract

A common way to calculate the glitch activity of a pulsar is an ordinary linear regression of the observed cumulative glitch history. This method however is likely to underestimate the errors on the activity, as it implicitly assumes a (long-term) linear dependence between glitch sizes and waiting times, as well as equal variance, i.e., homoscedasticity, in the fit residuals, both assumptions that are not well justified from pulsar data. In this paper, we review the extrapolation of the glitch activity parameter and explore two alternatives: the relaxation of the homoscedasticity hypothesis in the linear fit and the use of the bootstrap technique. We find a larger uncertainty in the activity with respect to that obtained by ordinary linear regression, especially for those objects in which it can be significantly affected by a single glitch. We discuss how this affects the theoretical upper bound on the moment of inertia associated with the region of a neutron star containing the superfluid reservoir of angular momentum released in a stationary sequence of glitches. We find that this upper bound is less tight if one considers the uncertainty on the activity estimated with the bootstrap method and allows for models in which the superfluid reservoir is entirely in the crust.

1. Introduction

To date, pulsar glitches are considered the most striking macroscopic manifestation of the presence of a neutron superfluid in the inner crust and outer core of neutron stars (see, e.g., [1] for a recent review). According to the current understanding, a rotating neutron star is comprised of two components, a normal one—strongly coupled to the magnetic field of the star and observed from Earth—and a superfluid one—which lags behind the normal one during the spin-down process [2,3]. Due to an unknown trigger, the two components can momentarily recouple (probably due to a mechanism known as vortex-mediated mutual friction [4,5]): the transfer of angular momentum from the neutron superfluid to the normal component results in a transient spin-up of the observable component, giving rise to a glitch.

Glitching behaviour can be very different from pulsar to pulsar [6,7,8], and its information can be encoded with a study of the glitch size and waiting time distributions (see, e.g., [9,10,11]); however, the identification of precise trends is difficult due to the scarcity of data for some objects: because of a combination of intrinsic physical properties and the different time spans of observations, some objects have presented only one glitch, while some others have displayed a statistically more relevant number of events (up to 45 in PSR J0537-6910 [12] and 23 in Vela [13]).

Some information about the structure of a glitching pulsar and the processes regulating this phenomenon can be obtained from its glitching behaviour, like the largest displayed glitch [14] or the short-time angular velocity evolution after a glitch [15,16,17,18]. Furthermore, it is possible also to obtain information about the neutron star structure (in particular, on the ratio between the moments of inertia of the normal and superfluid components [19,20,21,22,23]) from the observed activity parameter (i.e., the average spin-up due to glitches; see, e.g., [6]) of a pulsar. This parameter is particularly interesting as pulsar activity observations, together with theoretical modelling of the thermo-rotational evolution of a pulsar, have also been used to provide indirect mass estimates of isolated neutron stars [24,25,26].

Despite early models considered the superfluid neutrons in the core (see, e.g., [2,27]), after the seminal work of Anderson and Itoh [3], the superfluid participating in the glitch has been generally thought to be limited in the crust of the star, where vortex pinning is possible [28,29,30]. In fact, pulsar activity measured in the Vela pulsar seemed at first to be compatible with this idea of a crust-limited superfluid reservoir [19,20]. However, the introduction in the model of entrainment—a non-dissipative interaction that couples the two components [31]—and the calculation of this parameter in the neutron star crust [32] have posed serious issues to the modelling: entrainment coupling in the crust diminishes the effective angular momentum reservoir, making it difficult for a stellar model with a crust-limited reservoir to display a Vela-like activity.

Currently, to justify the observed activity of Vela, the neutron star crust must be sufficiently thick to store a significant amount of angular momentum, corresponding to a fraction of crust moment of inertia in a range going from 1.6% up to ∼10%, depending on the importance of the effect of crustal entrainment, which is currently under debate [33]. Alternatively, some models [17,24,26,34] also consider the possibility of a superfluid angular momentum that extends in the outer core: quantised vortex lines in the core superfluid could pin against the quantised flux lines of the proton superconductor ([35,36,37]; see also [38] for a review), so that it would be possible to store angular momentum in a region that is not just confined within the inner crust. However, there are some uncertainties about the nature of the proton superconductor [39], as well as the nature of neutron vortices in the outer core [40], so that pinning with flux tubes is quite uncertain to date.

Therefore, a reliable estimation of the glitch activity and its associated uncertainty is crucial to validate the crustal origin of pulsar glitches. The problem is particularly interesting for those pulsars that do not show a clear linear relation between the cumulative glitch size and the observational time. In this paper, we will deal with this problem, trying to find new ways to calculate glitch activity and in particular its uncertainty, stressing some subtleties regarding the latter value. Finally, we will employ the calculated activity parameter in a revised version of the original argument for the moment of inertia constraint found in [20,21,22].

2. Extracting the Activity Parameter from Observations

We consider a certain pulsar that has undergone $N_{\mathrm{gl}}$ glitches with size $\Delta {\Omega }_i$ (with $i=0,\dots ,N_{\mathrm{gl}}-1$) in a long observational time interval $T_{\mathrm{obs}}$. The absolute activity of a pulsar can be defined as:

$${\mathcal{A}}_a=\frac{1}{T_{\mathrm{obs}}}\sum _{i=0}^{N_{\mathrm{gl}}-1}\Delta {\Omega }_i.$$

(1)

Strictly speaking, this definition of ${\mathcal{A}}_a$ refers to a particular time window $T_{\mathrm{obs}}$. However, if the rate in (1) is almost constant when restricted to shorter time windows within $T_{\mathrm{obs}}$ (stationarity hypothesis), then a unique activity ${\mathcal{A}}_a$ can be defined for a long period, and (1) provides a reliable estimator for it (see, e.g., [41] for an analogous discussion on the mean seismic rate of earthquakes in a given time window). It is useful to introduce also the dimensionless activity parameter $\mathcal{G}$, defined as (see, e.g., [22]):

$$\mathcal{G}={|{\dot{\Omega }}_{\infty }|}^{-1}{\mathcal{A}}_a,$$

(2)

where $|{\dot{\Omega }}_{\infty }|$ is the absolute value of the secular spin-down rate of the pulsar, namely the average angular velocity derivative in the period $T_{\mathrm{obs}}$ containing several glitches. The variable $\mathcal{G}$ gives us an idea of the amount of spin-down reversed by glitches, and it allows for a better comparison of pulsars with different spin-down rates [26].

From the practical point of view, to use (1) with real data, the basic requirement is that the pulsar should regularly be monitored during the interval $T_{\mathrm{obs}}$, without missing any glitch. This is in general not the case, but while large glitches can be easily detected, very small glitches, easier to miss, should not contribute to the activity in a significant way (unless they are extremely frequent). Furthermore, it is not always clear whether the duration $T_{\mathrm{obs}}$ of the observational campaign has been long enough, so that ${\mathcal{A}}_a$ calculated via (1) really reflects the true activity of the pulsar under study. For this reason, we consider (1) as a theoretical definition of the true activity of a pulsar, in the limit of very long $T_{\mathrm{obs}}$. In the following, we will explore how to extract estimates of ${\mathcal{A}}_a$ from real data.

2.1. Ordinary Linear Regression on the Cumulative Glitch History

Let us assume that the information at our disposal consists only of a list of glitch dates $t_i$ and amplitudes $\Delta {\Omega }_i$. For simplicity, let us neglect the extra information that is possibly contained in the value of $T_{\mathrm{obs}}$, but note that $T_{\mathrm{obs}}>t_{N_{\mathrm{gl}}-1}-t_0$. Under this assumption, the absolute activity is usually calculated by fitting the cumulative glitch amplitude (see, e.g., [20,42]). In this case, the relationship between angular velocity and time is described by the equation:

$${\Omega }_i={\mathcal{A}}_a{t}_i+q+{\epsilon }_i,$$

(3)

where q is the vertical intercept and ${\epsilon }_i$ are independent random variables with zero expectation and the same variance (homoscedasticity). In the above formula, ${\Omega }_i$ and $t_i$ represent the angular velocity acquired by the star due to glitches and the time passed since the first glitch,

$${\Omega }_i=\sum _{j=1}^i\Delta {\Omega }_j{t}_i=\sum _{j=1}^i\Delta t_j,$$

(4)

where $\Delta t_j$ is the waiting time preceding the j-th glitch. This procedure sacrifices the information relative to the first glitch amplitude, $\Delta {\Omega }_0$, as the slope of the points in (4) does not change for vertical translations. One possibility to partially solve this issue is that of fitting the midpoints of the glitch steps drawn by the cumulative points $(t_i,{\Omega }_i)$, instead of the points themselves [26,42]. An example of the activity fit performed using this prescription is shown in Figure 1 for the six pulsars that have displayed the largest number of glitches $N_{gl}$ at the time of writing: the fit seems to capture the average slope of the glitch series, at least for J0537-6910 and Vela. All the glitch sizes and waiting times employed in this paper were retrieved from the Jodrell Bank Glitch Catalogue1 [8], while the spin-down rates of the stars from the ATNFPulsar Catalogue2 [43].

The central assumption behind this standard linear regression procedure is that the statistical properties of the processes underlying the glitch behaviour should not change if the window of observation is translated in time (i.e., the glitch series observed in a pulsar may be modelled as the outcome of a stationary stochastic process in the long run [44]). If this is the case, then the available data set should correspond to a stationary sequence of glitches where possible aftershocks and more quiet periods of activity are both present many times, intertwined in such a way to produce an overall stationary spin-up rate (eventually, also many very small and undetected glitches may be included, as their contribution to the cumulative glitch amplitude is negligible). In this way, the activity calculated on sub-intervals should fluctuate around the asymptotic value calculated by considering an interval $T_{\mathrm{obs}}$ containing several glitches [41].

For some pulsars, however, the observational window is so limited and the number of glitches detected so small that it may be unsafe to conclude that the observed value of ${\mathcal{A}}_a$ corresponds to the value that would be extrapolated by looking at a longer sequence of glitches. Practically, to be able to perform a standard linear regression on the cumulative data, we have to demand that the available data fulfil two practical requirements [26]. Firstly, the number of glitches should be significant, say $N_{\mathrm{gl}}>3$. Secondly, at least two glitches of size comparable to $\Delta {\Omega }_{\max }$, the maximum glitch amplitude in the sequence, should be present: a linear regression would poorly fit the set of data if the largest glitch ${\Omega }_{\max }$ is significantly larger than all the others. This last property can be quantified by demanding that the parameter $N_{\max }$, defined as [25,26]:

$$N_{\max }=\Delta {\Omega }_{\max }^{-1}\sum _{i=0}^{N_{\mathrm{gl}}-1}\Delta {\Omega }_i>1\mathrm{where}\Delta {\Omega }_{\max }=\underset{i=0,\dots ,N_{\mathrm{gl}}-1}{max}\Delta {\Omega }_i,$$

(5)

be larger than ∼2 [26]. However, in view of the hypothesis of stationarity, even if the glitches respect the above conditions, the ordinary least-squares linear regression may not be a well justified method, at least from the theoretical point of view. In fact, the technical point of how to extract ${\mathcal{A}}_a$ and the associated uncertainty, especially in the case of small $N_{\mathrm{gl}}$ or $N_{\max }$, have not been discussed in detail yet.

2.2. Linear Regression on Heteroscedastic Data

Besides all the issues related to linear regression mentioned in the previous section, there is one more linked to the fact that the data points in (4) are not independent, as they arise from a cumulative construction. Hence, the resulting residuals (which have to be minimised in the standard regression procedure) may not be independent and identically distributed: the fit residuals will not have the same variance, so it is no longer possible to assume homoscedasticity, which is a basic assumption of ordinary linear regression.

For this reason, we present here an alternative way to calculate the activity, but this time relaxing the hypothesis of homoscedasticity, by following a procedure for fitting cumulative data discussed by Mandel [45]. Let us assume that waiting times and glitch sizes are related by:

$$\Delta {\Omega }_i={\mathcal{A}}_a\Delta t_i+{\epsilon }_i.$$

(6)

Note that the above equation may not be true in general, as there seems not to be a correlation between glitch sizes and waiting times [46]. Using Equation (4), the cumulative times and sizes follow a relation:

$${\Omega }_i={\mathcal{A}}_a{t}_i+\sum _{j=1}^i{\epsilon }_j,$$

(7)

which justifies our assumption of heteroscedasticity, as the deviations ${\sum }_{j=1}^i{\epsilon }_j$ have different variance for each i. It can be shown that the best unbiased linear estimator is given by [45]:

$${\mathcal{A}}_a=\frac{{\sum }_i\Delta {\Omega }_i}{{\sum }_i\Delta t_i}=\frac{{\sum }_i\Delta {\Omega }_i}{t_{N_{\mathrm{gl}}-1}-t_0}i=1,...,N_{\mathrm{gl}}-1$$

(8)

which is similar to the definition in (1). The variance of this value is given by [45]:

$$\mathrm{Var}\left({\mathcal{A}}_a\right)=\frac{1}{(N_{\mathrm{gl}}-1)(t_{N_{\mathrm{gl}}-1}-t_0)}\sum _i\frac{{(\Delta {\Omega }_i-{\mathcal{A}}_a\Delta t_i)}^2}{\Delta t_i}i=1,...,N_{\mathrm{gl}}-1$$

(9)

We present in

Table 1. Dimensionless activities and their standard deviations, calculated for the six pulsars with the largest number of glitches, with a least-squares linear fit on the cumulative midpoints assuming homoscedasticity (${\mathcal{G}}_{\mathrm{hom}}$), with a least-squares linear fit assuming heteroscedasticity (${\mathcal{G}}_{\mathrm{het}}$), with a bootstrap on the size and waiting time samples separately (${\mathcal{G}}_{\mathrm{rand}}$), and on the pairs’ size, preceding waiting time (${\mathcal{G}}_{\mathrm{pre}}$), and size, following waiting time (${\mathcal{G}}_{\mathrm{post}}$).

Pulsar	${\mathcal{G}}_{\mathbf{hom}}(\%)$	${\mathcal{G}}_{\mathbf{het}}(\%)$	${\mathit{G}}_{\mathbf{rand}}(\%)$	${\mathit{G}}_{\mathbf{pre}}(\%)$	${\mathcal{G}}_{\mathbf{post}}(\%)$
0534+2200	0.0079 ± 0.0007	0.008 ± 0.006	0.008 ± 0.005	0.008 ± 0.005	0.008 ± 0.005
0537-6910	0.874 ± 0.003	0.85 ± 0.15	0.89 ± 0.11	0.86 ± 0.11	0.88 ± 0.03
0631+1036	1.77 ± 0.18	2.03 ± 1.95	2.11 ± 1.67	2.29 ± 1.80	1.80 ± 0.85
0835-4510	1.62 ± 0.02	1.6 ± 0.2	1.65 ± 0.3	1.6 ± 0.2	1.6 ± 0.2
1341-6220	1.52 ± 0.10	1.9 ± 0.6	2.0 ± 0.6	1.9 ± 0.6	1.9 ± 0.5
1740-3015	1.22 ± 0.04	1.3 ± 0.7	1.3 ± 0.5	1.3 ± 0.5	1.2 ± 0.45

the best estimator of the dimensionless activity $\mathcal{G}$, with its standard deviation obtained with this method. This value was obtained by neglecting the first glitch size $\Delta {\Omega }_0$ in each pulsar sequence, in order to have the same number of glitch sizes and preceding waiting times. The most interesting feature we notice from the results is that the standard deviation when we assume homoscedasticity is almost one order of magnitude smaller than that with heteroscedasticity.

2.3. Extracting the Activity from Glitch Size and Waiting Time Distributions

One alternative way to solve the issue in the calculation of the activity with the linear fit requires employing the probability distributions for the waiting times and sizes of the glitches of a particular pulsar. In this way, it is possible to relax the hypothesis of linear dependence between waiting times and glitch sizes.

Probability distributions of glitch sizes and waiting times have been obtained in several previous works [9,10,11], which show that—as a general trend—glitch sizes seem to be consistent with a power-law distribution, while the waiting times are consistent with an exponential distribution. Vela and PSR J0537-6910 are somehow exceptional, as they seem be well described by a normal distribution in both size and waiting time.

Starting from the probability distribution of the waiting times $\Delta t$ and sizes $\Delta \Omega$, it is possible to infer some information about the probability distribution $P_{{\mathcal{A}}_N}$ for the activity parameter ${\mathcal{A}}_N$ after N glitches. Note that ${\mathcal{A}}_N$ and ${\mathcal{A}}_M$, for $N\ne M$, are different random variables, distributed according to different laws, i.e., $P_{{\mathcal{A}}_N}\ne P_{{\mathcal{A}}_M}$.

Given the definition of activity, its distribution can be obtained by considering the ratio of two random variables, the sum of sizes $\Delta {\tilde{\Omega }}_N$ and of waiting times $\Delta {\tilde{t}}_N$. The latter, in turn, are the sum of random variables themselves, i.e., the single glitch size $\Delta {\Omega }_i$ and the single waiting time $\Delta t_i$, so that their densities $P_{\Delta {\tilde{\Omega }}_N}$ and $P_{\Delta {\tilde{t}}_N}$ can be obtained by means of repeated convolutions,

$$\begin{array}{cc}\hfill \Delta {\tilde{\Omega }}_N=\sum _{i=1}^N\Delta {\Omega }_i\Rightarrow & \Delta {\tilde{\Omega }}_N\sim P_{\Delta {\tilde{\Omega }}_N}=\underset{N\mathrm{times}}{\underbrace{P_{\Delta \Omega }\ast \dots \ast P_{\Delta \Omega }}}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \Delta {\tilde{t}}_N=\sum _{i=1}^N\Delta t_i\Rightarrow & \Delta {\tilde{t}}_N\sim P_{\Delta {\tilde{t}}_N}=\underset{N\mathrm{times}}{\underbrace{P_{\Delta t}\ast \dots \ast P_{\Delta t}}}\hfill \end{array}$$

(11)

Since ${\mathcal{A}}_N=\Delta {\tilde{\Omega }}_N/\Delta {\tilde{t}}_N$, its distribution can be obtained through:

$$P_{{\mathcal{A}}_N}\left(a\right)={\int }_{-\infty }^{\infty }\mathrm{d}x\left|x\right|P_{\Delta {\tilde{t}}_N}\left(x\right)P_{\Delta {\tilde{\Omega }}_N}\left(xa\right).$$

(12)

The main advantage of this method is that, with only the assumption of independence between sizes and waiting times in the stationary regime, it is possible to obtain the whole probability distribution of ${\mathcal{A}}_N$.

Although this is a possible method to extract ${\mathcal{A}}_N$, it is troublesome to numerically obtain this distribution, as a convolution of N probability distributions, albeit identical, starts to be infeasible when N becomes large. The calculation can be simplified by obtaining the first two moments of $P_{{\mathcal{A}}_N}$, the mean and the variance, by employing a generalised version of the central limit theorem, the delta method (see Appendix A). This method allows us to calculate the mean and the variance of $P_{{\mathcal{A}}_N}$ starting from the mean and variance of the two distributions $P_{\Delta \Omega }$ and $P_{\Delta t}$. This procedure, however, is problematic if $P_{\Delta t}$ or $P_{\Delta \Omega }$ have a non-well-defined variance.

2.4. Estimating the Uncertainty of Activity: The Bootstrap Method

An attractive alternative to the methods described above, which does not build on the linear assumption in (6), is the so-called “bootstrap method” [47]. The idea is that of resampling with replacement the original data in order to calculate some statistics, as, e.g., the mean and standard deviation of the calculated activity. In our case, the samples are two: the list of the waiting times $\Delta t_i$ (of length $N_{\mathrm{gl}}-1$) and the list of sizes $\Delta {\Omega }_i$ (of length $N_{\mathrm{gl}}$). Of course, we have to draw the same number ($N_{\mathrm{gl}}-1$) of waiting time-size couples in order to have a fair estimation of the activity and its standard deviation. To avoid the homoscedasticity problem, the pulsar activity is calculated by employing the definition in Equation (8) on each set of resampled data. We can also take into account the possibility of a dependency between a glitch size and the preceding or the subsequent waiting time, so it is useful to also bootstrap on the other two samples: the sample made up by ordered pairs ${\left\{(\Delta {\Omega }_i,\Delta t_i^{\mathrm{pre}})\right\}}_{i=1,\dots ,N_{\mathrm{gl}}-1}$, where $\Delta t_i^{\mathrm{pre}}$ is the waiting time preceding the glitch of size $\Delta {\Omega }_i$, and the sample comprised by ordered pairs ${\left\{(\Delta {\Omega }_i,\Delta t_i^{\mathrm{post}})\right\}}_{i=0,\dots ,N_{\mathrm{gl}}-2}$, where $\Delta t_i^{\mathrm{post}}$ is the waiting time following a glitch of size $\Delta {\Omega }_i$.

Figure 2 shows the histograms obtained by resampling the data ${10}^4$ times in all three cases described above. In the same plot, also the dimensionless activities obtained as a result of the ordinary linear regression on the cumulative glitch data are displayed.

We can see that the activity calculated by means of bootstrapping is compatible with the results obtained from an ordinary linear regression, but it generally has larger standard deviations (see also Table 1). It is interesting to notice the case of PSR J0537-6910, one of the few stars that presents a significant correlation between the glitch size and the following waiting time [48]. This correlation shows its effects also in Figure 2: the histogram for J0537-6910, in the particular case of the sample of size-following waiting time pairs, is much more peaked than the other two cases. At lower confidence, also PSR J0631+1036 shows a correlation between size and the following waiting time and Vela a correlation between size and the preceding waiting time [46]. These correlations show their effect in the histograms as well.

It is also interesting to notice the peculiar form of the PSR J0631+1036 activity distribution: it shows two clear peaks, one on very small values and one around $\mathcal{G}\approx 0.02$. This is probably because of the particular glitch sequence of this star (see Figure 1): it displays two very large glitches and many others with sizes several orders of magnitude smaller. Thus, it is likely that the peak on smaller values was generated by sampling the small glitches only, while the peak on larger values occurs when one or both large glitches were sampled. Moreover, as a consequence of the particular glitch sequence for this star, the value of $N_{\max }$ for this star is smaller than the ones of the other stars in the sample: this value may increase in the future, by observing large glitches, giving us a more reliable estimate of the activity [26].

In Figure 3, we try to give an idea of how much the activity changes when a new glitch occurs. The first point of each curve is the activity calculated with the linear fit on the cumulative midpoints using the first ten glitches, assuming homoscedastic data. Then, we update the activity value with the same method whenever a new glitch is displayed. We also plot the activity parameter calculated with the ordinary linear regression and all the glitches, along with its uncertainty. For PSR J0537-6910 and PSR J0631+1036, we present the bootstrap estimate ${\mathcal{G}}_{\mathrm{post}}$ and its standard deviation; for Vela, we show the same with ${\mathcal{G}}_{\mathrm{pre}}$; while for all the other pulsars, we present the case with uncorrelated glitch sizes and waiting times (${\mathcal{G}}_{\mathrm{rand}}$). The idea is that of taking into account the size-waiting time correlations, where present. We note that, except the very particular cases of Vela and PRS J0537-6910, the activity evolution of each star generally lies outside the error region for the linear fit with homoscedastic data for all the glitch history, except—of course—the latest glitch. The bootstrap uncertainty better describes the variance of the glitch history. A notable exception is that of PSR J1341-6220, which is well below the error bar for both the activity calculations, except for the last three glitches. This is because these glitches are three of the largest ones displayed by this pulsar (see also Figure 1). In general, however, it is interesting to notice how variable the activity parameter is. A single large glitch can change its value (see, e.g., the Crab pulsar after its November 2017 glitch [49]). This fact stresses the importance of having a much larger uncertainty on the activity, which is the result of not assuming homoscedasticity in the linear fit or linear dependence between glitch sizes and waiting times.

3. Moment of Inertia Constraint

The activity parameter allows extracting information on the moment of inertia fraction of the superfluid reservoir in a glitching pulsar [19,20]. In this section, we present a revised version of the constraint on the moment of inertia relative to the superfluid component derived by Link et al. [20]; see also Andersson et al. [21] and Chamel [22] for the inclusion of entrainment coupling between the normal and superfluid components. Our derivation is presented in Appendix B and takes into account the non-rigid rotation of the superfluid component, the stellar stratification and the non-uniform entrainment coupling between the components. The constraint is given by Equation (A16), which can be written in terms of the total moment of inertia I and of the moment of inertia of the superfluid component $I_v$ as:

$$\frac{I_v}{I-I_v}>\mathcal{G}.$$

(13)

Assuming that the superfluid reservoir extends in the layers between $R_c$ (the crust-core transition radius) and $R_d$ (the neutron drip radius), the value of $I_v$ in the slow-rotation approximation is:

$$I_v=\frac{8\pi }{3}{\int }_{R_c}^{R_d}dr{r}^4{e}^{\Lambda \left(r\right)-\Phi \left(r\right)}\left[\rho \left(r\right)+P\left(r\right)\right]\frac{y_n\left(r\right)}{1-{\epsilon }_n\left(r\right)}\frac{{\Omega }_p-\omega \left(r\right)}{{\Omega }_p},$$

(14)

while the total moment of inertia I is the usual one in the slow rotation framework [14,50,51], given in Equation (A19). In the above equation, $y_n$ is the superfluid neutron baryon density (limited to the region where pinning is possible) divided by the total baryon density, i.e., $y_n\left(r\right)$ is different from zero only where the superfluid can pin to inhomogeneities and maintain its state of motion while the normal component spins down: in this case, it is limited to the crust. The other quantities appearing in (14) are introduced in Appendix B, but it is important to remark here that the frame drag $\omega$ should contain a dependence on the angular velocities of both components [14,51,52], which was neglected in the derivation of (14).

The first measurements of the activity parameter of the Vela pulsar and the moment of inertia fraction estimates for different equations of state (EoSs) seemed to be in accordance with the constraint (13) with ${\epsilon }_n=0$ [20]. Only later, the entrainment parameter ${\epsilon }_n$ in the crust of a neutron star was calculated in [32], by estimating the effects of Bragg scattering on the superflow due to the presence of the crustal lattice. These calculations yield a negative entrainment parameter ${\epsilon }_n\sim -10$ in a substantial portion of the inner crust, which implies a severely hindered motion of the superfluid component. This would reduce the amount of extra angular momentum stored in the crust between two glitches—and thus of $I_v$—making the requirement in (13) more difficult to be met. As a result of that, the only way for the star to acquire enough angular momentum between glitches to explain the observed activity is to have a large region inside it to store angular momentum (larger than the crust of the star) or to have an unreasonably small mass, around $\sim 0.5M_{\odot }$. This problem has been highlighted in several papers [21,22,23,24,53].

If we assume the superfluid region to be limited in the crust of the star and we fix the microphysical parameters, namely the EoS and the entrainment parameter, then the moment of inertia fraction in (13) is a function of the mass of the star only.

In Figure 4, we plot the quantity $I_v/(I-I_v)$ appearing in (13) as a function of the stellar mass for some EoSs and for the entrainment parameter calculated by Chamel [32], assuming a superfluid reservoir limited to the crust.

As we can see, the constraint (13) imposes that the high $\mathcal{G}$ value of the Vela pulsar (PSR J0835-4510) can be explained only if the mass of Vela is small, ranging from $\approx 1.1M_{\odot }$ for the BSk20EoS to $\approx 0.8M_{\odot }$ for the SLy4EoS. Let us however consider also the $1\sigma$ uncertainty region, calculated with the bootstrap method described in Section 2, with random waiting times and glitch sizes. In this case for the BSk20 EoS, the Vela pulsar has an upper limit on the mass of about $1.2M_{\odot }$. Note that this value is slightly above the minimum mass of a neutron star estimated from the calculations of core-collapse supernovae (i.e., $1.17M_{\odot }$; see [58]) and the smallest mass measured in a neutron star ($1.174\pm 0.004M_{\odot }$, measured in PSR J0453+1559, [59]). If we consider the $3\sigma$ uncertainty range, then we obtain a limit of $\approx 1.5M_{\odot }$ for the same EoS. It is thus clear that a careful estimation of the activity parameter and the associated error is crucial if one is to set strict quantitative constraints on the moment of inertia involved in glitches and constrain the location of the superfluid reservoir.

A more careful estimate of the glitch activity parameter may thus play an important role in resolving the tension between strong entrainment and models with a crust-limited reservoir. Several other effects have, of course, been suggested and are likely to play a role in this problem, including a maximally stiff EoS [60], a Bayesian analysis of the EoS uncertainty [53] or an extension of the region where the neutron superfluid participates in the glitch beyond the crust-core transition, based on the assumption that only the superfluid in the ${}^1{S}_0$ state participates in the glitch phenomenon and on an analysis of the temperature of the star [24]. On the other hand, also different calculations of the entrainment parameter have been proposed (e.g., [33,61,62]), which yield milder entrainment effects in the crust. However, even if the high activity of Vela may be described in terms of a purely crustal reservoir by assuming a weaker entrainment, an analysis of the 2016 Vela glitch points to the need to invoke the neutron superfluid also in the core of the star [16,63,64], independently of the presence of entrainment in the model [18].

4. Conclusions

In this paper, we discussed different ways of calculating the activity parameter in glitching pulsars. The most commonly used one is that of performing an ordinary linear regression on the cumulative glitch data, which is justified if one considers the activity parameter as an intrinsic characteristic of a glitching pulsar, and thus inferable from a limited observation of its glitching behaviour (which, in general, may not be stationary). While this last statement might be true, it is also true that a strong autocorrelation or correlation between glitch sizes and waiting times has very rarely been observed [46,65]. Moreover, fitting the cumulative data may also affect the uncertainty calculated from a ordinary linear regression, as the hypothesis of the homoscedasticity of the data cannot be satisfied. Therefore, the main consequence of including these assumptions in the fitting procedure (namely, the linear dependence of glitch sizes and waiting times and the homoscedasticity of the cumulative data) is an underestimation of the activity uncertainty. This is an important point, as in some pulsars, an additional single glitch (or a sequence of a few glitches) can significantly affect its observed activity, especially for those objects with low $N_{\max }$ [26], like the Crab pulsar or PSR J0631+1036 (see Figure 3).

A first alternative to calculate the activity is the linear regression developed by Mandel [45], where the hypothesis of the homoscedasticity of the data was relaxed: this methodology is justified by the fact that cumulative data are not homoscedastic, since the variance of the data should increase as data are cumulated. As a result, the uncertainty of the activity parameter increases by about one order of magnitude with respect to that calculated with an ordinary linear regression.

As a second step, we also relaxed the hypothesis of linear dependence between glitch sizes and waiting times, by employing their probability distribution estimates. While employing the size and waiting time distributions is arguably the best way to include the information provided by observations of the pulsar in the activity parameter, it is also true that a study of the probability distribution of the activity parameter is computationally very challenging. Furthermore, an approximate version of this estimate, obtained by employing a more general version of the central limit theorem, leads to difficulties, due to the fat-tailed probability distribution of the glitch sizes.

We thus employed an alternative way to calculate the activity and its uncertainty, relaxing the hypothesis of linear dependence, by bootstrapping on the sizes and waiting time data sets and calculating the activity by summing sizes and waiting times and calculating the ratio. Much larger uncertainties are obtained, on the same order of magnitude of the linear fit on heteroscedastic data. A qualitative explanation of why a bootstrap estimate yields a larger uncertainty than the ordinary linear regression is the fact that the glitch activity is dominated by a few large events, especially for pulsars with low $N_{\max }/N_{\mathrm{gl}}\lesssim 0.1$, which is again the case for the Crab pulsar and J0631+1036; see Table 1.

We then used these results to study a revised version of the constraint on the moment of inertia of the superfluid component and on the mass of the star [20,21,22]. While the result obtained in this way is essentially the same for the mean value of the activity (i.e., the Vela activity does not allow a crust-limited reservoir and strong entrainment), the larger uncertainty allows for models in which the crust is enough, i.e., models in which the superfluid reservoir is located entirely in the crust for credible values of the star mass. For example, a model of a $1.2M_{\odot }$ neutron star described by the BSk20 EoS is within the $1\sigma$ uncertainty region of the activity that we calculate.