© The Authors 2026

1. Introduction

The X-ray pulsar PSR J0537−6910 in the Large Magellanic Cloud is exceptional in terms of its rotational properties. Its spin evolution is interrupted by large glitches, i.e. abrupt changes in the rotation and spin-down rates, followed by post-glitch relaxation (see reviews on pulsar glitches: Haskell & Melatos 2015; Zhou et al. 2022; Antonopoulou et al. 2022; Antonelli et al. 2022). This is the most frequently glitching known pulsar, with a quite regular rate of ≳3.2 glitches per year and glitch rotational frequency steps (Δν) among the largest observed in young pulsars (Middleditch et al. 2006; Antonopoulou et al. 2018; Ferdman et al. 2018; Abbott et al. 2021a; Ho et al. 2020, 2022).

The significant glitch activity makes the measurement of the braking index, $n=\nu \ddot{\nu }/{\dot{\nu }}^2$ (where ν, $\dot{\nu },$ and $\ddot{\nu }$ are the pulse frequency and its first and second time derivatives, respectively), quite ambiguous. A pulsar slowing down purely under magnetic dipole braking would have n = 3. Middleditch et al. (2006), Antonopoulou et al. (2018), and Ho et al. (2022) have determined an effective negative braking index from the overall increase in the spin-down rate (a more negative $\dot{\nu }$) in the long run, while Antonopoulou et al. (2018) and Ferdman et al. (2018) have measured high inter-glitch braking indices (n), typically over ∼7 and up to ∼100, from data between successive glitches. These anomalously high values are dominated by the internal torques that couple the crust to the neutron star superfluid interior; separating their effects in the timing data allows one to find the underlying braking index due to the external braking torque.

In Sect. 2 we overview the timing data and rotation of PSR J0537−6910. In Sect. 3 we discuss the post-glitch spin-down evolution of multi-component neutron stars with references to the vortex creep model. In Sect. 4 we present the relationships between the observable quantities and the hidden underlying braking index, with the results for the model parameters detailed in Sect. 5. Section 6 discusses a method for estimating the time to the subsequent glitch, which could be used to schedule an observation campaign covering the next glitch in PSR J0537−6910. We discuss our results and their implications in Sect. 7.

2. Timing behaviour of PSR J0537−6910

PSR J0537−6910, discovered with the Rossi X-ray Timing Explorer (RXTE; Marshall et al. 1998), is located inside the 1 to 5 kyr old supernova remnant N157B (Chen et al. 2006) in the Large Magellanic Cloud, 49.6 kpc away from Earth (Pietrzyński et al. 2019). It is the fastest rotating young pulsar, with a spin frequency ν = 62 Hz, and has a very high spin-down rate ($\dot{\nu }=-1.992\times {10}^{-10}$ Hz s⁻¹), resulting in the largest known spin-down luminosity Ė = 4.9 × 10³⁸ erg s⁻¹. The characteristic age ${\tau }_{\mathrm{c}}=\nu /2\dot{\nu }$ is about 4.9 kyr and the surface dipole magnetic field B_s is inferred to be ≅9.25 × 10¹¹ G. The source is observed primarily in the X-ray band, with some weak radio emission detected recently (Crawford 2024).

Glitches in PSR J0537−6910 were initially reported by Marshall et al. (1998, 2004). Later, all available RXTE observations were analysed by Antonopoulou et al. (2018), who reported a total of 45 glitches in the nearly 13-year period of data. Their findings were in agreement with Ferdman et al. (2018), who also reports no evidence of pulse flux and profile changes associated with the glitches. After the decommissioning of RXTE in 2012, PSR J0537−6910 was observed with the Neutron Star Interior Composition Explorer (NICER) that discovered 15 glitches (Ho et al. 2020; Abbott et al. 2021b; Ho et al. 2022). Five further glitches in the Jodrell Bank Glitch Catalogue¹ are not included in our analysis as the timing data are not yet published for these recent glitches.

All PSR J0537−6910 glitches have similar glitch steps in the rotation frequency (Δν) and the spin-down rate ($\mathrm{\Delta }\dot{\nu }$), and similar relaxation properties. The glitch activity parameter, calculated as $A_{\mathrm{g}}\cong {\displaystyle \sum _{j=1}^{N_{\mathrm{g}}}}{(\mathrm{\Delta }\nu /\nu )}_{\mathrm{j}}/T_{\mathrm{obs}}=9\times {10}^{-7}{\mathrm{yr}}^{-1}$ where T_obs is the total time span containing N_g observed glitches, each with fractional spin increases (Δν/ν)_j. This is the second highest glitch activity parameter after PSR J1023–5746, which has A_g = 14.5 × 10⁻⁷ yr⁻¹ (Gügercinoğlu et al. 2022).

The magnitude of the spin-down rate, $\left|\dot{\nu }\right|$, usually increases at glitches such that $\mathrm{\Delta }\dot{\nu }<0$. This step partially recovers with a nearly constant positive spin frequency second derivative $\ddot{\nu }$ until the time of the next glitch. This is typical inter-glitch timing behaviour in pulsars with large glitches. For PSR J0537−6910, the inter-glitch spin-down rate evolves with $\ddot{\nu }\sim {10}^{-20}$ Hz s⁻² on average, which implies an anomalously large average inter-glitch braking index n_ig ≅ 20. This value can be up to 5 times larger for some glitches. Exponentially relaxing components are not as prominent in PSR J0537−6910 glitches as for other pulsars, for example Vela. However, after the first glitch in the RXTE data, Antonopoulou et al. (2018) identified an exponentially decaying component with a timescale τ_d = 21 ± 4 days followed by a linear spin-down regime with an inter-glitch braking index 7.6 ± 0.1. Exponential recoveries have recently been detected for 11 additional glitches, most of which are followed by regimes with inter-glitch braking indices between 6 and 9 (Zubieta et al. 2026). Ferdman et al. (2018) applied Markov Chain Monte Carlo techniques to simultaneously fit all RXTE inter-glitch $\dot{\nu }$ data by aligning them at the glitch epoch. They measured a decay timescale ${\tau }_{\mathrm{d}}={27}_{-6}^{+7}$ d and an inter-glitch braking index ∼7.4. Using likewise superposed inter-glitch data in terms of $\ddot{\nu }(t)$ (Andersson et al. 2018) or in terms of the short-term braking index (Ho et al. 2020) yields timescales τ_d = 19 − 44 d for the post-glitch exponential relaxations. These exponential recovery timescales are consistent with the vortex creep response of the outer core of the neutron star (Gügercinoğlu 2017).

