Size-waiting-time Correlations in Pulsar Glitches

Glitches are impulsive, erratically occurring, spin-up events that interrupt the secular, electromagnetic spin-down of a rotation-powered pulsar. As the number of recorded events rises,³ there is growing evidence that glitching pulsars divide into two classes: Poisson-like glitchers, whose waiting times and sizes are described by exponential and power-law probability density functions (PDFs), respectively; and quasiperiodic glitchers, whose waiting times and sizes are distributed roughly normally around characteristic values (Wong et al. 2001; Melatos et al. 2008; Espinoza et al. 2011; Onuchukwu & Chukwude 2016; Yu & Liu 2017; Howitt et al. 2018). The physical mechanism that triggers glitch activity remains a mystery; see Haskell & Melatos (2015) for a recent review. Broadly speaking, however, it is believed that electromagnetic braking increases the elastic stress and differential rotation in the star, which then relax abruptly via some combination of starquakes and superfluid vortex avalanches, when a threshold is exceeded (Andersson et al. 2003; Middleditch et al. 2006; Glampedakis & Andersson 2009; Chugunov & Horowitz 2010; Warszawski & Melatos 2011).

Intuitively, one expects sizes and waiting times to correlate strongly in a threshold-driven, stress-release process, where "stress" refers to any disequilibrium variable including differential rotation. For example, after a larger glitch, one expects a longer delay until the next glitch, while the stress reservoir is replenished. That is, there should be a strong positive correlation between sizes and forward waiting times. Conversely, after a longer waiting time, one expects a larger glitch, because the stress reservoir is fuller. That is, there should also be a strong positive correlation between sizes and backward waiting times. The above intuition rests implicitly on the assumption that the stress reservoir is mostly emptied by each relaxation event.

In contrast, size-waiting-time correlations are rare in pulsar glitch data. A strong, linear correlation of 6.5 days μHz⁻¹ is observed between sizes and forward waiting times in PSR J0537−6910 (Middleditch et al. 2006; Antonopoulou et al. 2018; Ferdman et al. 2018), which can be exploited to reliably predict the epoch of the next glitch; see the "staircase plot" in Figure 8 in Middleditch et al. (2006). A similar claim has been made regarding PSR J1645−0317, where the slope of the correlation is measured to be 0.38 days pHz⁻¹ (Shabanova 2009). However, the glitches in PSR J1645−0317 rise gradually over ∼1 year and do not belong to the class of impulsive events studied in this paper. Beyond these two examples, there is scant evidence to date for statistically compelling correlations between sizes and forward waiting times (Yuan et al. 2010). Moreover, there is no evidence at all for a statistically significant correlation between sizes and backward waiting times in any object (Yuan et al. 2010; Fulgenzi et al. 2017) nor in microglitches (Onuchukwu & Chukwude 2016). Pulsar glitches are not unique in this regard. The absence of size-waiting-time correlations is mirrored in many threshold-driven, stick-slip, stress-release systems in nature, including sandpiles, earthquakes, solar flares, and flux tube avalanches in type II superconductors (Lu & Hamilton 1991; Field et al. 1995; Jensen 1998; Wheatland 2000; Sornette 2004). In these self-organized critical systems, only a small fraction of the stress reservoir empties at each relaxation event. The process randomly releases historical stress accumulated over an extended period covering many events, so that the stress released by any individual event is sometimes less than and sometimes greater than the stress added since the previous event (Jensen 1998; Melatos et al. 2008). Incomplete reservoir depletion is also inferred in some quasiperiodic glitchers, e.g., PSR J0537−6910 (Antonopoulou et al. 2018).

In this paper, we show quantitatively that sufficiently rapid electromagnetic braking produces a size-waiting-time correlation in certain pulsars, even when the microscopic dynamics of the underlying, self-organized critical process are uncorrelated, e.g., as in superfluid vortex avalanches (Warszawski & Melatos 2008, 2011; Melatos & Warszawski 2009). The paper is structured as follows. In Section 2, we review the statistical evidence for size-waiting-time correlations in the seven pulsars with the largest glitch samples, using the latest data from the Jodrell Bank Centre for Astrophysics and Australia Telescope National Facility catalogs (see footnote 3) (Manchester et al. 2005; Espinoza et al. 2011). In Sections 3 and 4, we interpret the data in terms of a state-dependent Poisson process (Cox 1955; Daly & Porporato 2007; Wheatland 2008; Fulgenzi et al. 2017) and show that there exists a critical spin-down rate, above which a strong correlation emerges between sizes and forward waiting times. The theory is quantitative and predictive and does not depend on the microphysics of the glitch trigger. We close in Section 5 by presenting a ranked list of targets predicted to display strong correlations, as a guide to designing the next generation of glitch monitoring campaigns at radio wavelengths, e.g., with phased arrays like LOFAR (Kramer & Stappers 2010), UTMOST (Caleb et al. 2016), and the Square Kilometer Array, as well as at other wavelengths, e.g., gamma-rays (Ray et al. 2011; Clark et al. 2017).

We emphasize that the state-dependent Poisson process analyzed here and by Fulgenzi et al. (2017) is a meta-model that is agnostic about the glitch microphysics. It applies to any threshold-based trigger mechanism, where the star is driven slowly away from equilibrium by electromagnetic spin-down and releases the cumulative stress impulsively, e.g., via superfluid vortex avalanches (differential rotation) or starquakes (elastic stresses). In this paper, we extend the framework in Fulgenzi et al. (2017) by developing size-waiting-time correlations as a new, quantitative, observational test of the model. Specifically, we apply the framework to existing glitch catalogs for the first time (Section 2), calculate correlation coefficients theoretically as functions of the spin-down rate and other variables (Section 4, Appendix), present a new recipe for inverting the correlation data to infer nuclear parameters, e.g., pinning strength (Sections 4, 5), and identify specific objects as targets for future correlation studies (Section 5).

Advances in pulsar timing methods, including multibeam surveys and multifrequency ephemerides, have expanded the total number of recorded glitches to 482 at the time of writing, with up to 42 in an individual object (PSR J0537−6910). Table 1 summarizes the size-waiting-time correlations observed in the seven objects with N ≥ 10 impulsive glitches, an arbitrary cutoff. The size of a glitch is defined by s = Δν/ν, where Δν is the instantaneous jump in pulse frequency ν. The forward (backward) waiting time from any given glitch to the next (previous) glitch is denoted by Δt₊ (Δt₋). For each object, the table displays the Pearson coefficients

$$\begin{eqnarray}&&{r}_{\pm }=\displaystyle \frac{\langle s{\rm{\Delta }}{t}_{\pm }\rangle -\langle s\rangle \langle {\rm{\Delta }}{t}_{\pm }\rangle }{{(\langle {s}^{2}\rangle -\langle s{\rangle }^{2})}^{1/2}{(\langle {\rm{\Delta }}{t}_{\pm }^{2}\rangle -\langle {\rm{\Delta }}{t}_{\pm }{\rangle }^{2})}^{1/2}}\end{eqnarray} \tag{ 1 }$$

