Gross–Pitaevskii model of pulsar glitches Free

Abstract

1 INTRODUCTION

Glitches are discrete, randomly timed jumps in the famously predictable angular velocity of a pulsar. The statistical distributions of glitches in individual pulsars (power-law sizes and exponential waiting times) point to an underlying threshold-triggered collective process (Melatos, Peralta & Wyithe 2008), akin to grain avalanches in sand piles (Wiesenfeld, Tang & Bak 1989; Morley & Schmidt 1996) or vortex avalanches in type II superconductors (Field et al. 1995). In pulsars, the ‘grains’ are quantized superfluid vortices. A comprehensive theory of why ≳10⁷ (out of ∼10¹⁸) vortices simultaneously unpin and self-reorganize to create a glitch remains elusive. A complete solution relies on connecting the microscopic (fm-scale) and macroscopic (km-scale) dynamics of the neutron superfluid in the pulsar interior. The challenge presented by this disparity of scales, especially in replicating collective phenomena, has proved formidable.

In this paper, we present the first comprehensive quantum mechanical simulations of pulsar glitches, using the Gross–Pitaevskii equation (GPE) to model the temporal evolution of the superfluid inside a pulsar. The GPE is a standard tool in condensed matter physics for investigating vortex dynamics in quantum condensates (Pethick & Smith 2002). Our simulations are unique in four ways: (1) the system is much larger than those routinely simulated in condensed matter applications (e.g. Tsubota, Kasamatsu & Ueda 2002; Penckwitt 2003; Tsubota, Kasamatsu & Kobayashi 2010), and therefore contains many more vortices; (2) in addition to the confining potential, we include a grid of pinning sites to model nuclear lattice defects in the pulsar’s crust (Jones 1998b; Donati & Pizzochero 2004, 2006; Avogadro et al. 2007, 2008; Barranco, Broglia & Vigezzi 2010); (3) the container and superfluid are coupled via a feedback torque and (4) we allow the angular velocity of the container to change non-adiabatically. Thus, the simulations describe scales ranging from individual vortex cores to the collective behaviour of several hundred vortices in the presence of a large grid of pinning sites. The results build on previous GPE studies of the general problem of spasmodic superfluid spin-down (Warszawski & Melatos 2010a) and the microphysics of individual vortex unpinning, including knock-on processes involving acoustic radiation and the proximity effect (Warszawski & Melatos 2010b).

The paper is structured as follows. Section 2 summarizes the vortex unpinning paradigm of pulsar glitches, the main features of the numerical model and the feedback mechanism. In Section 3 we describe our algorithm for extracting glitches from simulated spin-down curves. In Section 4 this algorithm is employed in a canonical example. In Section 5 we construct glitch size and waiting-time distributions from relatively small-scale simulations (∼30 vortices) for different pinning parameters. In Section 6, we vary the inertia of the star and the electromagnetic spin-down torque. Section 7 reports results from larger scale simulations involving ∼200 vortices, in which we look for evidence of collective dynamics. Parameter studies with the larger scale simulation are too expensive computationally at this stage. In Section 8 we synthesize the results and present our conclusions.

2 VORTEX UNPINNING PARADIGM

2.1 Overview

The superfluid in a pulsar is composed of neutron Cooper pairs, which interpenetrate the nuclear lattice in the crust (Migdal 1959). The superfluid attempts to match the rotation of the crust by forming quantized vortices, each carrying a quantum of circulation κ=h/(2m), whose areal density is proportional to the angular velocity of the superfluid. In the absence of pinning, vortices form a hexagonal Abrikosov lattice, causing the superfluid to rotate as a rigid body. In equilibrium, the mean areal density of vortices is n_F = 2Ω/κ. In reality, however, vortices pin at or between nuclear lattice sites and/or defects, thus distorting the equilibrium vortex configuration. The extent of distortion depends on the configuration and abundance of pinning centres supplied by the crustal lattice of nuclei. Pinning prevents the superfluid from decelerating smoothly with the crust. Glitches are believed to occur when a large number of vortices simultaneously unpin and move outward to new pinned positions. Indeed, GPE simulations reported in Warszawski & Melatos (2010a) exhibit stick-slip vortex dynamics and resulting glitch-like interruptions to the smooth spin-down of the rotating container. Similar spin-up events were observed in laboratory experiments using rotating containers filled with helium II (Tsakadze & Tsakadze 1980).

We construct a simple, two-dimensional pulsar model comprising a superfluid inside a rotating, decelerating crust, in the presence of a pinning grid that corotates with the crust. We include self-consistent coupling between the superfluid and crust by ensuring that angular momentum is conserved overall.

A popular explanation for pulsar glitches is that angular momentum is transferred from the superfluid to the crust in response to sudden unpinning, outward motion and repinning of vortices (Packard 1972; Anderson & Itoh 1975). Crust-superfluid differential rotation produces a fictitious Magnus force, F_M, that unpins vortices when it exceeds some threshold. F_M is transverse to the relative velocity between the vortex line and superfluid flow; it is similar in origin to the Bernoulli lift force in conventional hydrodynamics. For a vortex line moving with velocity v through a superfluid of density ρ and velocity v_s, the Magnus force is given by

1

where κ is directed along the vortex axis. For a pinned vortex, v is the crust velocity; for a vortex at r=b, one has ⁠. For N_v vortices in an Abrikosov lattice, the Feynman relation approximates (the result is exact in the infinite-vortex limit) the azimuthal component of the superfluid velocity at radius r as v_ϕ≈κN_v/(2πr).

The vortex unpinning paradigm presumes that a catalyst exists to spontaneously unpin between 10⁷ and 10¹⁵ vortices during a single glitch, with the range deduced from measurements of glitch sizes. To date, the catalyst has not been identified securely. If vortices unpin independently, the number of unpinning events in any time interval follows Poisson statistics (Hakonen, Avenel & Varoquaux 1998). On its own, a Poisson process does not account for the observed power-law distribution of glitch sizes. Hence, a mechanism is needed which correlates vortices, so that the unpinning rate of a particular vortex increases when another vortex (either nearby or far away) unpins; collective unpinning avalanches are essential to explaining pulsar glitches. In recent GPE simulations, Warszawski & Melatos (2010b) demonstrated two mechanisms – in addition to the growing global velocity shear as the crust spins down – for triggering and sustaining vortex avalanches: (i) acoustic radiation from moving vortices, and (ii) vortex proximity knock-on, as unpinned vortices move past their pinned neighbours.

Acoustic radiation results from the propagation of Kelvin waves on vortices (Vinen 2001), from vortex reconnection (Leadbeater et al. 2001) and due to dissipation by accelerating vortices (Parker et al. 2004; Barenghi et al. 2005). The disturbance in the density field as the sound waves propagate is capable of unpinning weakly pinned vortices. The amplitude of the disturbance is a function of vortex acceleration. Therefore, the strongest acoustic radiation results from the repinning phase, as a vortex spirals in to the new pinning site (Warszawski & Melatos 2010b). When the level of dissipation in the system is high, the acoustic disturbance generated by a repinning vortex is weaker than for a low-dissipation system, and hence it is harder to unpin vortices acoustically.

The second effect, proximity knock-on, occurs when two vortices approach sufficiently closely that the Magnus force on at least one of them exceeds the unpinning threshold. It was found that a ballistic unpinned vortex can unpin a pinned vortex at a larger radius via the proximity effect. The distance of closest approach between the two vortices is inversely proportional to pinning strength (Warszawski & Melatos 2010b).

2.2 Gross–Pitaevskii model

In this section we summarize the quantum mechanical pulsar model first introduced in Warszawski & Melatos (2010a). The interested reader who wishes to reproduce our results is referred to the latter paper for technical details. In the model, a two-dimensional, zero-temperature Bose–Einstein condensate represents an equatorial projection of the stellar superfluid. The time evolution of the superfluid density |ψ|², defined in terms of the order parameter ψ (x, t), in a potential V, with chemical potential μ, observed in a reference frame rotating with angular velocity Ω, is described by the dimensionless GPE:

2

In arriving at equation (2), we define the sound speed c_s = (n₀g/m)^1/2, the healing length (also the characteristic length-scale) ξ=ℎ/(mn₀g)^1/2, the characteristic time-scale ℎ/(n₀g) and the energy scale n₀g, where g is the self-interaction strength, m is the mass of each boson (twice the neutron mass) and n₀ is the mean particle density. The angular momentum operator is given by (with units ℎ/2), where ϕ is the azimuthal angle. V is the sum of the confining potential representing the crust and the pinning potential.