The high glitch activity of this source makes its braking index measurements ambiguous. In the long term, the evolution of the spin-down, $\dot{\nu }(t)$, has the characteristic sawtooth shape seen in pulsars with regular large glitches; however, the overall slope is negative. This translates into a long-term negative apparent braking index $n\prime \equiv {\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}{\nu }_0/{\dot{\nu }}_0^2\simeq -1.2$ (Antonopoulou et al. 2018; Ho et al. 2022), where ν₀ and ${\dot{\nu }}_0$ are the rotation frequency and spin-down rate, respectively, at some reference epoch and the long-term frequency second derivative ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}$ is measured for the entire data span, over many glitches. The long-term behaviour is likely the artefact of high glitch activity, as discussed in Antonopoulou et al. (2018), Ferdman et al. (2018), and Ho et al. (2022).

3. Spin-down of the multi-component neutron star with internal couplings

In steady-state, the crust and all interior components of the neutron star spin down at the same rate, dictated by the external (pulsar) braking torque and maintained by lags ω ≡ Ω_s − Ω_c > 0 between the rotation rates (Ω_s) of the various interior (superfluid) components and the rotation rate (Ω_c) of the outer crust (neutron star surface) and interior components co-rotating with the outer crust. In response to offsets from the steady-state lag introduced at each glitch (which decreases ω), the observed spin-down rate of the crust experiences step increases at each glitch. Subsequently, the lag will evolve back to the steady-state value under the effect of the external (pulsar) torque, resulting in relaxation of the observed crust spin-down rate (Alpar et al. 1984; Alpar 1989; Gügercinoğlu & Alpar 2020). This is the typical inter-glitch timing behaviour of pulsars with large glitches, of which the Vela pulsar is the prototype (Akbal et al. 2017; Gügercinoğlu et al. 2022). Internal contributions to the glitch-induced change in the spin-down rate will eventually fully relax back to the pre-glitch steady-state. Such features are common to multi-component systems driven by external torques or forces. In addition, as a result of possible structural changes that affect the coupling of the superfluid to normal components, there can be permanent, non-relaxing changes associated with glitches, similarly to what is seen as ‘persistent shifts’ in the spin-down rate of the Crab pulsar (Demianski & Proszynski 1983; Lyne et al. 1992, 2015).

If the coupling between some internal superfluid component and the crust is linear in the lag ω, then post-glitch relaxation will proceed exponentially. The return to steady-state on the timescale τ_d of the exponential decay can be discerned in the post-glitch timing data of many pulsars (e.g. Yu et al. 2013; Lower et al. 2021; Gügercinoğlu et al. 2022; Zubieta et al. 2023; Liu et al. 2024; Grover et al. 2025; Zubieta et al. 2026). On the other hand, if the coupling is non-linear, then the relaxation takes the form of either a power law, which has no scale, or some more complicated response function, which can involve ‘waiting times’ correlated with the offset introduced by the glitch, i.e. with the parameters of the glitch. Exponentials, power laws, and – occasionally – more complicated features are common in the inter-glitch behaviour of pulsars, indicating the presence of both linearly and non-linearly coupled regions. The most common manifestation of non-linear coupling is a roughly constant second derivative ($\ddot{\nu }$), called ‘anomalous’ because its value is too large to be associated with n ∼ 3 dipole radiation torque, as the contribution of the internal torque is larger, and in many cases much larger than that of the external pulsar torque (Alpar & Baykal 2006). Such signatures are common in glitching pulsars (Yu et al. 2013; Lower et al. 2021; Gügercinoğlu et al. 2022; Zubieta et al. 2024; Liu et al. 2024), and are studied most extensively in the Vela pulsar (Alpar et al. 1984; Akbal et al. 2017; Gügercinoğlu et al. 2022; Grover et al. 2025).

The scale-free power laws arising from the internal torque can dominate the spin-down rate for an entire inter-glitch interval, and must be accounted for in order to find the contribution of the external torque and the actual braking index. A dataset containing many glitches, as in the case of PSR J0537−6910, makes it possible to calculate inter-glitch fits that extend all the way until the next glitch (Antonopoulou et al. 2018; Ferdman et al. 2018; Ho et al. 2020, 2022; Abbott et al. 2021b). Ignoring for a moment the effect of the external torque and any permanent changes, the measured $\ddot{\nu }$ from such fits, together with the total increase in the spin-down rate observed at each glitch, gives a recovery timescale of $t_{\mathrm{g}}\cong |\mathrm{\Delta }\dot{\nu }|/\ddot{\nu }$. This timescale can be interpreted as the waiting time for return to steady-state inside the neutron star, when conditions will be suitable for another glitch to occur (Alpar et al. 1984).

The vortex creep model (Alpar et al. 1984) explains the inter-glitch constant second derivative term quite naturally: In the presence of local pinning sites that hinder the free motion of vortices, the radially outward motion of vortex lines required to spin the superfluid down in response to the pulsar torque is achieved by thermal activation of vortices over pinning barriers. The process is highly non-linear because of Boltzmann factors that depend exponentially on the lag ω = Ω_s − Ω_c between the superfluid rotation rate (Ω_s) and the crust rotation rate (Ω_c). In the steady state, vortex creep sustains a radially outward current of vortices and thereby a continuous internal torque, i.e. angular momentum transfer, from superfluid to crust.

