Generating quasiperiodic pulsar glitches using a state-dependent Poisson process

TL;DR

This paper demonstrates that a modified state-dependent Poisson process can replicate both Poisson-like and unimodal pulsar glitch behaviors by adjusting a single parameter, offering insights into pulsar microphysics.

Contribution

It introduces a modification to the existing Poisson process model that can generate both classes of pulsar glitch statistics with a single control parameter.

Findings

The model reproduces Poisson-like waiting times and sizes.

Adjusting the control parameter yields unimodal distributions.

Implications for understanding pulsar glitch mechanisms.

Abstract

Glitching pulsars fall broadly into two statistical classes: those with Poisson-like waiting times and power-law sizes, and those with unimodal waiting times and sizes. Previous glitch modeling based on a state-dependent Poisson process readily generates Poisson-like behaviour but struggles to produce unimodal waiting times or sizes. Here it is shown that, when some of the inputs to the model are modified, both classes of statistical behaviour can be reproduced by varying a single control parameter related to the spin-down rate. The implications for past and future glitch observations and the underlying microphysical mechanism are explored briefly.

Tables

Table 1. Table 1: Role of the conditional jump size distribution: functional forms tested and output PDFs (sizes and waiting times). The product of Heaviside functions, H ​ ( Δ ​ X ) ​ H ​ ( X p − Δ ​ X ) 𝐻 Δ 𝑋 𝐻 subscript 𝑋 p Δ 𝑋 H(\Delta X)H(X_{\textrm{p}}-\Delta X) , that imposes the restrictions outlined in Section 2 multiplies all the entries in the third column, along with a proportionality constant determined by normalization. The general trends noted in the final four columns are discussed in Sections 3 and 5 .

Shape	Fixed (F) or stretchable (S)	$\eta \left(\Delta X\mid X_{\text{p}}\right)$	$\text{p}\left(\Delta t\right)$, $\alpha =10$	$\text{p}\left(\Delta X\right)$, $\alpha =10$	$\text{p}\left(\Delta t\right)$, $\alpha =0.1$	$\text{p}\left(\Delta X\right)$, $\alpha =0.1$
Power law	N/A	$\Delta X^{-1.5}H\left(\Delta X-\beta X_{\text{p}}\right)$	Exponential	Power law	Power law	Power law
Uniform	N/A	$1$	Exponential	Monotonic decreasing	Uniform	Uniform
Increasing power law	N/A	$\Delta X^{1.5}$	Exponential	Monotonic decreasing	Monotonic increasing	Monotonic increasing
Exponential	S	$\exp \left(\frac{-{\beta }_{\text{E}}\Delta X}{X_{\text{p}}}\right)$	Exponential	Monotonic decreasing	Exponential	Monotonic decreasing
	F	$\exp \left(-{\beta }_{\text{E}}\Delta X\right)$
Gaussian	S	$\exp \left[\frac{-{\left(\Delta X-{\mu }_{\text{G}}{X}_{\text{p}}\right)}^2}{{\left({\sigma }_{\text{G}}{X}_{\text{p}}\right)}^2}\right]$	Exponential	Monotonic decreasing	Unimodal	Unimodal
	F	$\exp \left[\frac{-{\left(\Delta X-{\mu }_{\text{G}}\right)}^2}{{\left({\sigma }_{\text{G}}\right)}^2}\right]$
Trigonometric	F	${\sin }^2\left(3\pi \Delta X\right)$	Exponential	Monotonic decreasing	Non-monotonic	Unimodal
Narrow horseshoe	F	${\left\{\exp \left[\frac{{\left(\Delta X-{\mu }_{\text{H}}\right)}^2}{{\left({\sigma }_{\text{H}}\right)}^2}\right]-1\right\}}^2$	Exponential	Non-monotonic	Non-monotonic	Non-monotonic

Equations

$$X(t)=X(0)+\frac{N_{\text{em}}t}{I_{\text{c}}}-\sum _{i=1}^{N(t)}\Delta X^{(i)},$$

$$\Delta X^{(i)}=-\frac{2\pi (I_{\text{c}}+I_{\text{s}})\Delta {\nu }^{(i)}}{I_{\text{s}}},$$

$$\eta (\Delta X\mid X_{\text{p}})={\left[{\int }_0^{X_{\text{p}}}d\xi g(\xi ,X_{\text{p}})\right]}^{-1}g(\Delta X,X_{\text{p}}),$$

$$g(\xi ,X_{\text{p}})={\xi }^{-1.5}H(\xi -\beta X_{\text{p}})H(X_{\text{p}}-\xi ).$$

$$\text{p}\left[\Delta t\mid X(t)\right]=\lambda \left[X(t)+\frac{N_{\text{em}}\Delta t}{I_{\text{c}}}\right]\exp \left\{-{\int }_t^{t+\Delta t}d{t}^{\prime }\lambda \left[X(t^{\prime })\right]\right\}.$$

$$\lambda (X)={\lambda }_0{\left(1-X/X_{\text{cr}}\right)}^{-1}.$$

$$\lambda (\tilde{X})=\frac{\alpha }{1-\tilde{X}},$$