for the forward (s–Δt₊) and backward (s–Δt₋) correlations (where angular brackets denote an average), the standard errors

$$\begin{eqnarray}&&{\sigma }_{{r}_{\pm }}={\left(\displaystyle \frac{1-{r}_{\pm }^{2}}{N-3}\right)}^{1/2}\end{eqnarray} \tag{ 2 }$$

for the two correlations,⁴ and the epoch T₁ of the first glitch in each sample. For PSR J0534+2200 (Crab) and PSR J0835−4510 (Vela), the correlations are computed for the full historical data set and for subsets starting at Modified Julian Date (MJD) 46000, when nearly continuous, single telescope monitoring began using modern receivers and backends (Lyne et al. 2015).⁵ For PSR J0537−6910, the correlations are computed for the set of events that appear in at least two out of the three latest analyses (Middleditch et al. 2006; Antonopoulou et al. 2018; Ferdman et al. 2018). Several of the tabulated objects are likely to have experienced unpublished glitches in recent times; for example, PSR J0631+1036 glitched 15 times from MJD 50186 to MJD 55702, yet no glitches have been published since then. The data are plotted on a log–log scale in the form s versus Δt₊ and Δt₋ in Figures 1(a) and (b), respectively, for the seven largest samples.

Zoom In Zoom Out Reset image size

Table 1.  Size-waiting-time Correlations and Standard Errors for Actively Glitching Pulsars with N ≥ 10

PSR J	N	T1 (MJD)	r+	${\sigma }_{{r}_{+}}$	r−	${\sigma }_{{r}_{-}}$	ρ+	${\sigma }_{{\rho }_{+}}$	ρ−	${\sigma }_{{\rho }_{-}}$
0534+2200	27	40493	−0.075	0.204	0.328	0.193	0.024	0.204	0.464	0.181
	23	46664	−0.113	0.222	0.689	0.162	−0.095	0.223	0.506	0.193
0537−6910	42	51285	0.927	0.060	0.159	0.158	0.931	0.058	0.164	0.158
0631+1036	15	50186	0.701	0.206	−0.091	0.287	0.156	0.285	−0.145	0.286
0835−4510	21	40280	0.407	0.215	0.603	0.188	0.358	0.220	0.398	0.216
	14	46257	0.607	0.240	0.661	0.226	0.484	0.264	0.533	0.255
1341−6220	23	47989	0.293	0.214	−0.082	0.223	0.578	0.182	−0.145	0.221
1740−3015	35	47003	0.298	0.169	−0.065	0.176	0.264	0.171	−0.180	0.174
1801−2304	13	46907	0.764	0.204	−0.024	0.316	0.804	0.188	−0.042	0.316

Download table as:  ASCIITypeset image

The random variables ${r}_{\pm }/{\sigma }_{{r}_{\pm }}$ follow a Student's t-distribution with N − 3 degrees of freedom (see footnote 4) in the limit $N\to \infty$, if the null hypothesis (zero correlation) is true, and the underlying variables are drawn from a bivariate normal distribution. The asymptotic result holds approximately for moderate N, even if the underlying variables are not normally distributed. Thus, as a first pass, we can test the null hypothesis, that s and Δt_± are uncorrelated, by computing the corresponding p-value for the pulsars in Table 1, i.e., the probability that the measured $| {r}_{\pm }|$ or greater arises by chance, when the null hypothesis is true. At 99.7% confidence (3σ), using the full historical data set, only one pulsar exhibits a significant s–Δt₊ correlation, namely PSR J0537−6910, and zero pulsars exhibit a significant s–Δt₋ correlation. At 95.4% confidence (two sigma), using the full historical data set, PSR J0631+1036 and PSR J1801−2304 also exhibit significant s–Δt₊ correlations, and PSR J0835−4510 exhibits a significant s–Δt₋ correlation.⁶ Overall, the null hypothesis cannot be excluded for the majority of the objects with N ≥ 10, even though it is tempting to see correlations other than those above when inspecting Figure 1 visually.

It may be argued that the Pearson correlation is not optimal for glitch studies, because (i) s spans up to four decades in individual objects (e.g., PSR J0534+2200), biasing the covariance $\langle s{\rm{\Delta }}{t}_{\pm }\rangle$ unduly toward events with the highest s; and (ii) a nonlinear relation may exist between s and Δt_±, whereas the Pearson correlation tests for a linear relation. For safety, therefore, we also calculate the Spearman rank correlation, ρ_±, which tests for a monotonic relation, whether it is linear or not, and is less sensitive to outliers in the data. The results including standard errors ${\sigma }_{{\rho }_{\pm }}$ are quoted in the last four columns of Table 1. It is clear by inspection that r_± is broadly consistent with ρ_± within the standard errors, except possibly for the forward correlation in PSR J0631+1036, where we find $| {r}_{+}-{\rho }_{+}| =1.1({\sigma }_{{r}_{+}}+{\sigma }_{{\rho }_{+}})$. However, upon rechecking the p-values, we discover that the forward correlation for PSR J1801−2304 strengthens from 2σ to 3σ; the forward correlation for PSR J0631+1036 and the backward correlation for PSR J0835−4510 are no longer significant at 95.4% confidence; and new correlations arguably emerge for PSR J1341−6220 (forward; p-value 4.8 × 10⁻³) (Yuan et al. 2010) and PSR J0534+2200 (backward; p-value 1.7 × 10⁻²). Therefore, in what follows we adopt a conservative approach and only deem correlations to be significant, if they occur at the 3σ level in the full historical data set, i.e., s–Δt₊ for PSR J0537−6910 and PSR J1801−2304. Detailed Monte Carlo simulations with realistic underlying PDFs are deferred to a future paper, when more data become available and warrant a more detailed study.⁷

Quasiperiodic glitch activity is not always accompanied by a strong s–Δt₊ correlation. PSR J0537−6910 does glitch quasiperiodically, but so does PSR J0835−4510, whose s–Δt₊ correlation is weak, with Pearson and Spearman p-values >7.5 × 10⁻² in the full data set and >2.8 × 10⁻² in the truncated data set. This is interesting physically. Quasiperiodic glitches are thought to occur, when the stress reservoir empties almost completely at each event, whereupon a strong s–Δt₊ correlation is expected. Yet, a strong s–Δt₋ correlation is also expected under these circumstances, and no pulsar in Table 1 exhibits it at the 3σ level in the full data set.

It is sometimes argued that the glitch sizes in PSR J0835−4510 are bimodally distributed (Konar & Arjunwadkar 2014; Ashton et al. 2017). One can demonstrate, using kernel density estimator techniques, that the evidence for bimodality in glitch size PDFs is marginal for all the objects in Table 1 (Howitt et al. 2018). Nonetheless, for the sake of completeness, we repeat the analysis⁸ for PSR J0835−4510 after excluding the three smallest glitches, with s ≤ 10⁻⁷, defined consistently as "microglitches" by Palfreyman et al. (2016). We find r₊ = 0.421, r₋ = 0.500, ρ₊ = 0.291, and ρ₋ = 0.264. The latter coefficients are consistent with their counterparts in Table 1; the shifts lie well within the standard errors.⁹ At this juncture, the data yield no conclusive evidence for or against (i) the hypothesis of two independent glitch mechanisms in PSR J0835−4510, or (ii) a link between the s–Δt₊ correlation and quasiperiodicity. We look forward to these matters being clarified in the future, when more data become available.