In line with standard practice (Jin et al. 1996; Choi, Morgan & Burnett 1998; Penckwitt, Ballagh & Gardiner 2002; Tsubota et al. 2002; Kasamatsu, Tsubota & Ueda 2003), we include in equation (2) a phenomenological treatment of dissipation, modelled by the term −γ∂ψ/∂t. Its effect is to suppress sound waves, thus reducing the relaxation time-scale of the superfluid and ensuring that the wavefunction can adjust to non-adiabatic changes in Ω due to glitches. On the other hand, as reported in Warszawski & Melatos (2010b), sound waves play an important role in vortex unpinning avalanches, by facilitating knock on. Therefore, one must be careful not to raise γ to the point where it quells the latter effect. An accurate estimate of γ in a pulsar is not available at present.

The dissipationless GPE (γ = 0) models a superfluid at T = 0, which does not contain viscous excited-state components. In contrast, solutions to equation (2) for γ≠ 0 are by definition at non-zero temperature. Penckwitt et al. (2002) attributed dissipation in experiments to atom transfer between the viscous thermal cloud and the condensate. Experiments with Bose–Einstein condensates (Haljan et al. 2001) showed that the thermal cloud also influences vortex nucleation. A self-consistent description of this phenomenon requires additional terms in the GPE (Wouters & Carusotto 2007), which are not included in most studies, including ours. Laboratory experiments claim γ = 0.03 in helium II (Abo-Shaeer, Raman & Ketterle 2002); we use γ = 0.05. The dissipative term drives particle loss from the system, which we correct by adjusting the chemical potential, such that the number of particles remains constant (Kasamatsu et al. 2003; Warszawski & Melatos 2010a).

2.3 Container and pinning potential

Using a fourth-order Runge–Kutta algorithm in time, and a fourth-order finite difference scheme for the spatial derivatives, we solve equation (2) for a superfluid in a rotating potential V, in the frame corotating with the potential.

V is the sum of a confining potential, representing the edge of the pulsar inner crust, and the pinning potential, V_pin, representing the nuclear lattice in the crust, which corotates with the confining potential. V_pin comprises a square grid of equidistant spikes, each one described by

3

where R_i is the position of the centre of the pinning site, and V_i and Δ⁻¹ parametrize its strength and width, respectively. In all cases, the spikes in the pinning grid are equally spaced, with a mean areal density n_pin. In Section 5.3, we use the ratio n_pin/n_F as a measure of the relative abundance of pinning sites and vortices. In Section 5.5, we draw V_i from a top-hat distribution of pinning strengths with mean V₀ and width ΔV_i; in all other simulations, V_i is constant (V_i=V₀). A non-zero range of pinning strengths (ΔV_i/V₀∼ 1) is likely to play an important role in glitch dynamics. Each small unpinning event biases the distribution of pinning strengths of occupied pinning sites towards large V_i (Newman & Sneppen 1996; Melatos & Warszawski 2009), until a large global velocity shear or strong acoustic perturbation comes along and unpins many strongly pinned vortices simultaneously.

The parameters of vortex pinning within a pulsar have been studied theoretically in great detail (Blasio & Lazzari 1998; Jones 1998b, 2003; Hirasawa & Shibazaki 2001; Avogadro et al. 2008; Barranco et al. 2010). Of particular interest is the pinning strength, the density of pinning sites and the region within the star where pinning is strongest (the superfluid pairing state changes with depth). The strength of pinning is usually calculated from the energy difference between a vortex sitting atop a column of spherical defects and a vortex positioned interstitially. In general, vortices pin to nuclear sites to minimize the region from which superfluid is excluded (the vortex core is devoid of superfluid). Depending on the superfluid density, however, pinning can also occur interstitially (Pizzochero, Viverit & Broglia 1997; Donati & Pizzochero 2004, 2006; Avogadro et al. 2007; Pizzochero 2007). For a pulsar, the characteristic size of the pinning site (the width of a nucleus) is approximately equal to the superfluid healing length (∼1 fm), which is significantly smaller than the lattice (∼10 fm) and intervortex (∼10¹² fm) spacings.

2.4 Feedback torque

Vortex motion changes the superfluid angular momentum ⁠, the expectation value of the operator ⁠. The response of the crust to these changes is an essential element of any pulsar glitch theory. We employ a simple conservation argument, wherein changes in are communicated instantaneously to the crust. In reality, changes in the superfluid flow are communicated to the crust by sound and Kelvin waves along vortex lines (Jones 1998a), but these processes are fast. Given an electromagnetic (EM) spin-down torque, N_EM, the angular velocity of the crust, Ω, evolves according to

4

where I_c is the moment of inertia of the crust, and Ω is the same angular velocity that appears in the GPE. We emphasize that, even for smooth deceleration, N_EM/I_c does not equal the observed pulsar spin-down rate, as this rate is the combination of two effects: the external EM torque, and the spin-up response to gradual (non-glitch) spin-down of the superfluid (Warszawski & Melatos 2010a). We also note that if the vortex lattice is significantly distorted by pinning, the superfluid does not rotate rigidly. Hence it is not meaningful to assign to it a unique angular velocity or moment of inertia. We solve equation (2) simultaneously with equation (4) in the following sections.

3 AUTOMATED GLITCH FINDER

Our first task is to define exactly what we mean by a glitch. In practice, glitches in pulsar data are extracted ‘by hand’, by looking for discontinuities in the slope and/or curvature of the pulse phase. Recently, Chukwude & Urama (2010) employed an automated glitch-finding algorithm, based on a least-squares goodness-of-fit parameter, to identify microglitches. Since our simulations provide extremely high time resolution (ΩΔt∼ 10⁻⁴ revolutions time-step⁻¹ compared to ΩΔt∼ 10⁵ revolutions time-step⁻¹ for the observations), we resolve the glitch spin-up phase; we can identify glitches directly by looking for spin-up events, rather than resorting to polynomial fitting of the phase residuals. The difference between simulation and astrophysical time-scales is also reflected in which physical processes can be resolved. In our simulations, we track the motion of individual vortices between pinning sites, whereas observed glitches appear as unresolved spin-up events, which likely result from a multigeneration vortex cascade.

For the purposes of this investigation, we employ the following simple algorithm for identifying glitches in Ω(t) curves.

• Smooth Ω(t) with a top-hat window function w(t) of width Δt_sm, viz. Ω_sm(t) =∫^t₀dt′w(t−t′)Ω(t′).

• Calculate ΔΩ_sm(t_i) =Ω_sm(t_i) −Ω_sm(t_{i− 1}) for all discrete time points t_i in the simulated range.

• Find all t_s for which ΔΩ_sm(t_s) > 0. For each t_s, find the smallest t_g(>t_s) such that ΔΩ_sm(t_g) ≥ 0 and ΔΩ_sm(t_{g+ 1}) < 0 (points A and C in Fig. 1 meet this criterion), and record Ω_sm(t_g).

• Locate the global Ω_sm(t) minimum between t_{g− 1} and t_g, occurring at t_{g, min} (point B in Fig. 1). Define glitch size as ΔΩ/Ω=[Ω_sm(t_g) −Ω_sm(t_{g, min})]/Ω(t_g) (the vertical distance between the dashed lines in Fig. 1), and waiting time as Δt=t_g−t_{g− 1} (the horizontal distance between points A and C in Fig. 1).

Figure 1

Schematic of how the glitch-finding algorithm operates. A and B mark the positions of the two glitches in the interval 552 < t < 630 on the plot of angular velocity as a function of time Ω(t). The size of glitch C is the difference in Ω between B (the local minimum between A and C) and C (the vertical distance between the horizontal dashed lines). Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, n_pin/n_F = 0.97.

Open in new tabDownload slide

The principal weakness of this algorithm is the stipulation in step (iii) that ΔΩ_sm > 0, which does not capture decreases in that are smaller than N_EMΔt. Depending on the moment of inertia of the crust, ΔΩ_g due to real glitches may not be greater than the negative ΔΩ_sm due to the spin-down torque [(N_EM/I_c)Δt], and hence this algorithm is not suitable for observational data.

The arrows under the spin-down curve in the top panel of Fig. 2 indicate the positions of glitches found using the above algorithm for a curve smoothed with a top-hat window function of width Δt_sm = 0.1. The glitch sizes and waiting times are then used to construct probability density functions (pdfs). In Section 4 we discuss the two bottom panels of Fig. 2, which demonstrate how we subtract a linear fit to the pre-glitch spin-down curve (bottom left-hand panel), and the resulting phase residual curve (bottom right-hand panel).

Figure 2

Timing signature of a simulated glitch. Top: angular velocity of the crust as a function of time Ω(t) for Ω(t = 0) = 0.8. Ω(t) is smoothed with a top-hat function of width Δt_sm = 1.0. The arrows mark the positions of glitches identified using the automated glitch finder described in Section 3. Bottom left: close-up of a glitch at t = 736, bracketed by the dotted rectangle in the top panel. The dotted curve represents the linear model fitted to the pre-glitch spin-down curve. Bottom right: phase residual after removing a pre-glitch linear fit to the interval 710 < t < 730 [Ω(t = 710) = 0.942, ] from the interval 710 < t < 770 shown in the bottom left-hand panel. Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, n_pin/n_F = 0.97.

