Avalanche Dynamics of Radio Pulsar Glitches

ABSTRACT

We test statistically the hypothesis that radio pulsar glitches result from an avalanche process, in which angular momentum is transferred erratically from the flywheel-like superfluid in the star to the slowly decelerating, solid crust via spatially connected chains of local, impulsive, threshold-activated events, so that the system fluctuates around a self-organized critical state. Analysis of the glitch population (currently 285 events from 101 pulsars) demonstrates that the size distribution in individual pulsars is consistent with being scale invariant, as expected for an avalanche process. The measured power-law exponents fall in the range $-0.13\leq a\leq 2.4$ , with $a\approx 1.2$ for the youngest pulsars. The waiting-time distribution is consistent with being exponential in seven out of nine pulsars where it can be measured reliably, after adjusting for observational limits on the minimum waiting time, as for a constant-rate Poisson process. PSR J0537–6910 and PSR J0835–4510 are the exceptions; their waiting-time distributions show evidence of quasi-periodicity. In each object, stationarity requires that the rate λ equal $-\epsilon \dot{\nu }/ \langle \Delta \nu \rangle$ , where $\dot{\nu }$ is the angular acceleration of the crust, $\langle \Delta \nu \rangle$ is the mean glitch size, and $\epsilon \dot{\nu }$ is the relative angular acceleration of the crust and superfluid. Measurements yield $\epsilon \leq 7\times 10^{-5}$ for PSR J0358+5413 and $\epsilon \leq 1$ (trivially) for the other eight objects, which have a< 2

. There is no evidence that λ changes monotonically with spin-down age. The rate distribution itself is fitted reasonably well by an exponential for $\lambda \geq 0.25$ yr^–1, with $\langle \lambda \rangle =1.3^{+0.7}_{-0.6}$ yr^–1. For $\lambda < 0.25$ yr^–1 the exact form is unknown; the exponential overestimates the number of glitching pulsars observed at low λ, where the limited total observation time exercises a selection bias. In order to reproduce the aggregate waiting-time distribution of the glitch population as a whole, the fraction of pulsars with $\lambda > 0.25$ yr^–1 must exceed ∼70%.

1. INTRODUCTION

Glitches are tiny, impulsive, randomly timed increases in the spin frequency ν of a rotation-powered pulsar, sometimes accompanied by an impulsive change in the frequency derivative $\dot{\nu }$ . They are to be distinguished from timing noise, a type of rotational irregularity where pulse arrival times wander continuously, although there is evidence that timing noise is the cumulative result of frequent microglitches in certain pulsars (Cordes & Downs 1985; D’Alessandro et al. 1995).

At the time of writing, 285 glitches in total have been detected in 101 objects (∼6% of the known radio pulsar population), the majority in the last four years, as facilitated by the Parkes Multibeam Survey, refined multifrequency ephemerides, and better interference rejection algorithms (Hobbs 2002; Krawczyk et al. 2003; M. Kramer & A. Lyne 2005, private communication; D. Lewis 2005, private communication; Janssen & Stappers 2006). Efforts to analyze the data statistically have focused on the correlation of glitch activity with age (McKenna & Lyne 1990; Shemar & Lyne 1996; Urama & Okeke 1999; Lyne et al. 2000; Wang et al. 2000) and Reynolds number (Peralta 2006; Melatos & Peralta 2007), the postglitch relaxation timescale (Wang et al. 2000; Wong et al. 2001), the size distribution (Morley & García-Pelayo 1993a, 1993b; Peralta 2006), and the correlation between glitch sizes and waiting times (Wang et al. 2000; Wong et al. 2001; Middleditch et al. 2006; Peralta 2006). Hobbs (2002) reviewed the role of observational selection effects.

Most glitching pulsars (65%) have been seen to glitch once, but a minority glitch repeatedly; the current record holder is PSR J1740–3015, with 33 glitches. Of those objects which glitch repeatedly, most do so at unpredictable intervals, but two (PSR J0537–6910 and Vela) are quasi-periodic; Vela, in particular, has been likened to a relaxation oscillator (Lyne et al. 1996). The fractional increase in ν spans seven decades ( $3\times 10^{-11}\leq \Delta \nu / \nu \leq 2\times 10^{-4}$ ) across the glitch population and as many as four decades in a single object (e.g., $7\times 10^{-10}\leq \Delta \nu / \nu \leq 2\times 10^{-6}$ in PSR J1740–3015). The spin-down age $\tau _{c}=-\nu / ( 2\dot{\nu })$ of glitching pulsars spans four decades, from $1\times 10^{3}$ to $3\times 10^{7}$ yr. In many respects, therefore, the glitch phenomenon is scale invariant. This striking property invites physical interpretation.

Theories of pulsar glitches have focused mainly on the local microphysics of the superfluid in the stellar interior and its coupling to the solid crust, for example the strength of vortex pinning (Anderson & Itoh 1975; Jones 1998), the rate of vortex creep (Link & Epstein 1996), or the conditions for exciting superfluid turbulence (Peralta et al. 2005, 2006; Melatos & Peralta 2007; Andersson et al. 2007). Ultimately, however, the local microphysics must be synthesized with the global, collective dynamics in order to make full contact with observational data. (Likewise, a practical model of earthquakes must synthesize the microphysics of rock fracture with the macrodynamics of interacting tectonic plates.) For example, if approximately $10^{16}( \Delta \nu / 1\ \mathrm{Hz}\,)$ vortices unpin from crustal lattice sites in sympathy during a glitch, they must communicate rapidly across distances much greater than their separation. How? And why does the number that unpin fluctuate so dramatically (by up to 4 orders of magnitude) from glitch to glitch in a single pulsar, while always amounting to a small fraction ( $\Delta \nu / \nu$ ) of the total?

Such collective, scale-invariant behavior is a generic feature of a class of natural and synthetic far-from-equilibrium systems, called self-organized critical systems, that are discrete, interaction dominated, and slowly driven, and that adjust internally via erratic, spatially connected avalanches of local, impulsive, threshold-activated, relaxation events (Jensen 1998). Such systems fluctuate around a stationary state toward which they evolve spontaneously, in which global driving balances local relaxation on average over the long term. The archetype of a self-organized critical system is the sandpile (Bak et al. 1987).

In this paper, we study pulsar glitches as an avalanche process, as first proposed by Morley & García-Pelayo (1993a). After reviewing self-organized criticality in § 2, we define the statistical sample on which our study is based (§ 3) and analyze the observed distribution of glitch sizes (§ 4) and waiting times (§ 5). Some implications for glitch physics are explored in § 6. We only include radio pulsars in the sample, to preserve its homogeneity, even though glitches have now been observed in anomalous X-ray pulsars (magnetars) as well (Dall’Osso et al. 2003; Kaspi & Gavriil 2003).

2. AVALANCHE DYNAMICS

A system in a self-organized critical state exhibits the following distinguishing features (Jensen 1998).

It is composed of many discrete, mutually interacting elements, whose motions are dominated by local (e.g., nearest-neighbor) rather than global (e.g., mean field) forces.3
Each element moves when the local force exceeds a threshold (stick-slip motion). Hence stress accumulates sustainedly at certain random locations, while relaxing quickly elsewhere; at any instant, the system houses numerous metastable stress reservoirs, separated by relaxed zones.
An external force drives the system slowly, in the sense that elements adjust to local forces rapidly compared to the driver timescale. Combined with local thresholds, this ensures that the system evolves quasi-statically through a history-dependent sequence of metastable states (a huge number of which are available).
Transitions from one metastable state to the next occur via avalanches: spatially connected chains of local equilibration events, in which one element relaxes and redistributes some local stress to its neighbors, which in turn can exceed their thresholds and relax (knock-on effect). The duration of even the largest avalanches is short compared to the driving timescale (see previous point).
Avalanches have no preferred scale: they can involve a few (commonly) or all (rarely) of the elements in the system. Their sizes and lifetimes follow power-law distributions, whose exponents are related. The numerical values of the exponents depend on the spatial dimensionality of the system, the spatial symmetries of the local forces and redistributive channels, the strength of the local forces (Field et al. 1995), and the level of conservation (Olami et al. 1992).4
Over the long term, the system tends toward a critical state, which is stationary on average but not instantaneously. For example, on average, the power input by the external driver equals the energy per unit time released by avalanches. But there are fluctuations, because, at any instant, a random amount of energy is stored in metastable local reservoirs.