$$\alpha =\frac{I_{\text{c}}{X}_{\text{cr}}{\lambda }_0}{N_{\text{em}}}.$$

$$g(\xi ,X_{\text{p}})=f(\xi )H(X_{\text{p}}-\xi ),$$

$$\eta (\Delta X\mid X_{\text{p}})\propto \exp \left[\frac{-(\Delta X-{\mu }_{\text{G}}{)}^2}{{\sigma }_{\text{G}}^2}\right]H(X_{\text{p}}-\Delta X)H(\Delta X),$$

$$\eta (\Delta X\mid X_{\text{p}})\propto {\sin }^2(3\pi \Delta X)H(X_{\text{p}}-\Delta X)H(\Delta X).$$

$$\eta (\Delta X\mid X_{\text{p}})\propto \exp \left[\frac{-(\Delta X-{\mu }_{\text{G}}{X}_{\text{p}}{)}^2}{({\sigma }_{\text{G}}{X}_{\text{p}}{)}^2}\right]H(X_{\text{p}}-\Delta X)H(\Delta X).$$

$$\lambda \left[X(t)\right]=\frac{\alpha }{[1-X(t){]}^2}.$$

$$\lambda [X(t)]=\alpha X(t).$$

$$\lambda [X(t);a]=\frac{1}{a}\exp \left[\frac{X(t)-1}{a}\right],$$

Full text

††pubyear: 2018††pagerange: Generating quasiperiodic pulsar glitches using a state-dependent Poisson process–9

1 Introduction

Rotation-powered pulsars are observed to undergo impulsive spin-up events called “glitches” (Melatos et al., 2008; Espinoza et al., 2011), interrupting the monotonic spin down caused by electromagnetic braking (Taylor et al., 2000). While many statistical analysis of populations of glitching pulsars have been carried out (see Morley & Garcia-Pelayo, 1993; Shemar & Lyne, 1996; Lyne et al., 2000, Fuentes et al., 2017; Fuentes et al., 2018; among others), studies of the glitch statistics for individual pulsars have been stymied by small data sets; the five most prolific glitchers have between 20 and 42 glitches recorded up to 2018 May 28. Nevertheless, as these numbers grow, it is becoming possible to disaggregate the data and generate statistically meaningful probability density functions (PDFs) for the glitch sizes and waiting times (Ashton et al., 2017; Howitt et al., 2018). Most pulsars with adequate sample sizes exhibit Poisson-like (i.e. exponentially distributed) waiting times and scale-invariant (i.e. power-law distributed) event sizes. However, some pulsars show signs of quasiperiodicity in their glitch activity. Specifically, PSR J0835$-$4510 and PSR J0537$-$6910 have non-monotonic waiting time distributions (Melatos et al., 2008; Espinoza et al., 2011; Howitt et al., 2018).

Poisson-like and quasiperiodic behaviours are both consistent with a system existing in a state of self-organized criticality, in which stress accumulates under the action of a slow, global driver and is released sporadically by impulsive, stick-slip events (Jensen, 1998; Melatos et al., 2008). Bak et al. (1987) proposed a sandpile automaton as an idealized example of a self-organised critical system. Sandpiles exhibit power-law avalanche sizes and Poisson-like waiting times and fluctuate around a critical slope in a state of marginal stability. Experimental studies of real sandpiles often demonstrate large, system-spanning, quasiperiodic avalanches, and smaller avalanches in between (Bretz et al., 1992; Rosendahl et al., 1993). This distinction between large and small events also arises in seismology. Omori’s Law describes the phenomenon of small earthquakes occurring at an increased rate after large earthquakes (Omori, 1894; Utsu et al., 1995).

The exact physical mechanism that triggers glitches is unknown. Most microphysical models postulate that stress builds up in the system, as electromagnetic braking increases the elastic stress in the crust and/or crust-superfluid differential rotation. The stress is released spasmodically via superfluid vortex avalanches (Anderson & Itoh, 1975; Warszawski & Melatos, 2011; Warszawski et al., 2012), and/or starquakes (Larson & Link, 2002; Middleditch et al., 2006; Negi, 2011; Morley, 2018), among other possibilities; see Haskell & Melatos (2015) for a modern review. To simulate and understand stress-release mechanisms of this kind, Fulgenzi et al. (2017) modelled long-term glitch activity as a mean-field, state-dependent Poisson process, generalising work done by Daly & Porporato (2006) and Wheatland (2008), who developed versions of the model in the context of forest fires and solar flares respectively. This class of meta-model is agnostic about the microphysics; it applies equally to superfluid vortex avalanches and starquakes, for example. It makes microphysics-independent, falsifiable predictions regarding glitch observables, e.g. correlations between sizes and waiting times (Melatos et al., 2018).

Previous analyses of the state-dependent Poisson model have demonstrated that it generates power-law sizes and exponential waiting times naturally, when the spin-down rate is below a critical threshold (Fulgenzi et al., 2017; Melatos et al., 2018). Above the threshold, the model generates identical size and waiting time PDFs, as seen in quasiperiodic glitchers, but the functional form does not match observations; the model outputs power laws, whereas the data reveal Gaussian-like unimodal PDFs. It is therefore interesting to ask: can one modify the ingredients of the state-dependent Poisson model such that it reproduces both the observed Poisson-like and quasiperiodic behaviours by varying a single parameter — the spin-down rate, normalized as discussed by Fulgenzi et al. (2017) — below and above the threshold respectively? Answering this question systematically is the goal of this paper.