Open in new tabDownload slide

4 CANONICAL GLITCH

Let us now examine a canonical example, which demonstrates the ability of the model in Section 2 to generate astrophysically relevant results. Table 1 lists the set of canonical parameters for this investigation. In the first column, the simulation parameters are given in dimensionless units; their dimensional counterparts are provided in the second column. The third and fourth columns quote the range of values associated with a typical pulsar in dimensionless and dimensional units, respectively.

4.1 Set-up

The nuclear lattice of the neutron star inner crust is represented by an 11 × 11 square grid of equal-strength pinning sites (V₀ = 16.6, ΔV_i/V₀ = 0.0), within a box of side length 15, with crust radius R = 12.5 (in dimensionless units). The simulation grid is a 200 × 200 rectangular grid; the grey-scale plots appear circular because |ψ|² = 0 outside the crust. The time-step and grid spacing are Δt = 5 × 10⁻⁴ and Δx = 0.15, respectively.

The initial conditions are established by accelerating the crust from rest to Ω = 0.8 instantaneously and allowing the system to reach a steady state (defined as when the fractional change in total energy during one time-step drops below 10⁻⁵). Once a steady-state solution is obtained, the external spin-down torque N_EM( = 10⁻³I_c) is switched on, and the superfluid-crust system is allowed to evolve.

4.2 Spin-down curve

The spin-down observed after N_EM switches on is spasmodic. By watching movies of how the density field |ψ (x, t)|² evolves, we can easily see why pinning prevents vortices from moving smoothly towards the boundary as the crust decelerates; instead, vortices hop between pinning sites, resulting in step-like decreases in 〈L_z〉. A graph of versus time is presented in Fig. 3. The positions of glitches identified using the automated glitch-finding algorithm described in Section 3 are marked with arrows. Equation (4) tells us that each step down in is accompanied by a step up in Ω (so long as the change in exceeds N_EMΔt). We see this clearly by comparing the top panel of Figs 2 and 3.

Figure 3

Superfluid angular momentum as a function of time 〈L_z〉(t), in units of ℎ/2, for the canonical spin-down experiment described in Section 4. The arrows indicate the position of glitches identified using the glitch-finding algorithm defined in Section 3. Simulation parameters: η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, n_pin/n_F = 0.97.

Open in new tabDownload slide

A close-up of one of the glitches (bracketed by the dotted rectangle in the top panel of Fig. 2) is shown in the bottom left-hand panel of Fig. 2. Does this look like a pulsar glitch? An important difference is that the spin-down rate before and after the glitch does not change; the Ω(t) curve after the glitch parallels the linear timing model describing the pre-glitch spin-down. This contrasts with some (not all) observed glitches, which frequently produce a permanent increase in the spin-down rate (Wong, Backer & Lyne 2001). The simulation output also lacks a post-glitch recovery phase, during which spin-down is faster than the global average. Given that time-scales associated with the post-glitch recovery phase, i.e. ∼10²⁹ℎ/(n₀g) (typically lasting days to weeks; Wong et al. 2001), are associated with coupling of the superfluid to a viscous fluid component not included in these simulations, this result is unsurprising (van Eysden & Melatos 2010). In our simulations, the time for Ω(t) to return to its pre-glitch value is proportional to the size of the glitch and is comparable to the sound crossing of the system, R/c_s (Warszawski & Melatos 2010a).

Practical pulsar glitch detection involves comparing measured pulse times of arrival (TOAs) against TOAs predicted by a spin-down model. When a glitch occurs, pulses arrive more frequently than the model predicts, and the phase residuals (the difference between the number of revolutions predicted by the timing model and the number of revolutions inferred from the TOAs) turn negative and grow linearly with time. In the bottom right-hand panel of Fig. 2, we plot the difference between the simulation output and a timing model defined by

5

for the interval bracketed by the dotted rectangle in the top panel. The timing model, with slope ⁠, has been extracted by performing a least-squares fit of equation (5) to Ω(t) (710 ≤t≤ 730) and is graphed as a dotted curve in the bottom left-hand panel of Fig. 2. The phase residual Δϕ, in units of radians, is therefore

6

The shape of the residual curve, graphed in the bottom right-hand panel of Fig. 2, qualitatively resembles observations; see fig. 1 in Shemar & Lyne (1996), for example.

4.3 Statistics

We now construct pdfs of sizes and waiting times using the automated glitch finder. Motivated by the hypothesis that glitches are produced by avalanche dynamics, we look for evidence that the pdfs of sizes and waiting times are distributed as power laws [p(ΔΩ_g/Ω) ∝ (ΔΩ_g/Ω)^γ] and exponentials [p(Δt) =λ exp (−λΔt)], respectively. Evidence is building that pdfs of this form yield good fits to observational data (Melatos et al. 2008). To fit these functional forms, we find the parameter that minimizes the D statistic (the maximum distance between the cumulative distribution of the model and that of the data). The fit can be rejected with 1 −p_KS confidence, where p_KS is the Kolmogorov–Smirnov (K–S) probability (Press et al. 1992).

To begin with, we note that the smoothing width Δt_sm changes the glitch statistics, as summarized in Table 2. Glitch positions are marked on the spin-down curve in Fig. 4 using arrows for Δt_sm = 0.0, 1.0, 5.0 and 10.0 (solid, dotted, dashed and dot–dashed curves, respectively). The curves have been separated vertically for clarity. As Δt_sm increases, the number of glitches decreases from N_g = 105 for Δt_sm = 0.0 to N_g = 13 for Δt_sm = 10.0. The mean waiting time 〈Δt〉 and glitch size 〈ΔΩ/Ω〉 increase five- and ninefold, respectively. The trend in 〈ΔΩ/Ω〉 versus Δt_sm is evident in the size pdfs graphed on log–log axes in the left-hand panel of Fig. 5: the pdfs for longer Δt_sm are weighted further towards the large-(ΔΩ/Ω) end. Notably, the best-fitting power-law index, γ (we warn the reader that this is not the same quantity as the dissipation parameter γ that appears in equation 2), does not change dramatically with Δt_sm; we find γ=−0.994 for Δt_sm = 0 (γ can be rejected at the 25.5 per cent confidence level), and γ=−1.104 for Δt_sm = 10.0 (rejected at the 5.3 per cent confidence level).

Figure 4

Angular velocity as a function of time Ω(t) for the canonical example described in Section 4, smoothed with a top-hat window function of width Δt_sm = 0.0, 1.0, 5.0 and 10 (solid, dotted, dashed, dot–dashed curves, respectively; the curves are shifted vertically down by 0.0, 0.1, 0.2 and 0.3, respectively, to improve visibility). Arrows below each curve mark the positions of glitches found using the glitch-finding algorithm described in Section 3; as Δt_sm increases, the number of found glitches decreases. Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, n_pin/n_F = 0.97.

Open in new tabDownload slide

Figure 5

Glitch statistics as a function of the degree of smoothing for the spin-down experiment in Fig. 4. Left: pdf of fractional glitch sizes p(ΔΩ/Ω), binned logarithmically and graphed on log–log axes, for top-hat window functions of width Δt_sm = 0.0, 1.0, 5.0 and 10 (solid, dotted, dashed, dot–dashed curves, respectively). Right: cumulative distribution of glitch waiting times P(Δt) for top-hat window functions of width Δt_sm = 0.0, 1.0, 5.0 and 10 (solid, dotted, dashed, dot–dashed curves, respectively), graphed on linear axes. Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, n_pin/n_F = 0.97.

Open in new tabDownload slide

On the other hand, increasing Δt_sm changes markedly the shape of the cumulative waiting-time distribution (right-hand panel of Fig. 5, graphed on linear axes). For Δt_sm≤ 1.0, the cumulative probability is dominated by short waiting times (i.e. convex down). For Δt_sm≥ 5.0, waiting times are more evenly distributed. Exponential fits to the cumulative probability return mean glitch rates between λ = 0.011 (for Δt_sm = 10.0; rejected at the 87 per cent confidence level) and λ = 0.87 (for Δt_sm = 0.0; rejected with the near certainty). An exponential functional form for the cumulative distribution for Δt_sm = 1.0 is rejected with the least confidence (55.1 per cent).

Glitch positions and distributions depend on Δt_sm because smoothing Ω(t) affects the morphology of individual glitches. Understanding this tendency is essential to interpreting observational data properly using our algorithm. For example, increasing Δt_sm reduces the size of a glitch, and in some cases erases the glitch altogether. To visualize this, in Fig. 6 we graph an interval of Ω(t) containing an obvious bump (when viewed on the scale graphed in Fig. 4) for different smoothing widths. The arrows in the top panel mark the positions of glitches identified by the automated algorithm for different Δt_sm. In the unsmoothed (Δt_sm = 0.0, solid) and Δt_sm = 1.0 (dotted) curves, we see two obvious bumps at t = 371 and 367.5 (for Δt_sm = 0.0 the algorithm finds eight glitches that are not visible by eye on the scale plotted). However, for Δt_sm = 5.0 and 10.0, the two bumps merge, and the algorithm finds only one glitch. For Δt_sm = 10.0, the glitch is ∼16 per cent smaller than for Δt_sm = 5.0. Increased waiting times are a direct corollary of glitches being erased by smoothing.