We now show that the results in Table 1 and Figure 1 are consistent with an idealized yet general model of glitch activity as a state-dependent Poisson process (Fulgenzi et al. 2017). The model is not specific to a particular trigger mechanism. It describes the system in the mean-field approximation in terms of a global, random variable, x(t), which measures the spatially averaged differential rotation (vortex avalanche picture) or elastic stress (starquake picture) throughout the star as a function of time t. As the star spins down, x increases gradually, at a rate proportional to the electromagnetic torque N_em. When a glitch occurs, x drops discontinuously by a random percentage. Thus, x fluctuates around a mean value $\langle x\rangle$ over the long term. Glitch triggering is postulated to be a Poisson process, whose rate function λ(x) (the number of trigger events per unit time) increases monotonically with x. As x changes with time, so does λ(x).

Consider a rate function λ(x) of the form sketched in Figure 2, which diverges in the limit $x\to {x}_{\mathrm{cr}}$, where x_cr is the critical stress. The divergence is compatible with traditional glitch mechanisms (Haskell & Melatos 2015). In the vortex avalanche picture, x_cr is the critical crust-superfluid angular velocity lag, above which the Magnus force exceeds the pinning force throughout the star, and every vortex unpins (Link & Epstein 1991; Fulgenzi et al. 2017).¹⁰ In the starquake picture, x_cr is the critical elastic stress, above which the crustal lattice fails catastrophically (Middleditch et al. 2006; Horowitz & Kadau 2009; Chugunov & Horowitz 2010; Akbal & Alpar 2018). In the fluid instability picture, x_cr is the critical relative velocity between superfluid components, above which two-stream or Kelvin-wave instabilities are excited (Andersson et al. 2003; Mastrano & Melatos 2005; Peralta et al. 2006; Andersson et al. 2007; Glampedakis & Andersson 2009). Note that the rate in Figure 2 is small but nonzero in the limit $x\to 0$, e.g., due to thermal activation (Link & Epstein 1991).

Zoom In Zoom Out Reset image size

It is important to recognize that sizes and waiting times are likely to be uncorrelated in the avalanche microphysics underlying the vortex avalanche and starquake mechanisms. This is a well-known property of any self-organized critical system driven at a constant rate (Jensen 1998). For example, recent quantum mechanical, Gross–Pitaevskii simulations of vortex avalanches in a pinned, decelerating Bose–Einstein condensate show that the size of an avalanche is independent of the crust-superfluid angular velocity lag x immediately before the avalanche (Warszawski & Melatos 2011, 2013; Melatos et al. 2015), except that the avalanche size cannot exceed x, of course. One can approach x_cr closely yet trigger a tiny avalanche; counterintuitively, there is no tendency to trigger larger avalanches closer to the unpinning threshold.

Despite the absence of s–Δt_± correlations at the microscopic level, such correlations do emerge, when uncorrelated avalanches are combined with global spin-down. To see this, consider rapid spin-down first. The system climbs rapidly up the λ(x) curve in Figure 2 and almost reaches x_cr, before a glitch occurs. If s is relatively large, so is $| {\rm{\Delta }}x|$, the absolute value of the stress released by the glitch. Hence, x faces a relatively long climb $\propto | {\rm{\Delta }}x| /{N}_{\mathrm{em}}$ back to x_cr before the next glitch. On the other hand, if s is relatively small, so is $| {\rm{\Delta }}x|$, and the delay $\propto | {\rm{\Delta }}x| /{N}_{\mathrm{em}}$ until the next glitch is relatively short. This translates into a strong correlation between s and Δt₊. Note that the correlation emerges, even though there is zero correlation between $| {\rm{\Delta }}x|$ and the value of x ≈ x_cr just before the glitch. Also note that there is no significant s–Δt₋ correlation; the time taken by x to climb from its post-glitch starting point up to ≈x_cr due to spin-down has nothing to do with $| {\rm{\Delta }}x|$ (and hence s) at the next glitch.

Next consider slow spin-down. Now, the system does not reach x ≈ x_cr before every glitch; the avalanche is triggered at some intermediate value $x\approx \langle x\rangle$, with $0\lt \langle x\rangle \lt {x}_{\mathrm{cr}}$ and $\langle x\rangle \to 0$ as the spin-down rate decreases. Hence, the s–Δt₊ correlation in the previous paragraph almost vanishes. However, a weak s–Δt₋ correlation emerges instead. The physics of the avalanche process is such that the stress variable x cannot be negative, either for vortex avalanches or starquakes. If the waiting time before a glitch is relatively short, then x is relatively small just before the glitch, and so is the size of the avalanche, $| {\rm{\Delta }}x| \leqslant x;$ i.e., $| {\rm{\Delta }}x|$ is "capped," so that x remains positive. Conversely, if the waiting time is relatively long, $| {\rm{\Delta }}x|$ and hence s can be larger while always keeping x positive. This translates into a weak correlation between s and Δt₋.

To test these ideas, we investigate how r_± scales with spin-down rate for the objects in Table 1. Immediately, a question arises: what measure of spin-down rate is it best to use? The obvious candidate is $\dot{\nu }$, but $\dot{\nu }$ is clearly not the whole story; if glitches occur frequently, $\langle x\rangle$ can be much smaller than x_cr, even if $\dot{\nu }$ is large. Another possibility is $\dot{\nu }\langle {\rm{\Delta }}t\rangle /{x}_{\mathrm{cr}}$, which equals the mean stress accumulated between glitches normalized by the critical stress. The latter quantity has the advantages of being dimensionless and equaling the reciprocal of one of the control parameters in the quantitative theory presented in Section 4 (up to a factor of order unity; see Section 4.3 and Equation (6)). It has the disadvantage that x_cr is not observable. We therefore compromise and plot ρ₊ (red symbols) and ρ₋ (blue symbols) versus the dimensional yet observable quantity $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ in Figure 3 for the seven objects in Table 1. We discuss the implications of the compromise carefully in Section 4.3 from a theoretical perspective. The vertical error bars are given by ${\sigma }_{{\rho }_{\pm }}$, while the horizontal error bars are given by the standard error of the mean, ${(N-1)}^{-1/2}{\sigma }_{{\rm{\Delta }}t}$, where ${\sigma }_{{\rm{\Delta }}t}$ is the standard deviation of the measured waiting times.¹¹ The measurement uncertainty in $\dot{\nu }$ is negligible. Note that $\dot{\nu }$ is the long-term, average, spin-down rate after correcting for glitches and timing noise, as quoted in the Australia Telescope National Facility Pulsar Catalogue (Manchester et al. 2005).

Zoom In Zoom Out Reset image size