The paper is organized as follows. The state-dependent Poisson process model is described in detail in Section 2. In Section 3 we test several forms of the conditional jump size distribution, a key input into the model, in an effort to generate quasiperiodic glitches. In Section 4 we test several forms of the Poisson rate function, another key model input. The observational implications are discussed briefly in Section 5.

2 State-dependent Poisson process

2.1 Equation of motion

Fulgenzi et al. (2017) modelled long-term glitch activity as a mean-field, state-dependent Poisson process (Cox, 1955; Daly & Porporato, 2007; Wheatland, 2008). In this model, the instantaneous glitch rate, $\lambda(t)$, is governed by a single variable $X(t)$. In the vortex avalanche model, $X(t)$ is the spatially averaged lag between the angular velocities of the rigid crust and the superfluid interior. In the crustquake model $X(t)$ measures the build-up of crustal stress or strain. We frame the presentation in terms of vortex avalanches, but emphasize that the meta-model transfers to any stick-slip stress-release process (Melatos et al., 2018).

The angular velocity lag evolves according to

[TABLE]

where $N_{\text{em}}$ is the electromagnetic spin-down torque acting on the crust, $I_{\textrm{c}}$ is the moment of inertia of the crust, $X(0)$ is an arbitrary initial condition, $N(t)$ is the number of glitches that have occurred up to time $t$, and $\Delta X^{(i)}$ is the size of the $i$-th glitch. By conservation of angular momentum we can relate $\Delta X^{(i)}$ to $\Delta\nu^{(i)}$, the (observable) spin frequency gained by the crust at each glitch, viz.

[TABLE]

where $I_{\textrm{s}}$ is the moment of inertia of the superfluid interior.

Both $N(t)$ and $\Delta X^{(i)}$ are random variables. The avalanche sizes $\Delta X^{(i)}$, while random, depend on $X(t)$ before the glitch. $N(t)$ is a Poisson counting process governed by the waiting time distribution discussed in Section 2.3.

2.2 Sizes

The size of a glitch is governed statistically by $\eta(\Delta X\mid X_{\textrm{p}})$, a conditional jump distribution that depends on the lag immediately before the glitch, $X_{\textrm{p}}$. We posit that glitches always decrease the lag but cannot make it negative. The function $\eta(\Delta X\mid X_{\textrm{p}})$ is one of the ingredients we seek to modify in this paper in an effort to generate quasiperiodic behaviour in a natural way.

One choice for $\eta$, suggested by Gross-Pitaevskii simulations111The Gross-Pitaevskii equation describes the evolution of a zero temperature Bose-Einstein condensate, which is often used as an idealized model of a neutron star superfluid. See the recent review by Haskell & Melatos (2015) for details. (Warszawski & Melatos, 2011, 2013) and by analogy with other systems that exhibit signs of self-organised criticality (Jensen, 1998; Aschwanden et al., 2018), is a power law of the form

[TABLE]

where we define for convenience

[TABLE]

Equation (3) is automatically normalised in the domain $0\leq\Delta X\leq X_{\textrm{p}}$ by the pre-factor. The exponent $-1.5$ in $g(\xi,X_{\textrm{p}})$ is typical of avalanche processes and consistent with Gross-Pitaevskii simulations (Jensen, 1998). The factor $\beta<1$ in the Heaviside function $H$ ensures that the integral converges by setting a minimum glitch size as a fraction of $X_{\textrm{p}}$. It is obligatory for the power law in (4) but not for every other functional form. The impact of the functional form on the long-term statistics is explored in Section 3. There is freedom to choose $\eta(\Delta X\mid X_{\textrm{p}})$, as there is currently no possibility of observing it directly in a neutron star.

2.3 Waiting times

The instantaneous glitch rate is modelled as a time-dependent Poisson process, whose rate is a function of the instantaneous lag, i.e. $\lambda[X(t)]$. We assume the rate is a monotonically increasing function of $X(t)$ and diverges at a critical lag, $X_{\textrm{cr}}$, where a glitch is certain to occur. The function $\lambda(X)$, like $\eta(\Delta X\mid X_{\textrm{p}})$, is one of the ingredients we seek to modify in this paper to reproduce quasiperiodic glitch activity.

As $\lambda\left[X(t)\right]$ evolves deterministically between glitches, we can use the standard PDF for a variable-rate Poisson process to pick the waiting times, $\Delta t$, between glitches, conditional on the deterministic evolution of the lag in the system, $X(t)$, after the previous glitch:

[TABLE]

For the analysis in Section 3 we follow Fulgenzi et al. (2017) in choosing the following phenomenological functional form:

[TABLE]

Here $\lambda_{0}$ is a reference rate, e.g. $\lambda_{0}=\lambda(X_{\textrm{cr}}/2)/2$, and $X_{\textrm{cr}}$ is a complex function of the nuclear physics of vortex pinning. Fulgenzi et al. (2017) argued that the exact functional form of (6) does not affect the long-term statistics, as long as $\lambda(X)$ diverges at $X=X_{\textrm{cr}}$ and increases monotonically with $X$. We test this claim further in Section 4 by investigating the interaction between the rate function, the choice of conditional jump size distribution, and other control parameters.