Figure 6

Close-up of glitches identified using the glitch-finding algorithm described in Section 3, as a function of the degree of smoothing. The plots cover the interval 356 < t < 386 (taken from the Ω curve graphed in Fig. 4). The curves are smoothed with a top-hat window function of width Δt_sm = 0.0, 1.0, 5.0 and 10 (solid, dotted, dashed and dot–dashed curves, respectively). The number and position of glitches, shown as arrows in the top panel, changes with Δt_sm; for Δt_sm≥ 5, one glitch is detected. Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, n_pin/n_F = 0.97.

Open in new tabDownload slide

After watching movies of |ψ |²(x, t), we conclude that glitches too small to be discerned by eye on the scale depicted in Fig. 6 are associated with ‘jiggling’, not unpinning. By ‘jiggling’ we mean vortices wobble around their pinning site without unpinning. Multipeaked glitches arise from multistage vortex motion (i.e. one vortex moves, causing a second vortex to move a short time later), precisely the collective avalanche behaviour we are seeking. Observationally, we can only detect pulsar glitches resulting from the motion of many vortices (the smallest observed glitches correspond to ∼10⁷ unpinnings). For this reason, we henceforth adopt Δt_sm = 1.0, which also happens to produce the most exponential P(Δt). Thus we retain multipeaked glitches like the one in Fig. 6, whilst excluding glitches not associated with unpinning.

An important recent result from pulsar data analysis is the absence of a reservoir effect: glitch size is uncorrelated with waiting time (Wong et al. 2001; Melatos et al. 2008). In Table 2 we list the Pearson correlation coefficient, r_p, relating ΔΩ/Ω and Δt as a function of Δt_sm. For Δt_sm < 5.0, the two quantities are uncorrelated (r_p=−0.057 and 0.050 for Δt_sm = 0.0 and 1.0, respectively). For Δt_sm≥ 5.0, r_p differs significantly from zero (0.428 and 0.334 for Δt_sm = 5.0 and 10.0, respectively), suggesting that sizes and waiting times are correlated. The interpretation of this correlation depends on the glitch definition. Since the correlation emerges for smoothing scales that erase multipeaked glitches and ‘jiggling’, it is likely that this correlation exists if the glitch-finding algorithm does not distinguish between single- and multiple-vortex motion. Simulations large enough to house several isolated regions of unpinning activity should not exhibit such a correlation (Melatos et al. 2008).

As observational duty cycles and measurement precision improve, e.g. with the advent of multibeam telescopes, astronomers will begin to resolve the spin-up phase of a glitch, multipeaked glitches and post-glitch oscillations (Kramer & Stappers 2010). Our results suggest that improved measurement precision (equivalent to reducing Δt_sm) should lead to a large increase in the number of observed glitches, without altering the shape of the size distribution. For example, reducing Δt_sm fivefold doubles the number of glitches. Conversely, the mean glitch rate increases. More speculatively, we also predict that the low-ΔΩ/Ω end of the glitch size pdf is dominated by ongoing, non-collective vortex creep (see Alpar et al. 1984b, for example), likely characterized by a Gaussian cut-off.

5 PINNING PHYSICS

Systematic differences in the pinning strength from pulsar to pulsar may result from different cooling histories (Jones 2001), which alter the morphology and abundance of crystal defects. The role of the pinning strength V₀ in glitch physics is crucial. If pinning is weak, vortices spread out smoothly and homologously, allowing the superfluid to decelerate in sympathy with the crust; differential rotation does not build-up in stress reservoirs. On the other hand, if pinning is strong, it may stiffen the vortex lattice to the point where the crust cracks before the vortices are unpinned by the Magnus force (Horowitz & Kadau 2009). Intuitively, we expect that stronger pinning results in larger glitches; either vortices move further once they have unpinned, or more vortices unpin in order to reverse the crust-superfluid shear. In this section we use simulations to test these intuitive relations between glitch distributions and pinning strength.

5.1 Equal pinning potentials

We begin by examining how spin-down depends on V₀, in the situation where all the pinning sites are of equal strength (ΔV_i = 0). The four grey-scale plots of superfluid density at the top of Fig. 7 are the final states (t = 1025) for simulations with V₀ = 0, 8.3, 16.6 and 33.3 (left to right, respectively). Vortices are easily identified as darker spots. The lighter spots represent unoccupied pinning sites. All vortices, except for one in the north-east quadrant of the second panel, are pinned.

Figure 7

Vortex configurations and spin-down curves for different pinning strengths V₀. Top: grey-scale plots of the final (t = 1000) superfluid density for V_i = 0, 8.3, 16.6 and 33.7 (left to right, respectively). The density grey-scale runs from dark (low density) to light (high density). Bottom: angular velocity of the crust as a function of time Ω(t) for V_i = 0, 8.3, 16.6 and 33.7 (solid, dotted, dashed and dot–dashed, respectively). Simulation parameters: η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, n_pin/n_F = 0.97, Δt_sm = 1.0.

Open in new tabDownload slide

Without pinning, the vortex lattice expands smoothly and homologously, decreases smoothly and the superfluid exerts a continuous spin-up torque on the crust via equation (4). The spin-up torque lessens the gradient of Ω (the black curve in Fig. 7) when compared with no feedback, linear spin-down (the solid grey curve). Small bumps in the black curve in Fig. 7 are numerical noise.

We now switch on pinning, via an 11 × 11 equidistant, equal-strength grid. exhibits step-like decreases, accompanied by jumps in Ω, in the dotted, dashed and dot–dashed curves in Fig. 7. The first glitch (at t₁) occurs later when the pinning is stronger, viz. t₁ = 200 for V₀ = 8.3 compared to t₁ = 440 for V₀ = 33.3. The delay can be understood by considering the velocity shear between pinned vortices and the superfluid in the absence of feedback:

7

Whilst vortices remain pinned, v_s is approximately constant. Therefore, it takes longer for the Magnus force to unpin strongly pinned vortices.

Even in the strongest pinning case, before any unpinning has occurred, the spin-down curve is less steep than Ω(t) =tN_EM/I_c (the solid grey curve in Fig. 7). The reason is subtle; it is discussed in detail in Warszawski & Melatos (2010a). While remaining pinned, vortices respond to increasing differential rotation by migrating smoothly to the outer edge of their pinning site. Hence the crust experiences a gradual spin-up torque. Put another way, pinning is not a binary condition; it is weakened gradually as the shear grows.

At the end of the simulation, Ω is highest (Ω = 0.29) for the weakest pinning and lowest (Ω = 0.16) for the strongest pinning, because stronger pinning potentials can sustain greater levels of stress. [Note that this is not the reservoir effect (Link, Epstein & Lattimer 1999; Wong et al. 2001; Melatos et al. 2008), which predicts a correlation between glitch size and waiting time not found in observational data.] In other words, the areal density of remaining vortices after a fixed amount of time is proportional to the pinning strength.

At this point we draw attention to vast discrepancies between typical pulsar parameters and the physical scales modelled by our simulations. The extremely short simulation time (10⁻²¹ s in physical units), and the extremely high angular velocity (Ω = 10²⁴ Hz in pulsar units compared to the typical pulsar angular velocity Ω∼ 10 Hz) represent highly non-adiabatic conditions.

5.2 Statistics

A key goal of this paper is to understand the observed scale invariance of glitch sizes, which is believed to arise from collective vortex behaviour (Melatos et al. 2008). Global properties of the glitch distribution, such as the number of glitches, mean glitch size and mean waiting time for each V₀ studied in Section 5.1, are listed in Table 3; along with the activity parameter

8

used by pulsar astronomers to characterize the total glitch-induced spin-up of a pulsar (Lyne et al. 1998).

In Fig. 8 we graph logarithmically binned size pdfs on log–log axes. Our glitch-finding algorithm identifies 16, 48, 23 and 26 glitches for V₀ = 0, 8.3, 16.6 and 33.3, respectively. The pdfs for simulations with pinning (V₀ > 0) exhibit a bias towards larger glitches for larger V₀, which also shows up as a ∼fourfold increase in mean glitch size for V₀ = 33.3 compared to V₀ = 8.3. Power-law fits to the pdfs yield similar indices (listed with the K–S probability in parenthesis in the penultimate column of Table 3) for 8.3 ≤V₀≤ 33.3. However, the confidence level with which the fits can be rejected ranges from 23.8 per cent for V₀ = 16.6 to 70.7 per cent for V₀ = 33.3. We remind the reader that glitches detected in the V₀ = 0 case arise from numerical noise, rather than vortex unpinning; this assertion is supported by the small mean glitch size (10⁴〈ΔΩ/Ω〉 = 13.00).

