A Cellular Automaton Model of Pulsar Glitches

a r X i v :0806.4738v 1 [a s t r o -p h ] 29 J u n 2008

A cellular automaton model of pulsar glitches

L.Warszawski 1?and A.Melatos 1

?1

School of Physics,University of Melbourne,Parkville,VIC 3010,Australia

ABSTRACT

A cellular automaton model of pulsar glitches is described,based on the super?uid

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 e?cient means of modelling the collective behaviour of up to 1016vortices 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 ?tted power law exponents for the glitch sizes (a )and durations (b )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 ?nds ?4.3 a ?2.0and ?5.5 b ?2.2,and mean glitching rates in the range 0.0023 λ 0.13in units of inverse time.

Key words:pulsars:general —stars:interior —stars:neutron

1INTRODUCTION

It has recently been shown that glitches in individual pulsars obey avalanche statistics (Melatos et al.2007).Such statisti-cal 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 char-acterised 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 (10?11 δν/ν 10?4),and up to four decades in any individual pulsar.A natural way to model a scale invariant,avalanching system is to construct a cellular au-tomaton 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),granu-lar assemblies (Bak et al.1988;Wiesenfeld et al.1989;Bak 1996;Frette et al.1996;Pruessner &Jensen 2003),solar ?ares (Lu &Hamilton 1991),and superconducting vortices

?

lilaw@https://www.sodocs.net/doc/112634790.html,.au (L W)?amelatos@https://www.sodocs.net/doc/112634790.html,.au (AM)

(Jensen 1990;Field et al.1995;Bassler &Paczuski 1998;Linder et al.2004).

Theories of pulsar glitches centre around the mass unpinning model ?rst 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 super?uid in the stellar interior is threaded by many (～1016)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 super?uid.When the transverse Magnus force (di-rectly proportional to the lag)surpasses a threshold value (equal to the strength of the pinning force),a catastrophic unpinning of vortices occurs,transferring angular momen-tum to the crust.In order for this mechanism to generate glitches on the scale observed,it requires up to 1012vortices 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,b;De Blasio &Elgar?y 1999;Elgar?y &De Blasio

2L.Warszawski and A.Melatos

2001;Donati&Pizzochero2003;Avogadro et al.2007). In particular,the strength of the pinning force has re-cently 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 den-sity(Donati&Pizzochero2006).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 cen-tres available to the vortices,and their relative strengths (De Blasio&Lazzari1998).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 super?uid vor-tices and of the nuclear crustal lattice.

Several other theories of pulsar glitches,similarly fo-cused on the microphysics of the super?uid in the stel-lar interior and crust,are currently under consideration. These include thermally driven glitches arising from sud-den changes in the frictional coupling within the neu-tron star(Link&Epstein1996),super?uid turbulence (Peralta et al.2005,2006b,a;Melatos&Peralta2007),su-per?uid two-stream instabilities(Andersson et al.2004; Mastrano&Melatos2005),crust cracking(Ruderman1969; Bak1996;Alpar et al.1996;Middleditch et al.2006)and mass accretion(Morley&Schmidt1996;Morley1996).

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 de-scribe the observational data.In order to produce a realis-tic model of the global dynamics,some subtle aspects of the microphysics need to be sacri?ced.Cellular automa-ton models of complex systems aim to employ only the dominant aspects of the microphysics to generate replica-ble and e?cient automaton rules.Such an approach was attempted by Morley&Schmidt(1996)using a cellular au-tomaton platelet model of mass accretion onto the stellar crust.Although the microphysical details di?er from those treated in this paper,the underlying philosophy of using a simple model to describe large scale behaviour is iden-tical to ours.We are similarly encouraged by the marked parallels between super?uid vortices and those that pop-ulate hard,type II superconductors.The agreement be-tween cellular automaton models of?ux creep and vortex avalanches in superconductors and experimental data is ex-cellent(Linder et al.2004;Jensen1990;Field et al.1995; Bassler&Paczuski1998;Nicodemi&Jensen2001a,b),in particular with respect to the measurement of a power law over several decades in the avalanche(cf.glitch)size.

A general,?rst principles theory of self-organized criti-cality does not yet exist for any SOCS,let alone for pulsar glitches.In this paper we look for empirical agreement be-tween the output of our cellular automaton and the gross qualitative features in pulsar glitch data.A detailed quanti-tative comparison between a?rst 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 observa-tional pulsar glitch data in 2.3reviews the vortex unpin-ning paradigm and recent advances in the understanding of neutron star structure which are relevant to the paradigm. 4describes in detail the cellular automaton model and jus-ti?es physically the automaton rules.In5the 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 va-lidity and a summary of our?ndings are contained in 6. 2SELF-ORGANIZED CRITICAL SYSTEMS

2.1General 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 manip-ulation by an external agent.Critical implies that,like in phase transitions near the critical temperature,a small lo-cal perturbation can propagate throughout the entire system (Jensen1998).

For a system to be in a self-organized critical(SOC) state it must?rst satisfy the following three conditions (Jensen1998;Melatos et al.2007):

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

(2)Each element moves when the local force exceeds a threshold.In this way stress accumulates sustainedly at cer-tain 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 https://www.sodocs.net/doc/112634790.html,bined 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 proper-ties are generically observed(Jensen1998;Melatos et al. 2007):

(4)Transitions from one metastable state to the next oc-cur via avalanches:spatially connected chains of local equi-libration events,in which one element relaxes and redis-tributes some local stress to its neighbors,which in turn can exceed their thresholds and relax(knock-on e?ect).

(5)Avalanches have no preferred scale:they can involve

a few(commonly)or all(rarely)of the elements in the sys-tem.Their sizes and lifetimes follow power law distributions, whose exponents are related.Numerical values of the expo-nents 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&Jensen2002).

