GLITCH OR ANTI-GLITCH: A BAYESIAN VIEW

Recently, Archibald et al. (2013) discovered an unexpected anti-glitch phenomenon in magnetar 1E 2259+586. Unlike a normal glitch, which undergoes a sudden spin-up, this magnetar experienced a sudden spin-down. The mechanism that caused this phenomenon is still under discussion (e.g., Tong 2014; Lyutikov 2013; Katz 2014; Huang & Geng 2014; Ouyed et al. 2013), but to our knowledge no model explicitly predicted an anti-glitch prior to this discovery, although Thompson et al. (2000) predict a similar potential phenomenon in SGR 1900+14.

The data analysis performed by Archibald et al. (2013) shows that during the observation, 1E 2259+586 undergoes two timing events separated by 50–90 days. The first event is a certain anti-glitch, while the nature of the second event is less certain. If it is also an anti-glitch, we might need a qualitatively different physical model to explain its origin. Importantly, however, based on their analysis, Archibald et al. (2013) were unable to distinguish between these two types of glitch for the second event.

Since we know very little about the mechanism behind such a rare phenomenon, any information about it could be helpful to understand its physical cause. In this Letter we seek to use the data by employing Bayesian model selection to distinguish between two competing models, wherein the second event is either a glitch or anti-glitch, respectively. More specifically, we compute the ratio of the evidence for each model (as defined in Section 3) and investigate whether this ratio favors one model over the other.

The structure of this Letter is as follows. In Section 2 we briefly review relevant details of the model for the time of arrival (TOA) of pulses from the progenitor. Section 3 then describes our Bayesian inference method for carrying out model selection. Section 4 presents the results of our analysis, including a careful check of their robustness. Finally, Section 5 summarizes our conclusions.

Magnetar 1E 2259+586 was routinely observed by the X-Ray Telescope on board Swift every 2–3 weeks, with more frequent observations being made shortly after the discovery of the first anti-glitch event, reported in Archibald et al. (2013). The observations give the TOA of each X-ray pulse (which can be corrected to the solar system barycenter), which in turn gives the pulse phase of the magnetar. Together with each TOA, the X-ray flux is also recorded. An increase in the X-ray flux helps to pinpoint the epoch of the first glitch event, while for the second event no obvious flux change was detected, thus contributing to the confusion over the second event's type.

We model the magnetar's phase evolution, ϕ(t), with the standard Taylor expansion of frequency and frequency derivatives (e.g., Lorimer & Kramer 2005) using terms up to the second order. The effect of a (anti-)glitch on the phase timing model will be a sudden change in the frequency and frequency evolution after the event (e.g., Equation (11) of Hobbs et al. 2006). Therefore, the difference between the phase model with and without the (anti-)glitch (under the assumption that no decaying frequency increment is required due to none being fitted by Archibald et al. 2013) is

$$\begin{equation} \Delta \Phi _g(t) = \Delta \phi _g+\Delta \nu (t-T_g) + \frac{1}{2}\Delta \dot{\nu }(t-T_g)^2, \end{equation} \tag{ 1 }$$

where T_g is the time of the (anti-)glitch, Δϕ_g is an initial sudden phase change, and Δν and $\Delta \dot{\nu }$ are, respectively, the difference in frequency and its first derivative before and after the glitch.

We further define

$$\begin{equation*} R_i=\frac{\phi _i-N_i}{\nu } \end{equation*}$$

as the time residual after subtracting the model predictions from the data. Here, ϕ_i is the predicted pulse phase at the ith observation time t_i (i.e., ϕ(t_i)), N_i is the exact phase at the TOA, which by definition is an integer, and ν is the frequency according to the model. Together with the observed timing uncertainty, σ_i, we can form

$$\begin{equation} \chi ^2 = \sum _{i=1}^{N}\left(\frac{R_i}{\sigma _i}\right)^2. \end{equation} \tag{ 2 }$$