Figure 8

Glitch statistics for different pinning strengths corresponding to the Ω curves graphed in Fig. 7. Left: pdf of fractional glitch sizes p(ΔΩ/Ω) for V₀ = 0, 8.3, 16.6 and 33.7 (solid, dotted, dashed and dot–dashed, respectively) logarithmically binned and graphed on log–log axes. The dotted curve is a power-law least-squares fit to the dotted histogram with index −0.818. Right: cumulative probability of glitch waiting times P(Δt) graphed on log-linear axes for V₀ = 0, 8.3, 16.6 and 33.7 (solid, dotted, dashed and dot–dashed, respectively). Simulation parameters: η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, n_pin/n_F = 0.97, Δt_sm = 1.0.

Open in new tabDownload slide

In the right-hand panel of Fig. 8 we graph cumulative distributions of waiting times for the simulations shown in Fig. 7. Mean glitch rates, λ, derived from exponential fits, are listed in the final column of Table 3. In Melatos et al. (2008), it was found that waiting times for seven of the nine pulsars that have glitched more than six times are well represented by exponential distributions. Taking the results in Fig. 8 all together, we suggest that λ does not depend monotonically on V₀ because there is another degree of freedom, 〈ΔΩ/Ω〉, through which the additional stress can be released.

The activity parameter does not appear to scale monotonically with V₀. However, A almost doubles between V₀ = 8.3 and 33.3. We do not have a physical explanation for this behaviour at this time.

5.3 Pinning site abundance

If pinning occurs at individual nuclei, then ∼10⁹ pinning sites separates neighbouring vortices in a pulsar. Alternatively, if pinning occurs at grain boundaries or along macroscopic faults in the crustal lattice, the ratio of pinning sites to vortices is not necessarily large, and pinning sites may not be evenly spaced. In this section we explore the role of the pinning site density (n_pin sites per unit area) relative to vortex density (n_F = 2Ω/κ vortices per unit area). We present results for simulations with (i) few or no pinning sites, (ii) n_pin∼n_F and (iii) many more pinning sites than vortices. We parametrize the pinning site abundance in terms of the ratio n_pin/n_F in the initial state.

Four pinning scenarios are illustrated in the top panels of Fig. 9 as grey-scale plots of |ψ |²(x, t = 1000) for n_pin/n_F increasing from left to right. The middle panels both correspond to regime (ii). In every panel except the leftmost one (n_pin/n_F = 0.11), all vortices in the final state are pinned. For n_pin/n_F = 0.11, even in the final state, there are only enough pinning sites to pin nine out of 17 vortices (i.e. there are no unoccupied pinning sites). Notably, for n_pin/n_F = 0.11, the vortex lattice is regularly spaced, suggesting that pinning alters the orientation, but not the geometry, of the vortex lattice. For n_pin/n_F = 0.52 and n_pin/n_F = 0.97, vortices are irregularly spaced, as they are forced to adhere to the geometry of the pinning grid. When pinning sites significantly outnumber vortices (right-hand panel, n_pin/n_F = 1.90), the vortex lattice is once again evenly spaced, because enough pinning sites are available for a vortex to ‘choose’ a site near its equilibrium (Abrikosov) position. It should be noted that this equilibrium position changes as the crust spins down.

Figure 9

Results from numerical experiments for different pinning abundances n_pin/n_F. Top: grey-scale plots of the final superfluid density for n_pin/n_F = 0.11, 0.52, 0.97 and 1.90 (left to right, respectively). The density scale runs from dark (low density) to light (high density). The faint dots are unoccupied pinning sites; darker dots are pinned vortices. Bottom: angular velocity of the crust as a function of time Ω(t) for n_pin/n_F = 0.11, 0.52, 0.97 and 1.90 (solid, dotted, dashed and dot–dashed, respectively). Simulation parameters: V_i = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

The bottom panel of Fig. 9 graphs the spin-down curves for the four pinning abundances illustrated in the top panels. The example with few pinning sites (n_pin/n_F = 0.11, solid curve) is relatively smooth. Although our algorithm detects 25 glitches, they are small (〈ΔΩ/Ω〉 = 14.84 × 10⁻⁴; cf. 〈ΔΩ/Ω〉 = 13.00 × 10⁻⁴ for V₀ = 0) and well separated (〈Δt〉 = 37.10), suggesting that they arise in the rare event that an outwardly moving vortex runs over a pinning site. By contrast, for n_pin∼n_F (dotted, dashed and dash–dotted curves), the spin-down is roughly linear until vortices begin to unpin at t∼ 275, whereupon vortices move towards the edge in spasmodic jumps between pinning sites (vortices do not stop at all pinning sites in their path). For all n_pin/n_F > 0.11, the gradual outward migration of vortices across their pinning sites reduces ⁠, even before vortices have unpinned. For n_pin/n_F = 0.11, is considerably less than N_EM/I_c for t≳ 50, as the unpinned vortices move smoothly outward almost as soon as the spin-down torque is switched on.

It is physically interesting that the total spin-down of the crust over long time periods is almost identical for n_pin/n_F > 1 and n_pin/n_F∼ 1. The aggregate effect of glitches is insensitive to increases in n_pin/n_F over and above unity. One interpretation is to say that the average distance travelled by a vortex once it unpins is of order n^−1/2_F, rather than n^−1/2_pin, so that vortices maintain an approximate Abrikosov lattice, without developing the large-scale spatial inhomogeneities required in self-organized critical avalanche models (Warszawski & Melatos 2008). However, this hypothesis is challenged by three trends. First, the number of glitches and activity parameter increase for n_pin/n_F > 0.11; these qualities scale with n^−1/2_pin, even though n^−1/2_F is approximately fixed. Secondly, the mean waiting time is inversely proportional to n_pin/n_F; see Table 4 for details and the right-hand panel of Fig. 9 for cumulative distributions of Δt. Thirdly, the mean glitch size does not change monotonically with n_pin/n_F. Taken together, these results suggest that the typical distance travelled by a vortex in a glitch is a play off between the number of vortices that unpin, the distance to the next available pinning site and the distance the vortex must move to reverse the shear that unpinned it.

5.4 Maximum size

Glitch size is governed by the total change in and hence depends on the number of vortices that unpin, the distance that each vortex travels before repinning and the location from which it unpinned. From the pdfs of ΔΩ/Ω in Fig. 10, we see that the largest glitch is approximately the same for all n_pin≳n_F. If vortex motion is determined by the first encounter with a pinning site, rather than the imperative to reduce shear and maintain equal spacing, then we expect the mean glitch size to scale inversely with the density of pinning sites. This is not observed. By contrast, maximum glitch size does scale monotonically with V₀[max(10⁴ΔΩ/Ω) = 36, 538, 1057 and 1641 for V₀ = 0, 8.3, 16.6 and 33.3, respectively] because vortices need to move further to reduce the shear sustained by stronger pinning.

Figure 10

Glitch distributions for spin-down experiments with different pinning abundances n_pin/n_F for the Ω curves graphed in Fig. 9. Left: pdf of fractional glitch sizes p(ΔΩ/Ω) for n_pin/n_F = 0.11, 0.52, 0.97 and 1.90 (solid, dotted, dashed and dot–dashed, respectively) logarithmically binned and graphed on log–log axes. The solid grey curve is a power-law least-squares fit to the solid histogram with index −0.801. Right: cumulative probability of glitch waiting times p(Δt) graphed on log-linear axes for n_pin/n_F = 0.11, 0.52, 0.97 and 1.90 (solid, dotted, dashed and dot–dashed, respectively). Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

5.5 Unequal pinning potentials

Pinned vortices at a given radius experience an equal Magnus force. If every pinning site is of equal strength, then vortices at a given radius should unpin simultaneously, yielding periodic glitches of equal size, contrary to observations (Melatos et al. 2008; Melatos & Warszawski 2009). In this section, we investigate the effect of varying the width of the pinning strength distribution ΔV_i/V₀. The glitch statistics from these simulations are summarized in Table 5.

The spin-down curves in Fig. 11 demonstrate that varying ΔV_i/V₀ makes little difference to the total spin-down over the time-scale of the simulation; Ω(t = 1000) is approximately the same for all 0 ≤ΔV_i/V₀≤ 1. However, the timing of the first departure from linear spin-down scales as t_g1≈ 1250 + 200(1 −ΔV_i/V₀). In addition, the glitch statistics change systematically with ΔV_i/V₀. The pdfs of ΔΩ/Ω graphed in the left-hand panel of Fig. 12 confirm that the minimum glitch size depends on ΔV_i/V₀, via the minimum pinning strength. The correlation is not strict because the number of draws from the pinning strength distribution is relatively small (11 × 11 = 121), so the minimum V_i does not monotonically scale with ΔV_i/V₀. We detect no systematic change in the exponent of power-law fits to the size pdf. However, we note that the fits can be rejected with high confidence for ΔV_i/V₀ = 0 (78.8 per cent) and ΔV_i/V₀ = 1.0 (74.7 per cent).

