Statistical estimates of the pulsar glitch activity

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 assumption that are not well justified in 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. Our main finding is a much larger uncertainty on activity estimates, with respect to that obtained with an ordinary linear regression. We discuss how this affects the theoretical upper bound on the moment of inertia associated to 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.


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 on the glitch size and waiting time distributions (see, e.g., [9][10][11]), but the identification of precise trends is difficult due to the scarcity of data in some objects: because of a combination of intrinsic physical properties and different time span 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 the 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] has 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 the 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 in the nature of the proton superconductor [39], as well as on 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 which 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].

Extracting the activity parameter from observations
We consider a certain pulsar that has undergone N gl glitches with size ∆Ω i (with i = 0, . . . , N gl − 1) in a long observational time interval T obs . The absolute activity of a pulsar can be defined as Strictly speaking, this definition of A a refers to a particular time window T obs . However, if the rate in (1) is almost constant when restricted to shorter time windows within T obs (stationarity hypothesis), then a unique activity 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 G, defined as (see e.g. where |Ω ∞ | is the absolute value of the secular spin down rate of the pulsar, namely the average angular velocity derivative in the period T obs containing several glitches. The variable 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 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 obs of the observational campaign has been long enough, so that 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 obs . In the following we will explore how to extract estimates of A a from real data.

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 ∆Ω i . For simplicity, let us neglect the extra information that is possibly contained in the value of T obs , but note that T obs > t N 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 where q is the vertical intercept, and ε i are independent random variables with zero expectation and same variance (homoscedasticity). In the above formula, Ω i and t i represent the angular velocity acquired by the star due to glitches and the time passed since the first glitch, where ∆t j is the waiting time preceding the j-th glitch. This procedure sacrifices the information relative to the first glitch amplitude, ∆Ω 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 , Ω 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 which 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 the Vela. 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 on the long run [43]). 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 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 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   (4), for the six pulsars with the largest number of glitches N gl , as reported by the Jodrell Bank Glitch Catalogue [8]. The activity is calculated with a least-squares linear fit on the midpoints of the cumulative glitch sequence [see 42], giving the blue curve with the associated uncertainty on the slope. Following Montoli et al. [26], the plots are made by using the "nominal lag" t|Ω ∞ | on the horizontal axis (a rescaling of time t which gives a rough estimate of the typical angular velocity lag accumulated between the spinning down normal component and a pinned superfluid component). In this way, the slope of the blue curves is the dimensionless activity G.
standard linear regression on the cumulative data we have to demand that the available data fulfill two practical requirements [26]. Firstly, the number of glitches should be significant, say N gl > 3. Secondly, at least two glitches of size comparable to ∆Ω 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 Ω 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] is 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 A a and the associated uncertainty, especially in the case of small N gl or N max , has not been discussed in detail yet.

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 more 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 [44]. Let us assume that waiting times and glitch sizes are related by Note that the above equation may not be true in general, as there seems not to be correlation between glitch sizes and waiting times [45]. Using equation (4), the cumulative times and sizes follow a relation which justifies our assumption of heteroscedasticity, as the deviations ∑ i j=1 ε j have different variance for each i. It can be shown that the best unbiased linear estimator is given by [44] which is similar to the definition in (1). The variance on this value is given by [44] Var(A a ) = 1 We present in Table 1 the best estimator of the dimensionless activity G, with its standard deviation obtained with this method. This value has been obtained by neglecting the first glitch size ∆Ω 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.

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 to employ 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. The Vela and PSR J0537-6910, are somehow exceptional, as they seem be well described by normal distribution in both sizes and waiting times.
Starting from the probability distribution of the waiting times ∆t and sizes ∆Ω it is possible to infer some information about the probability distribution P A N for the activity parameter A N after N glitches. Note that A N and A M , for N = M, are different random variables, distributed according to different laws, i.e.
Given the definition of activity, its distribution can be obtained by considering the ratio of two random variables, the sum of sizes ∆Ω N and of waiting times ∆t N . The latter, in turn, are sum of random variables themselves, i.e. the single glitch size ∆Ω i and the single waiting time ∆t i , so that their densities P ∆Ω N and P ∆t N can be obtained by means of repeated convolutions, Since A N = ∆Ω N /∆t N , its distribution can be obtained through: The main advantage of this method is that, with the only assumption of independence between sizes and waiting times in the stationary regime, it is possible to obtain the whole probability distribution of A N . Although this is a possible method to extract 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 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 A N starting from the mean and variance of the two distributions P ∆Ω and P ∆t . This procedure, however, is problematic if P ∆t or P ∆Ω have a non-well-defined variance.

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" [46]. 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 ∆t i (of length N gl − 1) 1.52 ± 0.10 1.9 ± 0.6 2.0 ± 0.6 1.9 ± 0.6 1.9 ± 0.5 1740-3015 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 (G hom ), with a least-squares linear fit assuming heteroscedasticity (G het ), with a bootstrap on the size and waiting time samples separately (G rand ), and on the pairs size -preceding waiting time (G pre ) and sizefollowing waiting time (G post ). and the list of sizes ∆Ω i (of length N gl ). Of course, we have to draw the same number (N 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 other two samples: the sample made up by ordered pairs {(∆Ω i , ∆t where ∆t post i is the waiting time following a glitch of size ∆Ω i . Figure 2 shows the histograms obtained by resampling the data 10 4 times in all the 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 which presents a significant correlation between a glitch size and the following waiting time [47]. 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 [45]. 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 G ≈ 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 has been generated by sampling the small glitches only, while the peak on larger values occurs when one or both the large glitches have been 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 G post and its standard deviation, for Vela we show the same   with G pre , while for all the other pulsars we present the case with uncorrelated glitch sizes and waiting times (G 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 homeoscedastic 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, [48]). 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.