(6)Over the long term,the system tends to a critical state,which is stationary on average but?uctuates instan-taneously.

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

A cellular automaton model of pulsar glitches3 microscopic dynamics of a SOCS and hence is super?uous if

the system is conservative on microscopic scales.We claim

that our model is locally non-conservative and hence a driv-

ing 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 statis-

tically independent from other relaxation events in size and

when it occurs.Two or more simultaneous relaxation events

at di?erent locations cannot be distinguished by their global

e?ect(e.g.spin up of a neutron star),which is the cumu-

lative 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 con?rmed by models based on cellular

automata.

Alpar et al.(1995,1996)discussed the formation of de-

pleted 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 vor-

tex 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 lo-

cally into neighbouring regions,relaxing the capacitive re-

gion(stress reservoir).This qualitative picture bears the

hallmark of a SOCS.

2.2Pulsar glitch statistics

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?ve times.We consider here

the subset that have glitched more than six times(N g 6).

Of the eight pulsars satisfying this criterion,six are consis-

tent 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

P(δν)=(δν)1?a?(δν)1?a

min

N g X?t min exp(?λ?t min)?exp(?λ?t)

4L.Warszawski and A.

Melatos

Figure1.Cumulative distributions of glitch waiting times?t,for six pulsars with N g 6.The horizontal scaleλ?t is di?erent for each pulsar;λis the mean glitch rate that minimises the K-S statistic for that pulsar as listed in Table1.The solid curve represents the ideal,Poissonian,cumulative distribution,P(?t)=1?e?λ?t.The waiting time distribution in an ideal SOCS follows this Poissonian. Table1.Parameters of the glitch size and waiting time distributions for pulsars with N g 6(Melatos et al.2007).a andλare chosen to minimise the K-S statistic.

PSR Jλλminλmax a a min a max109(δν)min

F M=ρsκ×(v L?v s),(3)

whereρs is the super?uid density,κis a vector pointing

along the vortex,whose modulus|κ|=h/2m n is the circu-

lation around a single vortex(m n is the neutron mass),v L

is the vortex line velocity,and v s is the local bulk super?uid

velocity.If allowed to move freely,vortices organize into a

rectilinear array spaced evenly according to Feynman’s rule

(vortex area density n=4πν/κ).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～1016vortices in a neutron star,we group

vortices into bundles(in internal equilibrium)and regard

the bundles as discrete elements.The velocity?eld gener-

ated by each vortex is inversely proportional to the distance

from the vortex line,so local interactions dominate global

forces.

Condition(2)in 2.1exactly describes the pinning of

vortices to the pulsar crust,whether the pinning centres

are nuclear lattice sites,defects[monovacancies or impuri-

ties(De Blasio&Lazzari1998)],or grain boundaries.Con-

dition(2)does not rule out pinning of neutron and proton

vortices in the neutron star’s core.Vortices are unpinned by

the bulk super?uid?ow 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)calcu-

lated the monovacancy concentration needed to provide a

pinning force greater than the Magnus force,and claimed

that,forρs 3.4×1016kg m?3,the concentration is un-

physically large.

The external driver ful?lling condition(3)in 2.1is the

accumulating Magnus force generated by the build up of

rotational lag between the vortex lines and the super?uid

?ow due to the electromagnetic spin down of the star.The

A cellular automaton model of pulsar glitches5

e?ect 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; Jones1998a).Secondly,abrupt adjustments occur when lo-cally 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～1% 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 super?uid?ow,regardless of the direction of the unpinning force(Donnelly1991).Field et al.(1995)pointed out that, at least in the context of?ux creep in superconductors,such an assumption leads to‘stick-slip’dynamics that are more likely to resemble those of an ideal sandpile(the quintessen-tial SOCS),because inertial e?ects do not dominate the ef-fect 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 super?uid vor-tices are straight and parallel to the axis of rotation.Jones (1998b)argued that the high,but?nite,tension in vortex lines justi?es this assumption,as the displacements caused by the pinning centres are much smaller than the nuclear separation.Contrarily,it is the?nite tension in vortices that allows them to be drawn into the potential well of the pin-ning centres and pin at all(Jones1998b).Link&Cutler (2002)suggested pinning centres in the nuclear lattice ef-fectively span approximately thirty lattice sites,due to?-nite tension,where the internuclear spacing is～10?14m. Relative to the coarse graining scale of a practical cellu-lar automaton(linear cell dimension up to≈2×103m), this deviation from straight vortices is clearly insigni?cant. However,global hydrodynamic simulations suggest that the array breaks up into a vortex tangle above a critical ro-tational lag,driven by the Donnelly-Glaberson instability, with a net polarisation parallel to the rotation axis[see for example Peralta et al.(2005,2006a,b)].

4CELLULAR AUTOMATON

The cellular automaton model presented in this paper takes the basic physical features of the inner crust-super?uid sys-tem and interprets them as simple rules.Much of the mi-crophysics that has emerged from detailed studies of pin-ning strengths,crustal structure,and turbulent?ow(see references in1)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’(renor-malised)to the necessary coarse grain.In terms of reproduc-ing 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 compu-tational cost.Table2summarizes the nomenclature used in this and subsequent sections.4.1Grid

The simulation grid contains N2grid cells,covering a circular ‘star’of radius R,representing a projection onto the equa-torial 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 con?guration of super?uid vortices is a triangu-lar Abrikosov lattice.For the purposes of this model,we as-sume that the two con?gurations 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 in3there 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 calcu-lated using the Feynman relation

I v s·dl=κN v,(4)

where the integral is taken around a closed path(e.g.a circle of radius equal to the stellar radius),and N v is the total number of vortices enclosed by the path.In the limit of an in?nite vortex distribution,the equilibrium distribution is a rectilinear array.Although we have a?nite array of vortices, we assume that they form a rectilinear array,such that the total number of vortices is