Figure 11

Angular velocity of the crust as a function of time Ω(t) for spin-down experiments with ΔV_i/V₀ = 0.0, 0.25, 0.5, 0.75 and 1.0 (solid, dotted, dashed, dot–dashed and triple-dot–dashed, respectively). Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

Figure 12

Glitch distributions for spin-down experiments with different ΔV_i/V₀ for the Ω(t) curves graphed in Fig. 11. Left: pdf of fractional glitch sizes (ΔΩ/Ω) for ΔV_i/V₀ = 0.0, 0.25, 0.5, 0.75 and 1.0 (solid, dotted, dashed, dot–dashed and triple-dot–dashed, respectively) logarithmically binned and graphed on log–log axes. The solid grey curve is a power-law least-squares fit to the solid histogram with index −0.818. Right: cumulative probability of glitch waiting times p(Δt) graphed on log-linear axes for ΔV_i/V₀ = 0.0, 0.25, 0.5, 0.75 and 1.0 (solid, dotted, dashed, dot–dashed and triple-dot–dashed, respectively). Simulation parameters: V₀ = 16.6, η = 1, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

Cumulative waiting time distributions are graphed in the right-hand panel of Fig. 12. Exponential fits reveal a tendency for the mean glitch rate to increase with ΔV_i/V₀ (λ = 0.07 for ΔV_i/V₀ = 1.0 and 0.043 for ΔV_i/V₀ = 0.0). This finding disagrees with the uniform pinning scenarios discussed in Sections 5.1 and 5.2, in which λ does not scale inversely with V₀.

To interpret the above findings, we suggest that, even for ΔV_i/V₀ = 0, the effective potential at each pinning site is a function of its departure from an equilibrium position in the Abrikosov lattice. If a vortex is slightly inside or outside (radially) its Feynman position, the effective pinning energy is decreased or increased, respectively, and hence the distribution of pinning strengths can be treated as heterogeneous, even when ΔV_i = 0. In a pulsar (n_pin/n_F∼ 10¹⁴ if all nuclear sites can pin; see Table 1), the pinned vortices resemble an Abrikosov array. Hence we expect that ΔV_i/V₀ is more influential in a pulsar than in the simulations presented here.

6 STELLAR PARAMETERS

6.1 Crust moment of inertia

The inertia of the pulsar crust, relative to the inertia of the pinned superfluid, controls the crust’s response to changes in ⁠. Let η=I_c/I_S be the ratio of the crust’s moment of inertia, I_c, to the moment of inertia that the superfluid would have if it were a rigid body. If the crust is light (small η), small decreases in produce large crustal spin-up events. For sufficiently small η, even vortex creep (Alpar et al. 1984b; Link, Epstein & Baym 1993) should produce small but observable glitches. For larger η, on the other hand, crustal spin-up is smaller. In what follows, we adopt the view proffered by Link et al. (1999) that the ‘crust’ is composed of the rigid nuclear lattice, all charged fluid components that are coupled electromagnetically to the rigid lattice and the unpinned superfluid (which may be a large fraction of the total). Hydrodynamic models of glitch recovery suggest η∼ 1 (van Eysden & Melatos 2010). Results for the numerical experiments presented in this section are summarized in Table 6.

As η increases, the mean glitch size diminishes; a given change in results in a smaller change in Ω. For example, from the spin-down curve in Fig. 13, we see that doubling η halves Ω(t = 1000). Put another way, as η increases, unpinning and outward vortex motion do progressively less to reduce the crust-superfluid shear. Hence, the timing of the unpinning events changes with η. Cumulative waiting time distributions are graphed in the right-hand panel of Fig. 14. They show that the mean waiting time decreases as η increases. The waiting-time distribution is weighted to the low Δt end for η = 2.0 compared to η = 1.0 (e.g. 〈Δt〉 = 27.93, 25.19 and λ = 0.035, 0.065 for η = 1.0, 2.0).

Figure 13

Angular velocity of the crust as a function of time Ω(t) for spin-down experiments with lighter (η = 1) and heavier (η = 2) crusts (solid and dotted curves, respectively). Simulation parameters: V₀ = 16.6, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

Figure 14

Glitch distributions for spin-down experiments with different ΔV_i/V₀ for the Ω(t) curves graphed in Fig. 13. Left: pdf of fractional glitch sizes (ΔΩ/Ω) for η = 1 and 2 (solid and dotted curves, respectively) logarithmically binned and graphed on log–log axes. The dotted grey curve is a power-law least-squares fit to the dotted histogram with index −0.846. Right: cumulative pdf of glitch waiting times p(Δt) graphed on log-linear axes for η = 1 and 2 (solid and dotted curves, respectively). Simulation parameters: V_i = 16.6, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, N_EM/I_c = 10⁻³, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

The resistance of a heavy crust to spin-up also produces smaller glitches. This result is confirmed by considering the pdf of ΔΩ/Ω for two cases, η = 1.0 and 2.0. Looking at the left-hand panel in Fig. 14, we find that the heavier crust is accompanied by smaller glitches, e.g. 10⁴〈ΔΩ/Ω〉 = 188, 74.73 for η = 1.0, 2.0.

The vertical distance between the curves in Fig. 13 increases with time, demonstrating that the increased glitch rate is not sufficient to compensate for reduced glitch sizes. Over longer time periods than those simulated here, self-regulates so that the crust-superfluid shear exactly equals the value required to sustain the average unpinning rate leading to (given N_EM and η). The activity parameter A provides a means of comparing the relative changes in response to the spin-down torque in 〈Δt〉 and 〈ΔΩ/Ω〉; A decreases by ∼33 per cent as η doubles from η = 1.0 to 2.0, demonstrating again that the increased glitch rate does not compensate for the large decrease in glitch size.

6.2 Electromagnetic spin-down torque

The magnetic dipole torque acting on a pulsar ultimately drives the crust-superfluid shear essential to the (un)pinning model of pulsar glitches. Here, we quantify the torque in terms of an electromagnetic spin-down rate, N_EM/I_c, in the absence of feedback. The larger the torque, the more quickly Magnus stresses build up in the pinned vortex lattice. Results from simulations in which N_EM/I_c is varied are summarized in Table 7.

In Fig. 15 we graph spin-down curves for N_EM/I_c = 10^−2.5 (solid), 10⁻³ (dashed), 10^−3.5 (dotted) and 10⁻⁴ (dot–dash). For N_EM/I_c = 10^−4.0, we detect few glitches (N_g = 3), owing to the small shear that accumulates over the total lifetime of the simulation [t_totalN_EM/I_c = 0.4]. Importantly, the first glitch does not appear at the same value of Ω for different spin-down torques; Ω(t₁) = 0.45, 0.46, 0.61 and 0.44 for N_EM/I_c = 10^−4.0, 10^−3.5, 10⁻³ and 10^−2.5, respectively. This result reinforces the findings of hysteresis reported in Jackson & Barenghi (2006) and Warszawski & Melatos (2010a): the spin-down history, not merely the current Ω, determines the number and position of vortices.

Figure 15

Angular velocity of the crust as a function of time Ω(t) for spin-down experiments with different electromagnetic spin-down torques N_EM/I_c=−10^−2.5, −10^−3.0, −10^−3.5 and −10^−4.0 (solid, dotted, dashed and dot–dashed curves, respectively). Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

For a crust-superfluid system with η = 1.0, the observed crust spin-down rate is N_EM/(2I_c); see equation (4) with dL_z/dt=I_cdΩ/dt. For N_g≫ 1, the results in Fig. 15 demonstrate that the fraction of the torque-driven spin-down t_totalN_EM/I_c that is reversed by glitches decreases as N_EM/I_c increases. Quantitatively, we find I_c〈ΔΩ_g〉/(〈Δt〉N_EM) = 0.269, 0.245 and 0.137 for N_EM/I_c = 10^−3.5, 10^−3.0 and 10^−4.0, respectively. We interpret this to mean that weaker torques drive the system adiabatically, allowing the superfluid to match the crust’s deceleration. In a pulsar, the slow spin-down (⁠⁠) corresponds to the low-N_EM regime. Hence, we can expect on large time-scales the superfluid deceleration matches that of the crust, thus avoiding unfettered growth of differential motion.