Avalanche dynamics are generically observed in nature when conditions (1)-(3) are met, and properties (4)-(6) emerge irrespective of the detailed microphysics (Jensen 1998). Likewise, in this paper, we remain agnostic about the microphysics of pulsar glitches; the statistical analysis presented below makes no assumptions in this regard. Nevertheless, it is striking that the traditional glitch paradigm—a collective unpinning of quantized superfluid vortices interacting with an inhomogeneous, slowly decelerating crust—conforms closely with conditions (1)-(6) (Anderson & Itoh 1975; Alpar et al. 1996). So too does an alternative paradigm, based on crust fracture (Alpar et al. 1996; Middleditch et al. 2006), whose terrestrial counterpart (plate tectonics) is renowned as an archetype of self-organized criticality (Sornette et al. 1991). We elucidate the analogy briefly before continuing.

Consider a rectilinear array of quantized vortices, each carrying circulation κ, spaced evenly according to Feynman’s rule ( $4\pi \nu / \kappa$ vortices per unit area) in the neutron superfluid permeating the inner crust of a neutron star. A small percentage of the vortices are pinned to defects and/or nuclei at random locations in the crustal lattice, clustered to varying degrees (Alpar et al. 1996; Wong et al. 2001). As ν decreases gradually due to electromagnetic spin down, most vortices move apart, and the outermost ones are expelled. However, the pinned vortices stay (nearly) fixed, in metastable reservoirs separated by relaxed zones (see condition [2]), creeping slowly between adjacent pinning sites in response to thermal fluctuations (Link et al. 1993). The reservoirs are identical to the capacitive elements (vortex traps surrounded by vortex depletion regions) postulated by Alpar et al. (1996) and Wong et al. (2001). They may be seeded by starquakes, which create large numbers of fresh lattice dislocations with deep pinning potentials, or they may emerge spontaneously in the self-organized critical state, as successive generations of vortex avalanches traverse the crust. As the pinned vortices increasingly lag the regular, unpinned array, a gradient in vortex density is established, and the local Magnus force on a pinned vortex rises. When the pinning threshold is overcome, a pinned vortex unpins and moves abruptly away from the pinning site (stick-slip motion), disturbing the local superfluid velocity field (and hence the Magnus force) appreciably. Often, this is enough to push neighboring, barely subcritical, pinned vortices over their thresholds, triggering an avalanche. Most vortices in the avalanche rejoin the regular, unpinned array, and the crust spins up proportionately to compensate. The timescale for a vortex to adjust locally to the Magnus and pinning forces is much shorter than $\nu / \dot{\nu }$ , in keeping with condition (3).

Classic laboratory experiments on magnetic flux vortices in a type II superconductor (e.g., NbTi) immersed in a slowly changing magnetic field, an exactly analogous system, clearly exhibit properties (4)-(6) (Field et al. 1995). Vortices are expelled mostly in a continuous flow (cf. steady spin down) and occasionally in avalanches (cf. glitches). The distribution of avalanche sizes is measured to be a power law over several decades, whose exponent depends on the strength of the applied magnetic field (which controls the vortex spacing and hence the strength of the vortex-vortex interaction). The temporal fluctuation spectrum scales as an inverse power of frequency at high frequencies. After initial transients die away, the superconductor fluctuates around a self-organized critical state, called the Bean state, where the Lorentz force acting on each vortex is everywhere equal to the maximum pinning force.

If the pinning sites are sparsely distributed, so that global (mean field) forces dominate local forces between pinned vortex clusters, scale invariance breaks down (Jensen 1998). Avalanches still occur, but they are distributed narrowly around a characteristic size and lifetime, involving nearly all the vortices instead of small, independent subsets. In this regime, avalanches recur quasi-periodically, not stochastically. Similar behavior is observed when the external driver acts too rapidly, but this situation never arises in pulsars.

Scale-invariant avalanche dynamics and self-organized critical states are observed widely elsewhere, e.g., in sandpiles (Bak et al. 1987), earthquakes (Sornette et al. 1991), solar flares (Lu & Hamilton 1991; Wheatland 2000), and bursts from soft-gamma-ray repeaters (Göğüş et al. 2000). The analogy with pulsar glitches has been pointed out by Morley & García-Pelayo (1993a) and Carroll (1998), and modeled using a cellular automaton by Morley & Schmidt (1996).

3 Tectonic plates, or grains of sand in a pile, are terrestrial examples of interacting elements.

4 In this respect, far-from-equilibrium critical systems differ from equilibrium critical systems (e.g., second-order phase transition in a ferromagnet), whose exponents depend only on the dimensionality of the system and its order parameter(s).

3. DATA

Table 1 lists all 285 glitches discovered up to the time of writing and known to the authors. It is compiled from published sources (Shemar & Lyne 1996; Lyne et al. 2000; Wang et al. 2000; Hobbs 2002; Krawczyk et al. 2003; Janssen & Stappers 2006; Middleditch et al. 2006; Peralta 2006), the Australia Telescope National Facility Pulsar Catalogue (Manchester et al. 2005), which can be accessed online,5 and unpublished data communicated privately by M. Kramer, D. Lewis, and A. G. Lyne (2005). For each pulsar, the table lists its J2000.0 coordinates and the number of glitches detected ( $N_{g}$ ). The earliest and latest epochs observed ( $t_{\mathrm{min}\,}$ and $t_{\mathrm{max}\,}$ , respectively) are recorded separately in Table 2 for the nine pulsars with $N_{g}> 5$ . A footnote signifies that segmented data spans are not specified in the cited references; in this situation, $t_{\mathrm{min}\,}$ and $t_{\mathrm{max}\,}$ are estimated by eye from spin-down histories graphed in the cited references, where available, or else from the first and last glitches by default. For each glitch, Table 1 lists its epoch, the fractional increase in its spin frequency $\Delta \nu / \nu$ , and one or more bibliographic references. Uncertainties are quoted as a trailing integer in parentheses, corresponding to an absolute number of days for t [e.g., MJD 51141(248) means MJD $51141\pm 248$ ] and an uncertainty in the last significant digit for $\Delta \nu / \nu$ [e.g., 0.04(2) means $0.04\pm 0.02$ ]. For some newly discovered glitches, the information is incomplete. Epochs and sizes have been measured for 271 and 250 glitches, respectively. Other parameters, such as the healing fraction and postglitch relaxation timescale, are omitted, as they are not analyzed in this paper; please consult Peralta (2006) and references therein for a full catalog.

TABLE 1

Parameters of Pulsar Glitches

TABLE 2

Observing-Time Intervals for Pulsars with N_g > 5

5 See http://www.atnf.csiro.au/research/pulsar/psrcat.

4. SIZE DISTRIBUTION

4.1. Scale Invariance

If pulsar glitches are the result of an avalanche process, their size distribution should be scale invariant in any individual pulsar, with a probability density function

The exponent a is set by the dimensionality6 and symmetries of the local forces, which are likely to be universal, and the strength and level of conservation of these forces, which are functions of temperature and therefore not universal (see § 2). One therefore expects a to differ from pulsar to pulsar. As a corollary, the aggregate size distribution drawn from all pulsars is not expected to be a simple power law of the form of equation (1).

To test these ideas, we construct the observed cumulative size distributions of the nine known pulsars with $N_{g}> 5$ . The selection criterion $N_{g}> 5$ is arbitrary; it seeks to limit the impact of random errors while testing as many objects from Table 1 as possible. We then compare the data against the theoretical cumulative distribution