A glitch occurs when the difference between the local superfluid velocities and the rotation rate of the crustal normal matter becomes too large to be sustained by pinning forces. Unpinned vortices will move in the radially outward direction from the rotation axis. This results in a sudden transfer of angular momentum to the crust normal matter, observed as the glitch. The sudden decrease δΩ in the superfluid rotation rate and increase in the crust rotation rate reduces the lag, temporarily stopping the creep process. There is thus a sudden decrease in the internal torque (which is mediated by the creep process) on the crust, while the pulsar torque continues to spin the crust down. This results in an increase in the observed spin-down rate of the crust. The glitch-associated decrease (δΩ) in the superfluid rotation rate is proportional to the total number of unpinned vortices and relates to the size of the observed spin-up Δν. The total angular momentum transferred from the superfluid to the crust in a glitch is proportional to both δΩ and the moment of inertia of the distinct superfluid regions, including creep regions and vortex-depleted regions, through which the unpinned vortices have moved.

The vortex creep model also has a linear regime, prevailing in some superfluid regions whose response explains the prompt exponential relaxation and holds information on the temperature and pinning parameters. In PSR J0537−6910, in addition to the nearly constant $\ddot{\nu }$ behaviour and the exponentials seen in 12 glitches, ten inter-glitch intervals show more complicated features. These are not addressed in the present paper, but can provide additional insights about the non-linear creep process and further estimates of the pinning parameters and local temperature. These detailed applications of the vortex creep model to PSR J0537−6910 are deferred to a separate paper (Gügercinoğlu & Alpar, in prep.).

In addition to the recovering response of interior components to glitch-induced offsets, there can be permanent spin-down shifts. Persisting (but not necessarily permanent) shifts in the spin-down rate are seen for several pulsars. Using some of the Crab pulsar’s largest glitches, permanent changes with an average magnitude $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}/\dot{\nu }⟩=2.51(16)\times {10}^{-4}$ (Lyne et al. 2015), have been identified². It has already been noted by Antonopoulou et al. (2018) that the cumulative effect of persistent shifts associated with glitches could explain the low (under 3 in general) inferred long-term braking indices for several glitching pulsars (Espinoza et al. 2017), including the negative effective braking index of PSR J0537−6910 (Antonopoulou et al. 2018). In this work, we considered the existence of true permanent shifts, which remain once the part of the initial glitch step $\mathrm{\Delta }\dot{\nu }$ in spin-down rate has fully recovered due to both linear and non-linear internal torques (which lead to exponential and power-law-like recoveries, respectively). Note these permanent shifts differ from the measured persisting changes ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}$) often given in the literature, which typically refer to all components of $\mathrm{\Delta }\dot{\nu }$ that do not recover exponentially. In most pulsars, measured $|\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}|$ will thus contain a part that recovers under the rather stable $\ddot{\nu }$ due to the non-linear internal torques, and will be larger than any true permanent shift. To make the distinction we shall use the notation $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}$ to denote true permanent shifts, with $|\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}|\le |\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}|$.

The Crab pulsar’s glitch behaviour, extensively discussed by Lyne et al. (2015) is unique. Post-glitch recoveries show exponential components but do not display the high $\ddot{\nu }$ inter-glitch behaviour observed in PSR J0537−6910, the Vela and other pulsars. Long-term growing exponentials (in $\dot{\nu }(t)$), with τ = 320 days and initial amplitude $0.54\times |\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}|>0$, settling to a permanent shift $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}<0$ were seen following the 1975, 1989, and 2011 glitches, as described by Lyne et al. (2015, (in their Eq. (6), and displayed in their Fig. 3). They infer persistent shifts for ten Crab glitches in their Table 3 (there is a mistake in Eq. (6), which should contain the term $|\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}|$ rather than $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}$ to reflect the behaviour shown in Fig. 3). Persistent shifts were already noted in the 1975 glitch (Demianski & Proszynski 1983) and in the 1989 glitch by Lyne et al. (1992) who also found the ‘wrong sign’ long-term exponential. As the next glitch often arrives long after the longest exponential component with τ ≅ 320 days has settled to its asymptotic value, the return to steady-state seems to be achieved before the arrival of the next glitch. Thus, the $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}$ in the Crab appear to be permanent. In terms of the vortex creep model, the non-linear response of the internal torque in the Crab pulsar has a more complicated signature than the constant $\ddot{\nu }$ response because of the high internal temperature ACCP96. The observed behaviour led to the idea that permanent shifts occur simultaneously with Crab glitches, possibly associated with crust breaking and motion of crustal plates, leading to formation of new vortex traps or vortex depleted regions; in this scenario the ‘wrong sign’ exponential is explained qualitatively by plate motion towards the rotation axis, inducing inward vortex motion (Alpar et al. 1994, 1996; Gügercinoğlu & Alpar 2019).

Since PSR J0537−6910 is a young pulsar with similarly high spin and spin-down rate as the Crab, the temperature profile inside these neutron stars and the crustal and superfluid pinning stresses could also be similar for the two pulsars. It is, therefore, plausible that PSR J0537−6910 glitches can also result in permanent shifts like those observed from the Crab pulsar, likely to be due to sustained structural changes inside the neutron star that permanently affect the rotational dynamics of some superfluid components. This could be common if the vortex unpinning events take place in conjunction with crust breaking and formation of pinned vortex traps in crustal superfluid components, and such models have been employed for PSR J1119−6127 (Akbal et al. 2015), PSR J1048−5832 (Liu et al. 2025) and the Crab pulsar (Gügercinoğlu & Alpar 2019). For both the Crab and Vela pulsars, model-dependent estimates of the time to the next glitch improve significantly when permanent shifts of magnitude $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}/\dot{\nu }\gtrsim {10}^{-4}$ are assumed (Akbal et al. 2017; Gügercinoğlu & Alpar 2020). Note that whilst the inferred $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}/\dot{\nu }$ is an order of magnitude larger for the Vela pulsar than the corresponding fraction in the Crab pulsar, the permanent $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}$ shifts are similar since the Crab slows down about 25 times faster than Vela.

