A cellular automaton model of pulsar glitches
Abstract
A cellular automaton model of pulsar glitches is described, based on the superfluid vortex unpinning paradigm. Recent analyses of pulsar glitch data suggest that glitches result from scaleinvariant avalanches (Melatos
et al., 2007), which are consistent with a selforganized critical system (SOCS). A cellular automaton provides a computationally efficient means of modelling the collective behaviour of up to vortices in the pulsar interior, whilst ensuring that the dominant aspects of the microphysics are not lost. The automaton generates avalanche distributions that are qualitatively consistent with a SOCS and with glitch data. The probability density functions of glitch sizes and durations are power laws, and the probability density function of waiting times between successive glitches is Poissonian, consistent with statistically independent events. The output of the model depends on the physical and computational paramters used. The fitted power law exponents for the glitch sizes () and durations () decreases as the strength of the vortex pinning increases. Similarly the exponents increase as the fraction of vortices that are pinned decreases. For the physical and computational parameters considered, one finds and , and mean glitching rates in the range in units of inverse time.
keywords:
pulsars: general — stars: interior — stars: neutron1 Introduction
It has recently been shown that glitches in individual pulsars obey avalanche statistics (Melatos et al., 2007). Such statistical behaviour, in conjunction with the wide dynamical range of glitches, suggests that the underlying physics is that of a selforganized critical system (SOCS) (Pines & Alpar, 1985; Bak, 1996; Jensen, 1998; Melatos et al., 2007). SOCS are characterised by transitions between metastable states via scale invariant avalanches (Jensen, 1998). The fractional increase in spin frequency spans seven decades across the entire glitch population (), and up to four decades in any individual pulsar. A natural way to model a scale invariant, avalanching system is to construct a cellular automaton that is driven slowly to a threshold of instability (Field et al., 1995). Such models have been explored in detail in the context of earthquakes (Sornette et al., 1991), granular assemblies (Bak et al., 1988; Wiesenfeld et al., 1989; Bak, 1996; Frette et al., 1996; Pruessner & Jensen, 2003), solar flares (Lu & Hamilton, 1991), and superconducting vortices (Jensen, 1990; Field et al., 1995; Bassler & Paczuski, 1998; Linder et al., 2004).
Theories of pulsar glitches centre around the mass unpinning model first proposed by Anderson & Itoh (1975) and extended by many others (Alpar et al., 1984; Pines & Alpar, 1985; Alpar et al., 1989; Link & Cutler, 2002). The neutron superfluid in the stellar interior is threaded by many () vortices, approximately one percent of which are pinned to the stellar crust at grain boundaries and/or nuclear lattice sites (Alpar et al., 1989). As the pulsar crust spins down electromagnetically, a lag builds up between the velocity of the pinned vortex lines (corotating with the crust) and the superfluid. When the transverse Magnus force (directly proportional to the lag) surpasses a threshold value (equal to the strength of the pinning force), a catastrophic unpinning of vortices occurs, transferring angular momentum to the crust. In order for this mechanism to generate glitches on the scale observed, it requires up to vortices to unpin simultaneously, exhibiting a high level of collective, nonlocal behaviour.
The mass unpinning model spawned much activity in quantifying the microphysics that would lead to such macroscopic behaviour (Pines & Alpar, 1985; Jones, 1997, 1998a, 1998b; De Blasio & Elgarøy, 1999; Elgarøy & De Blasio, 2001; Donati & Pizzochero, 2003; Avogadro et al., 2007). In particular, the strength of the pinning force has recently been discussed in detail using a selfconsistent, semianalytical model for the vortexnucleus interaction (Avogadro et al., 2007), yielding the unexpected result that the strongest pinning occurs in regions of the crust that are of high and low density, rather than intermediate density (Donati & Pizzochero, 2006). An interest has also been taken in the evolution and morphology of the neutron star crust, with particular attention paid to impurities and defects, which determine the density of pinning centres available to the vortices, and their relative strengths (De Blasio & Lazzari, 1998). Link & Cutler (2002) looked at the role that precession, if present, plays in the unpinning process, and discussed the dependence of the pinnning force on various microphysical parameters of the superfluid vortices and of the nuclear crustal lattice.
Several other theories of pulsar glitches, similarly focused on the microphysics of the superfluid in the stellar interior and crust, are currently under consideration. These include thermally driven glitches arising from sudden changes in the frictional coupling within the neutron star (Link & Epstein, 1996), superfluid turbulence (Peralta et al., 2005, 2006b, 2006a; Melatos & Peralta, 2007), superfluid twostream instabilities (Andersson et al., 2004; Mastrano & Melatos, 2005), crust cracking (Ruderman, 1969; Bak, 1996; Alpar et al., 1996; Middleditch et al., 2006) and mass accretion (Morley & Schmidt, 1996; Morley, 1996).
Regardless of the model one chooses to examine, there is a common need to synthesise the microphysics with the global, collective dynamics in order to accurately describe the observational data. In order to produce a realistic model of the global dynamics, some subtle aspects of the microphysics need to be sacrificed. Cellular automaton models of complex systems aim to employ only the dominant aspects of the microphysics to generate replicable and efficient automaton rules. Such an approach was attempted by Morley & Schmidt (1996) using a cellular automaton platelet model of mass accretion onto the stellar crust. Although the microphysical details differ from those treated in this paper, the underlying philosophy of using a simple model to describe large scale behaviour is identical to ours. We are similarly encouraged by the marked parallels between superfluid vortices and those that populate hard, type II superconductors. The agreement between cellular automaton models of flux creep and vortex avalanches in superconductors and experimental data is excellent (Linder et al., 2004; Jensen, 1990; Field et al., 1995; Bassler & Paczuski, 1998; Nicodemi & Jensen, 2001a, b), in particular with respect to the measurement of a power law over several decades in the avalanche (cf. glitch) size.
A general, first principles theory of selforganized criticality does not yet exist for any SOCS, let alone for pulsar glitches. In this paper we look for empirical agreement between the output of our cellular automaton and the gross qualitative features in pulsar glitch data. A detailed quantitative comparison between a first principles theory and data is not possible at this stage.
In this paper we present a cellular automaton model of pulsar glitches based on the mass unpinning model. After describing the key features of a SOCS we revisit the observational pulsar glitch data in §2. §3 reviews the vortex unpinning paradigm and recent advances in the understanding of neutron star structure which are relevant to the paradigm. §4 describes in detail the cellular automaton model and justifies physically the automaton rules. In §5 the statistical behaviour of the model is explored with respect to driver parameters, internal physical parameters, and variations in the automaton rules. Further discussion of the model’s validity and a summary of our findings are contained in §6.
2 Selforganized critical systems
2.1 General properties
A system is described as a SOCS based on two underlying behavioural patterns. Selforganized implies the ability to develop structures and patterns in the absence of manipulation by an external agent. Critical implies that, like in phase transitions near the critical temperature, a small local perturbation can propagate throughout the entire system (Jensen, 1998).
For a system to be in a selforganized critical (SOC) state it must first satisfy the following three conditions (Jensen, 1998; Melatos et al., 2007):

It is composed of many discrete, mutually interacting elements, whose motions are dominated by local (e.g. nearestneighbour) rather than global (e.g. meanfield) forces.

Each element moves when the local force exceeds a threshold. In this way stress accumulates sustainedly at certain random locations (metastable stress reservoirs) while relaxing quickly elsewhere.

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 quasistatically through a historydependent sequence of metastable states.
If the above conditions are met, the following properties are generically observed (Jensen, 1998; Melatos et al., 2007): 
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 (knockon effect).

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. Numerical values of the exponents depend on the spatial dimensionality of the system, the symmetry and strength of the local forces (Field et al., 1995), and the level of conservation (Vespignani & Zapperi, 1998; Pruessner & Jensen, 2002).