derived from equation (1). The theoretical distribution is normalized after restricting it to the domain $( \Delta \nu / \nu ) _{\mathrm{min}\,}\leq \Delta \nu / \nu \leq ( \Delta \nu / \nu ) _{\mathrm{max}\,}$ , where $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ are the smallest and largest glitches observed in that pulsar respectively, as quoted in Table 1. There are more sophisticated ways to choose $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ , which we consider further below, but this is a conservative starting point.

For each object, we choose a to minimize the Kolmogorov-Smirnov (K-S) statistic D, i.e., the maximum unsigned distance between the curves. The numerical results are recorded in Table 3, while the measured and theoretical cumulative distributions are plotted together in Figure 1. (Cumulative distributions are free of binning bias.) The goodness of the fit at the optimal value of a is characterized by $P_{\mathrm{K}\,- \mathrm{S}\,}$ , which is defined such that $1-P_{\mathrm{K}\,- \mathrm{S}\,}$ equals the probability that the K-S null hypothesis (that the two data sets are drawn from the same underlying distribution) is false.7 The 1 σ lower and upper bounds, $a_{-}$ and $a_{+}$ , mark the range of a where the null hypothesis is rejected with less than 68% confidence. Note that the interval $[ a_{-},a_{+}]$ is asymmetric about the optimal a and widest for the best fits.

TABLE 3

Power-Law Size Distribution Parameters for Pulsars with N_g > 5

Figure 1.

Cumulative distributions of fractional glitch sizes $\Delta \nu / \nu$ in the nine pulsars that have glitched more than five times. The observational data (histograms) are plotted together with the best power-law fits (solid curves) given by eq. (2), where $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ are taken from Table 1 and a is taken from Table 3.

The results in Table 3 confirm what is apparent by eye from Figure 1: the null hypothesis that the size distribution is described by a power law for all nine pulsars with $N_{g}> 5$ is not ruled out at the 1 σ level of confidence. In turn, this is consistent with the avalanche hypothesis. However, in two objects, namely PSR J0537–6910 and PSR J0835–4510, the agreement is marginal. Interestingly, these two objects are also the only ones discovered so far that are believed to glitch quasi-periodically (Lyne et al. 1996; Middleditch et al. 2006).

4.2. PSR J0537–6910 and PSR J0835–4510

Quasi-periodicity is a natural feature of avalanche dynamics when mean field forces overwhelm local interactions, as described in § 2. We explore its manifestation in glitch waiting times in § 5. With respect to glitch sizes, we note that avalanches in the quasi-periodic regime tend to be distributed narrowly around a characteristic size (Jensen 1998). This can be modeled crudely by adding a term proportional to $\delta \left[\Delta \nu / \nu -( \Delta \nu / \nu ) _{c}\right]$ to equation (1), viz.,

where $( \Delta \nu / \nu ) _{c}$ denotes the characteristic size, and the scale-invariant and quasi-periodic components are weighted by the constants $C_{s}$ and $C_{p}$ , respectively. Normalization fixes $C_{s}$ in terms of $C_{p}$ (or vice versa), with $C_{s}+C_{p}\neq 1$ in general. The associated cumulative size distribution is given by

where $H( \cdot )$ denotes the Heaviside step function.

Parameters determined by fitting equation (4) to the data are recorded in Table 4, while the corresponding measured and theoretical cumulative distributions are plotted together in Figure 2. The fits are much improved, with $C_{p}\approx 0.2$ in both objects, although, to be fair, the delta-distributed component is not strictly required, at least not at the 1 σ level. Importantly, the delta-distributed component contains only ∼20% of the glitches, not all of them. This is consistent with the historical interpretation of the pulsar data (Lyne et al. 1995; Marshall et al. 2004). It is also seen in self-organized critical systems such as sandpiles, where large, system-spanning, quasi-periodic avalanches of a characteristic size are interspersed with small, randomly timed avalanches, which are power-law distributed (Rosendahl et al. 1993; Jensen 1998).

TABLE 4

Two-Component Size Distribution Parameters for Quasi-periodic Glitchers

Figure 2.

Cumulative distributions of fractional glitch sizes $\Delta \nu / \nu$ for the two pulsars that have a quasi-periodic component. The observational data (histograms) are plotted together with the best two-component fits (solid curves) given by eq. (4), where $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ are taken from Table 1, and a, $C_{p}$ , and $( \Delta \nu / \nu ) _{c}$ are taken from Table 4.

4.3. Upper and Lower Cutoffs

Strictly speaking, it is incorrect to normalize $P( \Delta \nu / \nu )$ by choosing $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ to be the smallest and largest glitches observed in a pulsar, respectively. A better choice of $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ is the actual resolution of the timing experiment, which varies with object and epoch. Janssen & Stappers (2006) simulated detection of a microglitch in a noisy time series and obtained $( \Delta \nu / \nu ) _{\mathrm{min}\,}=1\times 10^{-11}$ . Usually, this information is not provided explicitly and must be estimated from the size uncertainties quoted for detected glitches. All the same, the smallest glitch observed is likely to be a reasonable estimate of $\Delta \nu _{\mathrm{min}\,}$ , because the occurrence probability increases steeply as $\Delta \nu$ decreases, according to Table 3. On the other hand, $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ is limited by the total observing time. Its true value exceeds the largest glitch observed, but not by much, because the occurrence probability decreases steeply as $\Delta \nu$ increases.

To quantify these effects, we allow $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ to vary between 0.5 and 1.0 times the smallest glitch size observed, $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ to vary between 1.0 and 2.0 times the largest glitch size observed, and fit equation (4) again to the data. For every object, $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ shift only slightly, and a stays within the range $[ a_{-},a_{+}]$ in Tables 3 and 4. This confirms that the smallest and largest glitches provide reasonable estimates of $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ . At the 1 σ level, the constrained and unconstrained fits are both consistent with the data.

4.4. Aggregate Distribution

Figure 3 displays the cumulative size distribution for the glitch population in aggregate, together with the best power-law fit of the form of equation (2). The fit is poor. When all 250 glitches with measured sizes are included, the best fit corresponds to a=0.96

, $( \Delta \nu / \nu ) _{\mathrm{min}\,}=9.5\times 10^{-12}$ , $( \Delta \nu / \nu ) _{\mathrm{max}\,}=2.0\times 10^{-5}$ , and $P_{\mathrm{K}\,- \mathrm{S}\,}=7.1\times 10^{-4}$ . When the glitches from PSR J0537–6910 and PSR J0835–4510 are excluded, the best fit corresponds to a=0.98

and $P_{\mathrm{K}\,- \mathrm{S}\,}=3.2\times 10^{-4}$ , with $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ as before. Either way, we can state confidently that the theoretical and observed data are drawn from different underlying distributions. This is not surprising; the results in Figure 1 and Table 3 demonstrate clearly that the size distribution in individual pulsars is consistent with being scale invariant, but that a differs from object to object. Hence the aggregate distribution is expected to be a weighted sum of power laws, not a pure power law. Accordingly, the size distribution in individual pulsars is a more direct probe of glitch physics than the aggregate distribution (Lyne et al. 2000). The aggregate distribution can be inverted, in principle, to determine how a is distributed across the pulsar population. We defer this exercise until better historical estimates of $( \Delta \nu / \nu ) _{\mathrm{min}\,}$ and $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ , as well as more data, become available.

Figure 3.

Aggregate cumulative distributions of fractional glitch sizes $\Delta \nu / \nu$ for all glitching pulsars (left) and for all glitching pulsars except PSR J0537–6910 and PSR J0835–4510, which have a quasi-periodic component (right). The observational data (histograms) are plotted together with the best power-law fits (solid curves) given by eq. (2), for the best-fit parameters in § 4.4.

Janssen & Stappers (2006) claimed that the glitches in PSR J1740–3015 are drawn from a flat size distribution in $\mathrm{log}\,( \Delta \nu / \nu )$ , i.e., a=1