2.4 Dimensionless variables

The model is expressed usefully in dimensionless variables by setting $\tilde{X}=X/X_{\textrm{cr}}$ and $\tilde{t}=N_{\text{em}}t/(X_{\textrm{cr}}I_{\textrm{c}})$. The main control parameter is $\alpha$, which is introduced through the dimensionless version of (6),

[TABLE]

with

[TABLE]

Fulgenzi et al. (2017) showed that the model output is classified into two regimes: $\alpha\gtrsim\alpha_{c}(\beta)\approx\beta^{-1/2}$, called the slow spin-down regime, which generates exponentially distributed waiting times and power-law sizes; and $\alpha\lesssim\alpha_{c}(\beta)$, where the sizes and waiting times are identically distributed as power laws; see Section 4 in Fulgenzi et al. (2017).

2.5 Monte Carlo simulations

The evolution of $X$ can be modelled through a simple five-step Monte Carlo automaton:222Tildes are dropped here and henceforth in this paper for clarity.

Pick a random $\Delta t$ from (5) given the current lag $X$.

• 2.

Update the lag to $X+\Delta t$ to account for the deterministic evolution up to the glitch.

• 3.

Pick a random $\Delta X$ from (3) given the lag just prior to the glitch.

• 4.

Subtract $\Delta X$ from the lag.

• 5.

Repeat from step 1.

Random numbers are picked using a rejection method. The method is useful when the PDF has a finite upper bound. It handles functional forms in (3) and (5) that are not easily integrable or invertible, as required by the standard inverse cumulative algorithm (Press et al., 2007).

3 Jump size distribution $\eta(\Delta X\mid X_{\textrm{p}})$

The meta-model accommodates any physically plausible distribution for the jump sizes; the power law in (4) is not the only possible functional form. In the following analysis we keep the following two fundamental restrictions: the lag is never negative and it always decreases at a glitch, i.e.

1. $\eta(\Delta X\mid X_{\textrm{p}})=0$for $\Delta X>X_{\textrm{p}}$ and

2. $\eta(\Delta X\mid X_{\textrm{p}})=0$for $\Delta X\leq 0$.

3.1 Functional form

With the above restrictions in mind, we posit a suite of alternative jump size distributions, broadly categorized by two traits: “monotonic” or “unimodal”, and “fixed” or “stretchable” (for those that are not scale-invariant). The categorization attempts to provide a rough census of the function space covered by $\eta(\Delta X\mid X_{\textrm{p}})$. For completeness we also test some “quirky” functions, e.g. trigonometric and horseshoe-shaped. The functional forms of all jump size distributions tested in this paper are tabulated in Table 1.

A monotonic distribution is one that decreases with $\Delta X$. In contrast, a unimodal distribution is one with a single, well-defined peak, like a Gaussian. Two monotonic distributions are shown in the top row of panels in Figure 1: the red curve is the power law defined in (4); the blue curve is a stretchable exponential. Both curves “stretch” with $X_{\textrm{p}}$, i.e. their shape is the same, regardless of $X_{\textrm{p}}$.

We classify a distribution as “fixed” if the shape does not shift with $X_{\textrm{p}}$, i.e. if we have

[TABLE]

for some arbitrary function $f$, where the Heaviside function ensures restriction (i) above is maintained. A “stretchable” distribution on the other hand does shift with $X_{\textrm{p}}$, such that its shape is maintained regardless of $X_{\textrm{p}}$. This is seen clearly in the middle row of panels in Figure 1, where the green curves are stretchable, while the purple curves are fixed.

The bottom row of panels Figure 1 shows two “quirky” functional forms: a fixed trigonometric function (orange curves), and a fixed horseshoe-shaped function (grey curves).

Many of the functional forms require additional parameters to specify their shape as specified in Table 1, e.g. $\mu_{\textrm{G}}$ and $\sigma_{\textrm{G}}$ for the Gaussian, and $\beta_{\textrm{E}}$ for the exponential. In general, the specific choice of these parameters does not affect the overall shape of the long-term statistics produced by the model, i.e. the observable PDFs of waiting times, $\textit{p}(\Delta t)$, and sizes, $\textit{p}(\Delta X)$.

3.2 Monotonic versus unimodal

The effect of $\eta(\Delta X,X_{\textrm{p}})$ on $\textit{p}(\Delta t)$ and $\textit{p}(\Delta X)$ is explored in Figures 2, 3, and 4. We emphasise that $\textit{p}(\Delta t)$ and $\textit{p}(\Delta X)$ are not the same as equations (5) and (3) respectively, as the latter are instantaneous PDFs, while the former are generated by Monte Carlo simulations, during which $X(t)$ and $X_{\textrm{p}}$ fluctuate stochastically. Figures 2–4 are all constructed using the rate function (7). The effect of modifying $\lambda(X)$ is studied in Section 4.