One result stands out from Figure 3: the strongest s–Δt₊ correlation found in the sample is associated with PSR J0537−6910, which has the second-highest $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ among the plotted objects. This is consistent with the behavior predicted above for a state-dependent Poisson process with λ(x) qualitatively of the form sketched in Figure 2. Moreover, PSR J0537−6910 exhibits no statistically significant s–Δt₋ correlation, which also matches the predicted behavior of a state-dependent Poisson process. Beyond that, the picture is cloudy. PSR J0534+2200 has the largest $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ in the sample by ≈1 dex, yet it exhibits no significant s–Δt₊ correlation (Wong et al. 2001; Espinoza et al. 2014; Shaw et al. 2018), and, if anything, exhibits a 2σ s–Δt₋ correlation according to the Spearman test. PSR J1801−2304 does exhibit a 3σ s–Δt₊ correlation, yet it has the third-lowest $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ in the sample. It is hard to know what to make of these results without dividing $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ by x_cr, but x_cr is unknown and varies in general from pulsar to pulsar. We therefore postpone discussion of the less statistically significant features of Figure 3, until more data become available, and a better understanding of x_cr in specific objects develops.

To prepare for the arrival of more data, we predict the size-waiting-time correlation theoretically in this section. The calculation follows directly from the theory of a state-dependent Poisson process developed by Fulgenzi et al. (2017) for glitches triggered by superfluid vortex avalanches. It applies equally to starquakes for the reasons expressed in Section 3.

4.1. Equations of Motion

In general, the stress variable x(t) obeys a stochastic equation of motion of the form (Fulgenzi et al. 2017)

$$\begin{eqnarray}&&x(t)=x(0)+t-\displaystyle \sum _{i=1}^{N(t)}{\rm{\Delta }}{x}^{(i)},\end{eqnarray} \tag{ 3 }$$

where x(0) is an astrophysically irrelevant initial stress, the second term on the right side describes the secular increase in x, as the star spins down, N(t) is the number of glitches having occurred up to time t, and ${\rm{\Delta }}{x}^{(i)}$ is the absolute value of the step decrease in x due to the i-th glitch. Equation (3) is written in dimensionless form, with x and ${\rm{\Delta }}{x}^{(i)}$ expressed in units of x_cr, and t expressed in units of ${x}_{\mathrm{cr}}{I}_{{\rm{c}}}/{N}_{\mathrm{em}}$, where I_c is the moment of inertia of the stellar crust; see Section 3.4 of Fulgenzi et al. (2017) for details.¹²

Random processes like (3) are called doubly stochastic (Cox 1955), because both N(t) and ${\rm{\Delta }}{x}^{(i)}$ are random variables. In between two glitches, in the interval ${t}_{{\rm{g}}}\leqslant t^{\prime} \leqslant {t}_{{\rm{g}}}+{\rm{\Delta }}t$, the dimensionless stress evolves deterministically according to $x(t^{\prime} )=x({t}_{{\rm{g}}}^{+})+t^{\prime} -{t}_{{\rm{g}}}$ (i.e., spin-down), where $x({t}_{{\rm{g}}}^{+})$ denotes the stress immediately after the first glitch. The PDF of the waiting time Δt obeys the classic formula for a time-dependent Poisson process, viz.

$$\begin{eqnarray}&&p[{\rm{\Delta }}t| x({t}_{{\rm{g}}}^{+})]=\lambda [x({t}_{{\rm{g}}}^{+})+{\rm{\Delta }}t]\exp \left\{-{\int }_{{t}_{{\rm{g}}}}^{{t}_{{\rm{g}}}+{\rm{\Delta }}t}{dt}^{\prime} \,\lambda [x(t^{\prime} )]\right\}.\end{eqnarray} \tag{ 4 }$$

Following Fulgenzi et al. (2017), we work with the rate function

$$\begin{eqnarray}&&\lambda (x)=\displaystyle \frac{\alpha }{1-x},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{I}_{{\rm{c}}}{x}_{\mathrm{cr}}{\lambda }_{0}}{{N}_{\mathrm{em}}}\end{eqnarray} \tag{ 6 }$$

is a dimensionless control parameter proportional to the microscopic avalanche trigger (e.g., vortex unpinning) rate $2{\lambda }_{0}$ at the reference stress $x={x}_{\mathrm{cr}}/2$.¹³ Equation (6) embodies the properties discussed in Section 3 and sketched in Figure 2. Its specific, hyperbolic functional form is arbitrary; the results do not depend sensitively on it, e.g., $\lambda (x)=2\alpha \tan (\pi x/2)$ works just as well (Fulgenzi et al. 2017). The PDF of the jump sizes ${\rm{\Delta }}{x}^{(i)}\gt 0$ is given by the conditional jump probability

$$\begin{eqnarray}&&\eta (x| y)\propto {(y-x)}^{-1.5}H(y-x-\beta y),\end{eqnarray} \tag{ 7 }$$

where $\eta (x| y){dx}$ equals the probability of jumping from y to a stress value in the interval $(x,x+{dx})$, with $y-{\rm{\Delta }}{x}^{(i)}=x$ at the i-th glitch. Every glitch reduces the stress, the Heaviside function $H(...)$ in (7) ensures that no glitch makes x negative, and the minimum stress release is $\beta y$ ($0\lt \beta \lt 1$) (required for normalization). The power-law form and exponent of (7) are chosen to be consistent with the avalanche size PDFs seen universally in self-organized critical systems like sandpiles, earthquakes, and solar flares (Jensen 1998; Sornette 2004; Aschwanden et al. 2016) and specifically in Gross–Pitaevskii simulations of superfluid vortex avalanches in the neutron star context (Warszawski & Melatos 2011; Melatos et al. 2015). Monte Carlo simulations confirm that, as with λ(x), the output of the model does not depend sensitively on the specific functional form of $\eta (x| y)$ (Fulgenzi et al. 2017). There is no way at present of measuring $\eta (x| y)$ observationally or deriving it theoretically from first principles. A new generation of Gross–Pitaevskii simulations containing many more vortices than have been analyzed to date would be required, a challenging computational task.

4.2. Critical Spin-down Rate

The behavior of the model (3)–(7) was studied thoroughly as a function of the control parameters α and β by Fulgenzi et al. (2017) using Monte Carlo simulations and analytic theory. The behavior divides into two distinct regimes: large α ≳ α_c(β) (slow spin-down) and small α ≲ α_c(β) (fast spin-down), with

$$\begin{eqnarray}&&{\alpha }_{{\rm{c}}}\approx {\beta }^{-1/2}.\end{eqnarray} \tag{ 8 }$$

In the large-α regime, the simulations produce power-law and exponential PDFs for s and Δt, respectively (see Figures 6 and 8, respectively in Fulgenzi et al. 2017), consistent with observations of many pulsars (Melatos et al. 2008; Espinoza et al. 2011; Ashton et al. 2017; Howitt et al. 2018). In the small-α regime, p(s) and $p({\rm{\Delta }}t)$ have nearly the same functional form, i.e., $p(s)\approx p({\rm{\Delta }}t)$ in terms of dimensionless variables, which is consistent with observations of quasiperiodic objects, except that the functional form is a power law instead of a Gaussian for the specific jump distribution (7).

4.3. s–Δt_± Correlations