, with $P_{\mathrm{K}\,- \mathrm{S}\,}=0.902$ . This agrees with the results in Table 3.

Lyne et al. (2000) noted some evidence for an excess of large glitches, which is corroborated to some extent by Figure 3a. However, the excess largely disappears when the quasi-periodic glitchers are excluded, as in Figure 3b. Large glitches do not originate preferentially from any particular class of object. While it is true that the most active objects (e.g., PSR J0537–6910 and PSR J0835–4510) experience relatively large and narrowly distributed glitches, with $\Delta \nu / \nu > 10^{-7}$ , other active objects (e.g., PSR J1740–3015) experience a mix of small and large events, and there are several objects (e.g., PSR J1806–4212) that have only glitched once, with $\Delta \nu / \nu > 10^{-5}$ for that single glitch. Furthermore, although PSR J0534+2200 is sometimes portrayed as unusual for not experiencing large glitches, its size distribution is actually relatively flat ( $a\approx 1.2$ ). There is every reason to expect that it will experience large glitches in the future, but these events will be slow in coming, because PSR J0534+2200 builds up differential rotation between the crust and superfluid at a relatively slow rate, as we show in § 5.8

6 The effective dimension need not equal 3. For example, it may equal 2 in a rectilinear vortex array or faulting crust, where the local forces act in the transverse plane, and 3 in a turbulent vortex tangle (Peralta et al. 2005, 2006).

7 The K-S test is most sensitive to discrepancies near the median. An alternative test, based on Kuiper’s statistic, mitigates this bias (Press et al. 1986). It will be worth implementing when more data become available. The K-S test is also inefficient if the underlying probability density function contains a narrow notch, where the probability vanishes. Again, there are insufficient data at hand to look for such a notch; it is difficult to find in a cumulative distribution, and the probability density function is biased by binning when $N_{g}$ is small.

8 Wong et al. (2001) argued that the relative angular acceleration of the crust and superfluid in the Crab, inferred from the activity parameter, is much smaller than expected, given the large $\Delta \dot{\nu }/ \dot{\nu }\sim 10^{-4}$ observed during glitches. This paradox is resolved if most of the differential rotation is being stored temporarily, in advance of a large glitch in the future.

5. WAITING-TIME DISTRIBUTION

5.1. Poisson Process

If pulsar glitches are the result of an avalanche process, they should be statistically independent events. To understand why, recall that a system in a self-organized critical state configures itself into many metastable stress reservoirs insulated by relaxed zones (§ 2). Every avalanche relaxes one reservoir, typically occupying a small fraction of the system, and the next avalanche occurs at random, typically far from its predecessor. There is essentially no interference between successive avalanches; this is verified empirically in tests with cellular automata (Jensen 1998). Avalanches in the tail of the size distribution, which relax the whole system, are an (extremely rare) exception.

Given statistical independence, and assuming that the system is driven at a constant (mean) rate, the avalanche model predicts that the time between successive glitches, $\Delta t$ , termed the waiting time, obeys Poisson statistics (Jensen 1998; Wheatland 2000). Hence, in any individual pulsar, the waiting-time distribution is exponential, with a probability density function

The mean glitching rate λ is different for every pulsar. It depends on the rate at which differential rotation builds up between the superfluid and the crust (∝ $\dot{\nu }$ ), as well as the capacity to store the differential rotation (e.g., strength of pinning, rate of vortex creep, shear modulus of the crust). The storage capacity is presumably controlled by thermodynamic variables such as temperature, as well as the inhomogeneous nuclear structure of the crust. We do not expect λ to change appreciably during four decades of pulsar timing. In principle, however, as more data are collected in future, this claim can be tested by using a Bayesian blocks algorithm to divide the time series into a sequence of Poisson processes with piecewise-constant rates (Scargle 1998; Wheatland 2000; Connors & Carramiñana 2003).9

The avalanche model makes a further powerful prediction. Suppose the system tends toward a stationary, self-organized critical state, in which global driving balances local release in a time-averaged sense (property [6], § 2); i.e., there is no secular accumulation or leakage of stress. Stationarity implies that the mean waiting time $\langle \Delta t\rangle =\lambda ^{-1}$ , multiplied by the rate at which crust-superfluid differential rotation builds up ( $\epsilon \dot{\nu }$ ), equals the mean glitch size $\langle \Delta \nu \rangle$ , i.e.,

Here, $2\pi \epsilon \dot{\nu }$ is the relative angular acceleration of the crust and superfluid. Importantly, glitch data allow ϵ to be measured directly in principle (Wong et al. 2001). However, there is a serious question as to whether stationarity is achieved in practice, during the 40 yr that a typical pulsar has been observed. We discuss this issue further below.

To test the above ideas, we compare the observed cumulative waiting-time distributions of the same nine pulsars as in § 4, with $N_{g}> 5$ , against the theoretical cumulative distribution. In order to make the comparison fairly, we must first adjust for observational selection effects. Any given observation can detect waiting times in a range $\Delta t_{\mathrm{min}\,}\leq \Delta t\leq \Delta t_{\mathrm{max}\,}$ . The upper limit $\Delta t_{\mathrm{max}\,}$ is set by the total data span available for that pulsar, i.e., $\Delta t_{\mathrm{max}\,}=t_{\mathrm{max}\,}-t_{\mathrm{min}\,}$ . The lower limit $\Delta t_{\mathrm{min}\,}$ is different at different epochs. It is set by the gap between data spans in which a glitch is localized. For small glitches, the glitch epoch is determined by requiring continuity of pulse phase across the glitch. For larger glitches, where the phase winding number is ambiguous, the epoch is taken to be halfway between the bounding observations (Wang et al. 2000). Either way, $\Delta t_{\mathrm{min}\,}$ is different for each glitch, and is twice the absolute value of the epochal uncertainty quoted in Table 1 (Lyne et al. 2000; Wang et al. 2000; Janssen & Stappers 2006; Middleditch et al. 2006). Let $f( \Delta t_{\mathrm{min}\,}) d( \Delta t_{\mathrm{min}\,})$ be the observing-time-weighted probability that, when a glitch occurs, $\Delta t_{\mathrm{min}\,}$ lies in the range $\left[\Delta t_{\mathrm{min}\,},\Delta t_{\mathrm{min}\,}+d( \Delta t_{\mathrm{min}\,}) \right]$ , and let the smallest and largest values of $\Delta t_{\mathrm{min}\,}$ be $\Delta t^{( < ) }_{\mathrm{min}\,}$ and $\Delta t^{( > ) }_{\mathrm{min}\,}$ , respectively, as recorded in Table 5 for the nine pulsars with $N_{g}> 5$ . Then the cumulative waiting-time distribution is given by

where each glitch is weighted equally in the sum in equation (8) as a first approximation.

TABLE 5

Poissonian Waiting-Time Distribution Parameters for Pulsars with N_g > 5

In Figure 4, we plot as cumulative histograms the measured waiting-time distributions of the nine pulsars in Figure 1. The theoretical curves (eq. [8]) are overlaid, with $\Delta t_{\mathrm{min}\,}$ and $\Delta t_{\mathrm{max}\,}$ chosen according to the second through fourth columns in Table 5. For each object, we choose λ to minimize the K-S statistic D. The fitting parameters are displayed in the fifth through eighth columns in Table 5. As in § 4, the goodness of the fit is characterized by the K-S probability $P_{\mathrm{K}\,- \mathrm{S}\,}$ , with $P_{\mathrm{K}\,- \mathrm{S}\,}> 0.32$ in the interval $\left[\lambda _{-},\lambda _{+}\right]$ .

Figure 4.

Cumulative distribution of glitch waiting times $\Delta t$ (measured in yr) for the nine pulsars that have glitched more than five times. The observational data (histograms) are plotted together with the best Poisson fits (solid curves) given by eq. (8), where $\Delta t_{\mathrm{min}\,}$ is taken from Table 1 (twice the epochal uncertainty), $\Delta t_{\mathrm{max}\,}$ from Table 2, and λ from Table 6.