Over the long term, the system tends to a critical state, which is stationary on average but fluctuates instantaneously.
A driving force is not essential to the operation of a SOCS. Pruessner & Jensen (2002) assert that an external driving force compensates for a lack of conservation in the microscopic dynamics of a SOCS and hence is superfluous if the system is conservative on microscopic scales. We claim that our model is locally nonconservative and hence a driving force is essential to maintaining criticality (see §4.7).
Avalanches result from the relaxation of isolated metastable stress reservoirs. The storage capacity of these reservoirs is scale invariant. Each relaxation event is statistically independent from other relaxation events in size and when it occurs. Two or more simultaneous relaxation events at different locations cannot be distinguished by their global effect (e.g. spin up of a neutron star), which is the cumulative result of multiple relaxations. Likewise, relaxation of two or more smaller reservoirs causes an avalanche which is indistinguishable from that produced by the relaxation of a single large reservoir. Statistical independence implies that the waiting times between avalanches follow Poisson statistics, and this is confirmed by models based on cellular automata.
Alpar et al. (1995); Alpar et al. (1996) discussed the formation of depleted and capacitive regions of vortices in a neutron star interior, generated by steep gradients in the vortex pinning potential. Wherever the pinning potential is greater, the vortex density is also greater, resulting in a stronger local intervortex interaction (Magnus force). When the Magnus force exceeds the pinning threshold, vortices are redistributed locally into neighbouring regions, relaxing the capacitive region (stress reservoir). This qualitative picture bears the hallmark of a SOCS.
2.2 Pulsar glitch statistics
PSR J  