In the next section we will define the likelihood as proportional to exp (− χ²/2), and use it to evaluate the evidence for each model.

In Bayesian inference, given prior background information ${\cal I}$ and observed data D, the posterior probability $p(\theta |{\cal M},{\cal I})$ for a certain theoretical parameter set θ describing a model ${\cal M}$ follows from Bayes' theorem:

$$\begin{equation} p(\theta |D,{\cal M}, {\cal I})=\frac{p(D|\theta, {\cal M}, {\cal I})p(\theta |{\cal M},{\cal I})}{p(D|{\cal M},{\cal I})}. \end{equation} \tag{ 3 }$$

Here, $p(D|\theta,{\cal M}, {\cal I})$ is the likelihood, i.e., the probability of obtaining the observed data D given a particular set of parameters θ for model ${\cal M}$ and prior information ${\cal I}$. $p(\theta |{\cal M}, {\cal I})$ is the prior probability for the parameters θ of model ${\cal M}$, and $p(D|{\cal M}, {\cal I})$ is the evidence, i.e., the probability of obtaining the observed data given that model ${\cal M}$ is true. The evidence is defined as

$$\begin{equation} p(D|\mathcal {M,I}) = \int _{{\bf \Theta }}p(\theta |\mathcal {M},{\cal I})p(D|\mathcal {M},\theta,{\cal I})d\theta. \end{equation} \tag{ 4 }$$

In other words, the evidence is calculated by marginalizing the likelihood over the space of the model parameters. It is then obvious that the calculated evidence depends on the choice of prior range for the parameters θ.

Suppose there are several competing models $\mathcal {M}_i$ that can explain the data. The probability for each model, given the observed data D, is $p(\mathcal {M}_i|D,{\cal I})$. Applying Bayes' theorem, we obtain

$$\begin{equation} p(\mathcal {M}_i|D,{\cal I})=\frac{p(D|\mathcal {M}_i,{\cal I})p(\mathcal {M}_i|I)}{p(D|{\cal I})}. \end{equation} \tag{ 5 }$$

The odds ratio for two different models $\mathcal {M}_i$ and $\mathcal {M}_j$ can then be constructed as

$$\begin{equation} O_{ij}=\frac{p(\mathcal {M}_i|D,{\cal I})}{p(\mathcal {M}_j|D,{\cal I})}=\frac{p(D|\mathcal {M}_i,{\cal I})p(\mathcal {M}_i|{\cal I})}{p(D|\mathcal {M}_j,{\cal I})p(\mathcal {M}_j|{\cal I})}. \end{equation} \tag{ 6 }$$

Normally, in the absence of further information, the prior probabilities for the two models are taken to be equal. In this case, the odds ratio reduces to the ratio of the marginalized likelihoods, or evidences, which is known as the Bayes factor. Usually, a Bayes factor of 10 is already strong enough to favor one model over the other, while a Bayes factor of 100 is regarded as decisive (Jeffreys 1961).

The calculation of the evidence can be a time-consuming process for moderate to large dimensional parameter spaces. However, the nested sampling method (Skilling 2006) makes it practicable to calculate evidence efficiently. In this work, we follow closely the method proposed in Veitch & Vecchio (2010) to calculate the Bayes factor using nested sampling.

In our analysis, the two models under consideration only differ in the sign of the frequency change for the second event. In order to avoid undue influence of the prior range on the Bayes factor, we assign identical prior ranges to all common parameters in both models.

4.1. Setting the Priors

4.2. Comparing the Models

To calculate the evidence, a nested sampling code was applied with a stopping criterion set to equal e⁻⁵—i.e., when new live points made an additional contribution to the evidence that was smaller than a fraction e⁻⁵ of the total, the nested sampling code was stopped. The value of n used for setting our priors was initially taken to be 10—i.e., far beyond the 5σ region. The two models were found to have evidence values of ∼e⁻³³ and ∼e⁻²⁹, respectively, which yields a Bayes factor of 42.5 ± 3.4 in favor of the successive anti-glitch model over the anti-/normal glitch pair model. According to the definition of Jeffreys (1961), a Bayes factor larger than 40 is already very strong evidence in favor of the successive anti-glitch model.