For all nine pulsars in Figure 4 and Table 5, the null hypothesis that the waiting-time distribution is described by Poisson statistics is not ruled out at the 1 σ level. The data are therefore consistent with an avalanche process. An exponential waiting-time distribution was first postulated by Wong et al. (2001) for PSR J0534+2200, based on timing data up to and including the glitch on MJD 51,452. These authors obtained $\lambda =0.53$ yr^–1 and $P_{\mathrm{K}\,- \mathrm{S}\,}=0.7$ , which is marginally outside the 1 σ range in Table 5. The data analyzed here confirm that waiting times are consistent with Poisson statistics in several glitching pulsars, affording a key insight into the physics of the glitch mechanism. The implications of this result are discussed further in § 6.

5.2. Quasi-Periodicity

For seven pulsars in Figure 4, the Poisson distribution affords an excellent fit, both formally and by eye. However, for PSR J0537–6910 and PSR J0835–4510, the fits are marginal at the 1 σ level. (Indeed, in an earlier analysis of PSR J0835–4510, Wong et al. [2001] excluded a Poisson distribution with 96% confidence, on the basis of fewer data.) These are the same two objects whose size distributions are exceptional, and which are observed to glitch quasi-periodically.

Taking the same approach as in § 4, we model the quasi-periodicity crudely by adding a periodic component to equation (5), viz.,

In equation (9), $C^{\prime }_{s}$ and $C^{\prime }_{p}$ are the normalized relative weights of the Poisson and periodic components, respectively, and $\Delta t_{c}$ is the characteristic period. The associated cumulative distribution, weighted by $\Delta t_{\mathrm{min}\,}$ , is obtained by substituting equation (9) into equation (7), yielding

The two-component model (eq. [9]) yields improved fits to the data, with $C^{\prime }_{p}\approx 0.25$ for both objects. The fits are graphed with the data in Figure 5, and the best-fit parameters are recorded in Table 6. We obtain $\Delta t_{c}=( 0.3\pm 0.1)$ yr and $\Delta t_{c}=( 2.8\pm 0.1)$ yr for PSR J0537–6910 and PSR J0835–4510, respectively, in accord with previous authors (Lyne et al. 1996; Middleditch et al. 2006). Significantly, the data imply $C_{p}\approx C^{\prime }_{p}$ . In other words, the delta-distributed component accounts for the same fraction of the size and waiting-time distributions, even though the sizes and waiting times are statistically independent. This raises confidence in the model, and suggests that a quasi-periodic component is indeed present and distinct. It also suggests that the quasi-periodic component coexists with the Poisson component, instead of completely displacing it. Vela, for example, is likely to possess an extensively connected network of capacitive elements, but it also contains smaller subnetworks that are disconnected from the main network; cf. Alpar et al. (1996). This is natural for an avalanche process, as noted in § 4.2 (Rosendahl et al. 1993; Jensen 1998).

Figure 5.

Cumulative distribution of glitch waiting times $\Delta t$ (measured in yr) for the two pulsars that have a quasi-periodic component. The observational data (histograms) are plotted together with the best two-component fit (solid curves), given by eq. (10), where $\Delta t_{\mathrm{min}\,}$ is taken from Table 1 (twice the epochal uncertainty), $\Delta t_{\mathrm{max}\,}$ from Table 2, and λ, $C^{\prime }_{p}$ , and $\Delta t_{c}$ from Table 6.

TABLE 6

Two-Component Waiting-Time Distribution Parameters for Quasi-periodic Glitchers

5.3. Mean Rate

Stationarity of the avalanche model over long time intervals implies a relation between the Poisson rate, driving rate, and mean glitch size, given by equation (6). Unfortunately, for a< 2

, $\langle \Delta \nu \rangle$ is dominated by large glitches near the upper cutoff of $P( \Delta \nu / \nu )$ :10

In equation (11), $\Delta \nu _{\mathrm{lower}\,}$ and $\Delta \nu _{\mathrm{upper}\,}$ are the physical lower and upper cutoffs of the probability distribution function (eq. [1]). As large glitches are rare, stationarity is not achieved during the 40 yr interval over which a typical pulsar is observed; the largest observed size, $( \Delta \nu / \nu ) _{\mathrm{max}\,}$ , cannot be equated reliably with the maximum size allowed physically. Likewise, $\langle \Delta \nu \rangle$ is approximated poorly by the mean of the observed glitches. In practice, therefore, it is impossible to estimate $\langle \Delta \nu \rangle$ without much longer monitoring.

Physically, $\epsilon \dot{\nu }$ is the rate at which differential rotation builds up between the crust and superfluid. Hence, in the vortex-unpinning model, ϵ gives the time-averaged fraction of pinned vortices or capacitive elements. We can use equation (11) to place limits on ϵ, at least in principle.11 For example, the inequalities $\Delta \nu _{\mathrm{lower}\,}< \Delta \nu _{\mathrm{min}\,}$ and $\Delta \nu _{\mathrm{max}\,}< \Delta \nu _{\mathrm{upper}\,}< \nu$ must always be satisfied. For PSR J0358+5413, assuming that a=2.4

, we find $\langle \Delta \nu \rangle / \nu \leq 1.1\times 10^{-10}$ , and hence $\epsilon \leq 7\times 10^{-5}$ . This is lower than previous estimates of the pinned fraction for objects of that age, but in line with previous estimates of the pinned fraction in young objects such as the Crab (Lyne et al. 2000; Wong et al. 2001). For the remaining eight objects, the above inequalities lead to upper limits on ϵ that are greater than unity, and hence not useful. As a crude experiment, we check the result of setting $\Delta \nu _{\mathrm{upper}\,}=2\times 10^{-4}$ , the largest glitch observed in any pulsar over the last 40 years, for every pulsar. We obtain five slightly more useful upper limits ( $\epsilon \leq 0.04$ , 0.2, 0.8, 0.8, and 0.7 for PSR J0534+2200, PSR J0631+1036, PSR J0835–4510, PSR J1341–6220, and PSR J1740–3015, respectively). However, we emphasize that these values are still problematic, because there is no guarantee that a total effective observation interval of 40 yr × 101 pulsars is long enough in aggregate for stationarity to be observed. Moreover, these values are based on the assumption that all pulsars have the same physical $\Delta \nu _{\mathrm{upper}\,}$ , which is not necessarily the case.

Figure 6 displays the cumulative histogram of Poisson rates derived from the best-fit waiting-time distributions in Figures 4 and 5. Let $q( \lambda )$ denote the rate probability density function, such that $q( \lambda ) d\lambda$ is the probability that the mean rate lies in the interval $[ \lambda ,\lambda +d\lambda ]$ in a given pulsar. There is no obvious theoretical reason to prefer a particular analytic form of $q( \lambda )$ , which is controlled by the physics of the global driver, not the scale-invariant avalanche dynamics. In addition, the data in Figure 6 are insufficient to specify the analytic form of $q( \lambda )$ uniquely. However, motivated by the rate distribution observed in solar flares, which is measured reliably to be exponential (Wheatland 2000), we find that a distribution of the form

fits the data satisfactorily, with $\langle \lambda \rangle =1.3^{+0.7}_{-0.6}$ yr^–1 including quasi-periodic glitchers (Fig. 6, left) or $\langle \lambda \rangle =1.2^{+0.5}_{-0.4}$ yr^–1 excluding quasi-periodic glitchers (Fig. 6, right). Formally, the K-S probabilities are $P_{\mathrm{K}\,- \mathrm{S}\,}=0.9946$ and 0.82, respectively.

Figure 6.

Cumulative distribution of mean glitching rate λ (measured in yr^–1) for the nine pulsars that have glitched more than five times (left), and excluding the two quasi-periodic glitchers (right), showing the observational data (histograms) and the best fits obtained from eq. (12) (solid curves), where $\langle \lambda \rangle =1.3^{+0.7}_{-0.6}$ yr^–1 (left) and $\langle \lambda \rangle =1.2^{+0.5}_{-0.4}$ yr^–1 (right).