N v(R)=4π2νR2/κ.(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(corre-sponding to one bundle).For any given pulsar we model the observed glitch range:the minimum glitch size G min(what constitutes a glitch is discussed in 4.5)occurs when one cell unpins,and the maximum G max 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 B 1bundles 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 B vortex bundles. We choose our grid size such that

N2grid=

G max

6L.Warszawski and A.Melatos

example:

Figure 2.Schematic of the cellular automaton.The left side of the ?gure 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 super?uid velocity ?eld 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 ?eld 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 ?gure describes how the force imbalance due to the nearest-neighbour vortices,F NN ,is calculated.Each vortex bundle generates a velocity ?eld according to a right-hand rule ,with the circulation vector pointing out of the page,which falls o?inversely with distance from the centre of the bundle.

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 oc-cur at all scales,within the bounds of the system.To resolve smaller (δν/ν)min in the model,we must fractionate into smaller bundles.In order to make comparisons with the ob-servational data,we preclude a priori these smaller glitches from occurring in our model.If the system is truly scale in-variant,this e?ective window function (coarse-grain scale)should not interfere with its overall statistical behaviour.In the ?ux 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 his-tory dictates.Bassler &Paczuski (1998)take a similar ap-proach,allowing many point pins in an extended cell.4.2

Initial conditions

Initially the vortex bundles are laid out at random 3such that the mean occupancy of a cell is B bundles.Whether or not it is realistic to have a vortex distribution that is inhomogoneous on the scale of the grid cells (～2×103m)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 inho-mogeneous con?guration of pinned vortices,and take this to

3

Vortex bundles are placed at a randomly selected cell,one at a time,until all N bun bundles have been distributed.

be the initial state of the automaton.We do not model di-rectly how this inhomogeneous state arises from an initially homogeneous state because the latter con?guration is inac-cessible 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 de?ne a critical state of the system,we allow the simulation to complete many (～104)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 stel-lar crust,F thresh ,plays an essential role in the mass unpin-ning model of pulsar glitches.Recent calculations imply that pinning to nuclear lattice sites is strongest in regions of low (～4.9×1015kg m ?3)and high density (～1.6×1017kg m ?3),with a maximum energy of E p ～3MeV (Avogadro et al.2007).We take the pinning force to be the same everywhere on the grid,with value E p /ξ,where ξ～10?14m is the su-per?uid coherence length (Elgar?y &De Blasio 2001);i.e.F thresh ～5×1015J m ?1.For the purposes of the model we convert this to a threshold on the velocity lag,?v pin =(v L ?v s )pin ～104m s ?1.

Vortex pinning at lattice defects is also possible.Link &Cutler (2002)argued that if the pinned (coupled)?uid makes up 1%of the total moment of inertia,then pin-ning must occur throughout the crust,suggesting that the crust should be treated as an amorphous solid with ran-dom pinning sites.The randomness ensures that the con-tributions to the force of attraction towards the nuclear sites on either side of the vortex line do not cancel (Jones

A cellular automaton model of pulsar glitches7

1998a).Another possibility is to adopt the random pin-ning potential used to model?ux creep in superconductors (Bassler&Paczuski1998).

Vortices unpin when the relative velocity of the super-?uid and the vortex lines is large enough to produce a Mag-nus force F M greater than the pinning threshold F thresh.If the Magnus force does not exceed F thresh,the pinning force adjusts to balance the Magnus force exactly.Given that the vortices are assumed to be inertialess,the excess force be-yond the pinning threshold is of no relevance,as once the vortices unpin they?ow with the ambient super?uid:

F pin= F M if F M

0otherwise.

(8)

In our model,the super?uid 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-?eld approximation);and the contribution from the eight nearest-neighbour cells on the grid,as depicted in Fig. 2:

v s=v s,pinned+v s,unpinned+v s,NN.(9) The?rst two contributions,hereafter named the mean-?eld terms,are calculated using Eq.(4),assuming that the ve-locity?eld is exactly tangential to the closed path(a good approximation which becomes exact in the thermodynamic limit N v→∞).

By assuming that the vortices are arranged in a regular array,we can evaluate(4)with tangential v s,without know-ing the location of each vortex(the validity of this approach is veri?ed in 5.1).The Magnus force is then calculated from Eq.(3),with v L equal to the veloicty of the stellar crust and v s,tot given by Eq.(9).Initially,the mean Magnus force is very small.As the star spins down,however,the pinned con-tribution does not diminish,while the unpinned contribution does,and the mean?v 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 vor-tices is much less than the number of pinning sites.

4.4Driver

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 vor-tices and the unpinned super?uid.Over short time-scales,we include thermally activated unpinning of the vortex bundles in random cells.4In both scenarios(short and long time-scales),small variations in?v=v L?v s due to the eight

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.nearest-neighbour cells trigger avalanches in the system.By treating the nearest-neighbour interactions with greater pre-cision than the longer range contributions,we are empha-sizing the property outlined in condition(1)of 2.

4.5Operational de?nition of glitches