For each model, the posterior distribution of the model parameters was resampled appropriately from our nested sampling results. In Figure 1 we show posteriors for the parameters of the second event in the double anti-glitch model. The contour lines correspond to 68.3%, 95.5%, and 99.7% credible intervals. The maximum posterior corresponds to the following best-fitting parameter values: epoch = MJD 56088.4, Δν = −8.2 × 10⁻⁸ Hz, $\Delta \dot{\nu } = 5.2 \times 10^{-15}$ Hz s⁻¹, and phase change =−0.012 cycles.

Zoom In Zoom Out Reset image size

4.3. Robustness Check

We tested the robustness of our results by changing the width of our uniform priors and re-running the nested sampling analysis. Table 1 shows the mean Bayes Factor obtained as n is changed from 10 to 5 and then 3, given 10 independent nested sampling runs for each n. We see that the Bayes Factor fluctuates around a value of ∼45, but in all cases our conclusions are consistent.

Table 1. Estimated Bayes Factors with Uncertainties

n Value	10	5	3
Bayes factor	42.4 ± 3.4	43.5 ± 3.1	48.6 ± 3.4

Notes. Results were obtained from different prior ranges on the model parameters, represented by different n values (for the definition of n, see the main text). A Bayes factor of around 45 is obtained in each case, indicating consistently strong evidence favoring a successive anti-glitch scenario over an anti-/normal glitch pair.

Download table as:  ASCIITypeset image

The timing residuals, after subtracting the best-fitting double anti-glitch model, are shown in the upper panel of Figure 2. This consistency confirms that two anti-glitch events can explain the observed data well.

Zoom In Zoom Out Reset image size

So far, there have been some physical models proposed in order to explain the putative anti-glitch event. Some authors (e.g., Tong 2014; Lyutikov 2013) have suggested that the second timing event is not consistent with these physical models, and moreover have questioned whether the observational evidence for the second event is strong enough in the first place. However, as shown in the lower panel of Figure 2, if we consider only the first event the timing residual will quickly diverge away from zero thereafter, thus showing strong support for the existence of a second timing event. Note that the timing residuals for an anti-/normal glitch pair model are also similar to the upper panel in Figure 2, further supporting the case for a second timing event (but emphasizing that distinguishing between an anti-/normal glitch pair and an anti-glitch pair is less straightforward).

We also applied model selection to the case of two timing events versus one anti-glitch, and the Bayes factor was found to be e²⁰⁸—i.e., overwhelmingly favoring the two event scenario. This result demonstrates how the Bayes factor can favor a more complicated model, notwithstanding that it may require additional parameters, when the data are of sufficient quality and a simpler model cannot give a satisfactory fit.

Finally, a batch of simulated glitch-free residual data was also generated, with each point drawn from a Gaussian distribution using means and standard deviations from TEMPO2 fits and TOA errors. Nested sampling was applied to this simulated data, and the Bayes factor was computed for the comparison of a successive anti-glitch model and an anti-/normal glitch pair model. We calculated 15 Bayes factor ratios based on 15 realizations of fake anti-/normal/glitch-free data. We found that the Bayes factors fluctuated around unity, i.e., exp (0.17 ± 0.33), showing that intrinsic randomness in glitch-free data will not cause a preference for one model over the other.

We have shown that a model with two successive anti-glitches better explains the observed pulsar data presented in Archibald et al. (2013) when compared with an anti-/normal glitch pair model. Our analysis was robust against variations in the prior ranges, with a Bayes factor consistently larger than 40 in favor of two anti-glitches. Meanwhile, the Bayes factor between two events and one event is very large (e²⁰⁸), showing conclusively that the two events scenario is favored over one event.