Moment of inertia constraint
The activity parameter allows to extract 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 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 while the total moment of inertia I is the usual one in the slow rotation framework [14,49,50], 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 (r) 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 has been limited to the crust. The other quantities appearing in (14) are introduced in Appendix B, but it's important to remark here that the frame drag ω should contain a dependence on the angular velocities of both components [14,50,51], which has been 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 ε n = 0 [20]. Only later, the entrainment parameter ε n in the crust of a neutron star has been 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 ε n ∼ −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   [53], BSk20 and BSk21 [54], and the DDME2 EoS [55], glued with a SLy4 crust [following the method described in 56]. The entrainment parameter is that calculated in [32], and the superfluid reservoir is limited in the crust of the star. The dimensionless activity parameter G -calculated with the bootstrap method described in Section 2 (the case with random glitch sizes and waiting times) -is also plotted for the Vela pulsar, along with the 1σ, 2σ and 3σ uncertainties. 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 ∼ 0.5 M . This problem has been highlighted in several papers [21][22][23][24]52].
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 G value of Vela pulsar (PSR J0835-4510) can be explained only if the mass of the Vela is small, ranging from ≈ 1.1M of the BSk20 EoS to ≈ 0.8M for the SLy4 EoS. Let us however consider also the 1σ 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 . Note that this value is slightly above the minimum mass of a neutron star estimated from calculations of core-collapse supernovae (i.e. 1.17 M , see [57]) and the smallest mass measured in a neutron star (1.174 ± 0.004 M , measured in PSR J0453+1559, [58]). If we consider the 3σ uncertainty range, then we obtain a limit of ≈ 1.5M 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.

Conclusions
In this paper, we have discussed different ways of calculating the activity parameter in glitching pulsars. The most commonly used one is that of performing an ordinary linear regression to 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 [45,59]. Moreover, fitting the cumulative data may also affect the uncertainty calculated from a ordinary linear regression, as the hypothesis of 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 it can be shown that an additional single glitch or a sequence of a few glitches in a pulsar, can seriously affect its observed activity.
A first alternative way to calculate the activity includes the study of a linear regression where the hypothesis of homoscedasticity of the data has been relaxed. This relaxation is justified by the fact that cumulative data is not homoscedastic, as the variance of the data increases as data is cumulated. The main consequence of this assumption is that 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. Also 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. 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 neutron star described by the BSk20 EoS is within the 1σ uncertainty region of the activity that we calculate.
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 [52] 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 the 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].

Appendix A. Activity calculation with the delta method
We use the delta-method (a generalisation of the central limit theorem) to extract the mean and standard deviation of the activity parameter A N after N glitches, given a probability distribution of the sizes P ∆Ω and waiting times P ∆t of a pulsar. The method works under the simplifying assumption that the sizes and the waiting times are independent and identically distributed random variables. Let us recall the two random variables ∆Ω N and ∆t N defined in equations (10) and (11). The random variable A N is defined as Its expectation value is given by (a ∈ R + denotes the value assumed by A N ) Independence of the variables is loosely justified by observing little correlation and autocorrelation in glitch sizes and waiting times [45,59], while being identically distributed is a working assumption. The above equation boils down to To calculate the variance of A N , let us first calculate, for the sizes ∆Ω N : Now, using the above results, a simple direct calculation gives Let us now assume that P ∆t N converges to a normal distribution with mean Nθ and variance Nσ 2 as N → ∞. In this case, we can employ the delta method, which tells us that any function g of the random variable ∆t N will be distributed normally with mean g(Nθ) with variance N[g (Nθ) σ] 2 . Note that this assumption is not true if glitch sizes are described by a power law distribution with non defined variance (but should hold for those pulsars like PSR J0537-6910). Given the definition of ∆t N , we have that θ = E[∆t] and σ 2 = Var[∆t], by considering the cases g(∆t N ) = (∆t N ) −1 and g(∆t N ) = (∆t N ) −2 we obtain:

Appendix B. Derivation of the moment of inertia constraint
Let us write the dynamics of the total angular momentum L of the pulsar (neglecting temporal variations in the total moment of inertia I) as where we assumed that the normal component is rigidly rotating with angular velocity Ω p , while ∆L is the angular momentum reservoir due to the (non-rigid) angular velocity lag Ω np between the superfluid and the normal component [14]. Clearly, (A11) is valid only if we are assuming that the rigid component and the fluid one share a common rotation axis. It is useful to formally divide Ω p and ∆L into the contributions due to the smooth relaxation (R) and an impulsive one from glitches (G), During glitches we have∆L G < 0 andΩ G p > 0, while for the rest of the time∆L R > 0 andΩ R p < 0. We can average (A11) over a long time interval T obs , to get We can simplify the equation above by making two observations. Firstly, due to the angular momentum conservation during a glitch we must have that (note that Ω G p = A by definition) Secondly, over long timescales the star spins down as a whole: the reservoir ∆L fluctuates but remains bounded (i.e. ∆ L = [∆L(T obs ) − ∆L(0)]/T obs → 0 for T obs → ∞), so that if T obs is long enough. Now, we do not know the details of the inter-glitch dynamics (R), but it is possible to set an upper bound to∆L R by considering the hypothetical perfect pinning limit (P) in which the vortex creep is completely suppressed: 0 <∆L R <∆L P andΩ P p <Ω R p < 0. In this way, we obtain the constraint Note that it is not important whether or not the perfect pinning is realized in a real pulsar: what we are interested in is to give an estimate of∆L P , orΩ P p , in order to set a limit on the real averaged dynamics in (A13).
To do this, we follow the analysis in [14] and assume that the superfluid component is non-rigidly rotating, so that the angular velocity lag is 1 Ω np (x, z, t). Since we are considering a spacetime with circular symmetry, the spacetime metric reads ds 2 = −e 2Φ dt 2 + e 2Λ dr 2 + r 2 dϑ 2 + r 2 sin 2 ϑ (dϕ − ωdt) 2 . (A17) Following the analysis in [14], the angular momentum reservoir ∆L is while we define the total moment of inertia as whereω = ω/Ω p , P is the pressure, ρ is the internal energy and y n is the fraction of neutrons in the region where pinning is possible. However, in (A11) we explicitly neglected the temporal variations of I, but this is in sharp contrast with the fact thatω =ω(r, θ, Ω p , Ω np ), see (A19). Therefore, in the following we will ignore the time variations of the rescaled frame dragω for simplicity, although they may play a relevant role [50,51]. This amounts to neglect the dependence ofω on the small lag Ω np . In this way, in the limit of slow rotation, we have that ω(r, Ω p ) =ω(r)Ω p , whereω(r) is a fixed radial function [49]. For our numerical estimates we calculateω(r) by following the slow rotation prescription of Hartle [49] (in particular, all the structural functions ρ, P, Λ, Φ, ε n , y n are radial functions that can be obtained by solving the TOV equations [14,65]).
Since we have to invoke the perfect pinning condition (P), it is convenient to proceed by using the lag between the normal component and the superfluid momentum, defined as [14,18,65] where ε n is the entrainment parameter [22,66]. All the relations obtained until now are valid also in the presence of entrainment, the only difference is that ∆L is written as a function of the rescaled lag Ω vp instead of Ω np , by using equation (A20). Let us also define the lag derivative ∂ t Ω P vp > 0, which sets an upper limit on the value of ∂ t Ω R vp = (1 − ε n )∂ t Ω R np and is realised when the vortex configuration is perfectly pinned. The time derivative of Ω P vp reads [14] ∂ t (Ω P vp + Ω P p − ω(Ω P vp , Ω P p )) = 0 ⇒ ∂ t Ω P vp = −Ω P p + ∂ t ω(Ω P vp , Ω P p ) ≈ −Ω P p (1 −ω) .
Clearly, equation (A11) must hold also if the perfect pinning condition is assumed in the inter-glitch time, namely ∆L[∂ t Ω P np ] + IΩ P p = −I|Ω ∞ | . (A22) 1 We define as (r, ϑ, ϕ) the spherical coordinates, with ϑ and ϕ being the polar and azimuthal angles, respectively; the cylindrical coordinates are defined as (x, z, ϕ), with z being the ϑ = 0 axis.
Employing equation (A21) in the above relation, we finḋ where I v = ∆L[1 −ω], which coincides with the formula given in (14). Finally, let us now come back to equation (A16): by using the above result for Ω P p , we finally obtain which is equivalent to the constraint in equation (13). Let us remark that the present result is based on the quasi-stationary approach used in [14,65], an approximation that can be justified by the fact that the glitch rise time is expected to be orders of magnitude larger than the hydrodynamical timescale, as discussed in [51]. Of course, an analogous result can be obtained in a completely rigorous way (since there is no need to find approximations for the frame drag ω(Ω p , Ω np ) and its temporal derivative) also in a Newtonian context: the form is still the one in (13) and I v is given by (14), but with ω = Λ = Φ = 0.