Size and waiting-time distributions for N_EM/I_c > 10^−4.0 are plotted as logarithmically binned pdfs and cumulative probability distributions in the left- and right-hand panels of Fig. 16, respectively. As N_EM/I_c increases, one obtains more frequent, smaller glitches. The results reported in Table 7 show that the slopes of power-law fits to the size pdf also increase monotonically with N_EM/I_c, as does the rate extracted from exponential fits to waiting-time pdfs (from λ = 0.0289 to 0.112 for N_EM/I_c = 10^−3.5 and 10^−2.5, respectively; since the N_EM/I_c = 10⁻⁴ curve contains only three glitches, we do not trust the fitted λ). Systematic increases in the activity parameter A (from A = 0.0958 to 73.98 for N_EM/I_c = 10^−4.0 and 10^−2.5, respectively) confirm that increasing N_EM/I_c does not trigger the same vortex movements as for smaller N_EM/I_c.

Figure 16

Glitch distributions for spin-down experiments with different electromagnetic spin-down torques N_EM/I_c for the Ω(t) curves graphed in Fig. 15. Left: pdf of fractional glitch sizes (ΔΩ/Ω) for N_EM/I_c=−10^−2.5, −10^−3.0, −10^−3.5 and −10^−4.0 (solid, dotted, dashed and dot–dashed curves, respectively) logarithmically binned and graphed on log–log axes. The solid grey curve is a power-law least-squares fit to the solid histogram with index −0.984. Right: cumulative probability of glitch waiting times p(Δt) graphed on log-linear axes for N_EM/I_c=−10^−2.5, −10^−3.0, −10^−3.5 and −10^−4.0 (solid, dotted, dashed and dot–dashed curves, respectively). Simulation parameters: V₀ = 16.6, η = 1, ΔV_i/V₀ = 0.0, R = 12.5, Δx = 0.15, Δt = 5 × 10⁻⁴, Ω(t = 0) = 0.8, Δt_sm = 1.0.

Open in new tabDownload slide

7 LARGE-SCALE SIMULATIONS

The scale-invariant statistics of pulsar glitches are expected to emerge most clearly in the many-vortex limit. In this section we increase the system size, and hence the vortex number, by an order of magnitude. Prevented by computational expense from undertaking a full parameter study, we present examples of two systems with uniform and variable pinning strengths (ΔV_i/V₀ = 0.0 and 0.8, respectively), each with ∼200 vortices in total. The simulations ran for ∼2 weeks on a single 2.5-GHz processor. Statistics from the simulations are summarized in Table 8.

To maintain the same spatial resolution as in the simulation, we increase the crust radius to R = 38 and impose a 66 × 66 pinning grid (with V₀ = 36). We also increase the number of particles in the system 10-fold (via n₀), to provide sufficient contrast to distinguish vortices from the background density. We initialize the experiment by instantaneously accelerating the crust from rest to Ω = 0.8 and letting the system reach a steady state. Once a steady state is achieved, the crust is subjected to a spin-down torque N_EM = 10⁻³I_c.

The grey-scale plot in the left-hand panel of Fig. 17 depicts the initial superfluid density for uniform pinning (ΔV_i/V₀ = 0.0). The two insets zoom in on regions where vortices sit interstitially. The smooth increase in |ψ|² with radius, indicated by the radial colour gradient, represents the standard parabolic density profile, analogous to the meniscus of a rotating bucket of water (Warszawski & Melatos 2010a). Initially the vortices (dark dots in Fig. 17; light dots are unoccupied pinning sites) are surrounded by a vortex-free corona in the region 4 ≤r≤ 38, highlighting the distortion and hysteresis wrought by the pinning grid (see section IVB of Warszawski & Melatos 2010a). By the end of the simulation, the vortices are approximately evenly distributed across the pinning grid; see the grey-scale plot in the right-hand panel of Fig. 17. The radial density gradient is shallower because Ω is lower than its initial value. We can see that vortices have moved from their initial positions; vortices present in the insets in the left-hand panel of Fig. 17 are absent from the right-hand panel.

Figure 17

Grey-scale plots of the initial (Ω = 0.8, left) and final (Ω = 0.2, right) superfluid density |ψ|² for a 600 × 600 simulation with R = 38 and V_i = 36. The insets zoom in on regions with dimensions 3 × 3. Vortices are recognizable as the darker dots, initially clustered towards the centre of the crust and evenly spaced at the end of the simulation. Unpinned vortices can be seen in the initial but not the final insets. Simulation parameters: V₀ = 36, η = 1, N_EM/I_c = 10⁻³, R = 38, Δx = 0.15, Δt = 5 × 10⁻⁴, Ω(t = 0) = 0.8, Δt_sm = 0.1.

Open in new tabDownload slide

The spin-down curves for ΔV_i/V₀ = 0.0 and 0.8 (solid and dotted curves in Fig. 18, respectively), smoothed with a top-hat window function of width Δt_sm = 0.1, reveal only small glitches (10⁴〈ΔΩ_g/Ω_g〉 = 0.752 and 0.591 for ΔV_i/V₀ = 0.0 and 0.8, respectively). This result is unsurprising. In larger systems, each vortex corresponds to a smaller fraction of the total angular momentum. To achieve glitches as large as those produced by a R = 12.5 system, ∼10 times more vortices must unpin simultaneously. Dissipation, which suppresses sound waves, renders knock-on events less likely; hence such large-scale collective unpinning is rare in the systems simulated here. The insets show close-ups of the largest glitch (top right and bottom left for ΔV_i/V₀ = 0.0 and 0.8, respectively). The spin-up phase of the largest glitch in the ΔV_i/V₀ = 0.0 curve is bumpy. In fact, the vertical distance between the bumps is of order the height of the largest glitch in the ΔV_i/V₀ = 0.8 curve, which does not exhibit the same ‘bumpiness’.

Figure 18

Angular velocity of the crust as a function of time Ω(t) for large-scale spin-down experiments with ΔV_i/V₀ = 0.0 and 0.8 (solid and dotted curves, respectively). Grey-scale plots of the initial and final superfluid density |ψ|² for ΔV_i/V₀ = 0.0 are shown in Fig. 17. Left inset: close-up of glitch at t = 804.3 from ΔV_i/V₀ = 0.8 curve. Right inset: close-up of glitch at t = 806.1 from ΔV_i/V₀ = 0.0 curve. Glitch positions are marked as arrows. Simulation parameters: V₀ = 36, η = 1, N_EM/I_c = 10⁻³, R = 38, Δx = 0.15, Δt = 5 × 10⁻⁴, Ω(t = 0) = 0.8, Δt_sm = 0.1.

Open in new tabDownload slide

Pdfs of ΔΩ/Ω and cumulative distributions of Δt are graphed in the left- and right-hand columns of Fig. 19, respectively; the solid curves are for ΔV_i/V₀ = 0.0 and the dotted curves are for ΔV_i/V₀ = 0.8. Since there is an obvious turnover in the size distribution at ΔΩ/Ω≈ 7 × 10⁻⁷, power-law fits are to p(ΔΩ/Ω≥ 7 × 10⁻⁷). The top panels refer to glitches detected in the unsmoothed Ω curve, whereas the bottom panels are smoothed with Δt_sm = 0.1. For both values of ΔV_i/V₀, smoothing flattens the size distribution. Comparison with Fig. 8 reveals that the smallest glitch in the larger system is six times smaller than the smallest glitches for the R = 12.5 systems discussed in previous sections, which agrees with the ∼sixfold increase in vortex number N_v. If a single vortex is responsible for the smallest glitch, then the size of the smallest glitch should scale inversely with N_v. For Δt_sm = 0.1, the fitted power law is marginally shallower for ΔV_i/V₀ = 0.8 (γ=−1.05) than for ΔV_i/V₀ = 0.0 (γ=−1.23).

Figure 19

Left: pdf of fractional glitch sizes (ΔΩ/Ω) for the Ω(t) curves graphed in Fig. 18 for ΔV_i/V₀ = 0.0 and 0.8 (solid and dotted curves, respectively) unsmoothed (top) and smoothed with top-hat window function of width Δt_sm = 0.1 (bottom), logarithmically binned and graphed on log–log axes. The solid grey curves are power-law least-squares fits to the solid histograms for ΔΩ/Ω≥ 7 × 10⁻⁷, with indices −1.26 and −1.23, respectively. Right: cumulative pdf of glitch waiting times p(Δt) graphed on log-linear axes with ΔV_i/V₀ = 0.0 and 0.8 (solid and dotted curves, respectively) unsmoothed (top) and smoothed with top-hat window function of width Δt_sm = 0.1 (bottom). Simulation parameters: V₀ = 36, η = 1, N_EM/I_c = 10⁻³, R = 38, Δx = 0.15, Δt = 5 × 10⁻⁴, Ω(t = 0) = 0.8, Δt_sm = 0.1.

Open in new tabDownload slide

