A cellular automaton model of pulsar glitches

A cellular automaton model of pulsar glitches

L. Warszawski and A. Melatos
School of Physics, University of Melbourne, Parkville, VIC 3010, Australia
lilaw@ph.unimelb.edu.au (LW)amelatos@unimelb.edu.au (AM)

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 scale-invariant avalanches (Melatos et al., 2007), which are consistent with a self-organized 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.

pulsars: general — stars: interior — stars: neutron
pagerange: A cellular automaton model of pulsar glitchesReferencespubyear: 2007

1 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 self-organized 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, non-local 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 self-consistent, semi-analytical model for the vortex-nucleus 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 two-stream 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 self-organized 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 Self-organized critical systems

2.1 General properties

A system is described as a SOCS based on two underlying behavioural patterns. Self-organized 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 self-organized critical (SOC) state it must first satisfy the following three conditions (Jensen, 1998; Melatos et al., 2007):

  1. It is composed of many discrete, mutually interacting elements, whose motions are dominated by local (e.g. nearest-neighbour) rather than global (e.g. mean-field) forces.

  2. 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.

  3. An external force drives the system slowly, in the sense that elements adjust to local forces rapidly compared to the driver time-scale. Combined with local thresholds, this ensures that the system evolves quasistatically through a history-dependent sequence of metastable states.

    If the above conditions are met, the following properties are generically observed (Jensen, 1998; Melatos et al., 2007):

  4. 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).

  5. 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).

  6. 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 non-conservative 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 inter-vortex 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

Figure 1: Cumulative distributions of glitch waiting times , for six pulsars with . The horizontal scale is different for each pulsar; is the mean glitch rate that minimises the K-S statistic for that pulsar as listed in Table 1. The solid curve represents the ideal, Poissonian, cumulative distribution, . The waiting time distribution in an ideal SOCS follows this Poissonian.
1825-0935 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
1801-2304 0.55 0.35 0.88 0.57 0.092 1.1 4
1341-6220 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
1740-3015 1.5 1.2 2.5 1.1 0.98 1.3 0.7
Table 1: Parameters of the glitch size and waiting time distributions for pulsars with (Melatos et al., 2007). and are chosen to minimise the K-S statistic.

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


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 nearest-neighbour 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 Kolmogorov-Smirnov (K-S) test. The range of is presented in Table 1. Table 1 also gives the exponent that minimises the K-S 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 K-S 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 K-S 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


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 K-S 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 force111The 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


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 two-fold. 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 ‘stick-slip’ 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 time-scale of the readjustment is much shorter than the driving time-scale 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 Donnelly-Glaberson instability, with a net polarisation parallel to the rotation axis [see for example Peralta et al. (2005, 2006a, 2006b)].

4 Cellular automaton

Figure 2: Schematic of the cellular automaton. The left side of the figure is a cartoon of the grid, illustrating how the force imbalance on a cell is calculated. At the cell centred on the tip of the white arrow, the superfluid velocity field is calculated as the sum of three components: the circular velocity due to the the pinned vortices residing in the grey shaded circle; the circular velocity due to the unpinned vortices in the grey region; and the summed velocity field due to vortex bundles in the eight nearest-neighbour and next-nearest-neighbour cells (bundles are represented as black dots). The right side of the figure describes how the force imbalance due to the nearest-neighbour vortices, , is calculated. Each vortex bundle generates a velocity field according to a right-hand rule, with the circulation vector pointing out of the page, which falls off inversely with distance from the centre of the bundle.

The cellular automaton model presented in this paper takes the basic physical features of the inner crust-superfluid 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 large-scale, 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


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


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


Pulsar timing data suggests, via (6), a range of grid sizes varying from to a maximum of  222 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


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 self-organized 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 (coarse-grain 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 random333Vortex 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:


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 star-centred 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. mean-field approximation); and the contribution from the eight nearest-neighbour cells on the grid, as depicted in Fig. 2:


The first two contributions, hereafter named the mean-field 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 time-scales (), 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 time-scales, we include thermally activated unpinning of the vortex bundles in random cells.444The 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 time-scales), small variations in due to the eight nearest-neighbour cells trigger avalanches in the system. By treating the nearest-neighbour 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 quasi-equilibrium state. The reason for this is found in the Onsager-Feynman 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 spin-up 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 spin-down rate. The automaton applies to this sort of unresolved spin-up 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 time-resolved 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:

  1. The grid is initialised by randomly placing bundles of vortices on the cells that lie within the star, as in §4.1 and §4.2.

  2. 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.

  3. 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.

  4. The number of vortex bundles unpinned is recorded cumulatively.

  5. Steps 14 are repeated until no cells are supercritical.

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

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

    2. 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.

  7. Steps 17 are repeated for a predetermined number of time steps, corresponding to either an arbitrary period of time [if 7a is used] or the number of time steps multiplied by (where is the length of each time step), if 7b is used.

Throughout this paper we refer to the completion of steps 14 as a small time step, and steps 15 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 inter-cell 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 time-scale, 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 non-conservative. 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
Table 2: Default physical and computational parameters used to analyse the model.

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 self-similar 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)


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, two-dimensional 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 self-organized 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 waiting-time 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 mean-field 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 mean-field approximation for a random distribution of vortex bundles (with a mean occupancy of one bundle per cell). The result is compared with an exact, ‘N-body’ 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 mean-field 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 mean-field contribution to the Magnus force approximately equals the pinning threshold (vertical dashed lines in right panel of Fig. 6), we see that the nearest-neighbour 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 mean-field contribution consistently exceeds the pinning threshold, such that the nearest-neighbour variations are insufficient to alter the supercritical state of the cell. Out here, vortex bundles never pin. Similarly, well within this radius, the nearest-neighbour 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.