The distribution is incompletely sampled below an effective minimum rate $\lambda _{\mathrm{min}\,}$ , which is set by $\Delta t_{\mathrm{max}\,}$ . To illustrate, if we proclaim five glitches (arbitrarily) to be the minimum number needed for a reliable determination of λ, we obtain $\lambda _{\mathrm{min}\,}=5\Delta t^{-1}_{\mathrm{max}\,}\sim 0.2$ yr^–1. Careful modeling of this observational bias is deferred to a future paper. We describe a first attempt in § 5.4.

5.4. Aggregate Distribution

The nine pulsars in Figure 4 account for only 108 out of a total of 285 observed glitches. Most glitching pulsars have only glitched once or twice, but they still contribute statistical information regarding waiting times, the former via lower limits on $\Delta t$ . While these data cannot usefully constrain $P( \lambda ,\Delta t)$ in individual pulsars, they feed into the aggregate waiting-time distribution, and hence constrain $q( \lambda )$ more tightly than in § 5.3.

In Figure 7, we present the aggregate waiting-time distribution $P_{\mathrm{agg}\,}( \Delta t)$ including (left panel) and excluding (right panel) the quasi-periodic glitchers PSR J0537–6910 and PSR J0835–4510. The histogram is constructed to include all 182 waiting times in those objects that have glitched more than once. The K-S test confirms that the aggregate distribution is fitted poorly by a single, constant-rate, Poisson distribution of the form of equation (5). Furthermore, when we weight equation (5) by the exponential rate distribution (eq. [12]),12 viz.,

the fit remains poor. For example, the dotted curves in Figure 7 are computed by evaluating equation (13) with the mean values $\langle \lambda \rangle =1.3^{+0.7}_{-0.6}$ yr^–1 (left panel) and $\langle \lambda \rangle =1.2^{+0.5}_{-0.4}$ yr^–1 (right panel) extracted from Figure 6. They yield $P_{\mathrm{K}\,- \mathrm{S}\,}=4.3\times 10^{-2}$ and $2.4\times 10^{-2}$ , respectively. If, instead, we adjust $\langle \lambda \rangle$ to maximize $P_{\mathrm{K}\,- \mathrm{S}\,}$ while fitting equations (12) and (13) to the observed $P_{\mathrm{agg}\,}( \Delta t)$ , as shown by the dashed curves in Figure 7, we obtain $\langle \lambda \rangle =1.1$ yr^–1, $P_{\mathrm{K}\,- \mathrm{S}\,}=0.18$ (left panel) and $\langle \lambda \rangle =0.92$ yr^–1, $P_{\mathrm{K}\,- \mathrm{S}\,}=0.31$ (right panel), respectively.

Figure 7.

Aggregate cumulative waiting-time distribution $P_{\mathrm{agg}\,}( \Delta t)$ for all pulsars with $N_{g}\geq 2$ , including (left) and excluding (right) the two quasi-periodic glitchers. The histograms represent the observational data. The dotted curves are obtained by evaluating eq. (13) with eq. (12), using $\langle \lambda \rangle =1.3^{+0.7}_{-0.6}$ yr^–1 (left) and $\langle \lambda \rangle =1.2^{+0.5}_{-0.4}$ yr^–1 (right) as extracted from Fig. 6. The dashed curves are obtained by evaluating eq. (13) with eq. (12) and adjusting $\langle \lambda \rangle$ to maximize $P_{\mathrm{K}\,- \mathrm{S}\,}$ when fitting $P_{\mathrm{agg}\,}( \Delta t)$ . The solid curves are obtained from eq. (15), with $\langle \lambda \rangle =0.54$ yr^–1 and $\lambda _{\mathrm{min}\,}=0.25$ yr^–1 (left), and $\langle \lambda \rangle =0.43$ yr^–1 and $\lambda _{\mathrm{min}\,}=0.25$ (right).

We can exploit the extra information in Figure 7 from objects with $2\leq N_{g}\leq 5$ to determine $q( \lambda )$ more accurately. To do this, we assume a rate probability density function of the form

normalize $P( \lambda ,\Delta t)$ over the range $[ \Delta t_{\mathrm{min}\,},\Delta t_{\mathrm{max}\,}]$ , and evaluate equation (13) to obtain

In equation (15), we neglect for simplicity the observational bias introduced by uncertainties in glitch epoch discussed in § 5.1; lacking fuller information, we take $\Delta t_{\mathrm{min}\,}=0$ and $\Delta t_{\mathrm{max}\,}=28.03$ yr for all pulsars.

Excellent fits are obtained using equation (15). We find $\langle \lambda \rangle =0.54$ yr^–1, $\lambda _{\mathrm{min}\,}=0.25$ yr^–1, and $P_{\mathrm{K}\,- \mathrm{S}\,}=0.76$ for all nine pulsars with $N_{g}> 5$ , and $\langle \lambda \rangle =0.43$ yr^–1, $\lambda _{\mathrm{min}\,}=0.25$ yr^–1, and $P_{\mathrm{K}\,- \mathrm{S}\,}=0.98$ when the quasi-periodic glitchers are excluded. The fits are plotted as solid curves in Figure 7 (left and right, respectively). We verify the results by referring back to the measured $q( \lambda )$ distribution. Substituting the fitted values of $\langle \lambda \rangle$ and $\lambda _{\mathrm{min}\,}$ into equation (14), we obtain the solid curves in Figure 8, with $P_{\mathrm{K}\,- \mathrm{S}\,}=0.52$ (left panel) and 0.25 (right panel), respectively. Alternatively, the previous fits $\langle \lambda \rangle =1.1$ yr^–1, $\lambda _{\mathrm{min}\,}=0$ , and $\langle \lambda \rangle =0.54$ yr^–1, $\lambda _{\mathrm{min}\,}=0$ yield $P_{\mathrm{K}\,- \mathrm{S}\,}=0.72$ (left panel, dashed curve) and 0.52 (right panel, dashed curve).

Figure 8.

Cumulative distribution of mean glitching rate λ (measured in yr^–1) for the nine pulsars that have glitched more than five times (left), and excluding the two quasi-periodic glitchers (right), showing the observational data (histograms) and the theoretical rate distribution (eq. [14]) corresponding to the dashed and solid curves in Fig. 7.

In all cases, $P_{\mathrm{agg}\,}( \Delta t)$ points to the existence of more low-rate objects than the $N_{g}> 5$ sample in Figure 6 predicts. Specifically, up to ∼30 % of the population of glitching pulsars can have $\lambda < 0.25$ yr^–1 while still reproducing $P_{\mathrm{agg}\,}( \Delta t)$ . We emphasize again that Figure 7 constrains $q( \lambda )$ more tightly than Figure 6, because it contains information about $\Delta t$ from 1.3 times as many glitches, as well as useful information from pulsars that have glitched twice.

Undetectable microglitches probably occur between detected glitches without our knowledge, given that $p( \Delta \nu / \nu )$ is scale invariant. This effect subtracts from the lower end of the $\Delta t$ distribution and adds to the upper end. We do not correct for it here, because it is hard to quantify without better statistics. On two occasions, a pair of glitches occurred on the same date: once in the same pulsar, and once in different pulsars (M. Kramer & A. Lyne 2005, private communication). We take $\Delta t=0$ for these pairs. Phase-connected timing mitigates duty cycle biases, but it does not eliminate them.

5.5. Fluctuation Spectrum

The temporal fluctuations in a stochastic signal x( t)

carry independent statistical information about the underlying physical process. The power spectrum S( f)

, where f denotes the Fourier frequency, is related to the temporal autocorrelation function $G( \tau ) =\langle x( t) x( t+\tau ) \rangle -\langle x( t) \rangle ^{2}$ (where the average $\langle \ldots \rangle$ is performed over t for a stationary process) through the cosine transform