Just as p(s) and $p({\rm{\Delta }}t)$ change character at $\alpha \approx {\alpha }_{{\rm{c}}}(\beta )$, so do the s–Δt_± correlations. Figure 4 displays r_± versus α for the jump distribution (7). The plot spans the full range from small to large α, with β = 10⁻² and hence α_c ≈ 10 in this example. The behavior exactly matches what is predicted by the qualitative discussion in Section 3. A strong forward correlation emerges when the spin-down rate is fast because the stress approaches x ≈ x_cr before every glitch. A weak backward correlation emerges when the spin-down rate is slow because the size of a glitch is capped to ensure x ≥ 0. The transition occurs at α ≈ α_c in Figure 4.

Zoom In Zoom Out Reset image size

Do the measured values of r_± agree with the theoretical prediction in Figure 4? To answer this question, we need to know α for the objects in Table 1. Unfortunately, Equation (6) expresses α in terms of the quantities I_c, x_cr, λ₀, and N_em, none of which can be measured directly. We can write $2\pi \dot{\nu }\approx {N}_{\mathrm{em}}/{I}_{{\rm{c}}}$ to a good approximation, because the electromagnetic braking torque dominates the superfluid back-reaction torque on the crust (Espinoza et al. 2011). However, x_cr and ${\lambda }_{0}$ cannot be easily related to non-glitch observables. We therefore turn to the theoretical analysis in the Appendix for guidance. Although it applies to the special case where $\eta (x| y)$ is separable, nevertheless it turns out to offer useful clues. From (24), we find that α can be related to the observable mean waiting time, $\langle {\rm{\Delta }}t\rangle$, via $\alpha \approx {x}_{\mathrm{cr}}/(2\pi \langle {\rm{\Delta }}t\rangle \dot{\nu })$ up to a proportionality factor of order unity, where we now restore the dimensions to $\langle {\rm{\Delta }}t\rangle$. Clearly λ₀ drops out of the expression, leaving x_cr. Suppose we then make the assumption that x_cr does not vary much from one pulsar to the next because it is set by the balance of the Magnus and pinning forces (vortex avalanche picture) or crustal breaking strain (starquake picture), which are nuclear in origin and independent of the rotational state (ν, $\dot{\nu }$). Then, α is inversely proportional to the observable product $\dot{\nu }\langle {\rm{\Delta }}t\rangle$ , and it is possible to use this product to compare r_± across different pulsars. This motivates the choice of normalization of the abscissae in Figure 3, as foreshadowed in Section 3.

The crude first success of Figure 3—that the object with the highest r₊ also happens to have the second-highest value of $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$, in line with the theory—is an encouraging sign that a stick-slip process described by (3) may be at work. However, it is nothing more than a first indication; much more data are needed before we can say anything definite. Certainly, the assumption in the previous paragraph, that x_cr does not vary much from one pulsar to the next, is unlikely to hold exactly. Variation in x_cr between objects is one natural way to explain why PSR J0534+2200 fails to exhibit a strong s–Δt₊ correlation, despite having the highest $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ in Table 1; x_cr may be larger than average in this pulsar. Likewise, PSR J1801−2304 does exhibit a strong s–Δt₊ correlation, even though it has the third-lowest $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ in Table 1; x_cr may be smaller than average in this pulsar.

The reader might wonder whether some other observables, e.g., $\langle s\rangle$ or $\langle {({\rm{\Delta }}t)}^{2}\rangle$, depend on α and x_cr in a different combination, allowing us to disentangle the values of α and x_cr. Unfortunately, the prospects are dim. One finds from (24) and (39) that the dimensionless ratio $\langle {({\rm{\Delta }}t)}^{2}\rangle /\langle {\rm{\Delta }}t{\rangle }^{2}$ depends primarily on α, but the dependence is weak, and α is poorly constrained given the measurement uncertainties.¹⁴ Likewise, long-term conservation of angular momentum implies $\langle s\rangle /\langle {\rm{\Delta }}t\rangle =1$ up to a factor involving the crust and superfluid moments of inertia, and again x_cr cannot be disentangled; we have $\langle s\rangle \propto {x}_{\mathrm{cr}}$ and $\langle {\rm{\Delta }}t\rangle \propto {x}_{\mathrm{cr}}$, and hence x_cr cancels out in the ratio.

We conclude using the results in Sections 3 and 4 to predict what pulsars are likely to display strong size-waiting-time correlations in the future, when more data become available.

In Figure 5 we present the cumulative distribution function of $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ for all pulsars known to glitch at the time of writing with N ≥ 4, so that ${\sigma }_{{r}_{\pm }}$ is well-defined. In Table 2 we name the pulsars with the five highest and five lowest $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ values. From the results in Sections 3 and 4, we venture to make two predictions. First, if r₊ is measured to be high in a pulsar, then that particular object is likely to lie toward the top end of the $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ distribution, depending on its x_cr value. By and large, therefore, the objects in the top half of Table 2 represent good r₊ targets (except for PSR J0534+2200; see Table 1). Second, we predict that no glitching pulsar will exhibit a strong s–Δt₋ correlation, either now or in the future. Equation (40) implies ${r}_{-}\leqslant 0.5$ for separable $\eta (x| y)$ and ${r}_{-}\lesssim 0.1$ for typical parameters. Among the low r₋ measurements, we predict that the highest will lie toward the bottom end of the $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ distribution, again depending on x_cr. By and large, the objects in the bottom half of Table 2 represent good r₋ targets. For the sake of completeness, we quote r_± in the last two columns of the table, as computed from existing data. However, we urge the reader not to draw any conclusions at this stage about objects other than PSR J0534+2200 and PSR J0537−6910 in Table 2; the samples are simply too small (4 ≤ N ≤ 7) to say anything with confidence.

Zoom In Zoom Out Reset image size

Table 2.  Glitching Pulsars with the Highest (Top) and Lowest (Bottom) Values of $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ and N ≥ 4 Proposed as Targets for Future Correlation Analyses, When More Data Are Gathered

PSR J	N	$-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ (Hz)	r+	r−
0205+6449	6	3.32 × 10−3	0.947	−0.447
0534+2200	27	2.23 × 10−2	−0.075	0.328
0537−6910	42	1.90 × 10−3	0.927	0.159
1119−6127	4	4.31 × 10−3	0.900	0.860
2229+6114	6	1.68 × 10−3	0.874	−0.305
0528+2200	4	9.81 × 10−7	−0.645	0.768
1814−1744	7	2.89 × 10−6	0.042	0.222
1902+0615	6	1.99 × 10−6	0.490	−0.314
1957+2831	4	4.18 × 10−6	0.667	0.613
2225+6535	5	5.01 × 10−6	0.998	−0.325

Note. The current values of the Pearson coefficients r_± are tabulated for completeness, but most of the samples are too small (4 ≤ N ≤ 7) to draw any statistically significant conclusions.

Download table as:  ASCIITypeset image

