Statistical estimates of the pulsar glitch activity
Alessandro Montoli, Marco Antonelli, Brynmor Haskell, Pierre Pizzochero
AArticle
Statistical estimates of the pulsar glitch activity
Alessandro Montoli *, Marco Antonelli *, Brynmor Haskell and Pierre Pizzochero Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy Istituto Nazionale di Fisica Nucleare, sezione di Milano, Via Celoria 16, 20133 Milano, Italy Nicolaus Copernicus Astronomical Center of the Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warszawa,Poland * Email: [email protected] (A.M.), [email protected] (M.A.)
Abstract:
A common way to calculate the glitch activity of a pulsar is an ordinary linear regression ofthe observed cumulative glitch history. This method however is likely to underestimate the errors on theactivity, 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 welljustified in pulsar data. In this paper, we review the extrapolation of the glitch activity parameter andexplore two alternatives: the relaxation of the homoscedasticity hypothesis in the linear fit and the useof the bootstrap technique. Our main finding is a much larger uncertainty on activity estimates, withrespect to that obtained with an ordinary linear regression. We discuss how this affects the theoreticalupper bound on the moment of inertia associated to the region of a neutron star containing the superfluidreservoir of angular momentum released in a stationary sequence of glitches. We find that this upperbound 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.
Keywords: neutron stars; pulsars; superfluidity; hydrodynamics; general relativity
1. Introduction
To date, pulsar glitches are considered the most striking macroscopic manifestation of the presence ofa 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 normalone – 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 twocomponents can momentarily recouple (probably due to a mechanism known as vortex-mediated mutualfriction [4,5]): the transfer of angular momentum from the neutron superfluid to the normal componentresults 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–8], and its information canbe encoded with a study on the glitch size and waiting time distributions (see, e.g., [9–11]), but theidentification of precise trends is difficult due to the scarcity of data in some objects: because of acombination of intrinsic physical properties and different time span of observations, some objects havepresented 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 thisphenomenon can be obtained from its glitching behaviour, like the largest displayed glitch [14] or theshort-time angular velocity evolution after a glitch [15–18]. Furthermore, it is possible also to obtaininformation about the neutron star structure (in particular, on the ratio between the moments of inertia ofthe normal and superfluid components [19–23]) from the observed activity parameter (i.e. the average a r X i v : . [ a s t r o - ph . H E ] D ec of 19 spin up due to glitches, see e.g. [6]) of a pulsar. This parameter is particularly interesting as pulsar activityobservations, together with theoretical modelling of the thermo-rotational evolution of a pulsar, have alsobeen used to provide indirect mass estimates of isolated neutron stars [24–26].Despite early models considered the superfluid neutrons in the core (see e.g. [2,27]), after the seminalwork of Anderson and Itoh [3] the superfluid participating in the glitch has been generally thought to belimited in the crust of the star, where vortex pinning is possible [28–30]. In fact, pulsar activity measured inthe 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 twocomponents [31] – and the calculation of this parameter in the neutron star crust [32] has posed seriousissues to the modelling: entrainment coupling in the crust diminishes the effective angular momentumreservoir, 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 thickto store a significant amount of angular momentum, corresponding to a fraction of crust moment of inertiain a range going from 1.6 % up to ∼
2. 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 − T obs . The absolute activity of a pulsar can be defined as A a = T obs N gl − ∑ i = ∆Ω i . (1)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 aunique 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 isuseful to introduce also the dimensionless activity parameter G , defined as (see e.g. [22]) G = | ˙ Ω ∞ | − A a , (2)where | ˙ Ω ∞ | is the absolute value of the secular spin down rate of the pulsar, namely the average angularvelocity derivative in the period T obs containing several glitches. The variable G gives us an idea of the of 19 amount of spin down reversed by glitches, and it allows for a better comparison of pulsars with differentspin down rates [26].From the practical point of view, to use (1) with real data, the basic requirement is that the pulsarshould regularly be monitored during the interval T obs , without missing any glitch. This is in general notthe case, but while large glitches can be easily detected, very small glitches, easier to miss, should notcontribute to the activity in a significant way (unless they are extremely frequent). Furthermore, it is notalways 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 wewill explore how to extract estimates of A a from real data. Let us assume that the information at our disposal consists only of a list of glitch dates t i andamplitudes ∆Ω i . For simplicity, let us neglect the extra information that is possibly contained in the valueof T obs , but note that T obs > t N gl − − t . Under this assumption, the absolute activity is usually calculatedby fitting the cumulative glitch amplitude [see, e.g., 20,42]. In this case, the relationship between angularvelocity and time is described by the equation Ω i = A a t i + q + ε i , (3)where q is the vertical intercept, and ε i are independent random variables with zero expectation and samevariance (homoscedasticity). In the above formula, Ω i and t i represent the angular velocity acquired bythe star due to glitches and the time passed since the first glitch, Ω i = i ∑ j = ∆Ω j t i = i ∑ j = ∆ t j , (4)where ∆ t j is the waiting time preceding the j -th glitch. This procedure sacrifices the information relative tothe first glitch amplitude, ∆Ω , 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 bythe cumulative points ( t i , Ω i ) , instead of the points themselves [26,42]. An example of the activity fitperformed using this prescription is shown in Figure 1 for the six pulsars which have displayed the largestnumber 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 statisticalproperties of the processes underlying the glitch behaviour should not change if the window of observationis translated in time (i.e. the glitch series observed in a pulsar may be modelled as the outcome of astationary stochastic process on the long run [43]). If this is the case, then the available data set shouldcorrespond to a stationary sequence of glitches where possible aftershocks and more quiet periods ofactivity are both present many times, intertwined in such a way to produce an overall stationary spinup rate (eventually also many very small and undetected glitches may be included, as their contributionto the cumulative glitch amplitude is negligible). In this way, the activity calculated on sub-intervalsshould fluctuate around the asymptotic value calculated by considering an interval T obs containing severalglitches [41].For some pulsars, however, the observational window is so limited and the number of glitchesdetected so small that it may be unsafe to conclude that the observed value of A a corresponds to the valuethat would be extrapolated by looking at a longer sequence of glitches. Practically, to be able to perform a of 19 C u m u l a t i v e g li t c h a m p li t u d e [ r a d / s ] PSR J0534+2200 N gl = 26 N max = 2.1 0 10 20 30 40 50010203040 PSR J0537-6910 N gl = 45 N max = 16.70.0 0.2 0.40.00.20.40.60.81.0 C u m u l a t i v e g li t c h a m p li t u d e [ r a d / s ] PSR J0631+1036 N gl = 17 N max = 1.6 0 5 10 150510152025 PSR J0835-4510 N gl = 20 N max = 11.70.0 0.5 1.0 1.5 2.0 2.5Nominal lag [10 rad/s]012345 C u m u l a t i v e g li t c h a m p li t u d e [ r a d / s ] PSR J1341-6220 N gl = 23 N max = 5.5 0.0 0.2 0.4 0.6 0.8Nominal lag [10 rad/s]0.00.20.40.60.81.0 PSR J1740-3015 N gl = 36 N max = 3.7 Figure 1.
The glitch steps drawn by the cumulative points ( | ˙ Ω ∞ | t i , Ω i ) , defined in (4), for the six pulsarswith the largest number of glitches N gl , as reported by the Jodrell Bank Glitch Catalogue [8]. The activity iscalculated with a least-squares linear fit on the midpoints of the cumulative glitch sequence [see 42], givingthe blue curve with the associated uncertainty on the slope. Following Montoli et al. [26], the plots aremade by using the “nominal lag” t | ˙ Ω ∞ | on the horizontal axis (a rescaling of time t which gives a roughestimate of the typical angular velocity lag accumulated between the spinning down normal componentand a pinned superfluid component). In this way, the slope of the blue curves is the dimensionless activity G . of 19 standard linear regression on the cumulative data we have to demand that the available data fulfill twopractical requirements [26]. Firstly, the number of glitches should be significant, say N gl >
3. Secondly, atleast two glitches of size comparable to ∆Ω max , the maximum glitch amplitude in the sequence, should bepresent: a linear regression would poorly fit the set of data if the largest glitch Ω max is significantly largerthan all the others. This last property can be quantified by demanding that the parameter N max , defined as[25,26] N max = ∆Ω − N gl − ∑ i = ∆Ω i > ∆Ω max = max i = N gl − ∆Ω i , (5)is larger than ∼ A a and the associateduncertainty, especially in the case of small N gl or N max , has not been discussed in detail yet. Besides all the issues related to linear regression mentioned in the previous section, there is onemore linked to the fact that the data points in (4) are not independent, as they arise from a cumulativeconstruction. Hence, the resulting residuals (which have to be minimised in the standard regressionprocedure) may not be independent and identically distributed: the fit residuals will not have the samevariance, so it is no more possible to assume homoscedasticity, which is a basic assumption of ordinarylinear regression.For this reason, we present here an alternative way to calculate the activity, but this time relaxing thehypothesis 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 ∆Ω i = A a ∆ t i + ε i . (6)Note that the above equation may not be true in general, as there seems not to be correlation betweenglitch sizes and waiting times [45]. Using equation (4), the cumulative times and sizes follow a relation Ω i = A a t i + i ∑ j = ε j , (7)which justifies our assumption of heteroscedasticity, as the deviations ∑ ij = ε j have different variance foreach i . It can be shown that the best unbiased linear estimator is given by [44] A a = ∑ i ∆Ω i ∑ i ∆ t i = ∑ i ∆Ω i t N gl − − t i =
1, ..., N gl − ( A a ) = ( N gl − )( t N gl − − t ) ∑ i ( ∆Ω i − A a ∆ t i ) ∆ t i i =
1, ..., N gl − G , with its standard deviationobtained with this method. This value has been obtained by neglecting the first glitch size ∆Ω ineach pulsar sequence, in order to have the same number of glitch sizes and preceding waiting times. of 19 The most interesting feature we notice from the results is that the standard deviation when we assumehomoscedasticity is almost one order of magnitude smaller than that with heteroscedasticity.
One alternative way to solve the issue in the calculation of the activity with the linear fit requires toemploy the probability distributions for the waiting times and sizes of the glitches of a particular pulsar. Inthis 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 previousworks [9–11], which show that – as a general trend – glitch sizes seem to be consistent with a power-lawdistribution, while the waiting times are consistent with an exponential distribution. The Vela and PSRJ0537-6910, are somehow exceptional, as they seem be well described by normal distribution in both sizesand waiting times.Starting from the probability distribution of the waiting times ∆ t and sizes ∆Ω it is possible to infersome 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 (cid:54) = M , are different random variables, distributed according to different laws,i.e. P A N (cid:54) = P A M .Given the definition of activity, its distribution can be obtained by considering the ratio of two randomvariables, the sum of sizes ∆ ˜ Ω N and of waiting times ∆ ˜ t N . The latter, in turn, are sum of random variablesthemselves, 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, ∆ ˜ Ω N = N ∑ i = ∆Ω i ⇒ ∆ ˜ Ω N ∼ P ∆ ˜ Ω N = P ∆Ω ∗ . . . ∗ P ∆Ω (cid:124) (cid:123)(cid:122) (cid:125) N times (10) ∆ ˜ t N = N ∑ i = ∆ t i ⇒ ∆ ˜ t N ∼ P ∆ ˜ t N = P ∆ t ∗ . . . ∗ P ∆ t (cid:124) (cid:123)(cid:122) (cid:125) N times (11)Since A N = ∆ ˜ Ω N / ∆ ˜ t N , its distribution can be obtained through: P A N ( a ) = (cid:90) ∞ − ∞ d x | x | P ∆ ˜ t N ( x ) P ∆ ˜ Ω N ( xa ) . (12)The main advantage of this method is that, with the only assumption of independence between sizes andwaiting 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 thisdistribution, 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 meanand the variance, by employing a generalised version of the central limit theorem, the delta method (seeAppendix A). This method allows us to calculate the mean and the variance of P A N starting from the meanand 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. An attractive alternative to the methods described above, which does not build on the linearassumption in (6), is the so-called “bootstrap method” [46]. The idea is that of resampling with replacementthe original data in order to calculate some statistics, as, e.g., the mean and standard deviation of thecalculated activity. In our case, the samples are two: the list of the waiting times ∆ t i (of length N gl − of 19 Pulsar G hom [ % ] G het [ % ] G rand [ % ] G pre [ % ] G post [ % ] ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± Table 1.
Dimensionless activities and their standard deviations, calculated for the six pulsars with the largestnumber 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 andwaiting time samples separately ( G rand ), and on the pairs size - preceding waiting time ( G pre ) and size -following 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 waitingtime-size couples in order to have a fair estimation of the activity and its standard deviation. To avoid thehomoscedasticity 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 aglitch size and the preceding or the subsequent waiting time, so it is useful to also bootstrap on other twosamples: the sample made up by ordered pairs { ( ∆Ω i , ∆ t pre i ) } i = N gl − , where ∆ t pre i is the waiting timepreceding the glitch of size ∆Ω i , and the sample comprised by ordered pairs { ( ∆Ω i , ∆ t post i ) } i = N gl − ,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 times in all the three casesdescribed above. In the same plot, also the dimensionless activities obtained as a result of the ordinarylinear regression on the cumulative glitch data are displayed. We can see that the activity calculated bymeans of bootstrapping is compatible with the results obtained from an ordinary linear regression, butit generally has larger standard deviations (see also Table 1). It is interesting to notice the case of PSRJ0537-6910, one of the few stars which presents a significant correlation between a glitch size and thefollowing 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 theother two cases. At lower confidence, also PSR J0631+1036 shows a correlation between size and thefollowing waiting time, and Vela a correlation between size and the preceding waiting time [45]. Thesecorrelations 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 showstwo clear peaks, one on very small values and one around G ≈ N max for this star is smaller than the ones of the other stars in the sample: this value mayincrease 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. Thefirst point of each curve is the activity calculated with the linear fit on the cumulative midpoints usingthe first ten glitches, assuming homoscedastic data. Then, we update the activity value with the samemethod whenever a new glitch is displayed. We also plot the activity parameter calculated with theordinary linear regression and all the glitches, along with its uncertainty. For PSR J0537-6910 and PSRJ0631+1036, we present the bootstrap estimate G post and its standard deviation, for Vela we show the same of 19 P ( G ) PSR J0534+2200 0.004 0.006 0.008 0.010 0.012 0.0140200400600800100012001400 PSR J0537-69100.00 0.02 0.04 0.06 0.08 0.10020406080 P ( G ) PSR J0631+1036 0.010 0.015 0.020 0.025 0.030050100150200 PSR J0835-45100.01 0.02 0.03 0.04 0.05G01020304050607080 P ( G ) PSR J1341-6220 0.00 0.01 0.02 0.03 0.04G020406080 PSR J1740-3015Linear fitRandomprepost
Figure 2.
Dimensionless pulsar activity G , calculated by sampling both the size and waiting time samplesrandomly (in blue), by sampling the pair ( ∆Ω , ∆ t pre ) (in orange) and the pair ( ∆Ω , ∆ t post ) (in green). Theresults of the linear fit of the cumulative glitch data is also plotted (in red, the shaded area is the 1 σ region). of 19 G [ % ] PSR J0534+2200Linear fitRandomprepost 53000 54000 550000.840.860.880.90 PSR J0537-691054000 55000 56000 57000 580001.001.251.501.752.002.252.50 G [ % ] PSR J0631+1036 48000 50000 52000 54000 56000 580001.451.501.551.601.651.701.751.80 PSR J0835-451051000 52000 53000 54000 55000MJD1.21.41.61.82.02.22.42.6 G [ % ] PSR J1341-6220 50000 52000 54000 56000 58000MJD0.81.01.21.41.61.8 PSR J1740-3015
Figure 3.
Evolution of the glitch activity over time. The first point of each curve is the activity calculatedemploying the linear fit on the cumulative midpoints of the first ten glitches, assuming homoscedasticity.The subsequent points are calculated by gradually adding all the glitches that pulsar displayed. For eachpulsar, the estimate of the activity via ordinary linear regression with all glitches (in yellow) is shown, alongwith its uncertainty (shaded). For PSR J0537-6910 and PSR J0631+1036, G post and its uncertainty are alsoshown (green), for PSR J0835-4510 (Vela) G pre (orange), while for all other stars G rand (blue). with G pre , while for all the other pulsars we present the case with uncorrelated glitch sizes and waitingtimes ( 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 eachstar generally lies outside the error region for the linear fit with homeoscedastic data for all the glitchhistory, except – of course – the latest glitch. The bootstrap uncertainty better describes the variance of theglitch history. A notable exception is that of PSR J1341-6220, which is well below the error bar for boththe activity calculations, except for the last three glitches. This is because these glitches are three of thelargest ones displayed by this pulsar (see also Figure 1). In general, however, it is interesting to notice howvariable the activity parameter is. A single large glitch can change its value (see, e.g., the Crab pulsar afterits November 2017 glitch, [48]). This fact stresses the importance of having a much larger uncertainty onthe activity, which is the result of not assuming homoscedasticity in the linear fit or linear dependencebetween glitch sizes and waiting times.
3. Moment of inertia constraint
The activity parameter allows to extract information on the moment of inertia fraction of the superfluidreservoir in a glitching pulsar [19,20]. In this Section, we present a revised version of the constraint onthe 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 superfluidcomponents. Our derivation is presented in Appendix B and takes into account the non-rigid rotation ofthe superfluid component, the stellar stratification and the non-uniform entrainment coupling between thecomponents. The constraint is given by equation (A16), which can be written in terms of the total momentof inertia I and of the moment of inertia of the superfluid component I v as I v I − I v > 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 = π (cid:90) R d R c dr r e Λ ( r ) − Φ ( r ) [ ρ ( r ) + P ( r )] y n ( r ) − ε n ( r ) Ω p − ω ( r ) Ω p , (14)while the total moment of inertia I is the usual one in the slow rotation framework [14,49,50], givenin equation (A19). In the above equation, y n is the superfluid neutron baryon density (limited to theregion where pinning is possible) divided by the total baryon density, i.e. y n ( r ) is different from zeroonly where the superfluid can pin to inhomogeneities and maintain its state of motion while the normalcomponent 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 containa dependence on the angular velocities of both components [14,50,51], which has been neglected in thederivation of (14).The first measurements of the activity parameter of the Vela pulsar and the moment of inertia fractionestimates for different equations of state (EoSs) seemed to be in accordance with the constraint (13)with ε n = ε n in the crust of a neutron star has been calculatedin [32], by estimating the effects of Bragg scattering on the superflow due to the presence of the crustallattice. These calculations yield a negative entrainment parameter ε n ∼ −
10 in a substantial portion of theinner crust, which implies a severely hindered motion of the superfluid component. This would reducethe amount of extra angular momentum stored in the crust between two glitches – and thus of I v – makingthe requirement in (13) more difficult to be met. As a result of that, the only way for the star to acquire I v / ( II v ) PSR J0835-4510123SLy4BSk20BSk21DDME2
Figure 4.
Activity constraint on the superfluid component moment of inertia plotted for some EoS: SLy4,[53], BSk20 and BSk21 [54], and the DDME2 EoS [55], glued with a SLy4 crust [following the methoddescribed in 56]. The entrainment parameter is that calculated in [32], and the superfluid reservoir is limitedin the crust of the star. The dimensionless activity parameter G – calculated with the bootstrap methoddescribed in Section 2 (the case with random glitch sizes and waiting times) – is also plotted for the Velapulsar, along with the 1 σ , 2 σ and 3 σ uncertainties. enough angular momentum between glitches to explain the observed activity is to have a large regioninside it to store angular momentum (larger than the crust of the star), or to have an unreasonably smallmass, around ∼ M (cid:12) . This problem has been highlighted in several papers [21–24,52].If we assume the superfluid region to be limited in the crust of the star, and we fix the microphysicalparameters, namely the EoS and the entrainment parameter, then the moment of inertia fraction in (13) is afunction 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 forsome EoSs and for the entrainment parameter calculated by Chamel [32], assuming a superfluid reservoirlimited to the crust. As we can see, the constraint (13) imposes that the high G value of Vela pulsar (PSRJ0835-4510) can be explained only if the mass of the Vela is small, ranging from ≈ M (cid:12) of the BSk20 EoSto ≈ M (cid:12) for the SLy4 EoS. Let us however consider also the 1 σ uncertainty region, calculated with thebootstrap method described in Section 2, with random waiting times and glitch sizes. In this case for theBSk20 EoS, the Vela pulsar has an upper limit on the mass of about 1.2 M (cid:12) . Note that this value is slightlyabove the minimum mass of a neutron star estimated from calculations of core-collapse supernovae (i.e.1.17 M (cid:12) , see [57]) and the smallest mass measured in a neutron star (1.174 ± M (cid:12) , measured in PSRJ0453+1559, [58]). If we consider the 3 σ uncertainty range, then we obtain a limit of ≈ M (cid:12) for the sameEoS. It is thus clear that a careful estimation of the activity parameter and the associated error is crucial ifone is to set strict quantitative constraints on the moment of inertia involved in glitches, and constrain thelocation of the superfluid reservoir.
4. Conclusions
In this paper, we have discussed different ways of calculating the activity parameter in glitchingpulsars. The most commonly used one is that of performing an ordinary linear regression to the cumulativeglitch data, which is justified if one considers the activity parameter as an intrinsic characteristic of aglitching 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 autocorrelationor 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 consequenceof including these assumptions in the fitting procedure (namely, the linear dependence of glitch sizesand waiting times and the homoscedasticity of the cumulative data) is an underestimation of the activityuncertainty. This is an important point, as it can be shown that an additional single glitch or a sequence ofa 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 thehypothesis of homoscedasticity of the data has been relaxed. This relaxation is justified by the fact thatcumulative data is not homoscedastic, as the variance of the data increases as data is cumulated. The mainconsequence of this assumption is that the uncertainty of the activity parameter increases by about oneorder 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 andwaiting times, by employing their probability distribution estimates. While employing the size andwaiting time distributions is arguably the best way to include the information provided by observations ofthe pulsar in the activity parameter, it is also true that a study of the probability distribution of the activityparameter is computationally very challenging. Also an approximate version of this estimate, obtained byemploying a more general version of the central limit theorem, leads to difficulties, due to the fat-tailedprobability distribution of the glitch sizes.We thus employed an alternative way to calculate the activity and its uncertainty, relaxing thehypothesis of linear dependence, by bootstrapping on the sizes and waiting time data sets and calculatingthe activity by summing sizes and waiting times and calculating the ratio. Much larger uncertainties areobtained, on the same order of magnitude of the linear fit on heteroscedastic data. We then used theseresults to study a revised version of the constraint on the moment of inertia of the superfluid componentand on the mass of the star [20–22]. While the result obtained in this way is essentially the same forthe mean value of the activity (i.e. the Vela activity does not allow a crust-limited reservoir and strongentrainment), the larger uncertainty allows for models in which the crust is enough, i.e. models in whichthe superfluid reservoir is located entirely in the crust for credible values of the star mass. For examplea model of a 1.2 M (cid:12) neutron star described by the BSk20 EoS is within the 1 σ uncertainty region of theactivity that we calculate.A more careful estimate of the glitch activity parameter may thus play an important role in resolvingthe tension between strong entrainment and models with a crust-limited reservoir. Several other effectshave, of course, been suggested, and are likely to play a role in this problem, including a maximallystiff EoS [60], a Bayesian analysis of the EoS uncertainty [52] or an extension of the region where theneutron superfluid participates in the glitch beyond the crust-core transition, based on the assumptionthat only the superfluid in the S state participates in the glitch phenomenon and on an analysis of thetemperature of the star [24]. On the other hand, also different calculations of the entrainment parameterhave been proposed (e.g. [33,61,62]) which yield milder entrainment effects in the crust. However, even ifthe 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 inthe core of the star [16,63,64], independently of the presence of entrainment in the model [18]. Acknowledgments:
Partial support comes from PHAROS, COST Action CA16214. M.A. and B.H. acknowledgesupport from the Polish National Science Centre grants SONATA BIS 2015/18/E/ST9/00577 and OPUS2019/33/B/ST9/00942. We thank Cristóbal M. Espinoza and Andrew Melatos for providing interesting and usefulcomments.
Appendix A. Activity calculation with the delta method
We use the delta-method (a generalisation of the central limit theorem) to extract the mean andstandard deviation of the activity parameter A N after N glitches, given a probability distribution of thesizes P ∆Ω and waiting times P ∆ t of a pulsar. The method works under the simplifying assumption that thesizes and the waiting times are independent and identically distributed random variables. Let us recallthe two random variables ∆ ˜ Ω N and ∆ ˜ t N defined in equations (10) and (11). The random variable A N isdefined as A N = ∆ ˜ Ω N / ∆ ˜ t N . (A1)Its expectation value is given by ( a ∈ R + denotes the value assumed by A N )E [ A N ] = (cid:90) a P A N ( a ) d a = (cid:90) ∆ ˜ Ω N ∆ ˜ t N N ∏ i = P ∆ t ( ∆ t i ) d ∆ t i N ∏ j = P ∆Ω ( ∆Ω j ) d ∆Ω j . (A2)Independence of the variables is loosely justified by observing little correlation and autocorrelation inglitch sizes and waiting times [45,59], while being identically distributed is a working assumption. Theabove equation boils down to E [ A N ] = E [ ∆ ˜ Ω N ] E (cid:20) ∆ ˜ t N (cid:21) . (A3)To calculate the variance of A N , let us first calculate, for the sizes ∆ ˜ Ω N :E ( ∆ ˜ Ω N ) = N E [ ∆Ω ] + N ( N − ) E [ ∆Ω ] (A4)E ( ∆ ˜ Ω N ) = N E [ ∆Ω ] (A5)Var [ ∆ ˜ Ω N ] = N ( E [ ∆Ω ] − E [ ∆Ω ] ) = N Var [ ∆Ω ] . (A6)Now, using the above results, a simple direct calculation givesVar [ A N ] = N Var [ ∆Ω ] E (cid:34) ∆ ˜ t N (cid:35) + N E [ ∆Ω ] Var (cid:20) ∆ ˜ t N (cid:21) . (A7)Let us now assume that P ∆ ˜ t N converges to a normal distribution with mean N θ and variance N σ as N → ∞ . In this case, we can employ the delta method, which tells us that any function g of the randomvariable ∆ ˜ t N will be distributed normally with mean g ( N θ ) with variance N [ g (cid:48) ( N θ ) σ ] . Note that thisassumption 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 σ = Var [ ∆ t ] , by considering the cases g ( ∆ ˜ t N ) = ( ∆ ˜ t N ) − and g ( ∆ ˜ t N ) = ( ∆ ˜ t N ) − we obtain:E (cid:20) ∆ ˜ t N (cid:21) ≈ N E [ ∆ t ] E (cid:34) ∆ ˜ t N (cid:35) ≈ N E [ ∆ t ] Var (cid:20) ∆ ˜ t N (cid:21) ≈ Var [ ∆ t ] N E [ ∆ t ] (A8) Given the first of the above relations, the expectation value of A N in (A3) can be approximated asE [ A N ] ≈ E [ ∆Ω ] E [ ∆ t ] , (A9)while the variance of A N in (A7) is given byVar [ A N ] ≈ N − (cid:20) Var [ ∆Ω ] E ( ∆ t ) + E [ ∆Ω ] E [ ∆ t ] Var [ ∆ t ] (cid:21) . (A10) Appendix B. Derivation of the moment of inertia constraint
Let us write the dynamics of the total angular momentum L of the pulsar (neglecting temporalvariations in the total moment of inertia I ) as ∂ t L [ Ω p , Ω np ] = ˙ ∆ L + I ˙ Ω p = − I | ˙ Ω ∞ | , (A11)where we assumed that the normal component is rigidly rotating with angular velocity Ω p , while ∆ L isthe angular momentum reservoir due to the (non-rigid) angular velocity lag Ω np between the superfluidand the normal component [14]. Clearly, (A11) is valid only if we are assuming that the rigid componentand 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 ) andan impulsive one from glitches ( G ),˙ Ω p = ˙ Ω Gp + ˙ Ω Rp ˙ ∆ L = ˙ ∆ L G + ˙ ∆ L R . (A12)During glitches we have ˙ ∆ L G < Ω Gp >
0, while for the rest of the time ˙ ∆ L R > Ω Rp <
0. We canaverage (A11) over a long time interval T obs , to get (cid:104) ˙ ∆ L G (cid:105) + I (cid:104) ˙ Ω Gp (cid:105) + (cid:104) ˙ ∆ L R (cid:105) + I (cid:104) ˙ Ω Rp (cid:105) = − I | ˙ Ω ∞ | . (A13)We can simplify the equation above by making two observations. Firstly, due to the angular momentumconservation during a glitch we must have that (note that (cid:104) ˙ Ω Gp (cid:105) = A by definition) (cid:104) ˙ ∆ L G (cid:105) + I (cid:104) ˙ Ω Gp (cid:105) = ⇒ (cid:104) ˙ ∆ L G (cid:105) = − I A . (A14)Secondly, over long timescales the star spins down as a whole: the reservoir ∆ L fluctuates but remainsbounded (i.e. (cid:104) ˙ ∆ L (cid:105) = [ ∆ L ( T obs ) − ∆ L ( )] / T obs → T obs → ∞ ), so that (cid:104) ˙ ∆ L (cid:105) = (cid:104) ˙ ∆ L G (cid:105) + (cid:104) ˙ ∆ L R (cid:105) ≈ T obs is long enough. Now, we do not know the details of the inter-glitch dynamics ( R ) , but it is possibleto set an upper bound to ˙ ∆ L R by considering the hypothetical perfect pinning limit ( P ) in which the vortexcreep is completely suppressed: 0 < ˙ ∆ L R < ˙ ∆ L P and ˙ Ω Pp < ˙ Ω Rp <
0. In this way, we obtain the constraint (cid:104) ˙ ∆ L P (cid:105) > I A or (cid:104) ˙ Ω Pp (cid:105) + A < −| ˙ Ω ∞ | . (A16) Note that it is not important whether or not the perfect pinning is realized in a real pulsar: what we areinterested in is to give an estimate of ˙ ∆ L P , or ˙ Ω Pp , 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-rigidlyrotating, so that the angular velocity lag is Ω np ( x , z , t ) . Since we are considering a spacetime with circularsymmetry, the spacetime metric readsd s = − e Φ d t + e Λ d r + r d ϑ + r sin ϑ ( d ϕ − ω d t ) . (A17)Following the analysis in [14], the angular momentum reservoir ∆ L is ∆ L [ Ω np ] = (cid:90) d x x e Λ − Φ ( ρ + P ) y n Ω np , (A18)while we define the total moment of inertia as I = (cid:90) d x x e Λ − Φ ( ρ + P ) ( − ˜ ω ) , (A19)where ˜ ω = ω / Ω p , P is the pressure, ρ is the internal energy and y n is the fraction of neutrons in the regionwhere pinning is possible. However, in (A11) we explicitly neglected the temporal variations of I , butthis is in sharp contrast with the fact that ˜ ω = ˜ ω ( r , θ , Ω p , Ω np ) , see (A19). Therefore, in the followingwe will ignore the time variations of the rescaled frame drag ˜ ω for simplicity, although they may play arelevant role [50,51]. This amounts to neglect the dependence of ˜ ω on the small lag Ω np . In this way, inthe limit of slow rotation, we have that ω ( r , Ω p ) = ˜ ω ( r ) Ω p , where ˜ ω ( r ) is a fixed radial function [49]. Forour numerical estimates we calculate ˜ ω ( r ) by following the slow rotation prescription of Hartle [49] (inparticular, all the structural functions ρ , P , Λ , Φ , ε n , y n are radial functions that can be obtained by solvingthe TOV equations [14,65]).Since we have to invoke the perfect pinning condition ( P ) , it is convenient to proceed by using thelag between the normal component and the superfluid momentum, defined as [14,18,65] Ω vp = ( − ε n ) Ω np ⇒ ∆ L [ Ω np ] = ∆ L [ Ω vp / ( − ε n )] , (A20)where ε n is the entrainment parameter [22,66]. All the relations obtained until now are valid also in thepresence 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 Ω Pvp >
0, which setsan upper limit on the value of ∂ t Ω Rvp = ( − ε n ) ∂ t Ω Rnp and is realised when the vortex configuration isperfectly pinned. The time derivative of Ω Pvp reads [14] ∂ t ( Ω Pvp + Ω Pp − ω ( Ω Pvp , Ω Pp )) = ⇒ ∂ t Ω Pvp = − ˙ Ω Pp + ∂ t ω ( Ω Pvp , Ω Pp ) ≈ − ˙ Ω Pp ( − ˜ ω ) . (A21)Clearly, equation (A11) must hold also if the perfect pinning condition is assumed in the inter-glitch time,namely (cid:104) ∆ L [ ∂ t Ω Pnp ] (cid:105) + I ˙ Ω Pp = − I | ˙ Ω ∞ | . (A22) We define as ( r , ϑ , ϕ ) the spherical coordinates, with ϑ and ϕ being the polar and azimuthal angles, respectively; the cylindricalcoordinates are defined as ( x , z , ϕ ) , with z being the ϑ = Employing equation (A21) in the above relation, we find˙ Ω Pp ( I − I v ) = (cid:104) ˙ Ω Pp (cid:105) ( I − I v ) = − I | ˙ Ω ∞ | , (A23)where I v = ∆ L [ − ˜ ω ] , which coincides with the formula given in (14). Finally, let us now come back toequation (A16): by using the above result for (cid:104) ˙ Ω Pp (cid:105) , we finally obtain A − II − I v | ˙ Ω ∞ | < −| ˙ Ω ∞ | , (A24)which is equivalent to the constraint in equation (13). Let us remark that the present result is based on thequasi-stationary approach used in [14,65], an approximation that can be justified by the fact that the glitchrise 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 needto find approximations for the frame drag ω ( Ω p , Ω np ) and its temporal derivative) also in a Newtoniancontext: the form is still the one in (13) and I v is given by (14), but with ω = Λ = Φ = References
1. Haskell, B.; Melatos, A. Models of pulsar glitches.
International Journal of Modern Physics D , , 1530008,[arXiv:astro-ph.SR/1502.07062]. doi:10.1142/S0218271815300086.2. Baym, G.; Pethick, C.; Pines, D.; Ruderman, M. Spin Up in Neutron Stars : The Future of the Vela Pulsar. Nature , , 872–874. doi:10.1038/224872a0.3. Anderson, P.W.; Itoh, N. Pulsar glitches and restlessness as a hard superfluidity phenomenon. Nature , , 25–27. doi:10.1038/256025a0.4. Andersson, N.; Sidery, T.; Comer, G.L. Mutual friction in superfluid neutron stars. Mon. Not. R. Astron. Soc. , , 162–170, [astro-ph/0510057]. doi:10.1111/j.1365-2966.2006.10147.x.5. Antonelli, M.; Haskell, B. Superfluid vortex-mediated mutual friction in non-homogeneous neutron starinteriors. Mon. Not. R. Astron. Soc. , [arXiv:astro-ph.HE/2007.11720]. doi:10.1093/mnras/staa3097.6. McKenna, J.; Lyne, A.G. PSR1737-30 and period discontinuities in young pulsars.
Nature , , 349–350.doi:10.1038/343349a0.7. Lyne, A.G.; Shemar, S.L.; Graham Smith, F. Statistical studies of pulsarglitches. Monthly Notices of the Royal Astronomical Society , , 534–542,[http://oup.prod.sis.lan/mnras/article-pdf/315/3/534/2935525/315-3-534.pdf].doi:10.1046/j.1365-8711.2000.03415.x.8. Espinoza, C.M.; Lyne, A.G.; Stappers, B.W.; Kramer, M. A study of 315 glitches in the rotationof 102 pulsars. Mon. Not. R. Astron. Soc. , , 1679–1704, [arXiv:astro-ph.HE/1102.1743].doi:10.1111/j.1365-2966.2011.18503.x.9. Melatos, A.; Peralta, C.; Wyithe, J.S.B. Avalanche Dynamics of Radio Pulsar Glitches. Astrophys. J. , , 1103–1118, [0710.1021]. doi:10.1086/523349.10. Howitt, G.; Melatos, A.; Delaigle, A. Nonparametric Estimation of the Size and Waiting Time Distributions ofPulsar Glitches. Astrophys. J. , , 60, [arXiv:astro-ph.HE/1809.06038]. doi:10.3847/1538-4357/aae20a.11. Fuentes, J.R.; Espinoza, C.M.; Reisenegger, A. Glitch time series and size distributions in eight prolific pulsars. Astron. Astrophys. , , A115, [arXiv:astro-ph.HE/1907.09887]. doi:10.1051/0004-6361/201935939.12. Antonopoulou, D.; Espinoza, C.M.; Kuiper, L.; Andersson, N. Pulsar spin-down: the glitch-dominatedrotation of PSR J0537-6910. Mon. Not. R. Astron. Soc. , , 1644–1655, [arXiv:astro-ph.HE/1708.09459].doi:10.1093/mnras/stx2429.
13. Espinoza, C.M.; Antonopoulou, D.; Dodson, R.; Stepanova, M.; Scherer, A. Small glitches and other rotationalirregularities of the Vela pulsar. arXiv e-prints , p. arXiv:2007.02921, [arXiv:astro-ph.HE/2007.02921].14. Antonelli, M.; Montoli, A.; Pizzochero, P.M. Effects of general relativity on glitch amplitudes and pulsarmass upper bounds.
Mon. Not. R. Astron. Soc. , , 5403–5416, [arXiv:astro-ph.HE/1710.05879].doi:10.1093/mnras/sty130.15. Ashton, G.; Lasky, P.D.; Graber, V.; Palfreyman, J. Rotational evolution of the Vela pulsar during the 2016 glitch. Nature Astronomy , p. 417, [arXiv:astro-ph.HE/1907.01124]. doi:10.1038/s41550-019-0844-6.16. Pizzochero, P.M.; Montoli, A.; Antonelli, M. Core and crust contributions in overshooting glitches:the Vela pulsar 2016 glitch.
Astron. Astrophys. , , A101, [arXiv:astro-ph.HE/1910.00066].doi:10.1051/0004-6361/201937019.17. Sourie, A.; Chamel, N. Vortex pinning in the superfluid core of neutron stars and the rise of pulsar glitches. Mon. Not. R. Astron. Soc. , , L98–L102, [arXiv:astro-ph.HE/2001.09668]. doi:10.1093/mnrasl/slaa015.18. Montoli, A.; Antonelli, M.; Magistrelli, F.; Pizzochero, P.M. Bayesian estimate of the superfluidmoments of inertia from the 2016 glitch in the Vela pulsar. Astron. Astrophys. , , A223,[arXiv:astro-ph.HE/2005.01594]. doi:10.1051/0004-6361/202038340.19. Datta, B.; Alpar, M.A. Implications of the crustal moment of inertia for neutron-star equations of state. Astron.Astrophys. , , 210–212.20. Link, B.; Epstein, R.I.; Lattimer, J.M. Pulsar Constraints on Neutron Star Structure and Equation of State. Phys. Rev. Lett. , , 3362–3365, [astro-ph/9909146]. doi:10.1103/PhysRevLett.83.3362.21. Andersson, N.; Glampedakis, K.; Ho, W.C.G.; Espinoza, C.M. Pulsar Glitches: The Crust is not Enough. Phys. Rev. Lett. , , 241103, [arXiv:astro-ph.SR/1207.0633]. doi:10.1103/PhysRevLett.109.241103.22. Chamel, N. Crustal Entrainment and Pulsar Glitches. Phys. Rev. Lett. , , 011101,[arXiv:astro-ph.HE/1210.8177]. doi:10.1103/PhysRevLett.110.011101.23. Delsate, T.; Chamel, N.; Gürlebeck, N.; Fantina, A.F.; Pearson, J.M.; Ducoin, C. Giant pulsar glitchesand the inertia of neutron star crusts. Phys. Rev. D , , 023008, [arXiv:astro-ph.HE/1606.00016].doi:10.1103/PhysRevD.94.023008.24. Ho, W.C.G.; Espinoza, C.M.; Antonopoulou, D.; Andersson, N. Pinning down the superfluid and measuringmasses using pulsar glitches. Science Advances , , e1500578–e1500578, [arXiv:astro-ph.SR/1510.00395].doi:10.1126/sciadv.1500578.25. Pizzochero, P.M.; Antonelli, M.; Haskell, B.; Seveso, S. Constraints on pulsar masses from the maximumobserved glitch. Nature Astronomy , , 0134, [arXiv:astro-ph.HE/1611.10223]. doi:10.1038/s41550-017-0134.26. Montoli, A.; Antonelli, M.; Pizzochero, P.M. The role of mass, equation of state, and superfluid reservoirin large pulsar glitches. Mon. Not. R. Astron. Soc. , , 4837–4846, [arXiv:astro-ph.HE/1809.07834].doi:10.1093/mnras/staa149.27. Packard, R.E. Pulsar Speedups Related to Metastability of the Superfluid Neutron-Star Core. Phys. Rev. Lett. , , 1080–1082. doi:10.1103/PhysRevLett.28.1080.28. Alpar, M.A. Pinning and Threading of Quantized Vortices in the Pulsar Crust Superfluid. Astrophys. J. , , 527–530. doi:10.1086/155183.29. Donati, P.; Pizzochero, P.M. Fully consistent semi-classical treatment of vortex-nucleus interaction in rotatingneutron stars. Nucl. Phys. A , , 363–379. doi:10.1016/j.nuclphysa.2004.07.002.30. Seveso, S.; Pizzochero, P.M.; Grill, F.; Haskell, B. Mesoscopic pinning forces in neutron star crusts. Mon. Not. R.Astron. Soc. , , 3952–3967. doi:10.1093/mnras/stv2579.31. Andreev, A.F.; Bashkin, E.P. Three-velocity hydrodynamics of superfluid solutions. Soviet Journal of Experimentaland Theoretical Physics , , 164.32. Chamel, N. Neutron conduction in the inner crust of a neutron star in the framework of the band theory ofsolids. Phys. Rev. C , , 035801. doi:10.1103/PhysRevC.85.035801.33. Sauls, J.A.; Chamel, N.; Alpar, M.A. Superfluidity in Disordered Neutron Stars Crusts. arXiv e-prints , p.arXiv:2001.09959, [arXiv:astro-ph.HE/2001.09959].
34. Gügercino ˘glu, E.; Alpar, M.A. Vortex Creep Against Toroidal Flux Lines, Crustal Entrainment, and PulsarGlitches.
Astrophys. J. Lett. , , L11. doi:10.1088/2041-8205/788/1/L11.35. Muslimov, A.G.; Tsygan, A.I. Vortex lines in neutron star superfluids and decay of pulsar magnetic fields. Ap&SS , , 43. doi:10.1007/BF00653825.36. Srinivasan, G.; Bhattacharya, D.; Muslimov, A.G.; Tsygan, A.J. A novel mechanism for the decay of neutronstar magnetic fields. Current Science , , 31–38.37. Drummond, L.V.; Melatos, A. Stability of interlinked neutron vortex and proton flux tube arraysin a neutron star: equilibrium configurations. Mon. Not. R. Astron. Soc. , , 4851–4869,[arXiv:astro-ph.HE/1709.02254]. doi:10.1093/mnras/stx2301.38. Alpar, M.A. Flux-Vortex Pinning and Neutron Star Evolution. Journal of Astrophysics and Astronomy , , 44,[arXiv:astro-ph.HE/1709.07912]. doi:10.1007/s12036-017-9473-6.39. Wood, T.S.; Graber, V.; Newton, W.G. Superconducting phases in a two-component microscale model ofneutron star cores. arXiv e-prints , p. arXiv:2011.02873, [arXiv:cond-mat.supr-con/2011.02873].40. Leinson, L.B. Vortex lattice in rotating neutron spin-triplet superfluid. Monthly Notices of the Royal AstronomicalSociety , , 304–309, [https://academic.oup.com/mnras/article-pdf/498/1/304/33705847/staa2475.pdf].doi:10.1093/mnras/staa2475.41. Corral, Á. Dependence of earthquake recurrence times and independence of magnitudes on seismicity history. Tectonophysics , , 177–193. doi:10.1016/j.tecto.2006.03.035.42. Wong, T.; Backer, D.C.; Lyne, A.G. Observations of a Series of Six Recent Glitches in the Crab Pulsar. Astrophys.J. , , 447–459, [arXiv:astro-ph/astro-ph/0010010]. doi:10.1086/318657.43. Fulgenzi, W.; Melatos, A.; Hughes, B.D. Radio pulsar glitches as a state-dependent Poisson process. Mon. Not.R. Astron. Soc. , , 4307–4329. doi:10.1093/mnras/stx1353.44. Mandel, J. Fitting a Straight Line to Certain Types of Cumulative Data. Journal of the American StatisticalAssociation , , 552–566.45. Melatos, A.; Howitt, G.; Fulgenzi, W. Size-waiting-time Correlations in Pulsar Glitches. Astrophys. J. , , 196, [arXiv:astro-ph.HE/1809.03064]. doi:10.3847/1538-4357/aad228.46. Efron, B. Bootstrap Methods: Another Look at the Jackknife. Ann. Statist. , , 1–26.doi:10.1214/aos/1176344552.47. Middleditch, J.; Marshall, F.E.; Wang, Q.D.; Gotthelf, E.V.; Zhang, W. Predicting the Starquakes in PSRJ0537-6910. Astrophys. J. , , 1531–1546, [arXiv:astro-ph/astro-ph/0605007]. doi:10.1086/508736.48. Shaw, B.; Lyne, A.G.; Stappers, B.W.; Weltevrede, P.; Bassa, C.G.; Lien, A.Y.; Mickaliger, M.B.; Breton, R.P.;Jordan, C.A.; Keith, M.J.; Krimm, H.A. The largest glitch observed in the Crab pulsar. Mon. Not. R. Astron. Soc. , , 3832–3840, [arXiv:astro-ph.HE/1805.05110]. doi:10.1093/mnras/sty1294.49. Hartle, J.B. Slowly Rotating Relativistic Stars. I. Equations of Structure. Astrophys. J. , , 1005.doi:10.1086/149400.50. Andersson, N.; Comer, G.L. Slowly rotating general relativistic superfluid neutron stars. Classical and QuantumGravity , , 969–1002, [arXiv:gr-qc/gr-qc/0009089]. doi:10.1088/0264-9381/18/6/302.51. Sourie, A.; Chamel, N.; Novak, J.; Oertel, M. Global numerical simulations of the rise of vortex-mediated pulsarglitches in full general relativity. Mon. Not. R. Astron. Soc. , , 4641–4657, [arXiv:astro-ph.HE/1607.08213].doi:10.1093/mnras/stw2613.52. Carreau, T.; Gulminelli, F.; Margueron, J. General predictions for the neutron star crustal moment of inertia. Phys. Rev. C , , 055803, [arXiv:nucl-th/1810.00719]. doi:10.1103/PhysRevC.100.055803.53. Douchin, F.; Haensel, P. A unified equation of state of dense matter and neutron star structure. Astron.Astrophys. , , 151–167, [astro-ph/0111092]. doi:10.1051/0004-6361:20011402.54. Goriely, S.; Chamel, N.; Pearson, J.M. Further explorations of Skyrme-Hartree-Fock-Bogoliubov mass formulas.XII. Stiffness and stability of neutron-star matter. Phys. Rev. C , , 035804, [arXiv:nucl-th/1009.3840].doi:10.1103/PhysRevC.82.035804.55. Lalazissis, G.A.; Nikši´c, T.; Vretenar, D.; Ring, P. New relativistic mean-field interaction with density-dependentmeson-nucleon couplings. Phys. Rev. C , , 024312. doi:10.1103/PhysRevC.71.024312.
56. Fortin, M.; Providência, C.; Raduta, A.R.; Gulminelli, F.; Zdunik, J.L.; Haensel, P.; Bejger, M. Neutronstar radii and crusts: Uncertainties and unified equations of state.
Phys. Rev. C , , 035804,[arXiv:astro-ph.SR/1604.01944]. doi:10.1103/PhysRevC.94.035804.57. Suwa, Y.; Yoshida, T.; Shibata, M.; Umeda, H.; Takahashi, K. On the minimum mass of neutron stars. Mon. Not.R. Astron. Soc. , , 3305–3312, [arXiv:astro-ph.HE/1808.02328]. doi:10.1093/mnras/sty2460.58. Martinez, J.G.; Stovall, K.; Freire, P.C.C.; Deneva, J.S.; Jenet, F.A.; McLaughlin, M.A.; Bagchi, M.; Bates, S.D.;Ridolfi, A. Pulsar J0453+1559: A Double Neutron Star System with a Large Mass Asymmetry. Astrophys. J. , , 143, [arXiv:astro-ph.HE/1509.08805]. doi:10.1088/0004-637X/812/2/143.59. Carlin, J.B.; Melatos, A. Autocorrelations in pulsar glitch waiting times and sizes. Mon. Not. R. Astron. Soc. , , 4890–4896, [arXiv:astro-ph.HE/1907.09143]. doi:10.1093/mnras/stz2014.60. Piekarewicz, J.; Fattoyev, F.J.; Horowitz, C.J. Pulsar glitches: The crust may be enough. Phys. Rev. C , , 015803, [arXiv:nucl-th/1404.2660]. doi:10.1103/PhysRevC.90.015803.61. Martin, N.; Urban, M. Superfluid hydrodynamics in the inner crust of neutron stars. Phys. Rev. C , , 065801, [arXiv:nucl-th/1606.01126]. doi:10.1103/PhysRevC.94.065801.62. Watanabe, G.; Pethick, C.J. Superfluid Density of Neutrons in the Inner Crust of Neutron Stars:New Life for Pulsar Glitch Models. Phys. Rev. Lett. , , 062701, [arXiv:nucl-th/1704.08859].doi:10.1103/PhysRevLett.119.062701.63. Graber, V.; Cumming, A.; Andersson, N. Glitch Rises as a Test for Rapid Superfluid Coupling in Neutron Stars. Astrophys. J. , , 23, [arXiv:astro-ph.HE/1804.02706]. doi:10.3847/1538-4357/aad776.64. Haskell, B.; Khomenko, V.; Antonelli, M.; Antonopoulou, D. Crust or core? Insights from the slow rise of largeglitches in the Crab pulsar. Mon. Not. R. Astron. Soc. , , L146–L150. doi:10.1093/mnrasl/sly175.65. Gavassino, L.; Antonelli, M.; Pizzochero, P.M.; Haskell, B. A universal formula for the relativistic correctionto the mutual friction coupling time-scale in neutron stars. Mon. Not. R. Astron. Soc. , , 3562–3580,[arXiv:astro-ph.HE/2001.08951]. doi:10.1093/mnras/staa886.66. Chamel, N. Entrainment in Superfluid Neutron-Star Crusts: Hydrodynamic Description and MicroscopicOrigin. Journal of Low Temperature Physics ,189