Figure 2 reproduces the results of Fulgenzi et al. (2017), when $\eta(\Delta X,X_{\textrm{p}})$ is a power law with exponent $-1.5$ and $\beta=10^{-2}$, as in (4). Exponential waiting times are seen for $\alpha\geq 3$, while $\textit{p}(\Delta t)$ is a power law for $\alpha\leq 1$. This supports the conclusion of Fulgenzi et al. (2017): there exist two regimes of activity, $\alpha\gtrsim\beta^{-1/2}$ (slow spin-down) and $\alpha\lesssim\beta^{-1/2}$ (fast spin-down). The size distribution is a power law over roughly two decades for all values of $\alpha$ tested. The lower cut-off is dictated by $\alpha$ and $\beta$. When $\alpha$ is low, the lag climbs higher on average before a glitch, magnifying the impact of the minimum glitch size, which is specified as a fraction of the lag prior to the glitch.

Figure 3 shows the long-term statistics generated by a unimodal jump size distribution of the form

[TABLE]

where the proportionality constant is set by normalization. This truncated Gaussian is drawn as the purple curves in the middle row of panels in Figure 1. The dimensionless mean, $\mu_{\textrm{G}}$, and scale, $\sigma_{\textrm{G}}$, of the Gaussian are fixed at $0.5$ and $0.125$ respectively. Exponential waiting times are produced for $\alpha=10$. The size PDF is not a power law for any $\alpha$. However it is monotonically decreasing and roughly scale invariant for $\alpha\gtrsim 3$.

As $\alpha$ decreases, $\textit{p}(\Delta t)$ becomes more uniform and ultimately non-monotonic for $\alpha\lesssim 1$. In the fast spin-down (low $\alpha$) regime its shape resembles $\eta(\Delta X\mid X_{\textrm{p}})$ in the limit $X_{\textrm{p}}\rightarrow 1$, because we have $X_{\textrm{p}}\approx 1$ just before every glitch. The waiting time roughly equals $\Delta X$ at the previous glitch, as the system recovers back to $X\approx 1$. The same behaviour is seen in the size PDF, which approaches a Gaussian as $\alpha$ decreases.

Broadly speaking, adjusting $\mu_{\textrm{G}}$ and $\sigma_{\textrm{G}}$ does not change the shape of the results. Instead it shifts the mean and variance of the resultant PDFs. In the low $\alpha$ regime we find $\langle\Delta t\rangle$, $\langle\Delta X\rangle\rightarrow 1$ for $\mu_{\textrm{G}}\rightarrow 1$. Increasing $\sigma_{\textrm{G}}$ such that $\eta(\Delta X\mid X_{\textrm{p}})$ is a broader distribution generates broader $\textit{p}(\Delta t)$ and $\textit{p}(\Delta X)$, as expected.

Figure 4 shows roughly similar results, when an oddly-shaped jump size distribution is chosen, e.g.

[TABLE]

There is no physical motivation behind this choice. Its purpose is to see how the model responds when pushed in nontraditional directions. In the fast spin-down regime the shape of $\eta(\Delta X\mid X_{\textrm{p}})$ is reflected in both the waiting time and size PDF. In the slow spin-down regime exponential waiting times are seen, while the size PDF decreases monotonically, with a sharp upper cut-off at $\Delta X\approx$2\text{\times}{10}^{-1}$$.

3.3 Fixed versus stretchable

Both the Gaussian and trigonometric jump size distributions discussed in Section 3.2 are “fixed”, in the sense that the parameters that describe their shape do not scale with $X_{\textrm{p}}$. What happens if the parameters are stretchable, meaning that the shape of the distribution stretches with $X_{\textrm{p}}$, e.g. the blue and green curves in Figure 1? Figure 5 compares $\textit{p}(\Delta t)$ and $\textit{p}(\Delta X)$ at $\alpha=3$ for the Gaussian defined in (10) and the stretchable alternative

[TABLE]

The parameters $\mu_{\textrm{G}}$ and $\sigma_{\textrm{G}}$ are set to $0.5$ and $0.125$ respectively for both (10) and (12). The waiting time PDFs are similar for the two alternatives. The size PDF for the stretchable distribution is much narrower than for the fixed distribution. Most of the probability mass in the stretchable distribution lies halfway between 0 and $X_{\textrm{p}}$ making small glitches unlikely. The different behaviour between fixed and stretchable jump size distributions is broadly replicated for small and large $\alpha$ and various functional forms of $\eta(\Delta X\mid X_{\textrm{p}})$.

3.4 Summary

Based on Figures 2–5 and additional simulations not plotted here, the general behaviour can be broadly categorized against the criteria in Table 1.

1. In the high-$\alpha$ regime, all of the tested jump size distributions produce exponential waiting times and all have monotonically decreasing glitch size distributions, except for the “quirky” narrow horseshoe.

2. At high $\alpha$, none of the tested jump size distributions produce power-law distributed glitch sizes, except for (4). This includes monotonic $\eta$, such as the stretchable exponential (the blue curves in Figure 1) which qualitatively resembles a power law.

3. At low $\alpha$, we find $\textit{p}(\Delta t)\approx\textit{p}(\Delta X)\approx\eta(\Delta X\mid X_{\textrm{p}}\rightarrow 1)$ for all functional forms.