18250935  0.91  0.48  1.8  0.36  0.3  1.0  0.2 
0534+2200  0.91  0.57  1.3  1.2  1.1  1.4  1.1 
18012304  0.55  0.35  0.88  0.57  0.092  1.1  4 
13416220  1.8  1.2  5.6  1.4  1.2  2.1  10 
0631+1036  0.95  0.55  1.9  1.8  1.2  2.7  1.33 
17403015  1.5  1.2  2.5  1.1  0.98  1.3  0.7 
Melatos et al. (2007) compiled all of the available glitch data for pulsars, analysing in detail the nine pulsars that have been observed to glitch at least five times. We consider here the subset that have glitched more than six times (). Of the eight pulsars satisfying this criterion, six are consistent with a power law probability density function in the fractional glitch size , where is the spin frequency of the pulsar, implying a cumulative probability of the form
(1) 
where is the maximum observed glitch size, and is the mimimum observed glitch size. Melatos et al. (2007) found that takes a different value for each pulsar; when the data from all nine pulsars with are considered, Eq. (1) with a single value of is not a good description of the aggregate distribution. Two pulsars also exhibit quasiperiodic behaviour, which is also a natural manifestation of avalanche dynamics; when the global forces (driven externally) overwhelm the nearestneighbour interactions, quasiperiodicity in the waiting time distribution is expected (Jensen, 1998).
For each of the six pulsars analysed here, we find all values of the exponent , for which the data are not inconsistent with the analytic distribution at the confidence level according to the KolmogorovSmirnov (KS) test. The range of is presented in Table 1. Table 1 also gives the exponent that minimises the KS statistic, along with the smallest glitch , for each pulsar (Melatos et al., 2007). However, it should be noted that all values of within the range (bounded by and ) are equally likely; the exponent that minimises the KS statistic is not a best fit. To appreciate why, suppose that the true underlying distribution of is known for some pulsar. Each time the pulsar glitches it samples this underlying distribution. The glitch data for this pulsar represent one possible realisation of the distribution containing glitches. There is no reason to attach greater significance to the value that minimises the KS statistic with respect to the particular realisation than to any other in the range.
Melatos et al. (2007) also investigated the distribution of waiting times , between successive glitches. Based on the statistical independence of each avalanche and the assumption that the system is driven at a constant rate, it is expected that the probability density function for is Poissonian (Jensen, 1998), implying a cumulative probability distribution of the form
(2) 
The sum in (2) is taken over the minimum observable waiting time , which is set by the gap between data spans in which a glitch is localised. is a function of the observing schedule and hence of epoch; effectively, therefore, it takes a different value for each glitch (Melatos et al., 2007). is the total interval over which a pulsar is observed.
Figure 1 shows cumulative histograms of waiting times , for the same six pulsars listed in Table 1, plotted agaist , where is the Poisson rate for each pulsar in units of . Once again, is chosen to minimise the KS statistic for definiteness, but any value in the range would do just as well at the confidence level. The observational data is overlaid with the ideal cumulative Poissonian for comparison, taking and ; that is, we do not correct for data gaps and a finite total observation for simplicity, unlike in Melatos et al. (2007). The agreement with the universal Poissonian is good for these six objects, with . Observational confirmation that pulsar glitch waiting times obey Poisson statistics upholds the prediction that the glitches arise from the relaxation of metastable stress reservoirs isolated by relaxed zones.
The observed distribution of glitch durations cannot be discussed at this time, as very few glitches have been resolved temporally (McCulloch et al., 1990).
3 Vortex unpinning in a pulsar
The mass unpinning paradigm of pulsar glitches assumes that the inner crust is a dense nuclear lattice permeated by a superfluid that is threaded by many quantised vortices. The vortices interact repulsively via the Magnus force^{1}^{1}1The Magnus force is in fact a fictitious force, similar in nature to the Coriolis or centrifugal forces. A free vortex, in the absence of external forces, moves with the local superfluid flow. In order for a vortex to move with respect to the local flow, a force must be applied such that the total force per unit length of vortex line, , is sufficient to sustain the motion of the line with respect to the local flow. (per unit length of the vortex line) given by
(3) 
where is the superfluid density, is a vector pointing along the vortex, whose modulus is the circulation around a single vortex ( is the neutron mass), is the vortex line velocity, and is the local bulk superfluid velocity. If allowed to move freely, vortices organize into a rectilinear array spaced evenly according to Feynman’s rule (vortex area density ). In keeping with condition 1 in §2.1, each vortex is considered a discrete element of a SOCS. For the purpose of practical modelling, given that we are dealing with vortices in a neutron star, we group vortices into bundles (in internal equilibrium) and regard the bundles as discrete elements. The velocity field generated by each vortex is inversely proportional to the distance from the vortex line, so local interactions dominate global forces.
Condition 2 in §2.1 exactly describes the pinning of vortices to the pulsar crust, whether the pinning centres are nuclear lattice sites, defects [monovacancies or impurities (De Blasio & Lazzari, 1998)], or grain boundaries. Condition 2 does not rule out pinning of neutron and proton vortices in the neutron star’s core. Vortices are unpinned by the bulk superfluid flow when the Magnus force exceeds the pinning force. Jones (1997) suggested that vortex pinning is too weak to occur in the pulsar crust, and that strong magnetic interactions between neutron and proton vortices are a feasible alternative. In the event that strong pinning is provided by monovacancies alone, Jones (1998a, b) calculated the monovacancy concentration needed to provide a pinning force greater than the Magnus force, and claimed that, for , the concentration is unphysically large.
The external driver fulfilling condition 3 in §2.1 is the accumulating Magnus force generated by the build up of rotational lag between the vortex lines and the superfluid flow due to the electromagnetic spin down of the star. The effect of this slow build up is twofold. Firstly, the vortices tend to slowly move radially outward, as the system strives to maintain the Feynman vortex density (Alpar et al., 1984; Jones, 1998a). Secondly, abrupt adjustments occur when locally the pinning force threshold is exceeded. Once a vortex is unpinned, it readjusts rapidly, either by returning to the ‘soup’ of unpinned vortices (remembering that only of vortices are pinned at any one time), or by repinning at a nearby lattice or defect site. Most mesoscopic vortex models assume that the inertia of vortex lines is negligible: once unpinned, vortices move immediately with the bulk superfluid flow, regardless of the direction of the unpinning force (Donnelly, 1991). Field et al. (1995) pointed out that, at least in the context of flux creep in superconductors, such an assumption leads to ‘stickslip’ dynamics that are more likely to resemble those of an ideal sandpile (the quintessential SOCS), because inertial effects do not dominate the effect of the threshold. The timescale of the readjustment is much shorter than the driving timescale of the stellar spin down.
Throughout this paper we assume that superfluid vortices are straight and parallel to the axis of rotation. Jones (1998b) argued that the high, but finite, tension in vortex lines justifies this assumption, as the displacements caused by the pinning centres are much smaller than the nuclear separation. Contrarily, it is the finite tension in vortices that allows them to be drawn into the potential well of the pinning centres and pin at all (Jones, 1998b). Link & Cutler (2002) suggested pinning centres in the nuclear lattice effectively span approximately thirty lattice sites, due to finite tension, where the internuclear spacing is . Relative to the coarse graining scale of a practical cellular automaton (linear cell dimension up to ), this deviation from straight vortices is clearly insignificant. However, global hydrodynamic simulations suggest that the array breaks up into a vortex tangle above a critical rotational lag, driven by the DonnellyGlaberson instability, with a net polarisation parallel to the rotation axis [see for example Peralta et al. (2005, 2006a, 2006b)].
4 Cellular automaton
The cellular automaton model presented in this paper takes the basic physical features of the inner crustsuperfluid system and interprets them as simple rules. Much of the microphysics that has emerged from detailed studies of pinning strengths, crustal structure, and turbulent flow (see references in §1) has not been taken into account directly when constructing the automaton rules, but it is relevant for understanding how the system can be ‘scaled up’ (renormalised) to the necessary coarse grain. In terms of reproducing the largescale, collective dynamics of the SOCS, these simple cellular automaton rules have many advantages. In particular, the simple rules allow for a large range of system sizes and parameter space to be modelled for long enough to obtain stationary statistics, without prohibitive computational cost. Table 2 summarizes the nomenclature used in this and subsequent sections.
4.1 Grid
The simulation grid contains cells, covering a circular ‘star’ of radius , representing a projection onto the equatorial plane. The grid points are arranged in a rectilinear (square) array in order to easily identify nearest neighbours and to ensure that all cells are equal in size. In reality, the equilibrium configuration of superfluid vortices is a triangular Abrikosov lattice. For the purposes of this model, we assume that the two configurations are indistinguishable when vortices are bundled in large groups. Vortex positions are restricted to the coordinates of the centre of each grid cell. As mentioned in §3 there is an inherent assumption that the vortices are straight and parallel to the axis of rotation, and that the vortex dynamics are independent of the vortex length. Indeed, the intervortex force is calculated as a force per unit length.
The number of vortices that occupy the grid is calculated using the Feynman relation
(4) 
where the integral is taken around a closed path (e.g. a circle of radius equal to the stellar radius), and is the total number of vortices enclosed by the path. In the limit of an infinite vortex distribution, the equilibrium distribution is a rectilinear array. Although we have a finite array of vortices, we assume that they form a rectilinear array, such that the total number of vortices is
(5) 
Once we know how many vortices should occupy the system as a whole, we bundle them such that, on average, there is one bundle per cell. Partly, this is done to reduce the number of discrete elements and keep the computation practical. However, it is also done because uncertainties in pulsar timing limit the smallest glitch we can observe (corresponding to one bundle). For any given pulsar we model the observed glitch range: the minimum glitch size (what constitutes a glitch is discussed in §4.5) occurs when one cell unpins, and the maximum occurs when all cells unpin. This restriction ensures that our model does not produce glitches outside the observationally imposed range. This is in contrast with the model of Morley & Schmidt (1996) in which the minimum glitch size is set by timing noise and the maximum by the total available energy. Alternatively, we can have an average of bundles per cell. In this scenario, the minimum glitch size is less than the minimum observed glitch size, and hence we can disregard all glitches resulting from the unpinning of fewer than vortex bundles. We choose our grid size such that
(6) 
Pulsar timing data suggests, via (6), a range of grid sizes varying from to a maximum of ^{2}^{2}2 is the total number of cells in the square grid. However, only cells that lie within the stellar radius are actually available.. The number of vortices per bundle is then
(7) 
In an astrophysical context, however, we cannot assume that glitches smaller than the observable minimum do not take place, and with equal or greater frequency. If the system is in a selforganized critical state, we expect glitches to occur at all scales, within the bounds of the system. To resolve smaller in the model, we must fractionate into smaller bundles. In order to make comparisons with the observational data, we preclude a priori these smaller glitches from occurring in our model. If the system is truly scale invariant, this effective window function (coarsegrain scale) should not interfere with its overall statistical behaviour. In the flux creep model for superconductors developed by Jensen (1990), each cell holds at most one particle. In our model each cell holds as many bundles as the avalanche history dictates. Bassler & Paczuski (1998) take a similar approach, allowing many point pins in an extended cell.
4.2 Initial conditions
Initially the vortex bundles are laid out at random^{3}^{3}3Vortex bundles are placed at a randomly selected cell, one at a time, until all bundles have been distributed. such that the mean occupancy of a cell is bundles. Whether or not it is realistic to have a vortex distribution that is inhomogoneous on the scale of the grid cells (m) depends upon the distribution of the pinning centres. This point is discussed in more detail in §6. From the glitch data analysed by Melatos et al. (2007), we hypothesize that our system is in a SOC state at any given epoch, with an inhomogeneous configuration of pinned vortices, and take this to be the initial state of the automaton. We do not model directly how this inhomogeneous state arises from an initially homogeneous state because the latter configuration is inaccessible to observations; it is disrupted by avalanches almost immediately after the neutron star is born. Given that the random initial conditions of the model do not necessarily define a critical state of the system, we allow the simulation to complete many () time steps before considering the output, to ensure that criticality has been established.
4.3 Pinning threshold
The strength with which vortices are pinned to the stellar crust, , plays an essential role in the mass unpinning model of pulsar glitches. Recent calculations imply that pinning to nuclear lattice sites is strongest in regions of low () and high density (), with a maximum energy of MeV (Avogadro et al., 2007). We take the pinning force to be the same everywhere on the grid, with value , where m is the superfluid coherence length (Elgarøy & De Blasio, 2001); i.e. . For the purposes of the model we convert this to a threshold on the velocity lag, .
Vortex pinning at lattice defects is also possible. Link & Cutler (2002) argued that if the pinned (coupled) fluid makes up of the total moment of inertia, then pinning must occur throughout the crust, suggesting that the crust should be treated as an amorphous solid with random pinning sites. The randomness ensures that the contributions to the force of attraction towards the nuclear sites on either side of the vortex line do not cancel (Jones, 1998a). Another possibility is to adopt the random pinning potential used to model flux creep in superconductors (Bassler & Paczuski, 1998).
Vortices unpin when the relative velocity of the superfluid and the vortex lines is large enough to produce a Magnus force greater than the pinning threshold . If the Magnus force does not exceed , the pinning force adjusts to balance the Magnus force exactly. Given that the vortices are assumed to be inertialess, the excess force beyond the pinning threshold is of no relevance, as once the vortices unpin they flow with the ambient superfluid:
(8) 
In our model, the superfluid velocity at a given point is the sum of three components: the contribution from all pinned vortices that lie on or within the starcentred circle passing through the point (these vortices are assumed to be in a regular Feynman array on average); the contribution from all of the unpinned vortices lying within the above circle (also assumed to be in a regular Feynman array; i.e. meanfield approximation); and the contribution from the eight nearestneighbour cells on the grid, as depicted in Fig. 2:
(9) 
The first two contributions, hereafter named the meanfield terms, are calculated using Eq. (4), assuming that the velocity field is exactly tangential to the closed path (a good approximation which becomes exact in the thermodynamic limit ).
By assuming that the vortices are arranged in a regular array, we can evaluate (4) with tangential , without knowing the location of each vortex (the validity of this approach is verified in §5.1). The Magnus force is then calculated from Eq. (3), with equal to the veloicty of the stellar crust and given by Eq. (9). Initially, the mean Magnus force is very small. As the star spins down, however, the pinned contribution does not diminish, while the unpinned contribution does, and the mean increases.
We emphasise that the automaton represents a locally averaged model where each cell covers many pinning sites. Most of these sites are unoccupied because, for pinning at lattice defects (cf. seismic faults), the total number of vortices is much less than the number of pinning sites.
4.4 Driver
In this paper we consider two separate driving forces. Over long timescales (), the subject of §5.7, the Magnus force ramps up as the lag builds up between the pinned vortices and the unpinned superfluid. Over short timescales, we include thermally activated unpinning of the vortex bundles in random cells.^{4}^{4}4The contents of a randomly selected cell are deposited in the eight neighbouring cells. A more realistic implementation may be to select the thermally unpinned cell only from those cells that are already close to the pinning threshold, and therefore most likely to unpin thermally. In both scenarios (short and long timescales), small variations in due to the eight nearestneighbour cells trigger avalanches in the system. By treating the nearestneighbour interactions with greater precision than the longer range contributions, we are emphasizing the property outlined in condition 1 of §2.
4.5 Operational definition of glitches
In the automaton model a glitch corresponds to the unpinning of vortex bundles. Its size is determined by the total number of bundles (and hence vortices) unpinned. The associated fractional spin up of the pulsar is then given by the fraction of the total moment of inertia carried by the recently unpinned fluid. For example, if each vortex bundle represents the unpinning of fluid carrying of the moment of inertia of the system, then an avalanche of ten bundles results in a glitch of size .
In practice, the angular momentum transferred to the crust by vortex unpinning depends on the product of the number of vortices unpinned (equivalently, the reduction in superfluid rotation rate) and the moment of inertia of the region through which they move before repinning, or until the avalanche stops, reaching a new quasiequilibrium state. The reason for this is found in the OnsagerFeynman relation (4). The superfluid velocity at some distance from the rotation axis is proportional to the number of vortices enclosed, so the further the unpinned vortices travel radially, the more the superfluid decelerates, and hence the more angular momentum is transferred to the crust. The definition of glitch size given in this paper is therefore an approximation, valid only when three conditions are met: (i) all vortices travel the same radial distance during an iteration of the automaton; (ii) avalanches peter out after travelling one or two grid cells radially; and (iii) there are no steep radial gradients in vortex density (trivially satisfied by our uniform automaton, but not necessarily realistic). We show in §5.4 that these conditions are fulfilled by our automaton primarily because the unpinning activity is restricted to a narrow annulus. However, a more sophisticated automaton which incorporates radial gradients is needed to address this issue properly.
To facilitate future comparisons between the model and data, we clarify that the automaton describes spinup events where the jump in spin frequency is unresolved by current timing experiments; the maximum rise time consistent with existing data is 40 s (Dodson et al., 2002). Such glitches are accompanied by several recovery time scales, the longest of which may extend to the next glitch, together with a simultaneous jump in spindown rate. The automaton applies to this sort of unresolved spinup event, although it does not seek to model the recovery physics. Within the database of observed glitches, there are some glitches that do not conform to all facets of the above definition. In general, these nonstandard glitches are observed in pulsars that have only glitched one or two times. They are not described by the automaton in its current form, but they are usually large and hence relatively rare, so their impact on the statistics is muted. Another kind of nonstandard glitch, which we make no attempt to model, is the timeresolved secondary spin up (‘aftershock’) seen in the Crab to days after a standard glitch (Wong et al., 2001).
4.6 Automaton rules
The rules of our cellular automaton encompass the key features of SOC outlined in §2: discrete elements dominated by local forces, a threshold above which the elements move, and a slow driving force. The algorithm that governs the dynamics of the model can be summarised as follows:

The force imbalance is calculated at each cell based on the scheme outlined in §4.3. Cells where the force imbalance is positive are flagged as supercritical.

Half the vortex bundles in a supercritical cell are repositioned on the adjacent cell that lies in the path of the bulk superfluid velocity vector, which may lie on the diagonal, leaving the supercritical cell with half its original contents. Steps 2 and 3 are performed simulataneously, so that no bias arises from the order in which the cells are considered.

The number of vortex bundles unpinned is recorded cumulatively.

The total number of vortex bundles unpinned is recorded as the glitch size for the current big time step (where ‘big’ is defined below).


Either the vortex bundles occupying a randomly selected cell are redistributed into the neighbouring cells; or

the spin frequency of the crust (and hence ) is decremented by , the number of unpinned vortices is proprotionally reduced, and the number of pinned vortices remains unchanged.

Throughout this paper we refer to the completion of steps 1–4 as a small time step, and steps 1–5 as a big time step. Roughly speaking, a small time step lasts for the time it takes a vortex to be advected from one cell to its neighbour by the local superfluid flow. Of course, this varies somewhat with position. A big time step corresponds roughly to the mean waiting time between glitches, but please note that an avalanche does not always occur at every big time step given the above rules. The total duration of the simulation depends on whether the driver is thermally activated vortex creep (), or electromagnetic spin down () (see §4.3).
4.7 Vortex motion
The rules for moving vortices to an adjacent cell following an unpinning event rely on the (incorrect) assumption that is the same in the vicinity of all supercritical cells, such that the time taken to travel from one pinning site to the next is the same throughout the grid. That is, the discrete nature of the automaton means that the small time step sets the characteristic intercell travel time. This approximation is justified in two ways in the model. First, as long as the time to move one cell is much shorter than the driver timescale, it is reasonable to set the small time step using the slowest vortices. If, as we claim in §5.1, most of the avalanche activity takes place in a thin annulus, this approach is accurate, because is roughly uniform over the annulus. Second, only half the vortices in a supercritical cell are redistributed to its neighbours. The travel time for vortices at the far side of the supercritical (departure) cell (relative to the destination cell) is greater than the travel time for vortices starting near the border. Crudely, if vortices are evenly spaced throughout the departure cell, about half of them have sufficient time to reach the destination cell during a small time step. As a corollary, the model contains a dissipative mechanism, which ensures that avalanches eventually terminate. This feature means that our system is nonconservative. Following Pruessner & Jensen (2002), we require an external driving force to maintain the system in a critical state.
5 Statistics
5.1 General features
Parameter  Symbol  Default value 