We emphasize that the above predictions implicitly assume that x_cr does not vary much from one object to the next (see Section 4), so that α and $\dot{\nu }\langle {\rm{\Delta }}t\rangle$ can be used interchangeably. This seems unlikely, when one considers the nuclear physics of the crust, and may well explain the existing misfits PSR J0534+2200 (high $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$, low r₊) and PSR J1801−2304 (low $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$, high r₊). The predictions also assume that β and hence α_c ≈ β^−1/2 do not vary much between pulsars, which is an open question physically. Therefore, the predictions should be seen as a first step toward falsifiable tests of the correlation mechanism, to be refined as our understanding of x_cr and β in specific objects improves.

As more data become available, it will be possible in principle to turn around the above predictions and use measured correlations (or their absence) to constrain x_cr and the stress-release physics in glitches. Consider r₊, for example. Every glitching pulsar that is measured to have ${r}_{+}\ll 1$ also has α ≳ α_c and hence ${\beta }^{1/2}{x}_{\mathrm{cr}}\gtrsim -2\pi \dot{\nu }\langle {\rm{\Delta }}t\rangle$ in the theory in Section 4, upon relating α to x_cr as before and using (8). To illustrate what is possible in the future, we note that we obtain minimum values of ${\beta }^{1/2}{x}_{\mathrm{cr}}/(2\pi )$ between 2.3 × 10⁻² and 3.4 × 10⁻⁵ Hz for the objects in Table 1 with ${r}_{+}\ll 1$ and between 2.3 × 10⁻² and 9.8 × 10⁻⁷ Hz for the objects in Table 2. These bounds are consistent with sensible values of x_cr and β in the vortex avalanche picture (Link & Epstein 1991; Warszawski & Melatos 2011), e.g.,

