Bayesian Insights into post-Glitch Dynamics: Model comparison and parameter constraint from decades long observation data of the Crab pulsar

The Crab Pulsar has been observed for over half a century, particularly through radio band observations( eg. Lyne et al. (1992); Espinoza et al. (2011a); Fuentes et al. (2017); Serafin Nadeau et al. (2024)), resulting in a substantial accumulation of data on its rotational evolution, with significant contributions from the Jordrell Bank Observatory Espinoza et al. (2011b); Lyne et al. (2015); Basu et al. (2022). Nearly 30 instances of a sudden increase in the pulsar’s spin frequency, known as glitches, have been recorded. These glitches are typically accompanied by a permanent shift in the spin-down rate of the pulsar(see Alpar et al. (1985a, 1996); Akbal et al. (2017); Wang & Zheng (2019); Ge et al. (2020)). Despite extensive observations, the underlying mechanisms driving these glitches remain unclear. One hypothesis suggests that glitches may result from starquakes Baym & Pines (1971); Pines et al. (1972); Alpar & Baykal (1994); Link et al. (1998); Larson & Link (2002); Keer & Jones (2015); Akbal & Alpar (2018); Lu et al. (2023); Sourie & Chamel (2020) on the neutron star, while another model attributes them to the superfluidity of the neutron star’s matte Anderson & Itoh (1975); Ruderman (1976); Alpar (1977); Anderson et al. (1981); Alpar et al. (1984, 1985b); Link et al. (1993); Link & Cutler (2002); Glampedakis & Andersson (2009); Link (2014); Wlazłowski et al. (2016), which may also account for the abrupt increase in spin frequency.

The Crab Pulsar is particularly active in exhibiting glitches, likely due to its relative youth. Detailed modeling of glitch data from this pulsar is therefore crucial for understanding the phenomenon of pulsar glitches Alpar et al. (1985b, 1996); Lyne et al. (2000); Crawford & Demiański (2003); Wang & Zheng (2019); Sourie & Chamel (2020). In this study, we focus primarily on the persistent post-glitch shifts in spin-down rate following glitches, with a particular emphasis on the relationship between these effects and glitch sizes. This relationship may provide a useful probe for distinguishing among various glitch models, as the recovery behavior following a glitch can be modeled in greater detail due to its extended timescale, thereby offering rich information about the glitch mechanisms.

To analyze the glitch data, we will employ modern Bayesian techniques, which are well-suited for extracting information from observational datasets and for conducting model comparisons Kass & Raftery (1995). Bayesian inference allows for the systematic comparison of different models on a consistent basis, and also facilitates the constraints of model parameters within various frameworks. The objective of this study is to bridge the gap between models and observational data through systematic statistical modeling of post-glitch activity in the Crab Pulsar.

This paper is organized as follows: In Section 2, we will discuss the datasets used and the Bayesian framework implemented, along with the modeling motivations for the different data-fitting models. Section 3 will present the posterior results of the Bayesian inference, comparing the Bayesian evidence across different models, and explore the inter-model numerical conversion relations, to emphasis the deep connection between each two models. In section 4, the discussion of recent observation from 2017 event will be presented

The Bayesian technique requires several ingredients: observation data-set, the theoretical models, prior setting of all the model parameters, likelihood set-up and nested sampling algorithm. With accumulative observations of several decades, the observation data of the Crab pulsar already very significant, In Lyne et al. (2015), people summarized the 45 years of glitch observation on the Crab pulsar, with detailed record of the glitch sizes and the persistent shifts of post-glitch spin-down rates. (See Table 3 of Lyne et al. (2015)). Our primary focus is the relation between persistent shift $\Delta\dot{\nu}_{p}$ and glitch size $\Delta\nu$. The data is listed in the table 1. Larger glitch size normally corresponds to a larger step, this general correlation, in Lyne et al. (2015) was fitted by a linear model. While the detail treatment of this correlation is crucial for distinguishing different models, and only the linear model will prohibit a very bad explanation of the observation data. To explore this correlation systematically by modern statistical approach motivates this study.

Different model choices result in varying Bayesian evidence, reflecting how well each model explains the data. By analyzing this Bayesian evidence, we can perform model comparisons and draw data-driven conclusions about which models are favored by the observational data. Table 2 lists the logarithms of the Bayesian evidence for each models.