In the automaton model a glitch corresponds to the unpin-ning of vortex bundles.Its size is determined by the total number of bundles(and hence vortices)unpinned.The as-sociated fractional spin up of the pulsar is then given by the fraction of the total moment of inertia carried by the recently unpinned?uid.For example,if each vortex bun-dle represents the unpinning of?uid carrying10?7of the moment of inertia of the system,then an avalanche of ten bundles results in a glitch of sizeδν/ν=10?6.

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 super?uid 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 rela-tion(4).The super?uid velocity at some distance from the rotation axis is proportional to the number of vortices en-closed,so the further the unpinned vortices travel radially, the more the super?uid decelerates,and hence the more an-gular momentum is transferred to the crust.The de?nition of glitch size given in this paper is therefore an approxi-mation,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 satis?ed by our uniform automaton,but not necessarily realistic).We show in 5.4that these conditions are ful?lled 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～40s(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 simul-taneous 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 con-form to all facets of the above de?nition.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 Crab20to40days after a standard glitch(Wong et al.2001).

8L.Warszawski and A.Melatos

4.6Automaton rules

The rules of our cellular automaton encompass the key fea-tures of SOC outlined in2: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.1and 4.2.

(2)The force imbalance F M?F thresh is calculated at each cell based on the scheme outlined in 4.3.Cells where the force imbalance is positive are?agged as supercritical. (3)Half the vortex bundles in a supercritical cell are repo-sitioned on the adjacent cell that lies in the path of the bulk super?uid 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(1)–(4)are repeated until no cells are supercrit-ical.

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

(7)(a)Either the vortex bundles occupying a ran-domly selected cell are redistributed into the neighbouring cells;or

(b)the spin frequency of the crust(and hence v L/(2πR))is decremented by˙ν?t ts,the number of un-pinned vortices is proprotionally reduced,and the number of pinned vortices remains unchanged.

(8)Steps(1)–(7)are repeated for a predetermined num-ber of time steps,corresponding to either an arbitrary period of time[if(7)(a)is used]or the number of time steps mul-tiplied by?t ts(where?t ts is the length of each time step), if(7)(b)is used.

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 super?uid?ow.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 vor-tex creep(P?t ts?ν/˙ν),or electromagnetic spin down (P?t ts ν/˙ν)(see 4.3).

4.7Vortex motion

The rules for moving vortices to an adjacent cell following an unpinning event rely on the(incorrect)assumption that|v s| 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 justi?ed in two ways in the model.First,as long as the time Table2.Default physical and computational parameters used to analyse the model.

Parameter Symbol Default value 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 activ-ity takes place in a thin annulus,this approach is accurate, because|v s|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 destina-tion cell)is greater than the travel time for vortices start-ing near the border.Crudely,if vortices are evenly spaced throughout the departure cell,about half of them have suf-?cient time to reach the destination cell during a small time step.As a corollary,the model contains a dissipative mech-anism,which ensures that avalanches eventually terminate. This feature means that our system is non-conservative.Fol-lowing Pruessner&Jensen(2002),we require an external driving force to maintain the system in a critical state.

5STATISTICS

5.1General features

Three quantities in particular are used to characterise our cellular automaton model of pulsar glitches:the distribution of glitch sizesδν/ν,durations t,and waiting times?t.The glitch duration is de?ned 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 al-lowed to run for105big time steps.In 5.2–5.6we consider time periods over whichνdoes not decrease signi?cantly.In 5.7we consider longer time periods over which spin down becomes important.We characterise the model by explor-ing its response to changes in the controlling parameters, both physical and computational.The physical parameters are the spin frequencyν,the pinning threshold?v pin,and the fraction of vortices that are pinned?.The compuational parameters are the numbers of bundles N bun≡BN2grid and grid cells N2grid.

Fig.3shows the raw output of the model in the form of a time series of glitch sizes for a standard set of parame-ters(de?ned Table2).Unless otherwise speci?ed,these are the parameters used to generate the model output.This particular set of parameters is chosen to balance the need for su?cient avalanches to generate good statistics versus

A cellular automaton model of pulsar glitches9 computational e?ciency.The vertical axis gives the size in

terms of the fraction of the total super?uid moment of in-

ertia that unpins during the avalanche.The inset zooms in

on a segment of the time series spanning103big time steps.

Although the simulation runs for a total of105time steps

(without spin down),the?rst103time 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 plot105or103big time steps.

The lower panel of Fig.3shows the probability den-

sity functions of glitch sizesδν/ν,durations t,and wait-

ing times?t,plotted as histograms,for the standard set

of parameters de?ned in the caption.Bothδν/νand t

obey power law statistics,with probability density functions

p(δν/ν)∝(δν/ν)?a and p(t)∝t?b respectively.This is a

characteristic feature of a SOCS(see 2.1).For the standard

parameters in Fig.3,we obtain a=2.9and b=4.8.More

generally,for the parameter ranges studied in this paper,we

?nd2.0 a 4.3and2.2 b 5.5(see 5.2–5.6).In a

SOCS,a and b are related by(Jensen1998)

γ2

b=1+(a?1)

10L.Warszawski and A.

Melatos

Figure 3.Top :Automaton output for 9×104big time steps.The radius of the star is R =10km,the unpinning threshold is |?v pin |=1.0×104ms ?1,the neutron super?uid density is ρs =1016kg m ?3,the initial spin frequency is ν0=30.25Hz,and the current spin frequency is ν=13.9Hz.The avalanche size is given as the fraction of the total ?uid moment of inertia that unpins during the avalanche.The fraction of pinned vortices is taken to be ?=0.01.The main panel displays a time series of avalanche sizes beginning after 104time 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 5×103time steps.Bottom :Probability density funtions of δν/ν(left ),durations t (centre )and ?t (right )plotted as histograms on a log-log (left and centre )and log-linear (right )scale,using 102logarithmically (left and centre )and linearly (right )spaced bins.The same parameters are used as for the top panel.The size and duration distributions are ?tted with power laws of the form p (δν/ν)=(δν/ν)?a and p (δν/ν)=(δν/ν)?b respectively.The waiting time distribution is ?tted with an exponential of the form p (?t )=e ?λ?t ,where λis in units of the reciprocal of one big time step.All ?tted functions are plotted as dashed lines.

5.2

Dynamic range