The waiting times tend to lengthen as ΔV_i/V₀ increases, as observed in smaller systems (see Fig. 12). For ΔV_i/V₀ = 0.0, smoothing broadens the waiting-time distribution significantly. The top right-hand panel of Fig. 12 reveals a strong periodicity at Δt≈ 0.05 in the waiting-time distribution for Δt_sm = 0.0. This may be evidence of the periodic glitches expected from a system in which stress increases over time by the same amount on vortices pinned with equal strength. Despite the overall slower spin-down for ΔV_i/V₀ = 0.8 (cf. Fig. 18), the longer waiting times are not matched by larger glitches. This is likely due to one of two reasons: (i) much of the crust-superfluid shear is reversed by smooth vortex motion, unhindered by weak pinning or (ii) glitches are too small to cause a bump (positive change) in Ω, and hence are not detected by our glitch-finding algorithm.

Two significant differences between this experiment and those described in Sections 4–6 prevent us from attributing the altered spin-down behaviour (smaller glitches) entirely to system size. First, the ratio of pinning sites to vortices is significantly higher (n_pin/n_F≈ 5) in this experiment than in the smaller systems, whilst the spacing of the pinning sites, relative to the vortex core size, is smaller. It may be that such a closely spaced pinning grid presents the vortices with a nearly homogeneous, elevated potential, effectively weakening the pinning strength. Secondly, since the Magnus force on pinned vortices is proportional to |v−v_s | ≈r(〈L_z〉/I_s−Ω) at the pinning site, the larger radii of many of the pinned vortices in the R = 38 system make for larger Magnus forces for smaller values of (〈L_z〉/I_s−Ω).

8 DISCUSSION

The numerical experiments described in this paper demonstrate from first principles that a superfluid coupled to its container via vortex pinning spins down spasmodically. The experiments are designed to explore how the distributions of pulsar glitch sizes and waiting times depend on pinning and stellar parameters. Computational limitations prevent us from probing the parameter regime most appropriate to the pulsar glitch problem. However, we do our best to mimic astrophysical conditions by adopting the same ordering of dimensionless parameters as in nature. Qualitatively, and with respect to the canonical parameters listed in Table 1, the pulsar problem involves weak (but still effective) pinning, slow electromagnetic spin-down, a huge number of vortices, an even larger number of pinning sites and a heavy crust.

In order to speed up the process of compiling glitch statistics and minimize its subjective element, we devise an automated glitch finder. Since the glitch-finding algorithm identifies all positive jumps in ΔΩ, it is necessary to smooth Ω so as to filter out post-glitch oscillations or numerical noise. In the future, improvements in timing and observational duty cycle should allow real pulsar spin-down curves to be analysed in the same way.

Our simulations lend new insights into aspects of the physics of pulsar glitches.

• Stronger pinning causes larger glitches, with 〈ΔΩ_g/Ω〉≈ 7.6V_i. This is not as obvious as it may sound. Prima facie it is equally likely that strong pinning leads to a large average shear with small, frequent fluctuations. In fact, this does not happen.

• When pinning sites and vortices are equally numerous, vortices conform to the geometry of the pinning grid, producing larger glitches than if pinning sites are more abundant than vortices. Waiting times are also longer for n_pin∼n_F than for n_pin > n_F.

• Glitch sizes are insensitive to increases in the pinning site abundance beyond parity (i.e. n_F∼n_pin), because vortices move across many pinning sites to negate the local shear, rather than stopping at an adjacent pinning site. The distance travelled by an unpinned vortex is a key input to heuristic glitch models (e.g. Warszawski & Melatos 2008; Melatos & Warszawski 2009).

• A broader range of V_i results in a broader range of glitch sizes and more frequent glitches.

• A heavier crust produces more large glitches and longer waiting times than does a lighter crust.

• A strong spin-down torque catalyzes smaller, more frequent glitches than a weak torque.

• Pulsar astronomers derive the strength of the pulsar magnetic field from the spin-down rate of the crust without allowing for superfluid feedback. The results in this paper demonstrate that this is inaccurate. Even when no vortices have yet unpinned (e.g. t≤ 430 in the dot–dashed curve in Fig. 7), is a function of both the external torque and the superfluid dynamics (via equation 4). As vortices migrate towards the edge of their pinning sites, the smooth, internal spin-up torque effectively reduces the crust’s spin-down rate. However, this effect is much smaller in a real pulsar than in our simulations. For example, if all vortices move radially outward by the width of a nucleus, ⁠. A larger, smooth spin-up torque results from homogeneous outward vortex creep.

Table 9 summarizes how measures of glitch activity, like the mean size and waiting time and the glitch activity parameter, depend on the underlying physical levers, like the pinning strength/abundance and the electromagnetic torque. A dot indicates where there is no/weak dependence; a plus (minus) sign indicates that the measure of activity increases (decreases) as the physical variable increases. Intuitively, we expect that a positive correlation with mean glitch size is accompanied by a positive correlation with mean waiting time. For V₀, for example, this is not the case; mean glitch size increases with increasing V₀, whereas mean waiting time does not change monotonically with V₀. One possible explanation is that the simulation has not run long enough to establish a critical velocity shear that ensures a net outward flux of vortices that matches the superfluid spin-down to that of the crust (Alpar et al. 1984a; Melatos & Warszawski 2009). The increasing vertical distance between the spin-down curves in Fig. 7 confirms that this critical shear has not yet been established.

If pinning is not strong at the densities where only a neutron superfluid is present (Donati & Pizzochero 2003 show that V₀ depends on density), then glitches may arise from the physics of the outer pulsar core, where the neutron superfluid co-exists with a superconducting proton fluid. Magnetic flux tubes in the superconducting fluid complicate the superfluid vortex dynamics, rendering the simple Magnus-force dynamics in this paper invalid. Hence, simulations presented in this paper are relevant only to the pulsar inner crust, whose density lies in the range 10¹¹≲ρ≲ 10¹⁴ g cm⁻³, where the majority of protons are confined within nuclei (Pethick et al. 2000).

The discrete Ω jumps in the spin-down curves graphed throughout this paper are analogous to the step changes in magnetic field seen in fig. 1 of Field et al. (1995), caused by magnetic flux tube avalanches in type II superconductors. Superconductor experiments exhibit power-law size distributions over two decades, with exponents ranging from −1.4 to −2.2. Several authors have successfully modelled the superconducting system as a self-organized critical system fluctuating about a Bean state, governed by collective unpinning and motion of flux tubes (Bassler & Paczuski 1998). A natural theoretical framework for a self-organized critical system with knock-on mechanisms is a branching process. The characteristic power-law index of such systems is −1.5 (Carreras et al. 2002), consistently steeper than the size distributions reported in this paper. We therefore conclude that either (i) we are not simulating a self-organized critical system or (ii) our definition of event size and/or identification algorithm is not correct.

Although our glitch size distributions are not inconsistent with power laws, we do not find conclusive evidence of collective vortex behaviour, e.g. many-generation cascades of unpinning events in movies of |ψ|². This is not surprising. Warszawski & Melatos (2010b) found that sound waves from moving vortices unpin nearby vortices only when dissipation is low. In this paper’s experiments, for computational expediency, we use γ = 0.05, whereas Warszawski & Melatos (2010b) required γ = 0.0025 to facilitate an acoustic knock-on event given the same characteristic pinning strengths as studied in this paper. In contrast, the proximity knock-on mechanism identified by Warszawski & Melatos (2010b) does occur in the simulations presented here. Hence collective dynamics are not excluded completely.

Collective effects emerge most fully in large systems. We therefore present preliminary results from large, computationally expensive simulations with ∼200 vortices. The large simulations yield a broader range of glitch sizes than smaller systems. Power-law fits to the entire size distribution can be rejected with high confidence. Imposing a minimum cut-off, to create a subset of larger glitches, yields a distribution that can be fitted by a power law that is rejected with less confidence. Physically, this may be equivalent to distinguishing between creep-like vortex motion and avalanche events.

In addition to requiring larger system sizes, reliable pulsar glitch simulations must also address the role of pinning along the length of a vortex. When a three-dimensional vortex moves radially outward, it either ‘unzips’ from pinning sites along its length or moves in small vertical segments that behave effectively as individual vortices. Three dimensions also introduce the possibility of turbulent superflow, i.e. a polarized vortex tangle (Schwarz 1988; Barenghi et al. 2001; Kobayashi & Tsubota 2007; Mongiovì & Jou 2007). Tangled vortices can form bundles, which behave as multiply quantized vortices. A full understanding requires further investigation using 3D simulations.

Finally, as well as expanding GPE simulations to larger systems and higher resolution, we recommend that progress in understanding pulsar glitches can be achieved via cellular automata (Morley & Schmidt 1996; Bassler & Paczuski 1998; Warszawski & Melatos 2008). The generic behaviour observed in GPE simulations can be used to devise microscopic rules for such automata.

We thank Dr Andrew Martin and Dr Natalia Berloff for enlightening discussions on many of the topics addressed in this paper. LW acknowledges the support of an Australian Postgraduate Award.

REFERENCES