$$\begin{eqnarray}{\beta }^{1/2}{x}_{\mathrm{cr}}/(2\pi ) & = & 1\times {10}^{-3}{\left(\displaystyle \frac{\beta }{{10}^{-2}}\right)}^{1/2}\left(\displaystyle \frac{{F}_{\max }}{\mathrm{keV}\,{\mathrm{fm}}^{-1}}\right)\\ & & \times \,{\left(\displaystyle \frac{\rho }{{10}^{13}{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\displaystyle \frac{l}{{10}^{2}\mathrm{fm}}\right)}^{-1}\,\mathrm{Hz},\end{eqnarray} \tag{ 9 }$$

where F_max is the maximum pinning force per site, ρ is the superfluid density, and l is the pinning site separation. An analogous expression in the starquake picture can be deduced from the models in Middleditch et al. (2006) and Akbal & Alpar (2018). In the vortex avalanche picture, especially, a lot of complicated physics goes into F_max, including the form of the nuclear pinning potential, vortex tension, single- versus multi-site breakaway, and collective avalanche knock-on; see Haskell & Melatos (2015) and references therein. One therefore expects F_max to vary from one pulsar to the next.

In this paper, we quantify systematically the size-waiting-time correlations observed in pulsar glitches using the Pearson and Spearman coefficients. We find that, at the 3σ level, no objects exhibit a significant s–Δt₋ correlation, and only two, PSR J0537−6910 and PSR J1801−2304, exhibit significant s–Δt₊ correlations. We show that these results can be understood theoretically in terms of a state-dependent Poisson process, whose rate diverges when the system stress approaches a critical threshold x_cr in both the vortex avalanche and starquake pictures. The state-dependent Poisson process predicts a strong s–Δt₊ correlation (${r}_{+}\approx 1$) for fast spin-down, i.e., for $-\dot{\nu }\langle {\rm{\Delta }}t\rangle /{x}_{\mathrm{cr}}$ greater than a critical value related to the minimum avalanche size. It also predicts a weak s–Δt₋ correlation (${r}_{-}\approx 0$) for fast and slow spin-down. Applying the theory to the list of known, glitching pulsars with N ≥ 4, ranked by $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$, we identify the objects that are likely to display strong s–Δt₊ correlations (and weak or nonexistent s–Δt₋ correlations), as more data are collected. The prediction relies to some extent on assuming that x_cr and the minimum avalanche size, which are unobservable, do not vary much from one pulsar to the next. If future data are in accord with this assumption, measurements of r_± versus $-\dot{\nu }\langle {\rm{\Delta }}t\rangle$ can be turned around to constrain x_cr and hence the nuclear pinning forces (vortex avalanche picture) or crustal breaking strain (starquake picture) in individual pulsars.

The results in this paper extend the theoretical framework developed by Fulgenzi et al. (2017) by focusing on size-waiting-time correlations as a quantitative observational test of the model. The new elements include: (i) a systematic, multi-object analysis of the Pearson and Spearman coefficients derived from data in the Jodrell Bank and Australia Telescope glitch catalogs (Section 2); (ii) intuitive explanations for the strong s–Δt₊ and weak s–Δt₋ correlations expected in a state-dependent Poisson process (Section 3); (iii) closed-form integral expressions for r_± (Appendix A.2); (iv) a recipe for relating the correlation data to essential nuclear physics parameters, e.g., maximum pinning force (Section 4.3 and Equation (9)); and (v) predictions for what specific pulsars are most likely to exhibit emerging s–Δt_± correlations, as more observations are made.

In closing, we again emphasize that the theoretical framework is not specific to a particular version of the glitch microphysics. The state-dependent Poisson process is a meta-model that encompasses all the glitch mechanisms contemplated in the literature to date, e.g., starquakes and superfluid vortex avalanches. It rests on two assumptions of a general nature: (i) the stress x increases gradually between glitches and relaxes discontinuously at a glitch; and (ii) the trigger rate λ(x) increases with x and diverges at x_cr. If the meta-model is falsified in the future, with the arrival of more data and a better understanding of x_cr in specific objects, a fresh approach to the glitch problem will be required.

In order to take full advantage of the opportunity for falsification, more glitches need to be found. Improved data analysis techniques will play an important role in this regard. Recent innovations include algorithms that harness the power of distributed volunteer computing (Clark et al. 2017), alternatives to least-squares fitting for nongaussian noise (Wang et al. 2017), and Bayesian model selection (Shannon et al. 2016).

The authors thank Julian Carlin for assistance with the preparation of Figures 2 and 4. This research was supported by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), grant number CE170100004. G.H. is the recipient of an Australian Postgraduate Award. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY–1607611.

In this appendix, we summarize certain useful results from an analytic theory developed by Fulgenzi et al. (2017) to predict the long-term glitch statistics generated by (3)–(7). The aims are to justify the theoretical relations between observables (e.g., r_± and $\dot{\nu }\langle {\rm{\Delta }}t\rangle$) discussed in Section 3 onward and motivate the axis choices made in Figure 3 onward.

For $t\geqslant 1-x(0)$, the system (3)–(7) exhibits stationary behavior: x(t) fluctuates about a constant mean, $0\lt \langle x\rangle \lt 1$, governed by the balance between the second and third terms on the right side of (3). The system is self-regulating, because λ(x), which determines N(t), increases monotonically with x; as the stress rises, glitches occur more frequently and relax the system. Under stationary conditions, the PDF p(x) of the stress variable x satisfies the time-independent master equation (Warszawski & Melatos 2013; Fulgenzi et al. 2017),

$$\begin{eqnarray}&&0=-\displaystyle \frac{{dp}(x)}{{dx}}-\lambda (x)p(x)+{\int }_{x}^{1}{dy}\,p(y)\lambda (y)\eta (x| y).\end{eqnarray} \tag{ 10 }$$

Equations (10) and (3) describe exactly the same dynamics and are expressed in terms of the same dimensionless variables. The first two terms on the right side of (10) describe the probability lost from the interval $(x,x+{dx})$ due to secular spin-down and discontinuous jumps (glitches) out of the interval, respectively. The third term describes the integrated probability gained in the interval $(x,x+{dx})$, when glitches take the system from another state y into $(x,x+{dx})$. Once p(x) is known after solving (10), it is possible to calculate the statistical distributions of other system variables, including observables like s and Δt_±.

Equations (10) and (5)–(7) form a closed system, which can be solved by the methods developed by Fulgenzi et al. (2017). Monte Carlo simulations confirm that the solution is insensitive to the particular choices of λ(x) and $\eta (x| y)$, as the latter reference demonstrates. If $\eta (x| y)$ is separable, the theory can even be solved analytically. In this appendix, we present the analytic solution for

$$\begin{eqnarray}&&\eta (x| y)=(\delta +1){x}^{\delta }{y}^{-(\delta +1)}.\end{eqnarray} \tag{ 11 }$$

This choice is illustrative only; the vortex or starquake avalanche dynamics inside a neutron star cannot be measured experimentally at present. However it is consistent with the output of Gross–Pitaevskii simulations, viz. Equation (7) (Warszawski & Melatos 2011), and correctly favors small avalanches over large ones for δ > 0, with δ ≈ 3 yielding event statistics broadly in accord with those generated by (7). It also leads to generic scalings between observables, which are reproduced by other sensible choices of $\eta (x| y)$ too, as confirmed by Monte Carlo simulations with nonseparable $\eta (x| y)$ performed by Fulgenzi et al. (2017).

A.1. Stress, Size, and Waiting-time PDFs

Solving (10) and (11) by separation of variables, as in Appendices C and D in Fulgenzi et al. (2017), we find

$$\begin{eqnarray}&&p(x)={{Cx}}^{\delta +1}{(1-x)}^{\alpha },\end{eqnarray} \tag{ 12 }$$

with

$$\begin{eqnarray}&&C=\displaystyle \frac{{\rm{\Gamma }}(\alpha +\delta +3)}{{\rm{\Gamma }}(\alpha +1){\rm{\Gamma }}(\delta +2)},\end{eqnarray} \tag{ 13 }$$

where ${\rm{\Gamma }}(...)$ symbolizes the gamma function. Equation (12) implies $0\lt \langle x\rangle ={(\delta +2)(\alpha +\delta +3)}^{-1}\lt 1$. We can also calculate the PDFs of x immediately before and after a glitch, called ${p}_{{\rm{e}}}(x)$ and ${p}_{{\rm{s}}}(x)$, respectively, by Fulgenzi et al. (2017) and given by (see Equations (B2) and (B3) of the latter reference)¹⁵

$$\begin{eqnarray}&&{p}_{{\rm{e}}}(x)=\displaystyle \frac{\lambda (x)p(x)}{\langle \lambda \rangle }\,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&=\,C\langle \lambda {\rangle }^{-1}\alpha {x}^{\delta +1}{(1-x)}^{\alpha -1}\,\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&{p}_{{\rm{s}}}(x)=\displaystyle \frac{1}{\langle \lambda \rangle }{\int }_{x}^{1}{dy}\,\lambda (y)\eta (x| y)p(y)\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&=\,C\langle \lambda {\rangle }^{-1}(\delta +1){x}^{\delta }{(1-x)}^{\alpha }\,\end{eqnarray} \tag{ 17 }$$

with

$$\begin{eqnarray}&&\langle \lambda \rangle ={\int }_{0}^{1}{dx}\,\lambda (x)p(x)\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\,=\,\alpha +\delta +2.\,\end{eqnarray} \tag{ 19 }$$

The PDFs of the observable waiting times and sizes follow directly from (14)–(19). The waiting time leading up to a glitch is the random value of Δt generated by a Poisson process, whose rate $\lambda [x(t)]$ since the previous glitch evolves deterministically due to spin-down, conditional on the stress immediately after the previous glitch. The size of a glitch is the random value of ${\rm{\Delta }}x=y-x$ generated by $\eta (x| y)$, conditional on the stress y immediately before the glitch. Hence, applying Equations (34) and (35) in Fulgenzi et al. (2017), we obtain

$$\begin{eqnarray}&&p({\rm{\Delta }}t)={\int }_{0}^{1-{\rm{\Delta }}t}{dy}\,{p}_{{\rm{s}}}(y)p({\rm{\Delta }}t| y)\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&=\,(\alpha +\delta +1){(1-{\rm{\Delta }}t)}^{\alpha +\delta },\,\end{eqnarray} \tag{ 21 }$$

with $p({\rm{\Delta }}t| y)$ given by (4), as well as

$$\begin{eqnarray}&&p({\rm{\Delta }}x)={\int }_{{\rm{\Delta }}x}^{1}{dy}\,{p}_{{\rm{e}}}(y)\eta (y-{\rm{\Delta }}x| y)\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&=\,(\alpha +\delta +1){(1-{\rm{\Delta }}x)}^{\alpha +\delta }.\,\end{eqnarray} \tag{ 23 }$$

The PDFs (21) and (23) qualitatively resemble those observed in the pulsars in Table 1 but they do not match the data in detail, because the separable form of $\eta (x| y)$ in (11) represents an approximation. The moments, however, and their scalings with α, are insensitive to the functional form of $\eta (x| y)$. In particular, the first moment of $p({\rm{\Delta }}t)$ evaluates to yield the important result

$$\begin{eqnarray}&&\langle {\rm{\Delta }}t\rangle ={(\alpha +\delta +2)}^{-1},\end{eqnarray} \tag{ 24 }$$

which is used heavily in Section 4; see also Appendix A in Fulgenzi et al. (2017).

A.2. Size-waiting-time Correlations

To calculate the correlation coefficients r_±, we must first evaluate the joint probability of measuring size-waiting-time pairs $({\rm{\Delta }}x,{\rm{\Delta }}{t}_{\pm })$. There are subtleties involved. Consider an arbitrarily selected sequence of three consecutive glitches labeled by G₁, G₂, and G₃. Suppose that G₂ has a size Δx and forward and backward waiting times Δt₊ and Δt₋, respectively. Let y_s be the stress immediately after G₁. Then, deterministic evolution during the interval G₁G₂ implies that the stress immediately before G₂ is ${y}_{{\rm{e}}}={y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-};$ the event G₂ reduces the stress to ${y}_{{\rm{e}}}-{\rm{\Delta }}x$ immediately after G₂; and deterministic evolution during the interval ${G}_{2}{G}_{3}$ implies that the stress immediately before G₃ is ${y}_{{\rm{e}}}-{\rm{\Delta }}x+{\rm{\Delta }}{t}_{+}$. Putting everything together, the probability density of simultaneously measuring Δx and Δt₋ given y_s equals the conditional joint PDF

$$\begin{eqnarray}&&{q}_{-}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})=p({\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})\eta ({y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-}-{\rm{\Delta }}x| {y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\,\propto {(1-{y}_{{\rm{s}}})}^{-\alpha }{(1-{y}_{{\rm{s}}}-{\rm{\Delta }}{t}_{-})}^{\alpha -1}\\ &&\,\,\times \,{({y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-})}^{-(\delta +1)}{({y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-}-{\rm{\Delta }}x)}^{\delta },\end{eqnarray} \tag{ 26 }$$

where $p({\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})$ is given by (4). Likewise, the probability density of simultaneously measuring Δx and Δt₊ given y_e equals the conditional joint PDF

$$\begin{eqnarray}&&{q}_{+}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{+}| {y}_{{\rm{e}}})=p({\rm{\Delta }}{t}_{+}| {y}_{{\rm{e}}}-{\rm{\Delta }}x)\eta ({y}_{{\rm{e}}}-{\rm{\Delta }}x| {y}_{{\rm{e}}}),\,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\,\propto {(1-{y}_{{\rm{e}}}+{\rm{\Delta }}x)}^{-\alpha }{(1-{y}_{{\rm{e}}}+{\rm{\Delta }}x-{\rm{\Delta }}{t}_{+})}^{\alpha -1}\\ &&\quad \times \,{y}_{{\rm{e}}}^{-(\delta +1)}{({y}_{{\rm{e}}}-{\rm{\Delta }}x)}^{\delta }\end{eqnarray} \tag{ 28 }$$

where the first factor on the right side of (27) is given again by (4). The conditional joint PDFs are normalized according to

$$\begin{eqnarray}&&1={\int }_{0}^{1-{y}_{{\rm{s}}}}d({\rm{\Delta }}{t}_{-}){\int }_{0}^{{y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-}}d({\rm{\Delta }}x)\,{q}_{-}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})\end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}&&1={\int }_{0}^{{y}_{{\rm{e}}}}d({\rm{\Delta }}x){\int }_{0}^{1-{y}_{{\rm{e}}}+{\rm{\Delta }}x}d({\rm{\Delta }}{t}_{+})\,{q}_{+}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{+}| {y}_{{\rm{e}}}).\end{eqnarray} \tag{ 30 }$$