In 4.1we motivate the choice of grid size by referring to the observed dynamic range of pulsar glitches.In any indi-vidual pulsar the maximum dynamic range of δν/νis 104(for PSR J1704–3015),and the minimum is 101(for PSR J0537–69105).Hence,from (6),we run our simulation on grids ranging from N grid =10to N grid =100(the physi-cal size remains unchanged at R =20km).By construction,the range of avalanches outputted by the model should not exceed G max /G min (i.e.102for the 10×10grid,and 104for

5

This is one of two pulsars that glitch quasiperiodically;the other is Vela.

N grid =100grid).Fig.7con?rms this expectation.The top panels are histograms of the avalanches produced by sys-tems with N grid =100,N grid =60and N grid =10over 105time steps.The dynamic range is greatest for the largest grid (10?6.1 δν/ν 10?4.2)and smallest for the smallest grid (10?4.2 δν/ν 10?3.2).It is not expected that these dynamic ranges span the full range G min δν/ν G max ,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.7is overlaid with a power law best ?t of the form p (δν/ν)∝(δν/ν)?a ,estimated us-ing a least-squares algorithm,with the ?tted exponent in-dicated on each plot.By the central limit theorem,the un-

A cellular automaton model of pulsar glitches11

Figure5.Fitted power law exponents for the distribution of glitch sizes(?)and durations( ).The error bars represent the1σuncertainty returned by the least-squares?tting algorithm.Top left:10 N grid 100.Top right:9.5×104 |?v pin| 10.4×104m s?1.Bottom left:0.0092 ? 0.010.Bottom right:12.73 ν 15.91Hz.Automaton parameters:N grid=100,B=1,105big time steps.Physical parameters:ρs=1016kg m?3,ν0=30.25Hz,R=20km.

Figure6.Left:A vector plot of the velocity imbalance?v=v L?v s,generated using the mean-?eld approximation described in Eq.(9).The inset zooms in on the central10×10grid cells.The solid black circle shows the position of the critical radius,where the mean-?eld contribution to?v is approximately equal to?v pin.Right:Pro?les of?v along a straight-line cut through the grid,as a function of distance from the origin,as calculated by the mean-?eld(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?v pin.

12L.Warszawski and A.

Melatos

Figure7.Probability density functions p(δν/ν)of avalanche sizesδν/ν,plotted as histograms on a log-log scale,using102logarithmically spaced bins,for105big time steps.The histograms are?tted with a power law of the form p(δν/ν)∝(δν/ν)?a(dashed lines).Top: Grids with N grid=100(left),60(centre)and10(right).All panels have bundle factor B=1.Middle:Bundle factors B=1(left), 103(centre)and105(right).All panels have N grid=100.Bottom:B=105and N grid=100(left),60(centre)and20(right).All panels have B=105.The vertical grey line denotes the crossover scale.Physical parameters:|?v pin|=104m s?1,ρs=1016kg m?3,ν0=30.25Hz,ν=13.9Hz,?=0.01,R=20km.

certainty in each bin is inversely proportional to the square

root of its contents.The?tted exponent is graphed as a

function of N grid in Fig.5,with error bars taken from the

least-squares algorithm.There is a trend in both a and b to

increase with increasing N grid,in keeping with the increas-

ing dynamic range with N grid.When the simulation is run

for2×105time steps,the dynamic range does,as expected,

increase;however,it still falls well short of G max/G min.The

number of bundles available to unpin is restricted to a small

fraction(～15%)of BN2grid that lie in the vicinity of the

critical circle.To better match the desired dynamic range,

we should select the grid size such that G max/G min cells lie

in the active zone near the critical circle.

The glitch sizes generated by the model should be

A cellular automaton model of pulsar glitches 13

Figure 8.Probability density functions p (t )of avalanche durations t ,for bundle factors B =1(left ),103(centre )and 105(right ),plotted as histograms on a log-log scale,using 102logarithmically spaced bins.The histograms are ?tted with a power law of the form p (t )∝t ?b (dashed lines ).Automaton parameters:N grid =100,105big time steps.Physical parameters:|?v pin |=104m s ?1,ρs =1016kg m ?3,ν0=30.25Hz,ν=13.9Hz,?=0.01,R =20km.

Figure 9.Probability density functions p (δν/ν)of avalanche sizes δν/ν,plotted as histograms on a log-log scale,using 102logarithmically spaced bins.The histograms are ?tted with a power law of the form p (δν/ν)∝(δν/ν)?a (dashed lines).Top :Current spin frequency ν=15.9(left ),14.7(middle )and 13.1Hz (right ).Bottom :Pinning threshold |?v pin |=9.5×103(left ),9.9×103(middle )and 10.4×103m s ?1(right ).Automaton parameters:N grid =100,105big time steps,B =100.Physical parameters:ρs =1016kg m ?3,ν0=30.25Hz,?=0.01,R =20km.

14L.Warszawski and A.

Melatos

Figure 10.Probability density functions p (?t )of waiting times ?t ,plotted as histograms on a log-linear scale,using 102linearly spaced bins,for 105big time steps.The histograms are ?tted with an exponential of the form p (?t )∝e ?λ?t (dashed lines).Top :Grid size N grid =100(left ),60(centre )and 10(right ).Bundle factor B =1.Bottom :Bundle factor B =1(left ),5×102(centre )and 105(right ).Grid size N grid =100.Physical parameters:|?v pin |=104m s ?1,ρs =1016kg m ?3,ν0=30.25Hz,ν=13.9Hz,?=0.01,R =20km.

treated as a indicative only.Evidently,several of the pul-sars analyzed by Melatos et al.(2007)exhibit glitches sev-eral orders of magnitude smaller than those produced by the automaton.It is easy to elicit such small glitches from the automaton by adjusting F pin or ν?ν0so that the ac-tive 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.5for 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,de-pending 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 consis-tent with unbroken power laws.

5.3Vortex 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 B .

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 e?ect of increasing the number of vortex bundles (equivalent to decreasing the number of vortices per bundle,for ?xed ν)on the distribution of avalanches.

The avalanche size distributions for di?erent values of the bundle factor B are shown in the middle panels of Fig.7.The histograms are overlaid with a power law best ?t.For B =5×102and B =1×105there is an obvious ‘bump’in the glitch size distribution,whose peak occurs at the same value of δν/νfor di?erent B [log 10(δν/ν)bump =?4].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 de-cay.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 in-

A cellular automaton model of pulsar glitches15

creasing function of the system size(Jensen1998).In the examples discussed by Jensen(1998),the system size corre-sponds 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 B and becomes increasingly prominent for large B suggests a system size e?ect.On the other hand,the crossover scale,de?ned by the bump,does not increase with system size(i.e.B)in Fig.7.Nor does system size explain why the bump is a local maximum in p(δν/ν)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?eld 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 B=1×105) as in Fig.7.Despite the bump not being present for smaller systems(N grid 40),there is a value ofδν/ν,for each N grid, above which p(δν/ν)falls o?faster than a power law,which is the expected behaviour beyond the crossover scale de?ned by Jensen(1998)(see Fig.7).Bumps like the ones in Fig. 7also arise in SOCS when one changes the morphology of the grains(e.g.rice versus sand)(Frette et al.1996;Jensen 1998;Pruessner&Jensen2003).

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 B=1(one bun-dle per cell,on average).By increasing B without changing N grid,we are renormalising the vortices on a smaller coarse grain,whilst maintaining the coarseness of the grid,which is not a justi?able physical scenario.Pruessner&Jensen (2002)claimed that inappropriately chosen parameters pro-duce a bump in the size distribution,which further suggests that B?1alters 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 pre-ferred 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 bun-dles 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 signi?cantly outnumbered by the vortex bundles,the scale invariance of the system is broken.In the absence of computational lim-itations,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 cir-cumvents the subtleties of vortex-vortex interactions and reconnections by renormalizing these e?ects into e?ective vortex-vortex forces on a coarsened grid.