Grid dimension  100  
Stellar radius  
Initial spin frequency  
Current spin frequency  
Spin down rate  
Total vortices in simulation  
Vortices per bundle  
Bundle factor (avg. bundles per cell)  1  
Pinning threshold  
Pinned fraction  0.01 
Three quantities in particular are used to characterise our cellular automaton model of pulsar glitches: the distribution of glitch sizes , durations , and waiting times . The glitch duration is defined as the number of consecutive time steps at which there exist supercritical cells.
The hypothesis that the automaton results in avalanche dynamics, driven slowly by electromagnetic spin down or random thermal unpinning, is now tested. The model is allowed to run for big time steps. In §5.2–§5.6 we consider time periods over which does not decrease significantly. In §5.7 we consider longer time periods over which spin down becomes important. We characterise the model by exploring its response to changes in the controlling parameters, both physical and computational. The physical parameters are the spin frequency , the pinning threshold , and the fraction of vortices that are pinned . The compuational parameters are the numbers of bundles and grid cells .
Fig. 3 shows the raw output of the model in the form of a time series of glitch sizes for a standard set of parameters (defined Table 2). Unless otherwise specified, these are the parameters used to generate the model output. This particular set of parameters is chosen to balance the need for sufficient avalanches to generate good statistics versus computational efficiency. The vertical axis gives the size in terms of the fraction of the total superfluid moment of inertia that unpins during the avalanche. The inset zooms in on a segment of the time series spanning big time steps. Although the simulation runs for a total of time steps (without spin down), the first time steps are discarded to allow the system to establish stationary statistics. Jensen (1998) asserted that fractal dynamics necessarily arise in a SOCS. Although we do not endeavour to quantify this mathematically, we are able to demonstrate, by considering time series such as the one shown in Fig. 3, that our system produces selfsimilar dynamics; qualitatively, the time series looks the same whether we plot or big time steps.
The lower panel of Fig. 3 shows the probability density functions of glitch sizes , durations , and waiting times , plotted as histograms, for the standard set of parameters defined in the caption. Both and obey power law statistics, with probability density functions and respectively. This is a characteristic feature of a SOCS (see §2.1). For the standard parameters in Fig. 3, we obtain and . More generally, for the parameter ranges studied in this paper, we find and (see §5.2–§5.6). In a SOCS, and are related by (Jensen, 1998)
(10) 
such that the size of an avalanche scales with its linear extent and and the duration scales as . The exponents in Fig. 3 imply , which is exactly what one expects for compact, twodimensional avalanches with and (see §5.7). Even more encouragingly, we find in many of our numerical experiments where criticality applies, even though and individually cover wide ranges.
Based on the data analyzed in Melatos et al. (2007), as well as the likeness between our model and other cellular automata that display selforganized critical behaviour, we expect the distribution of waiting times to be exponential, in keeping with a Poisson process (statistically independent avalanches). The bottom right panel in Fig. 3 confirms this expectation: the waitingtime probability density function is , with in units of an inverse big time step. To further verify the property of statistical independence, we calculate the linear Pearson correlation coefficient relating the size of an avalanche and the size of the avalanche immediately preceding it, as plotted in Fig. 4. We find that for avalanches, the correlation coefficient is , implying that there is almost no correlation between an avalanche and its predecessor.
Contrary to the general idea that a SOCS tends naturally to its critical state, we are forced to fine tune the parameters to elicit avalanches. The system is particularly sensitive to and . Similarly, the exponents of the best power law fits to the glitch size and duration distributions (see §5.2 and §5.6) depend on the same tuned parameters. This is also observed to be the case experimentally when superconducting vortices respond to a change in the ramp rate of the magnetic field (Field et al., 1995). In cellular automaton models of forest fires, the number of trees between two ignitions needs to be tuned according to system size, to avoid a cutoff or bump in the size distribution (Pruessner & Jensen, 2002).
Before analysing the results produced by the cellular automaton, we assess the accuracy of the meanfield approximation outlined in §4.3. Fig. 6 shows the velocity imbalance across the grid (as both a vector field and along a linear cut) generated using the meanfield approximation for a random distribution of vortex bundles (with a mean occupancy of one bundle per cell). The result is compared with an exact, ‘Nbody’ algorithm, in which at any cell is computed by summing the produced at that cell by every bundle on the grid. The exact algorithm is not used elsewhere in this paper due to computational cost. Discrepancies between the meanfield and exact calculations are greater when a line is taken vertically through the velocity field, rather than diagonally. This stems from the discrete, rectilinear nature of the grid. The two calculations are also in greater agreement further from the origin. By considering in the vicinity of the critical circle, where the meanfield contribution to the Magnus force approximately equals the pinning threshold (vertical dashed lines in right panel of Fig. 6), we see that the nearestneighbour and gridding fluctuations are comparable. Yet the effect of gridding on the dynamics turns out to be negligible, because the output of the model is a power law in . If the fluctuations were dominated by gridding effects, we would not expect power law statistics.
The bold circle overlaid on the velocity field in Fig. 6 is called the critical circle. Well beyond this circle, the meanfield contribution consistently exceeds the pinning threshold, such that the nearestneighbour variations are insufficient to alter the supercritical state of the cell. Out here, vortex bundles never pin. Similarly, well within this radius, the nearestneighbour variations are insufficient to alter the subcritical state of the cell. In here, vortex bundles never unpin. Near the circle, however, the nearest neighbour fluctuations are pertinent; indeed, they cause the avalanches we observe. The location of the ‘active’ annulus thereby defined depends on the physical and the computational parameters that define the model.
5.2 Dynamic range
In §4.1 we motivate the choice of grid size by referring to the observed dynamic range of pulsar glitches. In any individual pulsar the maximum dynamic range of is (for PSR J1704–3015), and the minimum is (for PSR J0537–6910 ^{5}^{5}5This is one of two pulsars that glitch quasiperiodically; the other is Vela.). Hence, from (6), we run our simulation on grids ranging from to (the physical size remains unchanged at ). By construction, the range of avalanches outputted by the model should not exceed (i.e. for the grid, and for grid). Fig. 7 confirms this expectation. The top panels are histograms of the avalanches produced by systems with , and over time steps. The dynamic range is greatest for the largest grid () and smallest for the smallest grid (). It is not expected that these dynamic ranges span the full range , as many of the grid cells lie inactive, well within the critical circle discussed in §5.1. It is, however, encouraging that the larger grid sizes do in fact produce a larger dynamic range.
Each of the histograms in Fig. 7 is overlaid with a power law best fit of the form , estimated using a leastsquares algorithm, with the fitted exponent indicated on each plot. By the central limit theorem, the uncertainty in each bin is inversely proportional to the square root of its contents. The fitted exponent is graphed as a function of in Fig. 5, with error bars taken from the leastsquares algorithm. There is a trend in both and to increase with increasing , in keeping with the increasing dynamic range with . When the simulation is run for time steps, the dynamic range does, as expected, increase; however, it still falls well short of . The number of bundles available to unpin is restricted to a small fraction () of that lie in the vicinity of the critical circle. To better match the desired dynamic range, we should select the grid size such that cells lie in the active zone near the critical circle.
The glitch sizes generated by the model should be treated as a indicative only. Evidently, several of the pulsars analyzed by Melatos et al. (2007) exhibit glitches several orders of magnitude smaller than those produced by the automaton. It is easy to elicit such small glitches from the automaton by adjusting or so that the active circle lies closer to the rotation axis, thereby enclosing fewer vortices. We have not explored these possibilities here as there is too much freedom in the choice of parameters at present. (See §5.5 for further discussion.) If pinning is strongest in two or more annuli within the crust, the model predicts that glitches fall into two or more size brackets, depending on where the annuli are located. It should be noted, however, that Melatos et al. (2007) found no evidence for a bimodal distribution when analyzing the latest statistics from individual pulsars, whose size distributions are consistent with unbroken power laws.
5.3 Vortex bundling
An underlying assumption of the model is that the scale on which the pinning centres and vortices are coarse grained does not alter the behaviour of the output qualitatively. In this section, we explore the impact of maintaining the grid size whilst varying the mean number of bundles per cell . The simulation is designed such that the smallest (largest) avalanche occurs when one (every) bundle unpins. In a scale invariant system, avalanches should occur below and above the observational limits. For this reason we probe the effect of increasing the number of vortex bundles (equivalent to decreasing the number of vortices per bundle, for fixed ) on the distribution of avalanches.
The avalanche size distributions for different values of the bundle factor are shown in the middle panels of Fig. 7. The histograms are overlaid with a power law best fit. For and there is an obvious ‘bump’ in the glitch size distribution, whose peak occurs at the same value of for different []. The bump is also present in the duration distribution (see Fig. 8), and there is a clear turnover in the slope of the waiting time distribution (see Fig. 10), which may be interpreted as a crossover from power law statistics to exponential decay. Importantly, crossover is evidence that the system is no longer in a SOC state, as there is a preferred scale for avalanches.
In SOCS generally, there is a systemsizedependent crossover scale above which the power law distribution of glitch sizes gives way to exponential decay (see §2.1). In a genuine critical state, the crossover scale should be an increasing function of the system size (Jensen, 1998). In the examples discussed by Jensen (1998), the system size corresponds to the number of particles in the system (rice or sand grains, for example), not the physical extent of the system. That the bump does not appear for small and becomes increasingly prominent for large suggests a system size effect. On the other hand, the crossover scale, defined by the bump, does not increase with system size (i.e. ) in Fig. 7. Nor does system size explain why the bump is a local maximum in instead of a turnover. In their experiments on superconductors Field et al. (1995) saw a crossover to steeper decay in the avalanche size distribution for large values of the magnetic field ramp rate, but not a local maximum. If instead we interpret the system size to be the number of cells, not bundles, then the bump becomes more prominent for larger systems (e.g. with ) as in Fig. 7. Despite the bump not being present for smaller systems (), there is a value of , for each , above which falls off faster than a power law, which is the expected behaviour beyond the crossover scale defined by Jensen (1998) (see Fig. 7). Bumps like the ones in Fig. 7 also arise in SOCS when one changes the morphology of the grains (e.g. rice versus sand) (Frette et al., 1996; Jensen, 1998; Pruessner & Jensen, 2003).
The violation of scale invariance is at least partly due to a mismatch in how the pinning centres (cells) and vortices (bundles) are renormalized. The most scale invariant results, closest in spirit to a SOCS, are obtained for (one bundle per cell, on average). By increasing without changing , we are renormalising the vortices on a smaller coarse grain, whilst maintaining the coarseness of the grid, which is not a justifiable physical scenario. Pruessner & Jensen (2002) claimed that inappropriately chosen parameters produce a bump in the size distribution, which further suggests that alters the fundamental statistical properties of the driven system. The implications for the microphysics of pinning in pulsars are explored in §6.
In summary, bundling in the model introduce a preferred scale into the system. (This scale greatly exceeds the vortex quantum , which is too small to be resolved by a practical simulation.) When the average number of bundles per cell is kept near one, the system exhibits scaleinvariant behaviour, as demonstrated by the power laws in glitch size. In this regime, the grid cells and bundles are well matched. If, however, the number of cells is significantly outnumbered by the vortex bundles, the scale invariance of the system is broken. In the absence of computational limitations, one should represent each pinning site and each vortex uniquely, eliminating the need for bundling. On the other hand, bundling confers a distinct advantage: it circumvents the subtleties of vortexvortex interactions and reconnections by renormalizing these effects into effective vortexvortex forces on a coarsened grid.
5.4 Pinning threshold
The pinning threshold determines the minimum lag between the superfluid and pinned vortices necessary to unpin vortex bundles. The pinning threshold can be reexpressed in terms of a minimum between the pinned vortices (comoving with the crust) and the unpinned superfluid.
As increases, the number of avalanches occurring within a fixed time ( big time steps) decreases, from for m s to for m s, as shown in Fig. 9. There are two reasons for this trend. First, greater imbalance is needed in the nearest neighbour cells to overcome the pinning force and trigger an avalanche. The greater is the required imbalance, the less probable it becomes. Second is the issue of fine tuning. For a given there is a critical circle where the meanfield contribution to the velocity imbalance matches the pinning threshold almost exactly. Well outside the circle, every cell is supercritical and the vortices never pin. Near the circle, however, we find a mixture of sub and supercritical cells because the nearest neighbour imbalance can act to reduce at a given cell, if it opposes the meanfield contribution. As increases, the critical circle moves towards the origin, and hence the number of cells in its vicinity decreases. Fig. 11 depicts images of the grid after big time steps. The colours encode the value of on each cell. Supercritical cells are white; all other cells are subcritical. The pinning threshold in the upper panel of Fig. 11 is m s. We observe a smattering of supercritical cells near the critical circle. Inside the circle, all cells are subcritical. In the lower panel of Fig. 11, with m s, we observe an annulus in which most cells are supercritical (where the vortices never pin). Just inside the annulus, there is, like in the upper panel, a mixture of sub and supercritical cells. It is in these mixed regions where we claim that selforganized criticality controls the unpinning dynamics. These results are evidence that the unpinning activity is restricted to a narrow annulus, and hence unaffected by the use of a full circle as the simulation grid, as discussed in § 4.5.
A power law describes the size distribution well for all three values of in Fig. 9. The best fit power law and its exponent are shown in each panel of Fig. 9. We find that decreases as increases, however there is a slight upturn in both and for the largest values of considered. For larger values of , avalanches involving many cells (and hence many vortex bundles) are less likely to occur, both because the active annulus shrinks, and because it is relatively unlikely that the nearestneighbour imbalance sends a cell supercritical. Hence decreases because the dynamic range decreases; the area under the histogram equals the number of avalanches.
Fig. 6 demonstrates that, for a finite, square vortex array, the profile of along a linear cut flattens out in the region where , compared to smaller radii. Hence if approaches the maximum on the grid, more cells may participate in an avalanche, in comparison to when the critical circle lies further in. This may account for the upturn in versus seen in Fig. 5.
5.5 Pinned fraction
The pinned fraction impacts on the dynamics in two main ways: it changes the size of each vortex bundle, and alters the relative contribution to the velocity imbalance from each term in Eq. (9). Increasing the size of the vortex bundles increases the minimum avalanche size, whilst maintaining the maximum avalanche size. Larger bundles also enhance the role of the nearest neighbour cells. As more vortices pin and unpin, the first and third terms in Eq. (9) decrease with respect to the contribution from the unpinned vortices. Given that the unpinned vortices are homogeneously distributed, they are the foremost contributors to the meanfield critical circle discussed in §5.4. Decreasing their relative contribution reduces the meanfield effect, emphasizing the nearestneighbour, stochastic dynamics of the model. The lower left panel of Fig. 5 supports this argument: the error in the fitted and values decreases with increasing . The lower left panel of Fig. 5 also shows that both and increase with increasing . For the distribution flattens out, whereas the distribution continues to increase.
For no avalanches are observed; the model is badly tuned. Moreover, for , the model slows down computationally, as a large number of cells regularly pin and unpin. This computational limitation arises because our algorithm does not excise regions in which vortices are permanently unpinned (see the lower panel of Fig. 11), instead calculating unneccessarily in these regions at each small time step. We intend to address this flaw in future work.
5.6 Spin frequency
Two values of the spin frequency enter the model: at birth (), and today (). As a first approximation, the total number of pinned vortices is taken to be unchanged over the lifetime of the star and hence is dictated by as in Eq. (5). The pinned portion of the superfluid does not spin down with the stellar crust, whereas the unpinned superfluid does; its circulation is dictated by . In the current model we hold and constant (cf. §5.7), with Hz, and .
By decreasing , we increase and hence the force imbalance at every cell. The impact is similar to decreasing because the critical circle shrinks. The two righthand panels of Fig. 5 display a similar trend for increasing and . The main difference is that also determines the relative size of the first two terms in Eq. (9), given that the number of pinned vortices remains constant while the number of unpinned vortices decreases with . These effects increase the range of avalanche sizes as decreases as shown in Fig. 9. For large values of , the dynamic range goes to zero, as the critical circle has moved outside the star. For small , a turnover emerges in the size distribution for small , which is also present in the size distributions for small values of in Fig. 9. This may happen because the critical circle lies well within the star, so there is a sizable region outside the critical circle where vortices are permanently unpinned. The automaton does not recognise that permanently unpinned vortices do not, in fact, participate in avalanches in real physical systems, so the output is biased towards large avalanches.
5.7 Longterm activity
To model the longterm glitch behaviour of a pulsar, we must allow to mimic the electromagnetic spin down of the stellar crust. Observed spin down rates differ by several orders of magnitude; in this paper we take by way of example. Observations of pulsar glitches span a period of approximately forty years (1970 to present), which equates to a total spin down of Hz. In our model, each big time step must therefore span s, in order to represent a total observing time of forty years. Such a large time step determines the scale on which we can resolve individual avalanches. Observationally, glitches can be resolved down to a s window (McCulloch et al., 1990; Dodson et al., 2002).
Spin down causes the critical circle to shrink towards the axis of the star, because the number of pinned vortices is fixed, whilst the number of unpinned vortices decreases with . The rate of spin down thus determines how quickly the glitch behavior changes as the critical circle moves through the strata of the star. Fig. 12 shows the size distribution () for time steps , and s. Fig. 13 shows the fitted power law exponents for the sizes and durations, for . Within the uncertainty of the fitting algorithm, and are the same for all values of tested. For the largest , the time elapsed after big time steps is s, which equates to Hz. This moves the critical circle towards the origin by a negligible .
The long term fate of a pulsar, with respect to its glitch behaviour, depends on the location of its pinning zones. If pinning can occur througout the pulsar, glitches continue to occur indefinitely, as the critical circle moves inwards. As time passes, however, the glitches become smaller, on average, as the number of cells in the vicinity of the critical circle diminishes ^{6}^{6}6Interestingly, there is no guarantee that glitches emanating from different critical circles are, in fact, drawn from the same statistical distribution (e.g. value of ). If so, we do not expect a pure power law in over long times.. If, as claimed by Donati & Pizzochero (2003), pinning occurs at intermediate densities, no glitches occur once the critical circle shrinks inside the radius at which these densities occur. Alternatively, if pinning occurs at several strata [high and low densities, for example (Avogadro et al., 2007)], glitches occur when the critical circle lies in the vicinity of the pinning zones, which could mean that episodes of glitching are punctuated by ‘quiet’ periods.
In the models of postglitch relaxation developed by Alpar et al. (1996), the superfluid spins down when vortices move radially outward. The outward motion is driven by the distribution of capacitive and depleted zones; outward motion is energetically favourable due to a radial bias in vortexvortex scattering. In our automaton, spin down is an input; it does not arise selfconsistently from vortex motion. By reducing the number of unpinned vortices with time (e.g. § 5.7) we attempt to recreate the scenario in which outward radial motion eventually makes vortices exit the pinning zones in the star. The automaton rules do inadvertently allow for outward radial motion thanks to the rectangular nature of the grid: azimuthal motion (i.e. in the direction of rigid body rotation) is never truly possible, as the vortices can only ever move horizontally, vertically, or diagonally. To incorporate outward motion and hence spin down selfconsistently, the automaton rules should add a radial component to the bulk superfluid flow when computing . (This does not remove the artificial radial motion arising from the rectangular grid. In this context, an annular grid topology is preferable.)
6 Discussion
In this paper we present a cellular automaton model for glitches in neutron stars, based on the vortex unpinning paradigm (Anderson & Itoh, 1975; Alpar et al., 1984). Despite the gross idealisations in such a discrete model, the resulting statistical behaviour agrees with recent analyses of radio timing data (Wang et al., 2000; Wong et al., 2001; Melatos et al., 2007). Of particular interest are the power law distributions of glitch sizes and durations produced by the model, providing further evidence that glitches are a scale invariant phenomenon. For this reason, as observations improve in resolution and lengthen in duration, we expect to discover smaller and larger glitches than have been seen so far. Moreover, the exponential waitingtime distribution generated by the model supports the idea that the each vortex avalanche is a statistically independent event. The values of the power law exponents and Poisson glitching rates depend on the particular parameter set. Our fits return size and duration exponents in the ranges and respectively, and mean glitching rates in the range . Although we do not have a firstprinciples theory against which to compare our results, we claim that our cellular automaton is consistent with a SOCS.
The operation of the cellular automaton relies heavily on the assumption that vortices are inhomogeneously distributed. Alpar et al. (1996) suggested that crust cracking creates capacitive regions where many vortices pin simultaneously. These regions do not permit a vortex current, so pinned vortices in these regions do not participate in outward vortex creep. Crust cracking is a nonuniform process, suggesting that the regions themselves are inhomogeneously distributed. Our model conforms qualitatively with this picture. Of particular interest is the claim by Alpar et al. (1995) that the neutron star crust is divided coarsely into depleted and capacitive regions. If verified, this claim suggests that the level of coarse graining necessary to render our cellular automaton computationally tractable matches the level of inhomogeneity that is realistically presesnt in a neutron star.
Of course, the opposite viewpoint is also viable: the crust forms a large, defectfree crystal, allowing the vortices to occupy a contiguous sequence of pinning positions and preventing the type of inhomogeneity represented by our initial conditions (Jones, 1998a). Importantly, such a scenario challenges any model based on collective unpinning of numerous vortices, not just our cellular automaton. De Blasio & Lazzari (1998) concluded that the crust is relatively free of microcrystalline structures and vacancies, with impurities, or point defects, the most likely disruptions to the crystal lattice (faults and dislocations were not considered). Alternatively, if pinning occurs mostly at grain boundaries and cracks, as originally postulated by Anderson & Itoh (1975), then inhomogeneity on macroscopic scales is possible.
Irrespective of the density and distribution of pinning sites in reality, the model assumes for simplicity that vortex pinning occurs throughout the pulsar. Several recent studies have shown that pinning forces sufficient to oppose the Magnus force develop only in specific zones of the crust, e.g. regions of high and low, rather than intermediate, denisty (Avogadro et al., 2007). Restricting the simulation to these optimal pinning zones obviates the need to allow for a range of pinning strengths if the zones are thin layers. A model in which pinning does not occur everywhere in the star leads to a different computational topology, e.g. an annulus, or several nested annuli. However, we choose to be conservative in this introductory paper by not imposing a grid topology a priori, in order to see if a preferred topology emerges by itself. In fact, it does; we find that a narrow active annulus, situated at the critical radius, participates in the avalanches. However, it is important to experiment with other topologies in future work, because the active annulus we find in this paper may still be an artifact. If the active region lies near the rotation axis, an automaton executed in an annulus is free of the shortcomings introduced by the ‘permanently unpinned’ regions mentioned in §5.5. To use the full circular grid in a meaningful way, we must include pinning of proton and neutron vortices in the core of the neutron star, together with pinning of neutron vortices in the inner crust. In this case, the additional physics of entrainment (Link & Epstein, 1996; Ruderman et al., 1998; Sedrakian & Sedrakian, 1995; Andersson et al., 2004) and induced electric currents implies a more detailed set of automaton rules, outside the scope of the model presented here. It should also be noted that if the superfluid vortices are tangled (Peralta et al., 2006b, a), rather than forming a regular Abrikosov array, these twodimensional annular and circular topologies are clearly too restrictive.
The locations of pinning zones also determine the longterm fate of a pulsar as it spins down. If pinning occurs everywhere, spin down increases the proportion of vortices with which are always supercritical (outside the critical circle) and do not participate in avalanche dynamics as they never pin. If pinning is restricted to one or more annuli, the glitch behaviour of the pulsar changes significantly when the critical circle moves into and out of the pinning zones.
Clearly, our cellular automaton does not include the response of the pulsar to glitches (ie. the spin up of the pulsar crust and the subsequent ‘relaxation’). Nor do we account for the finite timescale on which the unpinned superfluid couples to the crust, governed by the HallVinenBekarevichKhalatnikov equations (Peralta et al., 2005, 2006b). We make the approximation that this timescale is considerably less than the timescale of our big time steps, motivated by glitch timing data (where the postglitch recovery phase typically ends well before the next glitch; cf. Vela).
To elicit avalanches from our model, we require fine tuning in both the physical and computational parameters, such that . This condition ensures that there are enough vortex bundles that switch between becoming sub and supercritical as time passes. We achieve this by choosing , and to place the critical circle near the surface of the star. We emphasize that for a larger (smaller) star, we would resize the critical circle proportinally, for computational rather than physical reasons. Fine tuning is also required in the ratio of vortex bundles to grid cells . For , the automaton output is no longer consistent with a SOCS.
In conclusion, we present an empirical cellular automaton model of pulsar glitches based on the mass vortex unpinning paradigm. We find that for certain physical and computational parameters the model produces dynamics that are consistent with a SOCS and with radio timing data from pulsars. There exists no general, firstprinciples theory of SOC , let alone of pulsar glitches, so many of our results are empirical. In particular, ther is no way known at present to predict theoretically the size and duration distribution exponents, and the mean rate of the waitingtime distribution. We do demonstrate that the basic physical principles governing intervortex interactions can produce the type of collective behaviour necessary to explain pulsar glitches within the mass unpinning paradigm, especially the puzzle of how so many vortices can unpin in sympathy during a glitch, and why their number varies so much from glitch to glitch.
We thank Carlos Peralta for sharing his uptodate glitch catalogue, and Stuart Wyithe for advice on statistical methods. LW acknowledges the support of an Australian Postgraduate Award.
References
 Alpar et al. (1984) Alpar M. A., Anderson P. W., Pines D., Shaham J., 1984, ApJ, 276, 325
 Alpar et al. (1996) Alpar M. A., Chau H. F., Cheng K. S., Pines D., 1996, ApJ, 459, 706
 Alpar et al. (1989) Alpar M. A., Cheng K. S., Pines D., 1989, ApJ, 346, 823
 Alpar et al. (1995) Alpar M. A., Kiziloglu U., van Paradijs J., 1995, The lives of the neutron stars. Dordrecht ; Boston : Kluwer Academic, c1995.
 Anderson & Itoh (1975) Anderson P. W., Itoh N., 1975, Nature, 256, 25
 Andersson et al. (2004) Andersson N., Comer G. L., Prix R., 2004, MNRAS, 354, 101
 Avogadro et al. (2007) Avogadro P., Barranco F., Broglia R. A., Vigezzi E., 2007, Phys. Rev. C, 75
 Bak (1996) Bak P., 1996, How nature works: the science of selforganized criticality. Copernicus, New York, NY, USA
 Bak et al. (1988) Bak P., Tang C., Wiesenfeld K., 1988, Phys. Rev. A, 38, 364
 Bassler & Paczuski (1998) Bassler K. E., Paczuski M., 1998, Phys. Rev. Lett., 81, 3761
 De Blasio & Elgarøy (1999) De Blasio F. V., Elgarøy O., 1999, Phys. Rev. Lett., 82, 1815
 De Blasio & Lazzari (1998) De Blasio F. V., Lazzari G., 1998, Nucl. Phys. A, 633, 391
 Dodson et al. (2002) Dodson R. G., McCulloch P. M., Lewis D. R., 2002, ApJ, 564, L85
 Donati & Pizzochero (2003) Donati P., Pizzochero P. M., 2003, Phys. Rev. Lett., 90
 Donati & Pizzochero (2006) Donati P., Pizzochero P. M., 2006, Phys. Rev. B, 640, 74
 Donnelly (1991) Donnelly R. J., 1991, Quantized vortices in Helium II. Press Syndicate of the University of Cambridge, Cambridge, New York NY, USA
 Elgarøy & De Blasio (2001) Elgarøy O., De Blasio F. V., 2001, A&A, 370, 939
 Field et al. (1995) Field S., Witt J., Nori F., 1995, Phys. Rev. Lett., 74, 1206
 Frette et al. (1996) Frette V., Christensen K., MaltheSørensen A., Feder J., Jøssang T., Meakin P., 1996, Nature, 379, 49
 Jensen (1990) Jensen H. J., 1990, Phys. Rev. Lett., 64, 3103
 Jensen (1998) Jensen H. J., 1998, Selforganized criticality: Emergent complex behavior in physical and bilogical systems. Cambridge University Press, Cambridge UK
 Jones (1997) Jones P. B., 1997, Phys. Rev. Lett., 79, 792
 Jones (1998a) Jones P. B., 1998a, Phys. Rev. Lett., 81, 4560
 Jones (1998b) Jones P. B., 1998b, MNRAS, 296, 217
 Linder et al. (2004) Linder J. F., Hughes S. B., Miller D. J., Thomas B. C., Weisenfeld K., 2004, Int. J. Bifurcation and Chaos, 14, 1155
 Link & Cutler (2002) Link B., Cutler C., 2002, MNRAS, 336, 211
 Link & Epstein (1996) Link B., Epstein R. I., 1996, ApJ, 457, 844
 Lu & Hamilton (1991) Lu E. T., Hamilton R. J., 1991, ApJ, 380, L89
 Mastrano & Melatos (2005) Mastrano A., Melatos A., 2005, MNRAS, 361, 927
 McCulloch et al. (1990) McCulloch P. M., Hamilton P. A., McConnell D., King E. A., 1990, Nature, 346, 822
 Melatos & Peralta (2007) Melatos A., Peralta C., 2007, ApJ, 662, L99
 Melatos et al. (2007) Melatos A., Peralta C., Wyithe J. S. B., 2007, MNRAS
 Middleditch et al. (2006) Middleditch J., Marshall F. E., Wang Q. D., Gotthelf E. V., Zhang W., 2006, ApJ, 652, 1531
 Morley (1996) Morley P. D., 1996, A&AS, 313, 204
 Morley & Schmidt (1996) Morley P. D., Schmidt I., 1996, Europhys. Lett., 33, 105
 Nicodemi & Jensen (2001a) Nicodemi M., Jensen H. J., 2001a, J. Phys. A: Math. and Gen., 34, 8425
 Nicodemi & Jensen (2001b) Nicodemi M., Jensen H. J., 2001b, Phys. Rev. Lett., 86
 Peralta et al. (2005) Peralta C., Melatos A., Giacobello M., Ooi A., 2005, ApJ, 635, 1224
 Peralta et al. (2006a) Peralta C., Melatos A., Giacobello M., Ooi A., 2006a, ApJ, 635, L53
 Peralta et al. (2006b) Peralta C., Melatos A., Giacobello M., Ooi A., 2006b, ApJ, 651, 1079
 Pines & Alpar (1985) Pines D., Alpar M. A., 1985, Nature, 316, 27
 Pruessner & Jensen (2002) Pruessner G., Jensen H. J., 2002, Europhys. Lett., 58, 250
 Pruessner & Jensen (2003) Pruessner G., Jensen H. J., 2003, Phys. Rev. Lett., 91
 Ruderman et al. (1998) Ruderman M., Zhu T., Chen K., 1998, ApJ, 492, 267
 Ruderman (1969) Ruderman R., 1969, Nature, 223, 597
 Sedrakian & Sedrakian (1995) Sedrakian A. D., Sedrakian D. M., 1995, ApJ, 447, 305
 Sornette et al. (1991) Sornette D., Sornette A., Vanneste C., 1991, Selforganized criticality, earthquakes and plate tectonics. Vol. 392, Berlin Springer Verlag, Berlin/Heidelberg
 Vespignani & Zapperi (1998) Vespignani A., Zapperi S., 1998, Phys. Rev. E, 57, 6345
 Wang et al. (2000) Wang N., Manchester R. N., Pace R. T., Bailes M., Kaspi V. M., Stappers B. W., Lyne A. G., 2000, MNRAS, 317, 843
 Wiesenfeld et al. (1989) Wiesenfeld K., Tang C., Bak P., 1989, J. Stat. Phys., 54, 1441
 Wong et al. (2001) Wong T., Backer D. C., Lyne A. G., 2001, ApJ, 548, 447