Model M1, the linear model that probably matches the starquake hypothesis, is the least favored by the observation data, with a very low value of -280.4. In contrast, the non-linear model, M2, shows a significant improvement in Bayesian evidence, rooted in the fact that the power-law function fits the large glitch-accompanied persistent shift better. Introducing a higher-order correction to Model M2, that is Model M3, further improves the Bayesian evidence significantly.

The final model, M4, particularly noteworthy because its functional form is quite another, yet it demonstrates the highest explanatory power for the observational data. This logarithmic-like dependence might represent some missed correction to the power-law models or a different mechanism altogether. In the following section, we will statistically show the convergence of each model, aiding in understanding the potential origin of this purely mathematical model.

However, it is important to note that, across all four models, the Bayesian evidence remains quite low, with the maximum value being -156.9. Considering this is the $\log(Z)$, the exponential of this value is extremely small. This may arise from the fact that the observational data for small glitch sizes accompanied by persistent shifts are well-observed, resulting in very small uncertainties, yet these phenomena seem challenging to explain smoothly using any elementary functional form. This raises questions about the reliability of the observations or suggests that the small glitch-related persistent shifts may require a more sophisticated model for adequate explanation.

In Figure 1, the heatmap of Bayes’ factors for different model pairs is presented. The Bayes factor is a standard method for model selection (see Kass & Raftery (1995)). A Bayes factor greater than 1 but less than 3.2 for model A versus model B suggests that the observations moderately favor model A, though this is not definitive. If the Bayes factor exceeds 3.2, it provides strong evidence that the data prefer model A. In this heatmap, the natural logarithm of the Bayes’ factors is shown.

This map indicates that from M1 to M4, each successive model represents a significant improvement over the previous one. The conclusion is robust, as the Bayes’ factors are exceptionally large, demonstrating that the measurements are sufficiently precise to distinguish among different theoretical models in this context. This provides overwhelmingly strong evidence from Bayesian inference to rule out the validity of the linear model. The large Bayes’ factor suggests that the linear model fails to capture some dominant feature of the system.

When comparing M3 to M2, the large Bayes factor indicates that the higher-order correction to an idealized power-law model is crucial for explaining the data. Even though this correction is theoretically minor, the data are already significant enough to distinguish this model from the one without the correction. The improvement from M3 to M4 is relatively smaller but still extremely strong. This gives us additional motivation to explore the physical meaning behind this model, as the improvement is critical, even though the model is currently only mathematical in nature.

The detailed posterior contour plots are provided in the supplementary materials. While the parameter space itself may be less interesting in this study—since the parameters lack direct physical interpretation and are purely model parameters—they become crucial in future theoretical discussions. However, translating these posterior samples back to the $\Delta\nu-\Delta\dot{\nu}_{p}$ space is of interest, as it allows us to see the distinct regions in which these models reside. By comparing these regions with the original observational data, we can gain insights into the performance of the various models and correlate these findings with the Bayesian evidence discussed earlier.

All the parameters posteriors are summarized in Table 3. In Figure 2, we illustrate the projection of these posterior distributions of the four different models onto the $\Delta\nu-\Delta\dot{\nu}_{p}$ plane. Our analysis shows that for all theoretical models, when the glitch size is smaller than approximately 3 $\mu$Hz, their performance is comparable, as they tend to occupy similar regions. However, for larger glitches, particularly the last data point with a glitch size exceeding 6 $\mu$Hz, the observed persistent shift is much smaller than predicted by M1 and M2, with the observation falling outside their 84% credible interval.

Comparing the overall bandwidth of different models, we observe that the 84% credible interval bandwidth increases sequentially from M1 to M4. This reflects that, from M1 to M4, the models expand the available parameter space to better explain the observational data. This expansion is also reflected in the Bayesian evidence, which measures the volume of the posterior parameter space.

An interesting observation is the sequentially increasing trend in Bayesian evidence from M1 to M4, which aligns with the models’ progressively better performance in explaining the last large glitch observation. This suggests that the large glitch observation is the key factor in distinguishing among the models from a statistical perspective. Notably, this large data point cannot be well-reproduced without compromising the fit for the small glitch accompanied by persistent shift observations, highlighting the challenge in balancing model accuracy across different scales of glitch sizes.