In general, for an avalanche process, the power spectrum depends jointly on the size, waiting-time, and lifetime distributions of the avalanches (Jensen 1998). For glitches, however, the lifetimes are too short to measure with current technology (see § 6). If, furthermore, we restrict attention to the unit-impulse signal $x( t) =\sum_{i}\delta ( t-t_{i})$ , where $t_{i}$ denotes the epoch of the ith glitch, then the sizes drop out of the problem too. The power spectrum then carries exactly the same information as the waiting-time distributions $P( \lambda ,\Delta t)$ and $P_{\mathrm{agg}\,}( \Delta t)$ , with

for any individual pulsar, and

for the pulsar population in aggregate.

At high frequencies $f\gg \langle \lambda \rangle$ , equations (17) and (18) [with $q( \lambda )$ given by eq. (12)] scale as $f^{-2}$ , with $O( f^{-4})$ and $O( f^{-4}\mathrm{sin}\,2f)$ corrections. These scalings are modified if the delta function in x( t)

is replaced by a nonsingular window function that embodies the shape of the signal from an individual avalanche. It will be instructive to revisit this question when it becomes possible to resolve the lifetimes of individual avalanches, e.g., in single- or giant-pulse timing experiments.

9 It is tempting to assume that λ is constant over decades, because the thermodynamic variables that control storage capacity (e.g., temperature) are nearly constant on such a timescale. Yet the Sun provides a cautionary counterexample: the dynamics of subphotospheric turbulence, and hence the rate of solar flaring, vary with the 11 yr solar cycle (Wheatland 2000).

10 For

, $\langle \Delta \nu \rangle$ is dominated by the upper cutoff, but the normalization of $p( \Delta \nu / \nu )$ is dominated by the lower cutoff.

11 The error in $\Delta \nu$ is $d\langle \Delta \nu \rangle / da$ multiplied by the error in a.

12 We compute $P_{\mathrm{agg}\,}( \Delta t)$ theoretically as a weighted sum of independent Poisson processes. In the same way, the waiting-time distribution for decays observed from a mixture of radioactive isotopes is a sum of constant-rate Poisson distributions, one per isotope, weighted by isotopic abundance.

6. DISCUSSION

In this paper, we analyze the size and waiting-time distributions of pulsar glitches, taking advantage of the enlarged data set produced by the Parkes Multibeam Survey. We conclude that the data are consistent with the hypothesis that pulsar glitches arise from an avalanche process. In each of seven pulsars with $N_{g}> 5$ , the size distribution is consistent with being scale invariant across the observed range of $\Delta \nu$ (up to four decades), and the waiting-time distribution is consistent with being Poissonian. These features are natural if the system is driven globally at a constant rate (as the pulsar spins down), and each glitch corresponds to a locally collective, threshold-activated relaxation of one of the many spatially independent, metastable stress reservoirs in the system (e.g., via a vortex-unpinning or crust-cracking avalanche). In two pulsars, PSR J0537–6910 and PSR J0835–4510, the dynamics may include a second, quasi-periodic component, comprising ∼20% of the events. The size and waiting-time distributions of the quasi-periodic component are narrowly peaked, as expected for rare, system-spanning avalanches, which relax a large fraction of the total stress accumulated in the system. This two-component behavior is observed widely in self-organized critical systems, including experiments on magnetic flux vortices in type II superconductors, which are closely analogous to neutron star superfluids (Field et al. 1995).

The power-law exponent of the size probability density function differs from pulsar to pulsar, spanning the range $-0.13\leq a\leq 2.4$ . Calculating a theoretically from first principles is a deep problem that has not yet been solved for other self-organized critical systems, let alone glitching pulsars, although some progress has been made on two-dimensional sandpiles using renormalization group techniques (Pietronero et al. 1994; Jensen 1998). In the mean field approximation, which is exact in four or more dimensions, theoretical calculations on sandpiles (and other systems in their universality class) yield a=1.5

, whereas three-dimensional cellular automata output a=1.3

(Jensen 1998).

The size distribution transmits two important lessons concerning the microphysics of glitches. First, the fact that a differs from pulsar to pulsar implies that the strength and level of conservation of the local (e.g., pinning and intervortex) forces also differs (Olami et al. 1992; Field et al. 1995). By contrast, in equilibrium critical systems such as ferromagnets, a depends only on the dimensionality of the system and its order parameter, and is therefore universal. Second, except for the two pulsars that show quasi-periodicity, a appears to vary smoothly with spin-down age, with $a\approx 1.2$ for the youngest pulsars (e.g., the Crab). Figure 9 depicts the trend between a and $\tau _{c}$ . It is suggestive; after all, local pinning forces do depend on temperature, and hence on $\tau _{c}$ . Interestingly, however, there is no clear trend between a and ν, even though the mean vortex spacing (and hence intervortex force) is proportional to $\nu ^{1/ 2}$ . It will pay to study these trends more thoroughly as more glitch data are collected.

Figure 9.

Spin-down age $\tau _{c}$ (in kyr) vs. power-law exponent a for the glitch size distribution. The error bars indicate the 1 σ range of allowable fits according to the K-S test. Systematic differences between $\tau _{c}$ and true age are not quantified here. Triangles show quasi-periodic glitchers with $N_{g}> 5$ , to which we fit a two-component (filled triangles) and one-component (open triangles) $P( \Delta \nu / \nu )$ , as in Tables 4 and 3, respectively. Squares show aperiodic glitchers with $N_{g}> 5$ , to which we fit a power law $P( \Delta \nu / \nu )$ , as in Table 3.

An avalanche process predicts a specific relation between the distributions of glitch sizes $\Delta \nu$ and lifetimes T (as opposed to waiting times $\Delta t$ ). Specifically, in a self-organized critical state, the lifetime probability density function is also a power law, $p( T) \propto T^{-b}$ , with

The constants $\gamma _{2}$ and $\gamma _{3}$ are defined such that the cardinality of an avalanche scales with its linear extent (L) as $L^{\gamma _{2}}$ and its lifetime (i.e., duration) scales as $L^{\gamma _{3}}$ (Jensen 1998). Both $\gamma _{2}$ and $\gamma _{3}$ depend on the effective dimensionality of the local forces, and can be calculated numerically using a cellular automaton. In two dimensions, avalanches are compact, not fractal, and one has $\gamma _{2}=2$ ; in three dimensions, one has $2< \gamma _{2}< 3$ . At present, radio timing experiments cannot measure T; most glitches are detected as unresolved, discontinuous, spin-up events with T< 120

s (McCulloch et al. 1990).13 In the future, however, single- and/or giant-pulse timing experiments with more sensitive instruments (e.g., the Square Kilometer Array) will test this prediction. If confirmed, it will independently corroborate the avalanche hypothesis.

The mean glitching rates of the nine pulsars studied here are fairly narrowly distributed, spanning the range $0.35\leq \lambda \leq 2.6$ yr^–1. The probability density function for λ is adequately fitted by an exponential, as for solar flare avalanches (Wheatland 2000), with $\langle \lambda \rangle =1.3^{+0.7}_{-0.6}$ yr^–1, or by an exponential with a lower cutoff, at $\lambda _{\mathrm{min}\,}\approx 0.25$ yr^–1. A theoretical derivation of $\langle \lambda \rangle$ from first principles is currently lacking, although estimates of how long it takes to crack the crust locally predict reasonable rates, if the critical strain angle approaches that of imperfect terrestrial metals (Alpar et al. 1996; Middleditch et al. 2006).

Figure 10 plots λ versus $\tau _{c}$ for the nine pulsars examined individually in this paper. There is no significant trend. The data are consistent with the notion that old pulsars glitch less frequently than young pulsars (Shemar & Lyne 1996), but they are equally consistent with the notion that the glitching rate is independent of age.

Figure 10.