Figure 3: Top: Automaton output for big time steps. The radius of the star is  km, the unpinning threshold is  ms, the neutron superfluid density is  kg m, the initial spin frequency is  Hz, and the current spin frequency is  Hz. The avalanche size is given as the fraction of the total fluid moment of inertia that unpins during the avalanche. The fraction of pinned vortices is taken to be . The main panel displays a time series of avalanche sizes beginning after time steps elapse, in order to ensure that a self-organized critical state has time to emerge from the initially random conditions. The inset zooms in on time steps. Bottom: Probability density funtions of (left), durations (centre) and (right) plotted as histograms on a log-log (left and centre) and log-linear (right) scale, using logarithmically (left and centre) and linearly (right) spaced bins. The same parameters are used as for the top panel. The size and duration distributions are fitted with power laws of the form and respectively. The waiting time distribution is fitted with an exponential of the form , where is in units of the reciprocal of one big time step. All fitted functions are plotted as dashed lines.
Figure 4: Next avalanche size versus previous avalanche size for avalanches. is the linear Pearson correlation coefficient. The automaton shows no sign of ‘memory’ as .
Figure 5: Fitted power law exponents for the distribution of glitch sizes () and durations (). The error bars represent the uncertainty returned by the least-squares fitting algorithm. Top left: . Top right: . Bottom left: . Bottom right: . Automaton parameters: , , big time steps. Physical parameters: , , .
Figure 6: Left: A vector plot of the velocity imbalance , generated using the mean-field approximation described in Eq. (9). The inset zooms in on the central grid cells. The solid black circle shows the position of the critical radius, where the mean-field contribution to is approximately equal to . Right: Profiles of along a straight-line cut through the grid, as a function of distance from the origin, as calculated by the mean-field (black) and exact (N-body) (grey) algorithms. The solid and dashed curves correspond to the vertical and diagonal cuts drawn in the left panel. The two vertical dashed lines show the position of the critical circle. The horizontal dashed line marks .

5.2 Dynamic range

Figure 7: Probability density functions of avalanche sizes , plotted as histograms on a log-log scale, using logarithmically spaced bins, for big time steps. The histograms are fitted with a power law of the form (dashed lines). Top: Grids with (left), (centre) and (right). All panels have bundle factor . Middle: Bundle factors (left), (centre) and (right). All panels have . Bottom: and (left), (centre) and (right). All panels have . The vertical grey line denotes the crossover scale. Physical parameters: , , , , , .
Figure 8: Probability density functions of avalanche durations , for bundle factors (left), (centre) and (right), plotted as histograms on a log-log scale, using logarithmically spaced bins. The histograms are fitted with a power law of the form (dashed lines). Automaton parameters: , big time steps. Physical parameters: , , , , , .
Figure 9: Probability density functions of avalanche sizes , plotted as histograms on a log-log scale, using logarithmically spaced bins. The histograms are fitted with a power law of the form (dashed lines). Top: Current spin frequency (left), (middle) and (right). Bottom: Pinning threshold (left), (middle) and (right). Automaton parameters: , big time steps, . Physical parameters: , , , .
Figure 10: Probability density functions of waiting times , plotted as histograms on a log-linear scale, using linearly spaced bins, for big time steps. The histograms are fitted with an exponential of the form (dashed lines). Top: Grid size (left), (centre) and (right). Bundle factor . Bottom: Bundle factor (left), (centre) and (right). Grid size . Physical parameters: , , , , , .

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 555This 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 least-squares 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 least-squares 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 system-size-dependent 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 scale-invariant 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 vortex-vortex interactions and reconnections by renormalizing these effects into effective vortex-vortex forces on a coarsened grid.

5.4 Pinning threshold

Figure 11: Images of the velocity imbalance for (upper) and (lower). White cells are supercritical; the darker the cell, the smaller the .

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 mean-field 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 mean-field 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 self-organized 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 nearest-neighbour 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 mean-field critical circle discussed in §5.4. Decreasing their relative contribution reduces the mean-field effect, emphasizing the nearest-neighbour, 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 right-hand 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 Long-term activity

Figure 12: Probability density functions of avalanche sizes , for time steps lasting (left), (centre) and (right), plotted as histograms on a log-log scale, using logarithmically spaced bins. The histograms are fitted with a power law of the form (dashed lines). Automaton parameters: , big time steps. Physical parameters: , , , , , .
Figure 13: Fitted power law exponents for the distribution of glitch sizes () and durations (), for time steps in the range . The error bars represent the uncertainty returned by the least-squares fitting algorithm. Automaton parameters: , big time steps. Physical parameters: , , , , , .

To model the long-term 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 666Interestingly, 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 post-glitch 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 vortex-vortex scattering. In our automaton, spin down is an input; it does not arise self-consistently 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 self-consistently, 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 waiting-time 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 first-principles 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, defect-free 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 two-dimensional annular and circular topologies are clearly too restrictive.

The locations of pinning zones also determine the long-term 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 time-scale on which the unpinned superfluid couples to the crust, governed by the Hall-Vinen-Bekarevich-Khalatnikov equations (Peralta et al., 2005, 2006b). We make the approximation that this time-scale is considerably less than the time-scale of our big time steps, motivated by glitch timing data (where the post-glitch 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, first-principles 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 waiting-time distribution. We do demonstrate that the basic physical principles governing inter-vortex 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 up-to-date glitch catalogue, and Stuart Wyithe for advice on statistical methods. LW acknowledges the support of an Australian Postgraduate Award.


  • 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 self-organized 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., Malthe-Sø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, Self-organized 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, Self-organized 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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description