This increasing trend of approaching the last observation also suggests an inter-model numerical convergence relationship. Each model performs similarly at lower glitch sizes, but in the larger glitch region, from M1 to M4, the models sequentially approach the largest data point. In Figure 3, we show the posterior projection residuals between each pair of models, defined as $\Delta\dot{\nu}_{p,M_{i}}-\Delta\dot{\nu}_{p,M_{j}}$, where the i-th model’s predicted persistent shift is subtracted from the j-th model’s prediction. This allows us to clearly see that all model comparisons approach zero for smaller glitches. However, for larger glitch sizes, the correction from the i-th model to the j-th model predicts a lower $\Delta\dot{\nu}_{p}$.

Statistically, this indicates that the models with better performance and larger Bayesian evidence account for higher-order corrections that become significant only in the case of large glitches. For the M3-M2 pair, this higher-order correction is theoretically well understood—it stems from new physics beyond the power-law process. As we show, this correction is relatively small for small glitches but becomes dominant for large glitch sizes.

For the M4-M3 pair, this higher-order correction is still evident. This suggests that the reason M4 outperforms the other models in Bayesian inference is likely due to this higher-order correction, assuming both models share the same foundational framework. However, the physical meaning behind the higher-order correction that M4 accounts for remains unclear and will require further modeling work. Despite this, the statistical behavior of the correction is clear from our analysis.

In 2017, a new glitch event was recorded, accompanied by a persistent shift in the spin-down rate of the star. This glitch was obscured by the occurrence of three much smaller glitches in later about 600 days Shaw et al. (2018); Basu et al. (2020, 2022); Shaw et al. (2021) and hence they were incorported into an effective event as 16.43 $\mu$Hz, with a persistent shift of 434(8) in units of $\left(10^{-15}\mathrm{~{}s}^{-2}\right)$. Since this event is more than twice the glitch size of the previous largest event, it presents a significant test of the current models discussed in this study. As a result, a detailed Bayesian inference was conducted with this new observation added to the dataset.

Table 4 summarizes the updated Bayesian evidence for the four models inclusion of the 2017 event. The overall sequence of Bayesian evidence remains consistent: from M1 to M4, the Bayesian evidence improves, supporting the conclusion that higher-order corrections are incorporated into the models, making them more accurate. Notably, the Bayesian evidence for the linear model has dropped to an extremely small value, which allows us to confidently exclude this linear model as an explanation and provides strong motivation to propose new theoretical models to interpret this observation.

In Figure 4, we present an updated heatmap of the Bayes’ factors. All the Bayes’ factors in the lower-left corner have increased, reflecting the fact that the inclusion of this new data makes the distinction among models even stronger. The only exception is the M4/M3 comparison, where the logarithm of the Bayes’ factor decreases from 12.90 to 0.3 after incorporating the 2017 event. This is a particularly interesting observation, as it suggests that the 2017 event might be better understood by M3. This single observation also is significant enough to support the overall conclusion. The dataset as a whole still favors Model 4, as the Bayes’ factor remains larger than 1, providing potential evidence that M4 is still the best model to explain the entire dataset.

In Figure 5, we compare the Bayesian evidence of different models with and without the 2017 event included. The conclusion is that, with the inclusion of this new observation, all models more and less face challenges in explaining this large event, as the logarithms of Bayesian evidence have dropped, and the posterior space has become more constrained. Especially, M1 and M2 exhibit a significant drop in evidence, suggesting that these two models are strongly constrained and potentially excluded by the new observation.

For M3 and M4, the decrease in Bayesian evidence is not as pronounced, particularly when compared to M1 and M2. If we compute the Bayes’ factor for the same model with and without the 2017 event, the logarithm of the Bayes’ factor is relatively small. Interestingly, the Bayes’ factor for M4 is even smaller than for M3, consistent with our earlier analysis showing that only the M4/M3 Bayes’ factor decreased when comparing the scenarios with and without the 2017 event.

Most of the models considered in this study are primarily mathematical in nature, with the physical interpretation of their parameters remaining largely undefined. However, even when treating these models as purely mathematical constructs, we can still gain insights into their characteristics by performing numerical expansions to various orders. This approach allows us to explore the inter-model conversions based on these expansions. It is possible that truncating higher-order effects could reduce the more data-favored models to resemble the less-favored ones.