4. Model

Glitches in PSR J0537−6910 have typical magnitudes Δν/ν = a few × 10⁻⁷ and $\mathrm{\Delta }\dot{\nu }/\dot{\nu }=\text{a few}\times {10}^{-4}$ and happen approximately every ∼100 days. To obtain the contribution to the inter-glitch $\ddot{\nu }$ of the external braking torque from the pulsar’s magnetosphere (hereafter ${\ddot{\nu }}_0$) and the equivalent braking index, the glitch contributions (in the form of both recovering and permanent components) should be assessed and removed. Our approach involves two ingredients:

(i) Most inter-glitch intervals display a roughly constant ${\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}$ (the subscript j indicates this is for the inter-glitch interval after glitch j), which can be measured from observations. This trend can be decomposed in two components, corresponding to the external and internal torque contributions, as ${\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}={\ddot{\nu }}_0+{\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$, in which ${\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$ reflects the superfluid response.

(ii) The total spin-down rate change $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{j}}$ contains a part that recovers during the inter-glitch interval after glitch j, denoted $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$, and a persisting part that has not recovered when the next glitch occurs and might be permanent: $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}\equiv \mathrm{\Delta }{\dot{\nu }}_{\mathrm{j}}-\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$. The cumulative effect of frequent shifts that persist at the time of the next glitch could cause the effective negative braking index of PSR J0537−6910, as noted by Antonopoulou et al. (2018) and Ferdman et al. (2018).

The glitches and persistent changes do not appear associated with changes in the electromagnetic signature of the pulsar, such as profile shape or intensity changes (Ferdman et al. 2018). Therefore, any permanent shifts likely reflect structural changes in the interior of the neutron star, namely in the solid crust. The natural physical model for such sudden and permanent structural changes is the star-quake model (Baym & Pines 1971). However, the overall glitch activity cannot be understood in terms of crust-quakes alone (see e.g. Antonopoulou et al. 2022 and Zhou et al. 2022) and requires the involvement of the internal superfluid (see e.g. the explanation in terms of vortex unpinning in Alpar et al. 1993 and Alpar 1995). Crust breaking associated with permanent shifts in the spin-down rate at the time of a glitch may be what triggers the sudden vortex unpinning, which in turn results in the angular momentum transfer to the outer crust observed as the glitch.

In this work we made the assumption (discussed further in the following) that a glitch takes place when the relaxation from the previous glitch due to the internal torques is complete, i.e. the recovery under ${\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$ stops at the epoch of glitch j + 1 and any remaining shifts in $\dot{\nu }$ are truly permanent. Making use of the θ (or Heaviside) function, θ(t) = 0 for t < 0, and θ(t) = 1 for t ≥ 0, we can then express the spin-down rate evolution over a time span including a total of N_g glitches, with recovering and permanent components, occurring at times t_j:

$$\begin{array}{cc}\hfill \dot{\nu }(t)& ={\dot{\nu }}_0+{\ddot{\nu }}_{0}t+{\displaystyle \sum _{j=1}^{N_{\mathrm{g}}}}\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}\theta (t-t_j)\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^{N_{\mathrm{g}}}}(\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}+{\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}(t-t_{\mathrm{j}}))\theta (t-t_{\mathrm{j}})\theta (t_{\mathrm{j}+1}-t)\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^{N_{\mathrm{g}}}}[\mathrm{\Delta }{\dot{\nu }}_{\mathrm{d},\mathrm{j}}exp(-(t-t_{\mathrm{j}})/{\tau }_{\mathrm{d},\mathrm{j}})]\theta (t-t_{\mathrm{j}}),\hfill \end{array}$$(1)

where ${\dot{\nu }}_0$ is a fiducial value and ${\ddot{\nu }}_0$ its (not directly observable) associated derivative based on the external dipole braking model. The first two terms then describe the evolution due to the external torque, whilst the rest are the glitch-associated changes. The third term, with $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}<0$, represents the permanent shifts, whilst the two last terms are the recovering components due to the internal superfluid torques: $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}<0$ are the changes that recover with rate ${\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$ during the interval t_g, j ≡ t_j + 1 − t_j until the next glitch j + 1, and $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{d},\mathrm{j}}<0$ represent the exponentially recovering components.

The total step increase in the spin-down rate observed at glitch j is the sum of the persistent component and the recovering components:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\dot{\nu }}_{\mathrm{obs},\mathrm{j}}\equiv \dot{\nu }(t_{\mathrm{j}+})-\dot{\nu }(t_{\mathrm{j}})=\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}+\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}+[\mathrm{\Delta }{\dot{\nu }}_{\mathrm{d},\mathrm{j}}],\end{array}$$(2)

where the epoch t_j+ is immediately after the glitch at time t_j. The exponentially recovering components are not considered in this work, as they are not apparent in all glitches and, when measured, carry large uncertainties. Nonetheless, the inferred decay timescales (τ) are under 40 d, so the later inter-glitch evolution will be dominated by the constant $\ddot{\nu }$ terms. The approach can be generalised to include exponentials in cases where it might prove necessary.