Mean glitching rate λ (in yr^–1) vs. spin-down age $\tau _{c}$ (in kyr). The error bars indicate the 1 σ range of allowable fits according to the K-S test. Systematic differences between $\tau _{c}$ and true age are not quantified here. Triangles show quasi-periodic glitchers with $N_{g}> 5$ , to which we fit a two-component (solid triangles) and one-component (open triangles) $P( \lambda ,\Delta t)$ , as in Tables 6 and 5, respectively. Squares show aperiodic glitchers with $N_{g}> 5$ , to which we fit a Poissonian $P( \lambda ,\Delta t)$ as in Table 5.

Many authors have searched for a correlation between waiting time and the size of the next glitch. Such a correlation appears to be absent from the data, e.g., Figure 17 in Wang et al. (2000) and Figure 10 in Wong et al. (2001). At first blush, this is surprising: the vortex-unpinning and crust fracture paradigms, which are driven by the accumulation of differential rotation and mechanical stress, respectively, seem to be natural candidates for a “reservoir effect.” Avalanche dynamics resolves this apparent paradox. The reservoir effect does operate locally, but the star contains many reservoirs, insulated from each other by relaxed zones, whose storage capacities evolve stochastically in response to the slow driver and avalanche history. During a glitch, a single reservoir (often small, but sometimes large) relaxes at random via an avalanche, releasing its stored $\Delta \nu$ (and destabilizing neighboring reservoirs in preparation for the next glitch). Some of the $\Delta \nu$ has accumulated since the previous glitch, but the remainder is “borrowed” from earlier epochs, when some other reservoir relaxed instead. All self-organized critical systems share these dynamics; the waiting time is uncorrelated with the size of the next avalanche (Jensen 1998). The only exceptions are large, system-spanning avalanches, which always have roughly the same sizes and waiting times, and which account for ∼20% of the glitches in PSR J0537–6910 and PSR J0835–4510.

A corollary of the previous paragraph is that the total $\Delta \nu$ released in glitches up to some epoch is less than the total crust-superfluid differential rotation accumulated since that epoch, viz.,

where $\epsilon \dot{\nu }$ is the relative angular acceleration of the crust and superfluid due to electromagnetic spin down. The “staircase” described by equation (20) has been noted previously (Shemar & Lyne 1996; Lyne et al. 2000), both in quasi-periodic glitchers such as PSR J0537–6910 (e.g., Figure 8 in Middleditch et al. 2006), where the reservoir effect is obvious, and in Poisson glitchers such as PSR J0534+2200, (e.g., Figure 12 in Wong et al. 2001), where the trend is more subtle because it reverts to the mean over long times, not after every glitch. On dividing equation (20) by $N_{g}$ , and averaging over long times, the inequality becomes an equality (provided there is no secular accumulation of differential rotation in the system), and we recover equation (6).

It is fundamentally impossible to measure ϵ for individual pulsars with current data, because $\langle \Delta \nu \rangle$ is dominated by large (and therefore rare) glitches for a< 2

. It is therefore wrong to assume stationarity over a typical, 40 yr observation interval. Consequently, we are prompted to reassess the familiar correlation between activity and spin-down age (Shemar & Lyne 1996). Our definition of $\epsilon \dot{\nu }$ is identical to $\dot{\nu }_{\mathrm{glitch}\,}$ in Lyne et al. (2000) (but for individual objects, not in aggregate) and $A_{g}$ in Wong et al. (2001). It is closely related to the original activity parameter defined by McKenna & Lyne (1990), which equals $N_{g}\nu ^{-1}\epsilon \left|\dot{\nu }\right|$ . For PSR J0358+5413, we measure $\epsilon \leq 7\times 10^{-5}$ , lower than the aggregate value $0.017\pm 0.002$ measured by Lyne et al. (2000) for objects with $\tau _{c}> 10$ kyr (binned by semidecades in $\dot{\nu }$ ).14 Interpreted in terms of the vortex-unpinning model, this result suggests that 0.007%-2% of the angular momentum outflow during spin down may be stored in metastable reservoirs on average over time. On the other hand, five other objects have $0.04\leq \epsilon _{\mathrm{max}\,}\leq 0.8$ , under the questionable assumption that the maximum physical size is $\Delta \nu _{\mathrm{upper}\,}=2\times 10^{-4}$ in all pulsars. Our data are therefore inadequate to update usefully the value $A_{g}/ \left|\dot{\nu }\right|=1\times 10^{-5}$ measured by Wong et al. (2001) for PSR J0534+2200.

In the context of vortex unpinning, it has been argued that the aggregated ϵ measured by Lyne et al. (2000) partly corroborates the hypothesis that younger pulsars are still in the process of forming their capacitive elements, e.g., by creating pinning centers through crust fracture, while older pulsars have mostly completed the task (Alpar et al. 1996; Wong et al. 2001). However, the full picture is more complicated. Vela’s quasi-periodic avalanches point to a richly connected network of reservoirs (Alpar et al. 1996), yet its aggregated value $\dot{\nu }_{\mathrm{glitch}\,}$ is relatively low. On the other hand, the other quasi-periodic glitcher, PSR J0537–6910, is relatively young ( $\tau _{c}=4.9$ kyr); how did it form a richly connected reservoir network so quickly? And, if its network is so richly connected, why is its aggregated $\dot{\nu }_{\mathrm{glitch}\,}$ value so low? Likewise, PSR J0358+5413 is the oldest object in the sample ( $\tau _{c}=560$ kyr), yet its ϵ value arguably points to a dearth of capacitive elements, characteristic of a young object. There are no obvious grounds (e.g., quasi-periodicity) on which to treat PSR J0358+5413 as exceptional.

Do all pulsars glitch eventually? It has been speculated in the past that there is something special physically about the minority of pulsars that do glitch. While it is impossible to reject this hypothesis unequivocally with the data at hand, the results presented here suggest that all pulsars are capable of glitching. However, most do so infrequently (low λ) and hence have not been detected during the last four decades of timing experiments. We find that up to ∼30% of the pulsar population can glitch at rates lower than $\lambda _{\mathrm{min}\,}=0.25$ yr^–1 and still conform with the measured aggregate waiting-time distribution.

Once verified, the claimed Poissonian nature of the glitch mechanism can be invoked to exclude broad classes of glitch theories, e.g., those that rely on “defects” or “turbulence” at special locations (such as the pole), or that involve a pair of dependent events (A. Martin 2007, private communication). It is important to interpret aftershocks carefully in this light (Wong et al. 2001). In self-organized critical systems, the excess number of avalanches following a large avalanche (over and above the Poissonian baseline following a small avalanche) scales inversely with the time elapsed, a property known as Omori’s law for earthquakes (Jensen 1998).

In this paper, we do not analyze postglitch relaxation times and glitch-activated changes in $\dot{\nu }$ in the context of avalanche processes, e.g., the correlation between $\Delta \dot{\nu }$ and the transient component of $\Delta \nu$ (Wong et al. 2001). We also assume implicitly that the quantized superfluid vortices in the vortex-unpinning model are organized in a rectilinear array, even though recent work suggests that meridional circulation destabilizes the array and converts it into a turbulent tangle (Peralta et al. 2005, 2006). Further study of these matters is deferred to future work.

13 In the Crab, some spin-up events seem to be resolved, e.g., at epochs MJD 50,260 ( $T\approx 0.5$ days) and MJD 50,489 (secondary spin up, $T\approx 2$ days) (Wong et al. 2001). If these are rare but otherwise standard glitches originating from the long- T tail of the lifetime distribution, it is puzzling that other shorter but still resolved (and presumably more common) spin-up events are not observed, e.g., with $T\sim 0.1$ or 0.01 days. Alternatively, the events at MJD 50,620 and MJD 50,489 may have been triggered by a different physical mechanism.

14 The aggregate value $\dot{\nu }_{\mathrm{glitch}\,}$ (Lyne et al. 2000), binned over semidecades in $\dot{\nu }$ , effectively averages together different pulsars. While this approach reduces the formal error bar on $\dot{\nu }_{\mathrm{glitch}\,}$ , its physical interpretation is less straightforward, given the likelihood that ϵ is different in different pulsars.