In this analysis, we focus exclusively on the best-fitting parameter sets for the different models, specifically selecting the optimal parameter combinations listed in Table 3. The methodology is as follows: for each model pair, we subtract the less-favored model from the more data-favored model, and then numerically integrate this difference from zero to the final glitch observation point. This integration provides a global residue that quantifies the discrepancy between the models. Our primary focus is on the M3-M2 and M4-M3 model pair conversions. The M2-M1 pair is straightforward, as the first-order expansion of M2 can easily revert to the linear M1 model.

For the M3-M2 pair, see the left panel of Figure 6. The x-axis represents the expansion orders of M3, where the expansion is based on a Taylor series. We observe that the residue is minimized at the fourth-order expansion; beyond this point, increasing the expansion order causes the residue to rise again. This suggests that the fourth-order expansion of M3 is the best numerical fit for M2. The superior performance of M3 can be attributed to its incorporation of higher-order modifications from unknown physical sources, thereby numerically demonstrating the correlation between different models.

The residue for the M4-M3 pair is shown in the right panel of Figure 6. This comparison is more complex because the usual Taylor expansion of the M4 logarithmic model does not monotonically increase with additional terms—each positive term is followed by a negative one, making the expansion less effective at reproducing the less-favored model. To address this, we employ a special expansion of the logarithmic function:

This expansion holds true as long as the condition $\operatorname{Re}\left(\xi_{2}+\alpha_{2}\Delta\nu\right)>1/2$ is satisfied. This condition is automatically met across the entire $\Delta\nu$ range, given that $\xi_{2}$ is set to 1.26 and the minimum value of $\Delta\nu$ is zero. Using this expansion, and as shown in Figure 6, the residue stabilizes after 14 terms. Thus, we can assert that the M3 modified power-law model is equivalent to the 14th-order expansion of the M4 logarithmic model. The better performance of M4 on explaining the observation data can be attributed to its inclusion of higher-order corrections absent in the previous model.

This finding is particularly intriguing because, even without an explicit theoretical framework to derive the logarithmic model, we can numerically establish the connection between existing models and identify the orders of the source of the corrections. Interpreting the physical significance of these higher-order terms will require further theoretical modeling.

In this study, we conducted a detailed Bayesian analysis of the glitch-induced persistent increases in the spin-down rate of the Crab pulsar. By systematically comparing four models, ranging from a simple linear model (M1) to more sophisticated power-law and modified power-law models (M2 and M3), as well as a logarithmic model (M4). We hope to understand how well these models explain the observational data.

Our results indicate that higher-order corrections are crucial for accurately explaining the persistent shift accompanying larger glitches. From M1 to M4, the Bayesian evidence increases, showing that the models progressively improve in capturing the post-glitch behavior, especially when considering larger glitches. The linear model (M1) is significantly disfavored by the data. This highlights the need for more complex models that incorporate effects like superfluidity and its associated corrections, as represented by M2 and M3.

The recent 2017 glitch is a extremely large event. Its inclusion in our analysis further constrained the models, particularly M1 and M2, which were strongly disfavored. However, M3 and M4 remain competitive, with M4 still slightly outperforming M3 statistically and M3 maybe better capturing the nuances of this 2017 largest glitch event property. Further modeling work is needed to fully understand success of model M4.

The post-glitch behavior of the Vela is a classical relaxation process, associated with a response to vortex creep. Unlike this standard case, the Crab pulsar cannot relax toward the pre-glitch state after a glitch and has a persistent torque increase. The effect was generally attributed to a external torque variation triggered by starquakes Link et al. (1992); Link & Epstein (1997); Epstein & Link (2000). According to “plate tectonic” proposed by Ruderman (1991), the persistent shift has linear dependence of glitch size found by Zheng et al. (2024). However our analysis evidently disfavors the linear model explanation of observed data, which means that Bayesian evidence rules out starquake hypothesis. This convinces us that possible physics inside neutron stars should be taken seriously. Our future work will focus on finding real post-glitch physics of the Crab pulsar by the understandings of these models.

CompactObject package: An full-scope open-source Bayesian inference framework especially designed for Neutron star physics https://github.com/ChunHuangPhy/CompactOject, documentation: https://chunhuangphy.github.io/CompactOject/. UltraNest: Buchner (2021), https://github.com/JohannesBuchner/UltraNest.