It is reasonable to expect that the time of completion of the inter-glitch recovery is close to the time of occurrence of the next glitch on average. This is because – even though the actual glitch trigger can have a stochastic component – physical conditions in the neutron star should be similar at the start of each glitch, given the regularity and similarity of their signatures in this pulsar. Moreover, it was observed for a sample of several pulsars that their average glitch waiting times scale with (and are in fact, very close to) their average $|\mathrm{\Delta }\dot{\nu }|/\ddot{\nu }$, further supporting this idea (Lower et al. 2021; Liu et al. 2024), and Ho et al. (2020, 2022) reported a similar correlation between $|\mathrm{\Delta }\dot{\nu }(t_{\mathrm{j}})|/{\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}$ and the time t_g, j from glitch j to the next glitch, j + 1 for PSR J0537−6910 specifically. We therefore made the aforementioned simplifying assumption, that the inter-glitch recovery is completed exactly at the time when the next glitch arrives. The full recovery of $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$ in the time interval t_g, j ≡ t_j + 1 − t_j until the next glitch means that the contribution to $\ddot{\nu }$ of the internal torque is

$$\begin{array}{c}\hfill {\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}=\frac{|\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}|}{t_{\mathrm{g},\mathrm{j}}}.\end{array}$$(3)

Unfortunately, the recovering and permanent parts $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$ and $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}$ of the glitch step in the spin-down rate cannot be separately inferred from the data, even under the assumption that $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$ fully recovers exactly before the next glitch, because of the unknown external torque contribution. If, however, we assume a constant value of the permanent shift for all glitches, then we can estimate the relaxing part $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$ of the glitch step $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{obs},\mathrm{j}}$ in glitch j as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}=\mathrm{\Delta }{\dot{\nu }}_{\mathrm{obs},\mathrm{j}}-⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩,\end{array}$$(4)

and use it in Eq. (3) to approximate ${\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$. Subtracting this internally driven ${\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$ from the total observed ${\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}$, gives ${\ddot{\nu }}_{0,\mathrm{j}}$ due to the external pulsar torque:

$$\begin{array}{c}\hfill {\ddot{\nu }}_{0,\mathrm{j}}={\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}-{\ddot{\nu }}_{\mathrm{ig},\mathrm{j}},\end{array}$$(5)

from which the braking index can be found. Equivalently, assuming a value for the braking index, the necessary permanent changes can be calculated. The effects of the external and internal torques and glitch-induced step changes in $\dot{\nu }(\mathrm{t})$ are sketched in Fig. 1.

Fig. 1. Upper panel: Model for the inter-glitch evolution of $\dot{\nu }$ from glitch j to glitch j + 1. The solid black line depicts the observed inter-glitch evolution. Assuming that the recovery of the glitch offset in the internal torque is completed at exactly the time of arrival of the next glitch, the observed glitch step ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{j}}$) in the spin-down rate is the sum of the permanent part ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}\mathrm{j}}$) and the recovering part ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$) under ${\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$ due to the internal torque, as shown in blue. Subtracting the contribution of the internal torque from the total observed ${\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}$ gives ${\ddot{\nu }}_{0,\mathrm{j}}$, the slope of the purple line, due to the external torque. The permanent shift ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}$) is only partially counteracted by an amount ${\ddot{\nu }}_{0,\mathrm{j}}{t}_{\mathrm{g},\mathrm{j}}$, as shown by the purple line. The combined effect results in a remaining decrease in $\dot{\nu }$, with an inferred slope of ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term},\mathrm{j}}$ (orange line). Lower panel: Long-term evolution over many glitches, shown in orange. This implies an effective negative ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}$ and a negative effective braking index, given by the slope of the red line.

After many glitches, over the entire time span of the observations, Eq. (1) reduces to

$$\begin{array}{c}\hfill \dot{\nu }(t)={\dot{\nu }}_0+{\ddot{\nu }}_{0}t+{\displaystyle \sum _{j=1}^{N_{\mathrm{g}}}}\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}\theta (t-t_{\mathrm{j}}).\end{array}$$(6)

The second sum in Eq. (1), which describes the recovering part of the glitch-induced steps in the spin-down rate, is zero based on our assumption of complete recoveries. Thus, it is the sum of the permanent shifts $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}<0$ over all glitches that results in the negative long-term slope of $\dot{\nu }(t)$, despite a positive ${\ddot{\nu }}_0$ arising from the external (pulsar) torque.

5. Results

In the following we use a sample of 32 glitches for which ${\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}$ and $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{obs},\mathrm{j}}$ are relatively accurately measured. We chose to exclude glitches with uncertainties larger than 50% in either of these quantities. This criterion leaves 23 RXTE and 9 NICER glitches, whose properties are given in the first 3 columns of Table 1. The average inter-glitch interval for this sample is ⟨t_g⟩≅ 122 days, with an rms deviation of 12 days. For a nominal value of the underlying braking index, taken as n = 3 in this instance, we used the observed parameters and Eqs. (3), (4), and 5) of the model to calculate the values of the permanent shifts $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}$ and ${\ddot{\nu }}_{\mathrm{ig}}$ for these glitches, which are also presented in Table 1.

If the long-term apparent braking index n′ is due to permanent shifts, we can calculate how big – on average – they must be, to bring an actual braking index of n due to the external torque down to the observed negative n′ value. The difference in $\ddot{\nu }$ yields – from Eqs. (4), (3) and (5) and a given n – the average permanent shift needed:

$$\begin{array}{c}\hfill {⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{\mathrm{n}}=⟨({\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}-{\ddot{\nu }}_{\mathrm{n}})t_{\mathrm{g}}⟩=(n\prime -n)\frac{{\dot{\nu }}^2}{\nu }⟨t_{\mathrm{g}}⟩,\end{array}$$(7)

where ${\ddot{\nu }}_0$ is the expected value for a spin-down law with a constant braking index (n).

For the long-term n′, we used the value n′= − 1.234(9), corresponding to ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}=-7.92(6)\times {10}^{-22}$ Hz s⁻², measured by Ho et al. (2022) from a linear fit to the combined RXTE and NICER data. We, moreover, took the reference value n = 3, an average inter-glitch interval of ⟨t_g⟩≅ 122 days, and a typical ${\dot{\nu }}^2/\nu$, and find that the average permanent shifts must be about ${⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{n=3}\simeq -3\times {10}^{-14}$ Hz s⁻¹. This result is not very sensitive to the chosen value of the braking index n; for a range of 1 ≤ n ≤ 7 the corresponding required average permanent shifts vary from −1.5 to −5.5 ×10⁻¹⁴ Hz s⁻¹. In general, the required average $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩$ for the theoretically predicted range of magnetosphere-driven braking models (see e.g. Melatos 1997; Ruderman et al. 1998; Yan et al. 2012; Arzamasskiy et al. 2015; Ekşi et al. 2016; Zhang et al. 2022) are typically a modest fraction of the total change $\mathrm{\Delta }\dot{\nu }$ (quoted in Table 1). An upper limit for the permanent shifts can be obtained from the persistent steps: after the first glitch, for which the exponential decay can be accounted for, the persisting change is found to be $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p},0537}=-7.0\pm 0.4\times {10}^{-14}$ Hz s⁻¹ (Antonopoulou et al. 2018). In the limiting case where the average permanent shifts are as large as the observed persisting one, the underlying braking index would be 9.

For comparison with possible permanent shifts in other pulsars, the ratio $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩/\dot{\nu }$ may be the relevant quantity, as it gives an effective fractional moment of inertia of the component responsible. In the Crab pulsar, the average of the persistent steps reported by Lyne et al. (2015) for ten glitches is ${⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{\mathrm{Crab}}\sim -9\times {10}^{-14}$ Hz s⁻¹, which gives a fractional change of $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩/\dot{\nu }{|}_{\mathrm{Crab}}=(3\pm 2)\times {10}^{-4}$. This value is rather close to the value we infer for PSR J0537−6910 under an assumed n = 3: ${⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{n=3}/\dot{\nu }{|}_{0537}=1.44\times {10}^{-4}$. For the Vela pulsar, Akbal et al. (2017) calculated the required permanent shifts in a model-dependent way, by assuming the same glitch signature as in Eq. (1) and calculating the average $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩$ that minimises the difference between observed and predicted (assuming full recovery) glitch epochs. They obtained $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩/\dot{\nu }{|}_{\mathrm{Vela}}=2.32\times {10}^{-3}$, which is an order of magnitude larger than the Crab and PSR J0537−6910 results³. However, their inferred ${⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{\mathrm{Vela}}=(-3.64)\times {10}^{-14}$ Hz s⁻¹ is very similar to the ones we find for PSR J0537−6910. Similarly to the upper bound we obtain on $|\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}|$ from observations of the first PSR J0537−6910 glitch, a limit for the permanent changes in Vela can be found from the observed persisting shifts (once exponentials have been removed). Using the persisting changes measured by Yu et al. (2013) for four Vela glitches, we find the limit $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩/\dot{\nu }{{|}_{\mathrm{Vela}}<⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{p}}⟩/\dot{\nu }|}_{\mathrm{Vela}}=(5\pm 1)\times {10}^{-3}$. It is thus conceivable that permanent changes such as the ones we assumed here occur in young glitching pulsars.

These considerations provide proof of concept that the irregular rotational evolution of PSR J0537−6910 can be explained by a regular underlying braking index, for instance n ∼ 1 − 7, and the effects of its glitches alone. The high inter-glitch braking indices are a result of the internal superfluid torque that dominates until the next glitch, whilst the anomalous long-term negative n′ is the cumulative result of glitch-related, abrupt permanent shifts with physically plausible magnitudes that can be accommodated by the timing data.

6. When is the next glitch due?

Within our model, given a chosen $(⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩,n)$ set of values obeying Eq. (7), one can deduct the contribution of the internal torque ${\ddot{\nu }}_{\mathrm{ig}}={\ddot{\nu }}_{\mathrm{obs}}-n{\dot{\nu }}^2/\nu$. If the last glitch of the pulsar has a measured $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{obs},\mathrm{last}}$ and we assume it involved a permanent shift of order $⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩$, then the recovering part is estimated as $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{last}}\simeq \mathrm{\Delta }{\dot{\nu }}_{\mathrm{obs},\mathrm{last}}-{⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{\mathrm{n}}$. From these we can estimate the time it takes for a full recovery:

$$\begin{array}{c}\hfill t_{\mathrm{g},\mathrm{next},\mathrm{n}}\cong \frac{|\mathrm{\Delta }\dot{\nu }(t_{\mathrm{ig},\mathrm{last}})|}{{\ddot{\nu }}_{\mathrm{ig},\mathrm{last}}}\cong \frac{\mathrm{\Delta }{\dot{\nu }}_{\mathrm{obs},\mathrm{last}}-(n\prime -n)({\dot{\nu }}^2/\nu )⟨t_{\mathrm{g}}⟩}{{\ddot{\nu }}_{\mathrm{obs},\mathrm{last}}-n{\dot{\nu }}^2/\nu },\end{array}$$(8)

in which we substituted ${⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{\mathrm{n}}$ from Eq. (7). This timescale approximates the waiting time until the next glitch.

Another estimate of the time to the next glitch can be found using the vortex creep model. The sudden motion of unpinned vortices in the glitch transfers angular momentum from the superfluid to the crust. Equating the angular momentum transfer from the superfluid to the angular momentum change of the crust (manifested as the observed glitch), and relating the moments of inertia of the vortex creep and vortex-depleted superfluid regions with the step in the spin-down rate gives (Alpar et al. 1984)

$$\begin{array}{c}\hfill \mathrm{\Delta }{\nu }_{\mathrm{j}}=\frac{\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}}{\dot{\nu }}(1/2+{\beta }_{\mathrm{j}})\delta {\nu }_{\mathrm{j}}\cong \frac{\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}}{\dot{\nu }}(1/2+{\beta }_{\mathrm{j}})|\dot{\nu }|t_{\mathrm{g},\mathrm{j}},\end{array}$$(9)