4 Rate function $\lambda(X)$

In an effort to produce quasiperiodic glitch activity, or a unimodal waiting time PDF, one may contemplate modifying the rate function defined in (7) instead of $\eta(\Delta X\mid X_{\textrm{p}})$. Fulgenzi et al. (2017) found that a rate function of the form $\lambda[X(t)]=\alpha\tan\left[\pi X(t)/2\right]$, which also diverges as $X(t)\rightarrow 1$, does not qualitatively change the long term statistics, when $\eta(\Delta X\mid X_{\textrm{p}})$ is a power law. However, if the divergence at $X\rightarrow 1$ is steeper than in (7), it is reasonable to expect some quasiperiodicity in the long-term statistics, as the system drives itself back to $X\approx 1$ faster.

We test the above hypothesis by proposing a “fast” alternative to (7), i.e.

[TABLE]

The exponent of $-2$ is arbitrary; the results do not change much if the divergence is steeper (i.e. higher exponent). Figure 6 shows that, to leading order, there is no change to the waiting time PDF when (7) is replaced by (13). The result holds for $\alpha=0.1$, $\alpha=10$, Gaussian $\eta$, and power-law $\eta$. The top panel displays $\textit{p}(\Delta t)$ for the power-law $\eta$ defined in (4). The bottom panel displays $\textit{p}(\Delta t)$ for the Gaussian $\eta$ defined in (10). In the top panel, (13) generates a thinner tail than (7), but the general shape is the same. In the bottom panel, (13) produces more dispersion in $\Delta t$ than (7) at low $\alpha$ and again has a thinner tail at high $\alpha$.

The rate function, $\lambda[X(t)]$, is an input into the conditional waiting time PDF, $\textit{p}(\Delta t\mid X)$, defined in (5). For all $X$, $\textit{p}(\Delta t\mid X)$ is a monotonically increasing function of $\Delta t$ for $\alpha<1$ and a monotonically decreasing function of $\Delta t$ for $\alpha>1$. The changeover shifts to $\alpha=2$ for the faster rate function (13), independent of $\eta(\Delta X\mid X_{\textrm{p}})$.

It may be tempting to reverse-engineer quasiperiodic waiting times by choosing a specific rate function that makes $\textit{p}(\Delta t\mid X)$ non-monotonic. One simple way to do this is to relax the requirement that the rate diverges as $X\rightarrow 1$, e.g.

[TABLE]

Equation (14) leads to a skewed Gaussian for $\textit{p}(\Delta t\mid X)$ and an exponential for $\textit{p}(\Delta t)$ for many choices of $\alpha$ and $\eta(\Delta X\mid X_{\textrm{p}})$. This means that (14) does not naturally lead to two distinct classes of glitchers: ones with exponential waiting times and power law sizes, and those with unimodal waiting times and sizes. This is not to say that reverse-engineering quasiperiodicity using the rate function is impossible. A rate function that is sharply peaked at a certain lag, e.g. $\lambda[X(t)]=\delta[1-X(t)]$, where $\delta$ is the Dirac delta function, always triggers glitches at $X(t)=1$. Hence $\textit{p}(\Delta t)$ and $\textit{p}(\Delta X)$ are the same as $\eta(\Delta X\mid X_{\textrm{p}}\rightarrow 1)$, i.e. they are unimodal if $\eta(\Delta X\mid X_{\textrm{p}}\rightarrow 1)$ is unimodal. However the foregoing strategy comes with the drawback that it does not involve $\alpha$, so it removes the elegant possibility that the state-dependent model encompasses Poisson and power-law glitch activity in a single framework just by varying $\alpha$.

An alternative approach is to consider a family of functions, $\lambda(X;a)$ which tend to $\delta(X-1)$ as $a\rightarrow 0$, e.g.

[TABLE]

where $a$ is a constant, analogous to the dimensionless control parameter $\alpha$. Using (15), we find $\textit{p}(\Delta t)=\textit{p}(\Delta X)=\eta(\Delta X\mid X_{\textrm{p}})$ for values of $a\lesssim 0.05$, as expected. However, for high values of $a$, $\textit{p}(\Delta t)$ becomes uniform if $\eta(\Delta X\mid X_{\textrm{p}})$ is unimodal, as in (10). Thus, like (14), this rate function does not produce both exponential and quasiperiodic waiting times easily by varying a single control parameter.

5 Observational implications

5.1 Waiting time and size PDFs

A recent nonparametric analysis of the waiting time and size PDFs for the five most active glitchers confirms that there are two general classes of activity: exponential waiting times with monotonically decreasing glitch sizes (PSR J0534$+$2200 and PSR J1740$-$3015), and quasiperiodic (i.e. unimodal) waiting times with roughly Gaussian glitch sizes (PSR J0537$-$6910 and PSRJ0835$-$4510) (Howitt et al., 2018). One pulsar (PSR J1341$-$6220) falls somewhere between these two extremes. Due to the relatively small samples ($N\leq 42$), there is not much information in the nonparametric PDFs beyond this broad dichotomy. In particular we cannot say much about the exact functional form. For example, is a monotonically decreasing PDF truly a scale-invariant power law, or is it the tail of a broad PDF with a scale, e.g. a broad Gaussian? PSR J0534$+$2200 is a convincing power law over 2 dex, but other pulsars are less clear-cut.