Prior to the discovery of an anti-glitch, there were already several published papers presenting mechanisms that can cause enhanced spin-down, while after its discovery, a number of further mechanisms have been proposed seeking to explain its physical origin. Roughly speaking, we can divide the proposed mechanisms into three groups: internal, accretion, and magnetospherical.

The internal mechanism is related to that causing glitches in normal pulsars, which can often be satisfactorily explained by the coupling of the crust with an inner faster-rotating superfluid (Link & Epstein 1996), where for a normal pulsar the superfluid interior could not spin slower than the crust. However, as the observed object 1E 2259+586 is a magnetar, where the dominant source of free energy is magnetism instead of rotation, the spin evolution could be vastly different from that of normal pulsars. Thompson & Duncan (1995) suggested that a magnetar could drive differential rotation, which allows a lag in the rotation of the superfluid interior. A sudden rearrangement of the inner structure could induce the interior and crust to corotate again, which would be observed as a sudden spin-down, or anti-glitch (Archibald et al. 2013). Another possible explanation for the faster-rotating crust might be the twist of a crust patch. As the superfluid vortex is pinned to the crust, a plastic deformation for such a patch will lead to a slower-rotating superfluid. A rapid twist would correspond to a conventional spin-up glitch, similar to a normal pulsar counterpart. However, while a gradual twist would have little effect on the secular spin evolution, a rapid unpinning of the associated vortices would give a sudden spin-down, or anti-glitch (Duncan 2013; Thompson et al. 2000).

Accretion mechanisms suggest that the anti-glitch is caused by the accretion of retrograde material from either a Keplerian ring (Ouyed et al. 2013) or from an asteroid (Huang & Geng 2014). Besides retrograde accretion, Katz (2014) also proposes an enhanced propeller effect to explain the anti-glitch. Although most accretion models are able to explain both events during the observation, either with or without being accompanied by radiation, this mechanism does not fit the model of magnetars, which has already been supported by many observations (e.g., Kaspi & Gavriil 2003).

Magnetospheric models (e.g., Tong 2014; Lyutikov 2013; Thompson et al. 2000) explain the observed anti-glitch with either an enhanced particle wind or a twisting of the magnetic field lines. Although these models fit the observational data for magnetars, most magnetospherical explanations are accompanied by strong radiation and/or a change in pulse profile—neither of which were observed during the second timing event for 1E 2259+586. The magnetospheric mechanism is not favored since Figure 2 shows that our analysis strongly favors the existence of a second event.

Among these three mechanisms, our analysis shows that the internal mechanism is most favored. We note that a satisfactory model should be able to explain the two successive anti-glitches that happened within a relatively short period. If the sudden unpinning of the quantum vortex due to the twist of a crust patch is responsible for the anti-glitch, for example, then the gradual plastic deformation of the crust patch should be able to accumulate enough angular momentum within a timescale of several months. If the two anti-glitches are caused by the same mechanism, then the observations may put some constraints on that mechanism. Enhanced radiation, pulse profile changes, and enhanced spin-down were observed for the first event while none of these phenomena were observed for the second. Future observations of similar phenomena with higher timing accuracy and sampling frequency will be helpful in order to more fully understand the mechanism responsible.

The authors thank Robert Archibald, Victoria Kaspi, and Neil Gehrels for helping with the data and discussions. We are also grateful for further useful discussions with Jonathan Katz, Yongfeng Huang, and Hao Tong. Additionally, we give special thanks to Ian Jones for his very helpful feedback on our article. The authors also thank the anonymous reviewers for constructive comments that helped improve the quality of this manuscript. The authors gratefully acknowledge the support of the Science and Technology Facilities Council of the United Kingdom and the Scottish Universities Physics Alliance. Y.-M.H. is supported by the China Scholarship Council.