where δν_j is the decrease in the superfluid rotation frequency resulting from the sudden vortex motion driving the spin-up Δν_j. In the vortex creep model $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}/\dot{\nu }$ is interpreted as the fractional moment of inertia, denoted I_A/I, of the vortex creep regions where vortices ‘creep’ against background of pinning sites, and β_j is a model parameter representing the ratio of the moments of inertia of the vortex-depleted superfluid regions (whose fractional moment of inertia is denoted by I_B/I in the vortex creep literature) to those regions that contain creeping vortices (Cheng et al. 1988; Alpar & Baykal 2006).

The next glitch is expected to occur roughly when the offset lag between the superfluid and the crust-normal matter rotation rates is recovered by the spin-down of the crust-normal matter. This leads to the estimation $t_{\mathrm{g},\mathrm{j}}\cong \delta {\nu }_{\mathrm{j}}/|\dot{\nu }|$ and to the second, approximate, equality in Eq. (9). Using Eqs. (3), (5), and (7) and assuming that the parameter β_j = β, a constant representative value for all glitches, we obtain

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{\nu }_{\mathrm{j}}}{{\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}-n{\dot{\nu }}^2/\nu }\cong (1/2+\beta ){t_{\mathrm{g},\mathrm{j}}}^2.\end{array}$$(10)

In the case of PSR J0537−6910, assuming n = 3 and the glitches in Table 1, a fit of Eq. (10) returns β = 18.4. This is much larger than inferred β parameters for other pulsars, which was found to vary from as low as ≲0.5 for the Crab and two more pulsars (Gügercinoğlu & Alpar 2019; Gügercinoğlu et al. 2022) up to ≅5 for other sources (Liu et al. 2024, 2025). The difference in inferred β for this pulsar poses some question for the interpretation within the creep model, but in principle β could be used in Eq. (10) together with Δν_last and ${\ddot{\nu }}_{\mathrm{ig},\mathrm{last}}$ following the last glitch, to approximate t_g, next.

Even though none of the current observed (Ho et al. 2020; Lower et al. 2021; Liu et al. 2024) or theoretical correlations provides a very tight estimate of the time of a future glitch, the predictability of glitches in this pulsar can be used for targeted observations. The strong correlation between Δν_obs, j and t_g, j that is well established for this pulsar (Middleditch et al. 2006; Antonopoulou et al. 2018) gives a typical uncertainty window around 10–20 d for the time of the next glitch, assuming good coverage to accurately measure the Δν of the previous glitch (Ho et al. 2020, 2022).

7. Discussion and conclusions

We have studied the timing behaviour of the source PSR J0537−6910, which is unique in terms of its many glitches, using published timing fits (Antonopoulou et al. 2018; Ho et al. 2020, 2022; Abbott et al. 2021a) that measure a high, relatively constant, inter-glitch frequency second derivative $(\ddot{\nu })>0$ and an effective ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}<0$ inferred from the decrease in $\dot{\nu }$ in the long term. We propose that the apparent long-term negative ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}$ of PSR J0537−6910 arises from permanent, non-relaxing steps in the spin-down rate associated with each glitch, like the persistent shifts observed in the Crab pulsar (Lyne et al. 2015) and inferred for the Vela pulsar (Akbal et al. 2017) and for PSRs J1119-6127 (Akbal et al. 2015) and J1048–5832 (Liu et al. 2025). The separate effects of the external torque, the internal torque, and that of permanent $\dot{\nu }$ changes cannot be extracted from the timing evolution. Making a justified approximation, however, that a glitch occurs when the recovery (both exponential and power-law-like) from the previous glitch is completed, we find a relationship between the braking index (n) of the external torque and the averaged required permanent changes that will produce the effective ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}<0$ over the years. Adopting an underlying pulsar torque of braking index n ∼ 3 yields the permanent shift values $\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}\sim {10}^{-14}$ Hz s⁻¹ listed in Table 1, with an average value ${⟨\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}}⟩}_{n=3}=-2.86\times {10}^{-14}$ Hz s⁻¹, of a similar order of magnitude as those inferred for the Crab and Vela pulsars. The fractional permanent changes inferred for the Vela pulsar are an order of magnitude larger, which is interpreted as indicating a larger connected network of vortex traps in the Vela pulsar (Akbal et al. 2017).

Thus, we have shown that internal torques combined with permanent shifts and external pulsar torques with braking indices not far from the canonical value of 3 suffice to explain the timing behaviour of PSR J0537−6910. We have based our discussion of the time between glitches on the expectation that the next glitch occurs at about the time when internal torques have recovered to pre-glitch equilibrium conditions. However, the requirement of permanent shifts in the spin-down rate suggests a possible rearrangement of the crust geometry in relation to the glitches. Crust breaking could be initiating local vortex unpinning, which then evolves into a large vortex avalanche and observed glitch, when the internal conditions are appropriate. Middleditch et al. (2006) discuss possible crust-quakes in PSR J0537−6910 in connection to its glitches, and crust failure is often suggested as a possible glitch trigger to explain the scale-invariant glitch size distributions seen in some pulsars like the Crab. Conversely, additional stresses to the crust by the pinned vortices could initiate both a crust-quake (Ruderman 1991) and (subsequent) large-scale unpinning when the superfluid lag grows large. The breaking and moving of a plate has also been invoked to explain some glitch signatures (and associated emission changes) in the highly magnetised pulsar PSR J1119−6127 akbal15; (Akbal & Alpar 2018). Such local breaking requires a non-axisymmetric treatment of elastic deformations, along the lines of recent work by Gangwar & Jones (2024) on mountains on neutron stars. A detailed discussion of the stresses induced by vortex pinning, the building of a vortex avalanche starting with vortices unpinned from a specific broken crustal plate, the motion of unpinned vortices under dissipative interactions with electrons, and their encounters with vortex traps (clusters) to unpin more vortices will be the subject of a separate paper (Gügercinoğlu & Alpar, in prep.). Here we just note in passing that the size of a broken plate, like the maximum height of a mountain above the neutron star surface, determines a local, non-axisymmetric deformation, which can source gravitational wave emission (Chamel & Haensel 2008; Gangwar & Jones 2024).