5.4Pinning threshold

The pinning threshold F pin determines the minimum lag be-tween the super?uid and pinned vortices necessary to unpin vortex bundles.The pinning threshold can be reexpressed in terms of a minimum|?v pin|between the pinned vortices (comoving with the crust)and the unpinned super?uid.

As|?v pin|increases,the number of avalanches occur-ring within a?xed time(105big time steps)decreases,from 33549for|?v pin|=9.5×103m s?1to2227for|?v pin|= 10.4×103m s?1,as shown in Fig.9.There are two rea-sons 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 im-balance,the less probable it becomes.Second is the issue of?ne tuning.For a given|?v pin|there is a critical circle where the mean-?eld 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?nd a mixture of sub-and supercritical cells because the nearest neighbour imbalance can act to reduce|?v|at a given cell,if it opposes the mean-?eld contribution.As|?v pin|increases,the critical circle moves towards the origin,and hence the number of cells in its vicinity decreases.Fig.11depicts images of the grid after105big time steps.The colours encode the value of |?v|on each cell.Supercritical cells are white;all other cells are subcritical.The pinning threshold in the upper panel of Fig.11is|?v pin|=104m s?1.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 ?v pin=9×103m s?1,we observe an annulus in which most cells are supercritical(where the vortices never pin).Just in-side 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 un-pinning activity is restricted to a narrow annulus,and hence una?ected 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|?v pin|in Fig.9.The best?t power law and its exponent are shown in each panel of Fig.9.We?nd that a decreases as|?v pin|increases,however there is a slight upturn in both a and b for the largest values of|?v pin|con-sidered.For larger values of|?v pin|,avalanches involving many cells(and hence many vortex bundles)are less likely to occur,both because the active annulus shrinks,and be-cause it is relatively unlikely that the nearest-neighbour im-balance sends a cell supercritical.Hence a decreases because the dynamic range decreases;the area under the histogram equals the number of avalanches.

Fig.6demonstrates that,for a?nite,square vortex ar-ray,the pro?le of?v along a linear cut?attens out in the region where r≈R,compared to smaller radii.Hence if |?v pin|approaches the maximum?v 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 a versus|?v pin|seen in Fig.5.

5.5Pinned 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 in-creases the minimum avalanche size,whilst maintaining the maximum avalanche https://www.sodocs.net/doc/112634790.html,rger bundles also enhance the role of the nearest neighbour cells.As more vortices pin and unpin,the?rst and third terms in Eq.(9)decrease with re-

16L.Warszawski and A.

Melatos

Figure11.Images of the velocity imbalance|?v|for|?v pin|=104m s?1(upper)and9×103m s?1(lower).White cells are supercritical; the darker the cell,the smaller the?v.

spect to the contribution from the unpinned vortices.Given that the unpinned vortices are homogeneously distributed, they are the foremost contributors to the mean-?eld critical circle discussed in 5.4.Decreasing their relative contribu-tion reduces the mean-?eld e?ect,emphasizing the nearest-neighbour,stochastic dynamics of the model.The lower left panel of Fig.5supports this argument:the1σerror in the ?tted a and b values decreases with increasing?.The lower left panel of Fig.5also shows that both a and b increase with increasing?.For? 0.0097the a distribution?attens out,whereas the b distribution continues to increase.

For? 0.009no avalanches are observed;the model is badly tuned.Moreover,for? 0.010,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),in-stead calculating?v unneccessarily in these regions at each small time step.We intend to address this?aw in future work.

5.6Spin frequency

Two values of the spin frequency enter the model:at birth (ν0),and today(ν).As a?rst approximation,the total num-ber of pinned vortices is taken to be unchanged over the lifetime of the star and hence is dictated byν0as in Eq.(5). The pinned portion of the super?uid does not spin down with the stellar crust,whereas the unpinned super?uid does; its circulation is dictated byν.In the current model we holdνandν0constant(cf. 5.7),withν0=30.25Hz,and 12.73 ν 15.91Hz.