The terminals on (29) and (30) ensure that the stress always stays in the domain [0, 1].

The law of total covariance states

$$\begin{eqnarray}\mathrm{cov}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}) & = & E[\mathrm{cov}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})]\\ & & +\,\mathrm{cov}[E({\rm{\Delta }}x| {y}_{{\rm{s}}}),E({\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})],\end{eqnarray} \tag{ 31 }$$

where $E(...)=\int {{dy}}_{{\rm{s}}}\,{p}_{{\rm{s}}}({y}_{{\rm{s}}})\times (...)$ denotes the expectation value when marginalizing over y_s, and $E({\rm{\Delta }}x| {y}_{{\rm{s}}})$ and $E({\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})$ are random variables themselves. An analogous result applies to $\mathrm{cov}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{+})$, except that one marginalizes over y_e. It turns out that the y_s integrals in (31) can be done analytically, viz.

$$\begin{eqnarray}E({\rm{\Delta }}x| {y}_{{\rm{s}}}) & = & {\int }_{0}^{1-{y}_{{\rm{s}}}}d({\rm{\Delta }}{t}_{-})\\ & & \times \,{\int }_{0}^{{y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-}}d({\rm{\Delta }}x)\,{\rm{\Delta }}x\,{q}_{-}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}}),\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1+\alpha {y}_{{\rm{s}}}}{(\alpha +1)(\delta +2)},\,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}E({\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}}) & = & {\int }_{0}^{1-{y}_{{\rm{s}}}}d({\rm{\Delta }}{t}_{-})\\ & & \times {\int }_{0}^{{y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-}}d({\rm{\Delta }}x)\,{\rm{\Delta }}{t}_{-}\,{q}_{-}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}}),\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\,=\,\displaystyle \frac{1-{y}_{{\rm{s}}}}{\alpha +1},\,\end{eqnarray} \tag{ 35 }$$

and

$$\begin{eqnarray}&&\mathrm{cov}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})={\displaystyle \int }_{0}^{1-{y}_{{\rm{s}}}}d({\rm{\Delta }}{t}_{-}){\displaystyle \int }_{0}^{{y}_{{\rm{s}}}+{\rm{\Delta }}{t}_{-}}d({\rm{\Delta }}x)\\ &&\quad \times \,[{\rm{\Delta }}x{\rm{\Delta }}{t}_{-}-E({\rm{\Delta }}x| {y}_{{\rm{s}}})E({\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}})]{q}_{-}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-}| {y}_{{\rm{s}}}),\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{\alpha {(1-{y}_{{\rm{s}}})}^{2}}{{(\alpha +1)}^{2}(\alpha +2)(\delta +2)}.\,\end{eqnarray} \tag{ 37 }$$

Upon substituting (33), (35), and (37) into (31), we obtain

$$\begin{eqnarray}&&\mathrm{cov}({\rm{\Delta }}x,{\rm{\Delta }}{t}_{-})=\displaystyle \frac{\alpha }{(\delta +2){(\alpha +\delta +2)}^{2}(\alpha +\delta +3)}.\end{eqnarray} \tag{ 38 }$$

Similarly, the total variances evaluate to give

$$\begin{eqnarray}&&\mathrm{var}({\rm{\Delta }}x)=\displaystyle \frac{\alpha +\delta +1}{{(\alpha +\delta +2)}^{2}(\alpha +\delta +3)}\end{eqnarray} \tag{ 39 }$$

and $\mathrm{var}({\rm{\Delta }}{t}_{-})=\mathrm{var}({\rm{\Delta }}x)$. Hence, from (38) and (39) we arrive at

$$\begin{eqnarray}&&{r}_{-}=\displaystyle \frac{\alpha }{(\delta +2)(\alpha +\delta +1)}\end{eqnarray} \tag{ 40 }$$

for the correlation between sizes and backward waiting times. Equation (40) exhibits the same behavior seen in Monte Carlo simulations and plotted in Figure 4. The correlation increases with α but it asymptotes to a value ${r}_{-}\to {(\delta +2)}^{-1}\leqslant 1/2$, which decreases as δ increases, i.e., as small avalanches are favored more heavily.

The Pearson coefficient r₊ for the correlation between sizes and forward waiting times is hard to calculate analytically. Instead, one can evaluate the integrals in the counterpart of (31) numerically if required. The result exhibits the same behavior seen in Monte Carlo simulations and plotted in Figure 4, i.e., the correlation r₊ decreases with α. It turns out that the relevant integrals diverge for α < 0.5 for the specific form of $\eta (x| y)$ given by (11). The divergence can be fixed by cutting off the domain of integration at some physically appropriate scale, in the same way that a Cauchy PDF (for example) does not have a well-defined mean or variance, unless a cutoff is introduced.