Table 1 is reductive; it does not encapsulate all of the information encoded in the long term statistics for all possible $\eta(\Delta X\mid X_{\textrm{p}})$. Nonetheless it reveals some broad trends which can be compared to observed nonparametric PDFs. In particular, the Gaussian $\eta(\Delta X\mid X_{\textrm{p}})$ produces exponential waiting times and monotonically decreasing glitch sizes for high $\alpha$, and unimodal waiting times and glitch sizes for low $\alpha$. If the message of the analysis by Howitt et al. (2018) is that size PDFs for the “Poisson-like” pulsars are strict power laws and not just broad, monotonically decreasing functions, the state-dependent Poisson process struggles to produce both power-law and unimodal glitch sizes just by changing the control parameter $\alpha$. On the other hand, if the samples are too small to be definitive about the exact functional form of $\textit{p}(\Delta X)$, then the Gaussian $\eta(\Delta X\mid X_{\textrm{p}})$ does a fair job of capturing the two types of behaviour in the low-$\alpha$ and high-$\alpha$ regimes respectively.

The aim of this paper is limited to exploring the behaviour of the theory under various input assumptions. We do not seek to fit the data for individual objects, partly because the samples available are small, and partly because the theory is idealized. However, to illustrate the points above, Figure 7 shows $\textit{p}(\Delta t)$ and $\textit{p}(\Delta\nu/\nu)$, where $\Delta\nu/\nu$ is the fractional glitch size, for PSR J0534$+$2200 (the canonical Poisson glitcher) and PSR J0835$-$4510 (the canonical quasiperiodic glitcher), compared against the theory. The measured PDFs are constructed using the nonparametric kernel density estimator (Howitt et al., 2018). The theoretical PDFs are generated using (10) and (7) and are scaled such that the mean of the generated waiting times (sizes) equals the mean observed waiting time (size). The $\alpha$ value for each object is selected as described in Appendix A. A Kolmogorov-Smirnov two-sample test indicates that the observed and generated waiting times are consistent with being drawn from the same distribution (p-values of 0.94 and 0.76 for PSR J0534$+$2200 and PSR J0835$-$4510 respectively). However the observed and generated sizes are not (p-values of $4\times 10^{-6}$ and $7\times 10^{-2}$). It would be premature to draw any conclusions from this illustrative exercise except to note that the shapes of the measured and theoretical PDFs are not entirely dissimilar.

5.2 Cross-correlations

It is interesting to ask whether the size-waiting time cross-correlations can constrain the functional form of $\eta$, in the light of recent analysis of the observed correlation coefficients and their dependence on $\alpha$ (Melatos et al., 2018). We present a preliminary study of this issue in Appendix A while we await larger data sets. The main finding is that a significant backward cross-correlation between sizes and waiting times is expected for most unimodal functions $\eta$. No significant backward cross-correlations are observed to date. If this situation persists, as more data are collected, it would indicate either that we have $\alpha\lesssim 1$ in every pulsar, or that $\eta$ is not a unimodal function, in the context of the state-dependent Poisson process.

5.3 Core and crust

One goal of glitch modelling is to probe the relative roles played by the superfluid in the core and the crust of the neutron star. Recent calculations of entrainment between the neutron superfluid and nuclear lattice combined with a phenomenological two-fluid model imply that the angular momentum reservoir carried by the superfluid in the inner crust is not large enough to explain the “regular” glitches in pulsars such as PSR J0537$-$6910, if the stress reservoir is completely emptied at each glitch (Andersson et al., 2012; Chamel et al., 2012). Two-fluid models have been applied to glitch trigger mechanisms (Andersson et al., 2003; Glampedakis & Andersson, 2009), glitch rise times (Sidery et al., 2010), and measuring pulsar masses using glitch data (Ho et al., 2015).

In the state-dependent Poisson model applied to the superfluid vortex avalanches, the moment-of-inertia ratio of the crust and core, $I_{\textrm{c}}/I_{\textrm{s}}$, enters through the relation between $\Delta X^{(i)}$ and $\Delta\nu^{(i)}$, viz. equation (2). As it is impossible to directly observe $\Delta X^{(i)}$, measuring $I_{\textrm{c}}/I_{\textrm{s}}$ is difficult. Maximum likelihood estimation of $I_{\textrm{c}}/I_{\textrm{s}}$ is possible in principle but it involves scanning over seven parameters at a minimum and falls outside the scope of this paper. The meta-model does not distinguish between superfluid in the inner crust and core; $I_{\textrm{s}}$ is proportional to the total angular momentum in both regions. Interestingly we find that the stress reservoir is not depleted totally under typical conditions in the state-dependent Poisson model (see figure 1 of Fulgenzi et al., 2017), i.e. we find $\Delta X^{(i)}\neq X_{\textrm{p}}$ regardless of the choice of $\eta$ or $\alpha$. In other words, the stress reservoir does not need to empty completely to produce quasiperiodic glitches.

6 Conclusions