The large inter-glitch braking indices seen for PSR J0537−6910 (Antonopoulou et al. 2018; Ferdman et al. 2018) have also evoked the idea of possible emission of gravitational waves arising from internal fluid motions from r-modes (Andersson et al. 2018) or crustal deformations (Cutler 2002), which give predicted braking indices (n) of 7 and 5, respectively. Such contributions to the spin-down evolution are not ruled out by the current analysis. The expected gravitational wave strain amplitude from r-modes was estimated as h ∼ 2 − 3 × 10⁻²⁶ (Andersson et al. 2018). Initial searches with the third LIGO (Laser Interferometer Gravitational-Wave Observatory)/Virgo survey have failed to detect gravitational radiation associated with PSR J0537−6910 (Fesik & Papa 2020; Abbott et al. 2021a,b, 2022). The most recent fourth LIGO-Virgo-KAGRA (Kamioka Gravitational Wave Detector) observing run (Abac et al. 2025), which covered the time span containing recent glitches of PSR J0537−6910 reported by Ho et al. (2020), place a constraining upper limit of h ≲ 10⁻²⁶ for the gravitational wave strain amplitude. It remains imperative to discover if gravitational wave emission contributes to the braking of PSR J0537−6910, which – together with continuous monitoring in the X-rays – can significantly alter inferences for the internal and magnetospheric torques of this pulsar.

Frequent monitoring of PSR J0537−6910 is also vital to better characterising its glitches for use as input in theoretical calculations. Model-dependent analyses and simulations can help resolve the ambiguity in separating any permanent shifts in spin-down rate and the ongoing power-law recovery due to internal torques. In the framework of the vortex creep model, this can lead to estimates of the time to the next glitch; a complete analysis of the superfluid recovery contributing to the constant $\ddot{\nu }$ behaviour, as well as the prompt exponential post-glitch relaxations, can yield information on the moments of inertia of the interior components involved, on the temperature, and on superfluid vortex pinning. Likewise, the reported correlations between glitch size and the time interval to the next glitch (Antonopoulou et al. 2018; Ferdman et al. 2018; Ho et al. 2020; Liu et al. 2024) is used to estimate superfluid moments of inertia and the angular momentum exchange at glitches (Gügercinoğlu & Alpar, in prep.).

As a final remark we would like to call attention to possible glitch-induced variations in pulse profile or emission properties. A handful of pulsars exhibited pulse profile changes associated with glitches (see e.g. Table 1 in Zhou et al. 2023). In particular, the best glitch resolved to date, the 2016 Vela glitch, displayed clear polarisation level variations and a nulling pulse around the spin-up (Palfreyman et al. 2018). These ephemeral changes might also point to a glitch-accompanying crust-quake, which can move the footpoints of the magnetic field lines close to the polar cap (Akbal et al. 2015) and release seismic energy into the magnetosphere (Bransgrove et al. 2020). While these short-lived and ephemeral changes observed at radio wavelengths are expected to arise from perturbations in the emission regions just above the polar cap, the origin of the X-rays likely lies in the outer gaps at much higher altitudes in the magnetosphere (Philippov & Kramer 2022). Because PSR J0537−6910 is only seen in X-rays, detecting a glitch-associated pulse profile variation could be more difficult; however, the predictability of a forthcoming glitch in PSR J0537−6910 can be used to plan future targeted observations.

Acknowledgments

We acknowledge the Scientific and Technological Research Council of Turkey (TÜBİTAK) for support under the grant 117F330. EG is supported by the Doctor Foundation of Qingdao Binhai University (No. BJZA2025025) and National Natural Science Foundation of China (NSFC) programme 11988101 under the Foreign Talents Grant QN2023061004L. D.A. acknowledges support from an EPSRC/STFC fellowship (EP/T017325/1). We thank the anonymous referee for many constructive comments and suggestions, that helped to improve the clarity of this article.

References

All Tables

All Figures

Fig. 1. Upper panel: Model for the inter-glitch evolution of $\dot{\nu }$ from glitch j to glitch j + 1. The solid black line depicts the observed inter-glitch evolution. Assuming that the recovery of the glitch offset in the internal torque is completed at exactly the time of arrival of the next glitch, the observed glitch step ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{j}}$) in the spin-down rate is the sum of the permanent part ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per}\mathrm{j}}$) and the recovering part ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{ig},\mathrm{j}}$) under ${\ddot{\nu }}_{\mathrm{ig},\mathrm{j}}$ due to the internal torque, as shown in blue. Subtracting the contribution of the internal torque from the total observed ${\ddot{\nu }}_{\mathrm{obs},\mathrm{j}}$ gives ${\ddot{\nu }}_{0,\mathrm{j}}$, the slope of the purple line, due to the external torque. The permanent shift ($\mathrm{\Delta }{\dot{\nu }}_{\mathrm{per},\mathrm{j}}$) is only partially counteracted by an amount ${\ddot{\nu }}_{0,\mathrm{j}}{t}_{\mathrm{g},\mathrm{j}}$, as shown by the purple line. The combined effect results in a remaining decrease in $\dot{\nu }$, with an inferred slope of ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term},\mathrm{j}}$ (orange line). Lower panel: Long-term evolution over many glitches, shown in orange. This implies an effective negative ${\ddot{\nu }}_{\mathrm{long}-\mathrm{term}}$ and a negative effective braking index, given by the slope of the red line.

In the text