By decreasingν,we increase|?v|and hence the force imbalance at every cell.The impact is similar to decreasing |?v pin|because the critical circle shrinks.The two right-hand panels of Fig.5display a similar trend for increasing |?v pin|andν.The main di?erence is thatνalso determines the relative size of the?rst two terms in Eq.(9),given that the number of pinned vortices remains constant while the number of unpinned vortices decreases withν.These e?ects 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|?v pin|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 perma-nently 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 bi-ased towards large avalanches.

5.7Long-term activity

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 di?er by several orders of magnitude;in this paper we take˙ν=?7×10?13Hz s?1by way of example.Observations of pulsar glitches span a pe-riod of approximately forty years(～1970to present),which equates to a total spin down of～0.1Hz.In our model,each big time step must therefore span～106s,in order to rep-resent a total observing time of forty years.Such a large time step determines the scale on which we can resolve indi-vidual avalanches.Observationally,glitches can be resolved down to a40s 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 ?xed,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.12shows the size distribution (δν/ν)for time steps?t ts=10,105and5×105s.Fig.13 shows the?tted power law exponents for the sizes and dura-tions,for10 ?t ts 5×105.Within the uncertainty of the ?tting algorithm,a and b are the same for all values of?t ts

A cellular automaton model of pulsar glitches

17 Figure12.Probability density functions p(δν/ν)of avalanche sizesδν/ν,for time steps lasting?t ts=10(left),105(centre)and 5×105(right),plotted as histograms on a log-log scale,using102logarithmically spaced bins.The histograms are?tted with a power law of the form p(δν/ν)∝(δν/ν)?a(dashed lines).Automaton parameters:N grid=100,105big time steps.Physical parameters: |?v pin|=104m s?1,ρs=1016kg m?3,ν0=30.25Hz,ν=13.9Hz,?=0.01,R=20

km.

Figure13.Fitted power law exponents for the distribution of

glitch sizes(?)and durations( ),for time steps in the range

10 ?t ts 5×105.The error bars represent the1σuncer-

tainty returned by the least-squares?tting algorithm.Automa-

ton parameters:N grid=100,105big time steps.Physical pa-

rameters:|?v pin|=104m s?1,ρs=1016kg m?3,ν0=30.25Hz,

ν=13.9Hz,?=0.01,R=20km.

tested.For the largest?t ts,the time elapsed after105big

time steps is5×105s,which equates to?ν=?0.007Hz.

This moves the critical circle towards the origin by a negli-

gible～0.007%.

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 inde?nitely,as the critical circle moves inwards.As

time passes,however,the glitches become smaller,on av-

erage,as the number of cells in the vicinity of the critical

circle diminishes6.If,as claimed by Donati&Pizzochero

(2003),pinning occurs at intermediate densities,no glitches

6Interestingly,there is no guarantee that glitches emanating

from di?erent critical circles are,in fact,drawn from the same

occur once the critical circle shrinks inside the radius at

which these densities occur.Alternatively,if pinning oc-

curs at several strata[high and low densities,for example

(Avogadro et al.2007)],glitches occur when the critical cir-

cle 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 super?uid spins down when vortices

move radially outward.The outward motion is driven by the

distribution of capacitive and depleted zones;outward mo-

tion 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 re-

ducing 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

super?uid?ow when computing v L.(This does not remove

the arti?cial radial motion arising from the rectangular grid.

In this context,an annular grid topology is preferable.)

6DISCUSSION

In this paper we present a cellular automaton model for

glitches in neutron stars,based on the vortex unpinning

paradigm(Anderson&Itoh1975;Alpar et al.1984).De-

spite 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

statistical distribution(e.g.value of a).If so,we do not expect a

pure power law inδν/νover long times.

18L.Warszawski and A.Melatos

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 vor-tex avalanche is a statistically independent event.The values of the power law exponents and Poisson glitching rates de-pend on the particular parameter set.Our?ts return size and duration exponents in the ranges2.0 a 4.3and 2.2 b 5.5respectively,and mean glitching rates in the range0.0023 λ 0.13.Although we do not have a?rst-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 dis-tributed.Alpar et al.(1996)suggested that crust cracking creates capacitive regions where many vortices pin simulta-neously.These regions do not permit a vortex current,so pinned vortices in these regions do not participate in out-ward vortex creep.Crust cracking is a nonuniform process, suggesting that the regions themselves are inhomogeneously distributed.Our model conforms qualitatively with this pic-ture.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 veri?ed,this claim suggests that the level of coarse graining necessary to render our cellu-lar 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 vor-tices to occupy a contiguous sequence of pinning positions and preventing the type of inhomogeneity represented by our initial conditions(Jones1998a).Importantly,such a scenario challenges any model based on collective unpin-ning of numerous vortices,not just our cellular automa-ton.De Blasio&Lazzari(1998)concluded that the crust is relatively free of microcrystalline structures and vacan-cies,with impurities,or point defects,the most likely dis-ruptions 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 macro-scopic 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 stud-ies have shown that pinning forces su?cient to oppose the Magnus force develop only in speci?c 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 di?erent computational topology,e.g.an annulus,or sev-eral 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?nd that a narrow active annulus, situated at the critical radius,participates in the avalanches.However,it is important to experiment with other topolo-gies in future work,because the active annulus we?nd 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&Sedrakian1995; 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 super?uid 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|?v|>|?v pin|which are always supercritical(outside the critical circle)and do not participate in avalanche dy-namics as they never pin.If pinning is restricted to one or more annuli,the glitch behaviour of the pulsar changes sig-ni?cantly when the critical circle moves into and out of the pinning zones.

Clearly,our cellular automaton does not include the re-sponse of the pulsar to glitches(ie.the spin up of the pulsar crust and the subsequent‘relaxation’).Nor do we account for the?nite time-scale on which the unpinned super?uid cou-ples 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?ne tun-ing in both the physical and computational parameters,such that|?v max|≈|?v pin|.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?v pin 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 B. For B?1,the automaton output is no longer consistent with a SOCS.