Glitching pulsars broadly fall into two statistical classes: those with exponential waiting times and monotonically decreasing sizes, and those with unimodal waiting times and sizes. A microphysics-agnostic meta-model based on a state-dependent Poisson process can generate both quasiperiodic and exponential waiting times through varying the control parameter $\alpha$, when the conditional jump size distribution, $\eta(\Delta X\mid X_{\textrm{p}})$, is unimodal. Likewise, the size PDFs generated by the model capture the broad features of the data. The size PDF is not a strict power law in the large-$\alpha$ regime, but it is unclear whether the relatively small glitch samples available to date absolutely require a power law either. This is a step forward from previous analyses with $\eta(\Delta X\mid X_{\textrm{p}})$ of power-law form, which failed to produce unimodal $\textit{p}(\Delta t)$ and $\textit{p}(\Delta X)$ in the small-$\alpha$ regime.

Somewhat counterintuitively, altering the phenomenological rate function does not seem to be the best way to produce non-monotonic waiting times from the model. With the inputs tested in this paper the model struggles to generate both classes of glitching behaviour, unless the jump size distribution is changed from the power law used in previous work.

Physically, if $\eta(\Delta X\mid X_{\textrm{p}})$ is unimodal, there exists a characteristic size for the stress released at each glitch. There are ways to achieve this microphysically, of course, but matters are complicated by the spatially correlated nature of crustquakes and superfluid vortex avalanches, which leads naturally to scale invariance. Alternatively, one can interpret the results as evidence against one glitch mechanism governing all pulsars, or that the canonical view of glitches as a result of a marginally critical system is incomplete.

7 Acknowledgements

Parts of this research are supported by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav) (project number CE170100004).

Appendix A Size-waiting time cross-correlations

The state-dependent Poisson process provides a framework to predict the correlations between the sizes and waiting times to (since) the next (previous) glitch, henceforth called the forward (backward) cross-correlation. In the fast spin-down regime ($\alpha\lesssim 1$), one obtains significant forward cross-correlations. The stress, $X(t)$, is driven rapidly to $\approx X_{\textrm{cr}}$ before each glitch. If the glitch is large, $X(t)$ takes longer to recover to $X(t)\approx X_{\textrm{cr}}$, than if the glitch is small. There is no significant backward cross-correlation in this regime, because $\eta(\Delta X\mid X_{\textrm{p}}\approx 1)$ is independent of the backward waiting time. On the other hand, in the slow spin-down regime ($\alpha\gtrsim 1$), cross-correlations show more varied behaviour, as they depend on the particular choice of $\eta(\Delta X\mid X_{\textrm{p}})$.

Simulations bear out the trends described above, as we see in Figure 8. At low values of $\alpha$ there is a strong forward cross-correlation and no backward cross-correlation for monotonic and unimodal forms of $\eta(\Delta X\mid X_{\textrm{p}})$. At high values of $\alpha$ neither the power law nor the stretchable exponential show significant forward or backward cross-correlations. Small glitches are favored for both these functional forms, rendering the restriction $\Delta X\leq X_{\textrm{p}}$ largely irrelevant. In contrast, a functional form for $\eta(\Delta\mid X_{\textrm{p}})$ which favours larger glitches (e.g. the Gaussians in the bottom-left and bottom-right panels) produces significant backward cross-correlations for $\alpha\gtrsim 1$. In the slow spin-down regime, we have $X(t)\approx 0$, and the restriction $\Delta X\leq X_{\textrm{p}}$ plays an important role. More glitches return the reservoir close to zero, and a correlation between the size of the glitch and the time since the last glitch (i.e. the available reservoir) emerges.

Figure 8 indicates that we should expect to see some backward cross-correlations in pulsars with a unimodal $\eta$ but no forward cross-correlations in the regime $\alpha\gtrsim 1$. The fact that significant backward cross-correlations have not been observed so far may be a product of the small samples. Alternatively it may indicate that we have $\alpha\lesssim 1$ for all objects, or that $\eta$ is not unimodal.

Precisely inferring parameters, such as $\alpha$, from available glitch data is computationally expensive and falls outside the scope of this work. It is also premature given the small samples available. However we can make some broad statements by matching the output waiting time and size PDFs of the state-dependent Poisson model to observed PDFs. Figure 9 shows the nonparametric kernel density estimates for the waiting time and size PDFs for the five pulsars with the most recorded glitches. Overlaid on each plot are PDFs generated with the state-dependent Poisson model, using the Gaussian $\eta$ defined in (10) and the standard rate function defined in (7), each with a different value of $\alpha$. PDFs are generated for 20 logarithmically spaced $\alpha$ values, ranging from $\alpha=0.01$ to $\alpha=100$. Each set of $N=10^{6}$ glitches is scaled such that the mean waiting time equals that of the observed glitches (and equivalently for sizes). The generated PDF highlighted in orange is the one that has the smallest Kolmogorov-Smirnov two-sample statistic for the waiting time distribution. Of the values tested, the statistic is minimised for $\alpha=15$, $\alpha=100$, $\alpha=0.8$, $\alpha=0.8$ and $\alpha=0.5$ for PSR J0534$+$2200, PSR J1740$-$3015, PSR J1341$-$6220, PSR J0537$-$6910 and PSR J0835$-$4510 respectively.