{'fraction_non_alphanumeric': 0.032491555412578414, 'fraction_numerical': 0.014422819151788108, 'mean_word_length': 4.479435957696827, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The effect of superfluid hydrodynamics on pulsar glitch sizes and waiting times

Pulsar glitches, sudden jumps in frequency observed in many radio pulsars, may be the macroscopic manifestation of superfluid vortex avalanches on the microscopic scale. Small scale quantum mechanical simulations of vortex motion in a decelerating container have shown that such events are possible and predict power-law distributions for the size of the events, and exponential distributions for the waiting time. Despite a paucity of data, this prediction is consistent with the size and waiting time distributions of most glitching pulsars. Nevertheless a few object appear to glitch quasi-periodically, and exhibit many large glitches, while a recent study of the Crab pulsar has suggested a cut-off deviations from a power-law distribution for smaller glitches. In this paper we incorporate the results of quantum mechanical simulations in a macroscopic scale superfluid hydrodynamics simulation. We show that the effect of vortex coupling to the neutron and proton fluids in the neutron star naturally leads to deviations from power-law distributions for sizes and from exponential distributions for waiting times. In particular we predict a cut-off in the size distribution for small glitches.


INTRODUCTION
The rotation rate of radio pulsars, rotating magnetised Neutron Stars (NSs), is exquisitely stable, in some cases rivalling the stability of atomic clocks. However some pulsars undergo sudden jumps in frequency, known as 'glitches', that are instantaneous to the accuracy of the data. The origin of these events is still debated forty years after their first discovery, but it is generally thought that they are due to a large scale superfluid component in the NS interior that is only weakly coupled to the 'normal' component of the star which is tracked by the radio signal. On theoretical grounds neutrons are expected to be superfluid, as a mature NS is likely to be cold enough for most regions in the interior to be below the superfluid transition temperature (Baym et al. 1969). Superfluidity has a strong impact on the dynamics of the system, as a superfluid rotates by forming an array of quantised vortices, which mediate a dissipative interaction between the superfluid and the normal fluid; the mutual friction. Vortices, however, can also be strongly attracted, or 'pinned', to ions in the crust (Alpar 1977) or superconducting flux tubes in the core of the star (Link 2003). In this case the superfluid cannot expel vorticity and cannot spin down together with the normal component, thus lagging behind and storing angular momentum. The sudden re-coupling of the two components leads to an exchange of angular momentum and a glitch, as was suggested early on by Anderson & Itoh (1975).
Despite the success of this paradigm in explaining many qualitative features of pulsar glitches, and the bulk of work that has been devoted to studying the response of the star to a glitch (see Haskell & Melatos (2015a) for a review), it is still unclear what triggers such an event. Several mechanisms have been suggested, including starquakes (Ruderman 1969), superfluid instabilities (Andersson et al. 2003) and vortex avalanches (Cheng et al. 1988;Melatos et al. 2008). This last mechanism is based on the idea that vortices in a NS may form a Self Organised Critical (SOC) system, in which global stresses, due to classical drag forces acting on pinned vortices, are relieved locally via discrete avalanches, so that the system self-regulates and is always close to the critical threshold for unpinning (Warszawski & Melatos 2008;Warszawski et al. 2012). As a consequence of self organised criticality, the distribution of glitch sizes in a pulsar is expected to be a power law and the distribution of waiting times an exponential. This is approximately true for many glitching pulsars , although the small number of events does not allow for strong statistical conclusions. There are, however, at least three arXiv:1603.04304v1 [astro-ph.SR] 14 Mar 2016 pulsars, the best known example being the Vela pulsar, which appear to mostly have glitches of a typical size which occur quasi-periodically (Melatos et al. , 2015. This behaviour is more reminiscent of the so-called 'snowplow' mechanism (Pizzochero 2011), in which the global reservoir of angular momentum is depleted periodically once vortices can no longer be held in place by the pinning force. The periodicity and size of Vela glitches can be well reproduced in this framework, as the system has a natural length-scale and time-scale set by the height and position of the maximum of the microphysical pinning force. Additionally, a recent analysis of glitches in the Crab pulsar suggests that there may be a deviation from a power law distribution for the sizes, and that glitches may have a substantial minimum size (Espinoza et al. 2014).
In this paper we investigate how local, small scale, distributions of unpinning events are affected by large scale hydrodynamics and coupling between the fluids. In particular we will show that even if unpinning events are distributed as a power-laws on a microscopic level, the large scale glitch distribution can be substantially different, and exhibit a cutoff for small glitch sizes.

METHODS
We model the NS as a two-fluid system of superfluid neutrons and a charge neutral component consisting of the crust and electromagnetically bound protons and electrons. We thus do not consider vortex motion directly, but rather average over many vortices to consider the macroscopic motion of two dynamical degrees of freedom. Following  we can write conservation laws for each species x: where ρx is the density of constituent x (with x = n for superfluid neutrons and x = p for the locked proton,electron and crust fluid). The Euler equations are: where we assume sums over repeated spacial indices, εx is the entrainment parameter, Φ is the gravitational potential andμx = µx/mx is the chemical potential per unit mass (and we will assume mp = mn). Finally f x i is the mutual friction force, which for straight vortices takes the form: with Ω j the angular frequency of the neutron fluid and κ = h/2mn the quantum of circulation. For the vortex number per unit area nv one has the relation: wherer is the cylindrical radius and γ < 1 represents the fraction of vortices which are free (Jahan-Miri 2006), with the remaining vortices pinned. The strength of the mutual friction is parametrised by the dimensionless constants B and B . Following Haskell et al. (2012) we consider the two components to be rotating around the same axis defined byΩp, and assume that the proton fluid is rigidly rotating on the timescales of interest. This assumption is justified on longer post-glitch timescales, on which the crust can be considered as a solid and the magnetic field of the star will lock it to the protons in the core. On shorter timescales, however, several modes of oscillation of the fluids could be present (Sidery et al. 2010;van Eysden 2014), and are neglected in the current treatment. The equations of motion, averaged over vortex length (assuming straight vortices aligned with the z axis), take the form: where we have definedQ(r) = ρnγκnvB(Ωp − Ωn) andΩO is the contribution to the spin evolution from the external spin down torque acting on the star. In the following we take g(εn) = 0. This is not accurate in the crust, as entrainment will be strong and reduce the amount of angular momentum available for a glitch (Chamel 2012(Chamel , 2013?). For our purposes this is, however, simply a rescaling of the allowed size of glitches, and as we are interested in distributions and not in fitting absolute sizes, we will neglect this term for computational ease. We have also included the contribution due to pinning,Fp, which we define as: Fp = 0 for Ωn − Ωp ∆c (8) where ∆c is the critical lag for unpinning, and must be determined from microphysical calculations of pinning forces. We approximate the realistic results of Seveso et al. (2016) by taking a Gaussian profile of the form: (10) where rm is the location of the maximum, which we take to be at ρ = 7 × 10 13 g/cm 3 , and we take as typical values ∆M = 10 −2 , ∆min = 2 × 10 −4 and σ = 10 −2 (rm − Rc). We model the equation of state as an n = 1 polytrope. The timescale τ is the minimal coupling timescale that we resolve in our numerical formulation, which we fix as τ = 60s, to approximate the observational upper limits on the rise time of a glitch (Dodson et al. 2007). The external boundary of our simulations R nd is the neutron drip radius and Rc is the the radius of the crust-core interface (taken at ρ = 1.6×10 14 g/cm 3 ). Ri is the internal boundary, which we take to be the radius at which τMF = τ , with τMF = 1/(2Ωn < B >), with < B > the standard mutual friction coefficient for electronvortex scattering from , averaged over the vortex length (Haskell et al. 2012). The moment of inertia Ip in equations (5)-(6) is thus the combined moment of inertia of the outer crust, protons in the computational domain, and all components forr Ri. Finally, we fix the initial value of γB over the computational domain to γB =ΩO/2Ωn∆c, which ensures that the neutrons are spinning down together with the proton fluid, with a fixed lag close to the critical lag. This not only makes certain that the system is sub-critical from the start, but also allows us to circumvent the substantial uncertainties in microphysical estimates of B in the crust.

Random unpinning and avalanches
Evolving the system of equations in (5)-(6), as described in the previous section, will not lead to glitches. An unpinning trigger has to be added to initiate a glitch. Haskell et al. (2012) followed the so-called 'snowplow' model (Pizzochero 2011) and assumed that a vortex sheet forms close to the maximum of the critical lag ∆c. Once said lag is exceeded vortices are free to move out, leading to an increase in B (assumed to be due to Kelvin waves being excited by rapid vortex motion) and a glitch. This mechanism predicts well the sizes and waiting-times of the Vela and other pulsars that exhibit mainly giant glitches, however it cannot explain why the sizes of glitches in many other pulsars span several decades and are consistent with a power-law distribution ). On the other hand quantum mechanical Gross-Pitaevskii (GP) simulations of pinned vortices in a decelerating trap show that vortex avalanches can lead to glitches and naturally give rise to power-law distributions for their sizes (Warszawski & Melatos 2011). In this paper we take a first step towards reconciling these two approaches by taking the large scale two-fluid hydrodynamical NS model described in the previous section as a background over which to evolve small scale fluctuations in the number of pinned vortices, as predicted by GP simulations.
In practice we will follow the approach of Haskell & Antonopoulou (2014) and assume that vortex avalanches can randomly unpin vortices, increasing the unpinned fraction from γ to αγ in a region R nd − Rav < r < R nd , where both α and Rav are drawn from a power law distribution of the form where x = Rav, α; the power-law index is k, and xm = 1, xM = 10 6 . We also assume that waiting times t between events are exponentially distributed: with tw the mean waiting time between unpinning events, which we stress is not necessarily the mean waiting time between observed glitches. We thus draw a value α > 1 for the region R nd − Rav < r < R nd , and the increased coupling leads to a 'glitch', although for small values of α, the event is slow enough and weak enough that it will not appear as a sudden jump in frequency, but as a gradual increase of the spin-down rate.
Following the results of Haskell & Melatos (2015b) we assume that vortices can repin if the lag Ωn − Ωp 0.98∆C . In practice once the lag has been reduced below this limit, we set γ = 0, until the lag increases again to Ωn − Ωp = ∆C , when we set once again γ = 1.
Let us note that before running the hydrodynamical simulations we have verified that simply drawing Rav from the distribution in (11) and transferring the angular momentum of the superfluid in the region R nd − Rav < r < R nd to the crust leads to a power-law distribution for the glitches for The probability distribution function is estimated both with a histogram and a gaussian kernel density estimator. The distribution deviates significantly from a power law at lower sizes, with a notable absence of small glitches. This is due to the fact that very small events do not appear as sudden jumps, i.e. glitches, but rather as gradual changes in the spin-down rate (i.e. timing noise) and are not flagged by our glitch-finding algorithm.
which we recover the original microphysical index k, and the smallest glitch size we can resolve is ∆Ωp/Ωp = 8.7 × 10 −14 .

RESULTS
We have run a number of simulations, each comprising approximately 100 glitches, with varying ∆C , mass M , mean waiting time tw and power-law index k and analysed the results. After running the simulations we run a simple glitch finding algorithm, that identifies all events in which the spin frequency of the 'crust' Ωp rises, i.e.Ωp > 0, and measures the maximum size of the event from the start of the rise to the maximum of the frequency. The main conclusion is that, independently of the choice of parameters, the macroscopic glitch size and waiting time distributions that we extract from the simulation differ significantly from our microphysical inputs. For all cases size distributions deviate from power-law distributions at the lower end, with a strong drop off in the number of glitches for smaller sizes, as can be seen from the examples in figures 1 and 2. This is due to the fact that for small values of α, the effective mutual friction parameter γB is small, leading to a slow, more gradual event that also exchanges a small amount of angular momentum, due to the low value of Rav. This event is thus a gradual change in the spin-down rate, more similar to timing-noise, and is not identified as a glitch by our glitch finding algorithm. An example of such an event is shown in figure 3. This is consistent with observations of glitches in the Crab pulsar, which reveal deviations of the size distribution from a power-law with a lack of observed small glitches (Espinoza et al. 2014). Our results show that the effect of superfluid hydrodynamics are sizeable and lead to . Example of a spin down curve obtained from a simulation with k = −1.05, tw = 1 day (thick red line), compared to a linear fit to the spin down from previous data (dashed blue line). Two events are notable, event (a) which corresponds to a slow decrease of the spin down rate, and is not detected by the glitch finder, and event (b), which is a standard glitch.
deviations from power law distributions, even though on a microphysical scale vortex unpinning events may still be power-law distributed. Figures 1 and 2 also reveal that a steeper power-law with index k = −1.5, rather than k = −1.05, on the microphysical level leads to a macroscopic distribution of glitch sizes that appears more narrowly peaked around larger glitches sizes. This is due to the fact that a steeper power law leads to more small events that are not picked up by the glitch finder. Size distributions for k = −1.5 thus differ significantly from a power-law and show an abundance of larger glitches. Fitting a power law to the full distribution returns an index k > −1, and even cutting off the low end of the distribution returns a fit that is rejected by a Kolmogorov-Smirnov (KS) test. For k = −1.5 the waiting time distribution is also skewed towards longer waiting times compared to the k = −1.05 case, even if the microphysical choice of waiting times is the same. Furthermore for k = −1.5 larger glitches are associated with longer waiting times, as can be seen in figure 4, where we plot the waiting time versus size distribution for a microscopic waiting time of tw = 0.1 days for the two cases k = −1.5 and k = −1.05. However, we find no significant correlation between waiting times and sizes of glitches, which is consistent with the lack of any such correlation in the observed distributions of glitches in pulsars , with the sole exception of the pulsar J0537-6910, in which (Middleditch et al. 2006) suggest the existence of a correlation between the size of a glitch and the waiting time to the next event.
Mass is also a key factor in determining glitch sizes, with lower mass stars exhibiting generally larger glitches than high mass stars, as can be seen from figure (5), where we compare a M = 1.1 M and a M = 1.6 M neutron star. The waiting time distributions, on the other hand, show no strong dependance on mass.
We have also tested whether the waiting time distribution depends on changes the spin-down rateΩO ( or equivalently the height of the maximum of ∆C ). The snowplow model of Pizzochero (2011) predicts, in fact, that large glitches should occur quasi-periodically as the system builds up the maximum lag ∆ max C , on a timescale τg = ∆ max C /ΩO. We test this hypothesis by comparing the waiting time distribution for k = −1.5 and tw = 0.1 days foṙ ΩO = −6 × 10 −10 andΩO = −6 × 10 −11 . We find no strong evidence for a correlation between the spin-down rate and the glitch waiting time and we do not recover the relation predicted by Pizzochero (2011) and a KS test returns a 48 % probability of the two datasets being drawn from the same distribution. We caution the reader, however, that we have assumed that there is only a single cause for glitches, namely vortex avalanches. The 'snowplow' mechanism suggests that large glitches occur periodically when a vortex sheet forms in the crust where the pinning forces is strongest. If such a mechanism where active together with regularly occurring avalanches one may expect to see an additional periodicity and abundance of large glitches in the distributions. We intend to investigate this hypothesis in future work. In both cases the distributions differ significantly from power-laws for the sizes and exponentials for the waiting times. There is also no significant correlations between sizes and waiting times, as can be seen from the Pearson r and p coefficients. However, for k = −1.5 there is an apparent excess of larger glitches which is associated with a weak quasi-periodicity. The probability distribution function is obtained a gaussian kernel density estimator. We can see that for lower mass stars the distribution is more narrowly peaked around higher values for the glitch size.