In conclusion,we present an empirical cellular automa-ton model of pulsar glitches based on the mass vortex unpin-ning paradigm.We?nd that for certain physical and compu-tational parameters the model produces dynamics that are consistent with a SOCS and with radio timing data from pulsars.There exists no general,?rst-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 ex-ponents,and the mean rate of the waiting-time distribution. We do demonstrate that the basic physical principles gov-erning inter-vortex interactions can produce the type of col-lective behaviour necessary to explain pulsar glitches within

A cellular automaton model of pulsar glitches19

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 statisti-cal methods.L W acknowledges the support of an Australian Postgraduate Award.

REFERENCES

Alpar M.A.,Anderson P.W.,Pines D.,Shaham J.,1984, ApJ,276,325

Alpar M.A.,Chau H.F.,Cheng K.S.,Pines D.,1996, ApJ,459,706

Alpar M.A.,Cheng K.S.,Pines D.,1989,ApJ,346,823 Alpar M.A.,Kiziloglu U.,van Paradijs J.,1995,The lives of the neutron stars.Dordrecht;Boston:Kluwer Aca-demic,c1995.

Anderson P.W.,Itoh N.,1975,Nature,256,25 Andersson N.,Comer G.L.,Prix R.,2004,MNRAS,354, 101

Avogadro P.,Barranco F.,Broglia R.A.,Vigezzi E.,2007, Phys.Rev.C,75

Bak P.,1996,How nature works:the science of self-organized criticality.Copernicus,New York,NY,USA Bak P.,Tang C.,Wiesenfeld K.,1988,Phys.Rev.A,38, 364

Bassler K.E.,Paczuski M.,1998,Phys.Rev.Lett.,81,3761 De Blasio F.V.,Elgar?y O.,1999,Phys.Rev.Lett.,82, 1815

De Blasio F.V.,Lazzari G.,1998,Nucl.Phys.A,633,391 Dodson R.G.,McCulloch P.M.,Lewis D.R.,2002,ApJ, 564,L85

Donati P.,Pizzochero P.M.,2003,Phys.Rev.Lett.,90 Donati P.,Pizzochero P.M.,2006,Phys.Rev.B,640,74 Donnelly R.J.,1991,Quantized vortices in Helium II.Press Syndicate of the University of Cambridge,Cambridge, New York NY,USA

Elgar?y O.,De Blasio F.V.,2001,A&A,370,939

Field S.,Witt J.,Nori F.,1995,Phys.Rev.Lett.,74,1206 Frette V.,Christensen K.,Malthe-S?rensen A.,Feder J., J?ssang T.,Meakin P.,1996,Nature,379,49

Jensen H.J.,1990,Phys.Rev.Lett.,64,3103

Jensen H.J.,1998,Self-organized criticality:Emergent complex behavior in physical and bilogical systems.Cam-bridge University Press,Cambridge UK

Jones P.B.,1997,Phys.Rev.Lett.,79,792

Jones P.B.,1998a,Phys.Rev.Lett.,81,4560

Jones P.B.,1998b,MNRAS,296,217

Linder J.F.,Hughes S.B.,Miller D.J.,Thomas B.C., Weisenfeld K.,2004,Int.J.Bifurcation and Chaos,14, 1155

Link B.,Cutler C.,2002,MNRAS,336,211

Link B.,Epstein R.I.,1996,ApJ,457,844

Lu E.T.,Hamilton R.J.,1991,ApJ,380,L89

Mastrano A.,Melatos A.,2005,MNRAS,361,927 McCulloch P.M.,Hamilton P.A.,McConnell D.,King E.A.,1990,Nature,346,822

Melatos A.,Peralta C.,2007,ApJ,662,L99

Melatos A.,Peralta C.,Wyithe J.S.B.,2007,MNRAS Middleditch J.,Marshall F.E.,Wang Q.D.,Gotthelf E.V., Zhang W.,2006,ApJ,652,1531

Morley P.D.,1996,A&AS,313,204

Morley P.D.,Schmidt I.,1996,Europhys.Lett.,33,105 Nicodemi M.,Jensen H.J.,2001a,J.Phys.A:Math.and Gen.,34,8425

Nicodemi M.,Jensen H.J.,2001b,Phys.Rev.Lett.,86 Peralta C.,Melatos A.,Giacobello M.,Ooi A.,2005,ApJ, 635,1224

Peralta C.,Melatos A.,Giacobello M.,Ooi A.,2006a,ApJ, 635,L53

Peralta C.,Melatos A.,Giacobello M.,Ooi A.,2006b,ApJ, 651,1079

Pines D.,Alpar M.A.,1985,Nature,316,27

Pruessner G.,Jensen H.J.,2002,Europhys.Lett.,58,250 Pruessner G.,Jensen H.J.,2003,Phys.Rev.Lett.,91 Ruderman M.,Zhu T.,Chen K.,1998,ApJ,492,267 Ruderman R.,1969,Nature,223,597

Sedrakian A.D.,Sedrakian D.M.,1995,ApJ,447,305 Sornette D.,Sornette A.,Vanneste C.,1991,Self-organized criticality,earthquakes and plate tectonics.Vol.392, Berlin Springer Verlag,Berlin/Heidelberg

Vespignani A.,Zapperi S.,1998,Phys.Rev.E,57,6345 Wang N.,Manchester R.N.,Pace R.T.,Bailes M.,Kaspi V.M.,Stappers B.W.,Lyne A.G.,2000,MNRAS,317, 843

Wiesenfeld K.,Tang C.,Bak P.,1989,J.Stat.Phys.,54, 1441

Wong T.,Backer D.C.,Lyne A.G.,2001,ApJ,548,447