The TRAPUM Small Magellanic Cloud pulsar survey with MeerKAT –
II. Nine new radio timing solutions and glitches from young pulsars

We conducted three pseudo-logarithmically spaced timing campaigns with the full array of the MeerKAT telescope at L-band (central frequency of 1284 MHz, Lehmensiek & Theron, 2012). The data were recorded on the APSUSE on-site computing cluster (Accelerated Pulsar Searching User-Supplied Equipment, used by TRAPUM, Barr 2017, Padmanabh 2021) with a bandwidth of 856 MHz split into 4096 channels and a sampling time of 76 $\upmu$s, in sigproc filterbank format (Lorimer, 2008) with no coherent de-dispersion. The observations started with a short separation (half a day to a day) between observations to ensure phase connection. The separation between observations then increased in a pseudo-logarithmic fashion to an intended total time span of about nine months per campaign.

We used the Filterbanking Beamformer User Supplied Equipment (FBFUSE, used by TRAPUM, Barr, 2017; Padmanabh, 2021; Chen et al., 2021) to form several simultaneous MeerKAT closely packed coherent beams on and around the location of each timed pulsar within a pointing. We initially employed the seeKAT multi-beam localisation software (Bezuidenhout et al., 2023) to obtain a more precise position than had been possible from the discovery data. Then, we repeated the process by centring the next observation’s beam tiling on the new localisation. This increased the S/N of the pulsars (reducing the integration time from 30 minutes to 10-15 minutes) and provided a precise pre-timing position. When this precision was sufficient, pulsars were observed with a single beam. The accurate position obtained is crucial in reducing the timing baseline where a covariance between the spin period derivative and position error is present to about 6 to 9 months. Typically, with larger position errors, this covariance would not be eliminated for closer to a year.

The first campaign aimed to determine a timing solution for the first six Carli et al. (2024) discoveries and the two Titus et al. (2019) discoveries. A pseudo-logarithmic cadence was adopted from the third observation. The campaign started in February 2022 and ended in December 2022 with a total of 11 L-band observations of each of two pointing positions. Pointing ‘SMCTIMING1a’⁴⁴4These target names can be used to search the records of the SARAO web archive. comprised the pulsars J0040$-$7326, J0040$-$7335, J0040$-$7337, J0044$-$7314, J0048$-$7317, J0043$-$7319 and the bright known pulsar J0045$-$7319 (McConnell et al., 1991) for testing purposes. We used the multi-beam capability of the TRAPUM backends to time these pulsars simultaneously. Pointing ‘SMCTIMING1b’ recorded pulsars J0054$-$7228 and J0051$-$7204. As part of this campaign, a single 40-minute observation of ‘SMCTIMING1a’ and ‘SMCTIMING1b’ was completed with MeerKAT’s Ultra High Frequency receiver (UHF, 544-1088 MHz, Lehmensiek & Theron 2014) with a bandwidth of 544 MHz split into 4096 channels and a sampling time of 60 $\upmu$s to characterise eight pulsars in this band. The resulting pulse profiles are shown in subsection 4.2. The data were not used as part of the timing solution.

The second pseudo-logarithmically spaced MeerKAT campaign consisted of observations to obtain the timing solution for a later discovery, PSR J0105$-$7208. It is spatially separated from the other pulsars in this work, and thus was timed on its own, at the centre of the ‘SMCTIMING2a’ pointing location. The campaign started in December 2022 and ended in January 2023 with a total of 8 observations. One observation was affected by lightning but could still be used. The rest of the campaign was observed with Murriyang, which is described in the next section.

The third pseudo-logarithmically spaced MeerKAT campaign was performed to provide a phase-connected post-glitch timing solution for PSR J0040$-$7335 and PSR J0040$-$7337 as they proved difficult to time with Murriyang (described in the next section). The campaign started in June 2023 and ended in August 2023 with a total of 7 observations of one pointing position. Pointing ‘SMCGLITCHERS’ was aimed directly at PSR J0040$-$7335 to obtain maximal gain on this faint pulsar. Taking advantage of the multi-beam capacity of TRAPUM, we also timed all the nearby pulsars that were in the field of ‘SMCPOINTING1a’.

We reported in Paper I that, upon discovery of a pulsar, we attempted to obtain detections in archival Murriyang data. None of the pulsars could be detected. This prevented the extension of our timing solutions into the pre-discovery past. We also conducted our own observations with Murriyang, using the Ultra-Wide-band Low receiver (UWL, Hobbs et al., 2020), covering a frequency range from 0.7 to 4 GHz. Most pulsar discoveries of Paper I are too faint to be detected in one hour or less at Murriyang with the UWL.

We began observing our brightest discovery, PSR J0048$-$7317, with Murriyang early in the MeerKAT SMC survey. This was initially done using Director’s Time, then registered under Project P1054 (‘Follow-up of pulsar discoveries from MeerKAT searches’) from June 2021. PSR J0048$-$7317 then became part of the P1054 project from September 2021 until today. Initially, these (usually) 1-h long observations were performed with a close cadence to obtain a first timing solution, which revealed that the pulsar was young, then monthly to monitor the pulsar for possible glitches and timing noise. Extra observations were performed using Director’s Time between September 2022 and March 2023 under project code PX092⁵⁵5Except for the initial September 2022 Director’s Time observation that was recorded under P1054. to obtain a new coherent timing solution after the pulsar underwent a large glitch in May 2022. Search mode (64 $\upmu$s sampling, 3328 frequency channels, no polarisation) was initially used while there was no timing solution, then followed by search and fold mode from December 2021 to May 2022, and in fold mode with coherent de-dispersion (1024 pulse profile bins, 30 s sub-integrations, 3328 frequency channels, all Stokes parameters recorded) only afterwards when it was deemed reliable. This amounted to a total of 37 successful observations with only one non-detection.

After the first MeerKAT timing campaign determined that PSR J0040$-$7335 and PSR J0040$-$7337 were young pulsars, we began monitoring them with Murriyang Director’s Time to reveal possible glitches or timing noise. They were observed in the same beam, under the target names ‘PWN0040-7337’ or ‘J0040-7335’, in search mode (128 $\upmu$s sampling, 3328 frequency channels) to be able to fold both pulsars from the same observation. The observations began in January 2022 under P1054, then were attributed the Director’s Time project code PX096 from December 2022. The flux of the pulsars was very variable, and often too faint to be detected in short observation times. From March to May 2023 a denser campaign was attempted to obtain a coherent timing solution after both pulsars glitched (PSR J0040$-$7335 glitched in December 2022, PSR J0040$-$7337 glitched in February 2023). This was taken over at MeerKAT due to insufficient sensitivity. In total, 14 observations were performed successfully until December 2023, with integration times between 20 minutes and 4 hours.

Finally, after a successful test in November 2022, two 1.5 h PX096 search mode observations were conducted in March and May 2023 to extend the MeerKAT pseudo-logarithmically spaced campaign on PSR J0105$-$7208. They yielded weak detections, particularly affected by Radio Frequency Interference (RFI). In August 2023, a test 1.5 h fold mode observation was unsuccessful in detecting the pulsar, again due to RFI.

The only pulsar discoveries of Paper I that were not observed with Murriyang are J0040$-$7326 and J0044$-$7314. The two Titus et al. (2019) discoveries PSR J0051$-$7204 and PSR J0043$-$7319 were also only observed with MeerKAT. These four pulsars were timed as part of the multi-beam targets described in the previous section. A test Murriyang observation of PSR J0054$-$7228 was performed on 7 June 2021 but the pulsar was deemed too weak to observe regularly there.

For each MeerKAT observation, the RFI in the raw data were cleaned using pulsarX’s filtool (Men et al., 2023). We then created a high-resolution phase-folded archive with dspsr (Van Straten & Bailes, 2011), downsampling the raw data to 512 channels, 10-s sub-integrations, and 1024 phase bins. The initial phase folding parameters were the discovery localisation⁶⁶6As the MeerKAT timing observations progressed, the localisations were refined with seeKAT as detailed in subsection 2.1., DM and period of each pulsar (which are stated in Paper I). The folded data were further RFI-cleaned with clfd (Morello, 2023), and if required, additional manual cleaning was performed using psrchive (Hotan et al., 2004). Initially, before an ephemeris was obtained, we searched for the highest S/N total profile of each observation with psrchive’s pdmp tool if deemed necessary. We optimised the sum of the phase-folded data over a small range of periods and DMs. Each observation’s data were independently re-folded with these parameters if needed. Eventually, ephemerides were used to re-fold all the high-resolution archives with a single ephemeris. We thus obtained a single (frequency and time summed) pulse profile for each observation with psrchive’s pam tool.

A noiseless template profile was first fitted to the discovery pulse profile of each pulsar, and later a new template was made from a higher S/N detection. This was done using psrchive’s paas tool, manually adding von Mises function components (see subsection 4.2). The final template profile was built from adding up all⁷⁷7One or two were omitted for some pulsars. L-band observations in phase using the ephemerides stated in section 4 with psrchive’s psradd tool, and optimising the archive over 20 pc cm^-3 with pdmp to determine the final DM value of each pulsar exclusively timed at MeerKAT.

We then fitted a topocentric Time of Arrival (ToA) and ToA error to each phase-added observation with psrchive’s pat tool. We used the Fourier domain algorithm with Markov chain Monte Carlo to cross-correlate the noiseless template to each observation’s pulse profile. For the non-glitching pulsars, we fitted a timing model containing celestial coordinates and pulsar spin parameters to the ToAs with tempo2 (Hobbs et al., 2006; Edwards et al., 2006). tempo2 transformed our topocentric ToAs to the Solar System barycentre using the Jet Propulsion Laboratory’s DE405 Solar System ephemeris (Standish, 1998). We used the tempo2 model version number 5.00 and the TT(TAI) clock correction procedure. The model parameters and residuals are presented in subsection 4.1.

The pulse width of the pulsars was obtained from the final summed pulse profile of all observations (that was used to determine ToAs). The initial paas model components were refined and debased using psrsalsa’s fitvonMises tool (Weltevrede, 2016). Then, we added noise to the refined model profile using the off-pulse noise statistics of the summed profile. We again refined the model, fitting a pulse width at 50 per cent of the peak intensity to the model with added noise, and repeated this process 1000 times. The mean and standard deviation of the resulting width distribution are quoted in subsection 4.2.

The reduction of the Murriyang observations was performed similarly to the steps detailed in the previous section. The differences in processing are detailed here. filtool was not applied to the observations of PSR J0048$-$7317 and dspsr’s digifil was used to merge raw data time segments (‘search mode’), unless the data were already folded on-line (‘fold mode’). A static UWL channel mask was generally employed. The fold-mode data were flux calibrated using psrchive’s pac command and flux calibrators provided by the Murriyang observatory⁸⁸8https://www.parkes.atnf.csiro.au/observing/Calibration_and_Data_Processing_Files.html. Where polarisation products were recorded, the fold-mode data were polarisation calibrated using short noise diode recordings performed before each observation. The high-resolution stored archives have 416 channels. Due to an oversight, polarisation calibration was performed on archives with this frequency resolution. When adding observations with psradd, the phase alignment option was used if the observations were folded with an ephemeris that contained red noise terms. The Murriyang pulse template was rotated to have its peak at the same phase as MeerKAT’s, except PSR J0105$-$7208 for which our observations did not accumulate enough total S/N to produce a UWL pulse profile template. The Dispersion Measure obtained from optimising the Murriyang phase-added observations was used as the final value for pulsars timed with both MeerKAT and Murriyang, as the wide frequency range of the UWL allowed a more precise determination⁹⁹9This was not done for PSR J0040$-$7335 and PSR J0105$-$7208 due to their low S/N in Murriyang observations. of this value than MeerKAT’s L-band, though they were consistent. Finally, we fitted a temporal jump between the two observatories with tempo2.

The timing parameters from our best fit solutions for the pulsars that have not undergone a glitch are given in Table 1. The time span of the PSR J0105$-$7208 data is too short to fit a position, therefore we used our best seeKAT timing observation localisation (for which 3 $\upsigma$ errors are quoted but not used in the fit). We find that the apparent change in rotation period that would be caused by the seeKAT position error (if the period was compared at the near and far side of the Earth’s orbit) is negligible compared to the error on the first derivative of frequency for this pulsar. The DM errors quoted are obtained from a pdmp optimisation over 20 pc cm^-3 of the total added L-band observations, and are not used in fitting. We note that none of the pulsars required additional parameters from e.g. binary motion. Additional derived parameters calculated by tempo2: the characteristic age, surface magnetic field strength, and rotational energy loss are also given. The residuals from these best fit models are given in Figure 2. The glitching pulsars’ timing solutions are given in subsection 5.2.

We note that two of the pulsars timed only at MeerKAT have an ‘adolescent’ age of around 200 kyr. The first of these, PSR J0051$-$7204, shows no notable timing noise. However, we could not phase-connect the discovery observation of the second adolescent pulsar PSR J0105$-$7208, which occurred five months before the beginning of the radio timing campaign presented here. We would normally expect the timing solution to successfully extrapolate back to the discovery observation, as the phase-connected time span is of similar length as the extrapolation period, thus we cannot exclude a glitch or timing noise in that observation gap. Furthermore, PSR J0043$-$7319, which has a characteristic age of 1.6 Myr, requires a $\ddot{\nu}$ component to fit the residuals, which is about two orders of magnitude larger than the expected contribution from dipole braking for a braking index of 3. The origin of this is not known, but a glitch or timing noise are again not excluded.

We present the total summed pulse profiles of the pulsars in Figure 3. We noted no significant pulse profile shape change between the MeerKAT L-band/UHF and Murriyang UWL band during observations, except for a longer scattering tail for PSR J0040$-$7335 in the UHF band.

We also present the polarisation profile of PSR J0048$-$7317 over the UWL band in Figure 4. We first performed the steps described in subsection 3.2 to obtain a total flux and polarisation-calibrated phase-folded data archive with 1024 pulse profile bins and 416 frequency channels from 26 fold-mode Murriyang observations. We performed a search for Rotation Measure (RM) using psrsalsa’s rmsynth (Weltevrede, 2016) over $\pm 6\times 10^{4}$ rad m^-2 (in case of a large PWN contribution, see e.g. Piro & Gaensler 2018 and the high-RM Fast Radio Bursts with persistent radio sources like FRB 20121102A Michilli et al. 2018). At the minimum, central and highest frequencies of the UWL band, a RM of $\pm 381$ rad m^-2, $\pm 1.5\times 10^{4}$ rad m^-2, and $\pm 7.5\times 10^{4}$ rad m^-2 is required to cause a rotation of $\frac{\pi}{2}$ rad in a frequency channel respectively. The known range of RMs in the SMC spans $-57$ to $+128$ rad m^-2 (Johnston et al., 2021), and this can result in 0.7–3.6 rotations across the UWL band. We were not able to constrain the RM, possibly due to the polarised S/N being too low and/or the frequency resolution too coarse. The RM can reduce the linear polarisation fraction, thus we can only provide a lower limit on the linear polarisation fraction of the total pulse profile of 6 per cent, which we estimated using psrchive’s psrstat. The circular polarisation fraction is about 20 per cent. Cotton et al. (2024) weakly detected a circular polarisation fraction of 11 per cent at the 2–3 sigma level in their image of the head of the nebula. If the circular polarisation fraction being larger than the linear fraction was an intrinsic property of the emission, this would be the first such case among extragalactic pulsars (Johnston et al., 2021), and though frequency-dependent, it is also relatively uncommon among non-recycled galactic pulsars (e.g. Serylak et al., 2021; Oswald et al., 2023). We have not yet accumulated enough polarised S/N to measure the polarisation Position Angle (PA) variation across the pulse.

To model glitch and red noise parameters, we used the software run_enterprise (Keith et al., 2022; Ellis et al., 2020), a Bayesian timing parameter estimation toolkit, as expanded by Liu et al. (2024) to include glitch parameters. The fitting process is described in detail in Liu et al. (2024). We input initial solutions obtained as described in section 3 with an initial glitch size fitted with tempo2, as a starting point for the Bayesian parameter estimation. We tested two models on PSR J0048$-$7317 and PSR J0040$-$7337: a single glitch of permanent changes in frequency $\Delta\nu$ (Hz) and spin-down rate $\Delta\dot{\nu}$ (Hz s^-1), with or without a single exponential relaxation term (see Equation 1). We did not include in our analysis a glitch-induced change in the second order frequency derivative, $\Delta\ddot{\nu}$, as initial modelling showed that it was consistent with zero. Likewise, we did not include a second, longer-term exponential recovery (with a timescale of up to 1000 days). Indeed, the prior on the second decay time was either too large to be well sampled, or resulted in an unconstrained decay time and amplitudes consistent with zero for the second exponential recovery.

The glitch model is thus expressed as:

where $\Delta\phi$ is the difference in predicted phase after the glitch in radians, $\Delta t$ is the time since the glitch in seconds, $\tau_{\text{decay}}$ is the exponential relaxation timescale in seconds and $\Delta\nu_{\mathrm{decay}}$ its amplitude in Hertz.

In the case of PSR J0040$-$7337, we fitted the timing model parameters to a dataset of one ToA per observation. There are about 30 days between the last pre-glitch and the first post-glitch observation¹⁰¹⁰10The possibility that a glitch had occurred since the previous observation was established by inspection of the first post-glitch observation, which showed a drift in the phase of the folded pulse as a function of time.. As with most large glitches, we are unable to unambiguously track the rotation of the pulsar during this observational gap, as there is a degeneracy between the glitch epoch and the number of phase wraps occurring between these ToAs. Hence the glitch epoch cannot be precisely determined (see the Appendix of Basu et al., 2022 for a discussion of this), and we chose to set the number of rotations between the pre-glitch and post-glitch ToAs that resulted in the least difference in phase. This choice results in the latest possible glitch epoch as the initial input for the analysis (just before the first post-glitch ToA), and the glitch epoch prior was set as uniform within $10$ days of this initial value. The prior range on the exponential recovery timescale was set to be sampled uniformly in log-timescale between 0.1 to 100 days, and the remaining priors were set automatically by run_enterprise as described in Liu et al. (2024), section 3.4: the red noise amplitude prior was sampled uniformly in log-scale between $-16$ to $-8$ yr${}^{\frac{3}{2}}$, the red noise power-law index was sampled uniformly in linear space between $0$ and $10$, the exponential amplitude was sampled uniformly in linear space with a range of $\pm 0.8$ times the input glitch step change in frequency, and finally the priors in the glitch step change in frequency and frequency derivative spanned $\pm 0.8$ times the input values, sampled uniformly in linear space.

We initially fixed the DM, temporal jump between observatories, period and period derivatives using values obtained as described in section 3. The position was also fixed to the point source position (0.2–8 keV band) measured from a Chandra observation of the PWN, RA(J2000)$=$0^h40^m46$\aas@@fstack{s}$3850, Dec(J2000)$=-$73°37′07$\aas@@fstack{\prime\prime}$030 (see subsection 6.1 and Figure 11). After fitting the aforementioned glitch and red noise parameters with run_enterprise, we then input the maximum likelihood solution into a new run_enterprise MCMC process (see Liu et al. 2024) to now simultaneously fit for position, period, period derivatives and refine the red noise parameters, while the glitch solution remained fixed (as well as the DM and temporal jump between observatories). This produced the final ephemeris presented in Table 2.

In the case of PSR J0048$-$7317, the glitch occurred between observations at MeerKAT and Murriyang separated by only three days. The pre-glitch detection with MeerKAT had a S/N of 45.3, which we divided into 10 ToAs, and the post-glitch observation with Murriyang had a S/N of 6.5, which we divided into 2 ToAs. All other observations were included as single ToAs. The next observation was performed one month after the first post-glitch observation. Initial glitch modelling indicated that the first observation after the glitch was not well modelled by a simple step change in frequency and frequency derivative (see Figure 5), and hence a short-duration transient recovery was present.

We attempted to fit the transient recovery with an exponential, however since the exponential is largely determined by two residuals from a single observation, there is a wide range of valid solutions. In order to efficiently explore the parameter space, we applied a logarithm to the prior of the amplitude of the exponential, however we re-weighted the exponential amplitude posterior, multiplying it by the amplitude values. This results in a uniform space posterior distribution (revealing larger exponential recovery sizes) between $2.5\times 10^{-14}$ and $2.5\times 10^{-5}\,$Hz, though the amplitudes were sampled linearly in a logarithmic prior space. The exponential timescale prior was sampled uniformly in linear timescale between 0.01 and 100 days. With a short-duration exponential, there may be some correlation between the glitch epoch and the exponential parameters, and hence we want to allow solutions in the full range of glitch epochs. To achieve this, we fitted for an integer number ($\pm 10$) of phase wraps just prior to the first post-glitch ToA (MJD 59708.7).

We also ran a separate fit with the first two observations after the glitch removed, and the same priors and the two-stage fitting as described for PSR J0040$-$7337 (except that the position used was from the initial timing solution, see section 3). To visualise long-term changes in the frequency and frequency derivative using groups of ToAs, we also used the stride fitting method (see Shaw et al. 2018) on this dataset.

In the case of PSR J0040$-$7335, the glitch was small and phase connection appears to be preserved. As the ToAs after the glitch have quite large errors, the glitch epoch is poorly constrained. We thus set the glitch epoch to the last ToA that seemed to be consistent with the pre-glitch timing parameters and fitted only for the post-glitch change in frequency and frequency derivative with tempo2 (having first obtained a preliminary timing solution as described in section 3).

For PSR J0040$-$7337, an informative exponential relaxation fit cannot be performed to the sparse post-glitch data (there are about 30 days between the last pre-glitch point and the first post-glitch point, and again 20 days before the next). A model with no recovery is preferred by run_enterprise, with an odds ratio of $22\pm 5$.

In the case of PSR J0048$-$7317, models containing an exponential recovery term are preferred by run_enterprise when using the full dataset (i.e. they have a larger likelihood); however, the constraints on the exponential timescale and amplitude are very poor due to the low number of ToAs soon after the glitch (there is also a gap of 30 days between the first and second observation after the glitch). We show the posterior distribution of the model with exponential recovery fitted on the full dataset, as described in the previous section, in Figure 10. A strong covariance between the fitted exponential recovery timescale and amplitude is present: the larger the amplitude, the shorter the timescale – to vanishingly small timescales and large amplitude. They are both poorly constrained, thus we can only provide a 95 per cent upper limit for the recovery timescale of 22 days, and for the exponential decay amplitude the 95 per cent credible interval around the median is $1.2^{+4.1}_{-0.4}\times 10^{-7}$ Hz. As can be seen in the first column of Figure 10, a larger number of fitted phase wraps between the pre- and post-glitch observations results in an earlier fitted glitch epoch, with 5 likely glitch epoch solutions in this time range.

We therefore elected to remove the first two observations after the glitch from the dataset to ignore the strongest transient effect and thus get better constraints on the permanent glitch parameters. We model the data in the same way as PSR J0040$-$7337 (see subsection 5.1). A no-recovery model is preferred with an odds ratio of about 2.7$\pm$0.3. This is the model reported in Table 2, shown in Figure 8 and with stride fitting in Figure 6. This model’s residuals with the first two post-glitch observations included are shown in Figure 5.

We report the timing model parameters for the three glitching pulsars in Table 2. The DM errors quoted are obtained from a pdmp optimisation over 20 pc cm^-3 of the total added UWL (L-band for PSR J0040$-$7335) observations and are not used in fitting. The run_enterprise model residuals for PSR J0048$-$7317 and PSR J0040$-$7337 are shown in Figures 8 and 9, while tempo2’s PSR J0040$-$7335 model residuals are given in Figure 7. For J0048$-$7317, the stride fitting results are shown in Figure 6: the datapoints in the stride windows are in accordance with our model.

Two aspects of the aforementioned glitch activity in the SMC pulsars are worth noticing. First, the inferred glitch sizes $\Delta\nu$ for PSRs J0048$-$7317 and J0040$-$7337 are amongst the largest reported, falling at the high end of the known glitch size distribution (see for example Basu et al. 2022). Secondly, all three glitching pulsars displayed a glitch shortly after their discovery (just over a year later). It is therefore compelling to examine how their glitch activity compares to that of the total glitching pulsar population, which predominantly consists of Galactic neutron stars.

The two SMC sources with major glitches of $\Delta\nu\sim 10^{-5}$ Hz have relatively high spin-down rates and a similar characteristic age $\tau_{\rm c}$, around $10^{4}$ yr. On the other hand, PSR J0040$-$7335 has a lower $|\dot{\nu}|$, an age $\tau_{\rm c}\sim 10^{5}$ yr, and its glitch has a small inferred size. Although it is too early to draw any conclusions, the above is in general agreement with the finding that younger pulsars tend to present larger glitches (e.g. Basu et al., 2022). Glitches of small sizes like the one in PSR J0040$-$7335 are common in both young and old glitching pulsars (Basu et al. 2022; and we use theirs and Espinoza et al. 2011’s Jodrell Bank Centre for Astrophysics Glitch Catalogue in this analysis). Small glitches are harder to detect in the presence of timing noise and infrequent monitoring, thus their population is potentially underestimated.

Glitches are customarily regarded as rare events, with the vast majority of pulsars never seen to glitch at all. The situation changes, however, when we take into account the young characteristic age (and high $|\dot{\nu}|$) of the SMC glitching sources. In the population of pulsars with $5\times 10^{3}<\tau_{\rm c}\text{(yr)}<5\times 10^{5}$, over 30 per cent of sources have at least one reported glitch¹¹¹¹11Similarly, if we consider pulsars in the range $5\times 10^{-13}<|\dot{\nu}|\text{(Hz\,s}^{-1})<5\times 10^{-11}$ (which also encompasses all three SMC glitching pulsars), about 50 per cent have known glitches.. This demonstrates that glitches are not an uncommon feature for pulsars similar to J0040$-$7337, J0048$-$7317, and J0040$-$7335. Their red noise parameters are also typical (Parthasarathy et al., 2019). Moreover, of the glitching pulsars in the above age range, 52 per cent have had giant glitches, which are defined by a magnitude of $\Delta\nu/\nu$ greater than $10^{-6}$; in fact, over $1/4$ of the total detected glitches are ‘giant’ according to the Jodrell Bank Glitch Catalogue.

Earlier works on glitches noted an anti-correlation between the characteristic age $\tau_{\rm c}$ of a pulsar and its glitching rate (e.g. Espinoza et al., 2011; Fuentes et al., 2017). The two most recent studies of glitch rate, which use a large sample of pulsars, confirm this relationship and explore the scaling of glitch rate with other pulsar parameters, such as $\nu$ and $\dot{\nu}$, under different assumptions. The sample of the first study (Basu et al., 2022) consists of pulsars observed at the Jodrell Bank Observatory (JBO) that have at least one detected glitch. An average glitch rate is calculated for each pulsar simply as

with $N_{\rm g}$ the total number of detected glitches and $T_{\rm obs}$ the total monitoring time span (which could be accurately determined for JBO observations). A negative power-law scaling with $\tau_{\rm c}$ was assumed, of the form:

with $\tau_{\rm c}$ in kyr; and similarly, the trend for increasing rate with increasing $\dot{\nu}$ is described as

with $\alpha$, $\beta$ and the proportionality constants determined by the data. To somewhat compensate for the fact that glitch rates can be underestimated, only pulsars with $r>0.05\;{\rm yr^{-1}}$ were considered (currently the smallest inferred rate is about $0.02\;{\rm yr^{-1}}$).

The second study, Millhouse et al. (2022), assumes that glitches occur stochastically and describes glitch activity as a constant-rate Poisson process with the probability of observing $N_{\rm g}$ glitches over an observing timespan $T_{\rm obs}$ being

for a pulsar with a Poisson glitch rate $\lambda$. This assumption is to some degree supported by observations of several frequently glitching pulsars, for which the distribution of waiting times between consecutive glitches is consistent with an exponential probability density function. Not all pulsars are consistent with this description, however. For example, the Crab pulsar observations are inconsistent with a homogeneous Poisson process (e.g. Carlin & Melatos, 2019). Moreover, some neutron stars present a characteristic glitch waiting time – the Vela pulsar being a prominent example of this behaviour. Three pulsars for which the regularity in glitching patterns is well established were excluded from the study (namely, the Vela pulsar, the LMC pulsar J0537$-$6910, and PSR J1341$-$6220).

Two pulsar samples were used to calculate $\lambda$ and infer its scaling with pulsar parameters (power-law relations with $\nu$, $\dot{\nu}$, or a combination of the two were examined): one consisted of 174 glitching pulsars, the other additionally included 233 pulsars that have $N_{\rm g}=0$. The glitching pulsars sample significantly overlaps with the one used in Basu et al. (2022), but the used values of $T_{\rm obs}$ contained more uncertainties (see section 7 in Millhouse et al. (2022)). The preferred model (based on Bayes factor) was an anti-correlation of the Poisson glitch rate $\lambda$ with characteristic age, modelled as

where $\tau_{\rm ref}=1$ yr. Note that the optimal values for the parameter $A$ in Table 3 of Millhouse et al. (2022) are in units $\rm d^{-1}$ (Millhouse, private communication).

It is noteworthy that the power-law index ($\alpha$ for Basu et al. 2022 and $-\gamma$ for Millhouse et al. 2022) of the rate ($r_{\tau}$ and $\lambda$ respectively) to characteristic age relationship for glitching pulsars is consistent between the two studies. Using Table 2, we calculate the observed glitch rate as in Equation 2, as well as the predicted rate based on Equations 3, 4, and 6. The results are presented in Table 3. The rate $\lambda$ is denoted by index $1$ when the optimal parameters for the $N_{\rm g}\geq 1$ sample are used, and by index $0$ for the $N_{\rm g}\geq 0$ sample. Direct comparisons between $r_{\rm obs}$ and $\lambda$ cannot be made: ideally, the median $\lambda_{\rm obs}$ should be calculated for each pulsar individually from a posterior distribution of $\lambda$ based on Equation 5. Nonetheless, the predicted rates $\lambda$ based on the Millhouse et al. (2022) relationship are of the same order of magnitude as $r_{\rm obs}$, and do not lead to very low probabilities of a glitch occurring during our monitoring period.

With a single glitch per pulsar, the derived rate $r_{\rm obs}$ can only be a rough approximation of the true rate, under the assumption that a stationary mean rate can indeed be defined for these sources. We used the entire time interval of observations in Equation 2 to calculate the $r_{\rm obs}$ in Table 3. Upper limits can be derived if the maximum time interval between the glitch epoch and the endpoints of the total observing span is used instead, and are close to $1\;\rm{yr^{-1}}$ for all three pulsars. Given the results in Table 3 for either the $\dot{\nu}$ or $\tau_{\rm c}$ relationship of Basu et al. (2022), the observed rate in the SMC pulsars is in accordance with predictions based (almost exclusively) on Galactic pulsars. Therefore, whilst lucky, it cannot – at the moment – be considered extraordinary that a glitch was detected so soon after their discovery.

The above discussion on glitch rates does not take into account the fact that the observed glitches of PSRs J0048$-$7317 and J0040$-$7337 had amplitudes $\Delta\nu/\nu\gtrsim 10^{-6}$, at the higher end of the known distribution. Whilst the analysis of Basu et al. (2022) included all JBO-observed glitching pulsars, we reiterate that the Millhouse et al. (2022) analysis excludes ‘regularly’ glitching pulsars and assumes a homogeneous Poisson process, which would lead to exponential distributions of interglitch waiting times. Typically, however, pulsars for which such a waiting time distribution is an adequate fit to observational data present at the same time a power law-like distribution of their glitch sizes, with large events being far less common than small ones. On the other hand, large glitches are ordinary in regularly-glitching pulsars. The archetype of this glitch behaviour is the Vela pulsar (PSR B0833$-$45, with $\nu\simeq 11.2$, $\dot{\nu}\simeq-1.567\times 10^{-11}\,{\rm Hz/s}$ and $\tau_{\rm c}\sim 10^{4}$ yr), characterised by predominantly giant glitches ($\Delta\nu/\nu\gtrsim 10^{-6}$) that occur every few years. Both PSR J0040$-$7337 and PSR J0048$-$7317 present very strong similarities to Vela, not only due to their timing parameters and inferred age, but also because of their (preliminary) deduced glitch rate and their large spin-up size and spin-down rate change. For regularly-glitching neutron stars the waiting times between large glitches varies from approximately one year up to about a decade (depending on the pulsar), but typical recurrence times are in the range of 1000 to 3000 days¹²¹²12The LMC PSR J0537$-$6910 also has regular large glitches, very frequently at a rate of about 3 per year (calculated with glitches of any size taken into account). It is though younger ($\tau_{\rm c}\sim 5$ kyr) and has an exceptionally high $\dot{\nu}=-1.99\times 10^{-10}$.. The description of the post-glitch recovery in such pulsars often requires one or more exponentially-relaxing terms, which were not included in the models presented in subsection 5.2. Yet, for PSR J0048$-$7317, the presence of a strong transient component following the glitch, with a large decaying amplitude and a short relaxing timescale (see Figure 10 for their indicative magnitudes), further strengthens the resemblance to the Vela pulsar, for which decaying terms on similar timescales have been identified (Dodson et al., 2002). Another very distinct characteristic of post-glitch recoveries of regular glitchers is a high value of $\ddot{\nu}$ that – once any exponential terms decay – dominates the time interval between glitches and leads to ‘anomalous’ braking indices. As can be seen from Table 2, the apparent braking index $n=\nu\ddot{\nu}/\dot{\nu}^{2}$ for PSR J0048$-$7317 is about $152$, whilst for PSR J0040$-$7337 $n\simeq 91$. This is also suggestive of possible Vela-like glitching behaviour, though it cannot be excluded that such values of $n$ are sometimes an intrinsic property of pulsar spin-down (Parthasarathy et al., 2020). It might take a long time to confirm if that is the case as, by extrapolating from the known regularly-glitching pulsars (see for example the scaling between average waiting time and $\dot{\nu}$ presented in Haskell et al. (2012)), we estimate an average time interval between large glitches of about 3.5 years for PSR J0040$-$7337 and around 6 years for PSR J0048$-$7317. It remains to be seen whether these pulsars are truly Vela-like, they are certainly sources of interest and should be closely monitored in the future.

In paper I, PSR J0040$-$7337 was associated with the Pulsar Wind Nebula in the SNR DEM S5 (Haberl et al., 2000; Filipović et al., 2008; Alsaberi et al., 2019), due to a seeKAT multi-beam localisation using the discovery observation and the period of the pulsar being indicative of a young system. In Figure 11, we compare the radio timing position of PSR J0040$-$7337 with a new 97 ks Chandra ACIS-S image (OBSIDs 24633, 22704) centred on the PWN in SNR DEM S5. The point-like source with soft diffuse emission detected with XMM-Newton (Alsaberi et al., 2019) can be further resolved into a compact nebula and a point source in the Chandra image, which our radio timing position error region overlaps. A typical bow-shock morphology is detected in X-rays with harder emission at the position of the pulsar (as detected in the XMM-Newton image, Alsaberi et al., 2019), that is compatible with the radio morphology which extends further out given the longer lifetime of the radio-emitting electrons. Therefore, the position of the pulsar is consistent with the point source position and places it right behind the bow-shock.

Using the pulsar age, one can better constrain the PWN kinematics in the SNR environment, including the pulsar birth kick velocity. We note that the characteristic age from radio timing assumes a simple magnetic dipole radiation model (introduced in section 1) with a braking index, $n$ of 3 as well as a birth period much less than the current period, and thus may not reflect this pulsar’s true age. Alsaberi et al. (2019) proposed an age between 10 and 28 kyr for the pulsar-PWN system, and the characteristic age matches that upper value. For this higher age, Alsaberi et al. (2019) suggest a kick velocity of 700–800 km s^-1, based on a distance range of 60–67.5 kpc estimated from the SMC depth. They also find that this velocity is plausible based on the PWN and SNR morphology in the local environment. This is within the known distribution of transverse velocities of bow-shock nebula pulsars, which spans 60–2000 km s^-1 (Kargaltsev et al., 2017).

It is important to consider if the kick velocity is in fact different, due to the characteristic age being an incorrect evaluation of the true age. $\tau_{\rm c}$ is related to the true age $\tau$ as

where $P_{\rm birth}$ is the birth spin period of PSR J0040$-$7337. $n$ is extremely difficult to measure (e.g. Johnston & Galloway, 1999), though a range of 1.8–3 can be assumed (Melatos, 1997). The young Galactic pulsars, Vela and the Crab (PSRs J0835$-$4510 and J0534$+$2200) have had their braking indices measured to be 1.4$\pm$0.2 (Lyne et al., 1996) and 2.509$\pm$0.001 (Lyne et al., 1988, 1993) respectively, so we take a range of 1.4–3 for this calculation.

We plot the relationship between the true age, the birth period and the kick velocity¹³¹³13We used a distance of 60 kpc, if instead we use a distance of 67.5 kpc for the Western part of the SMC (Scowcroft et al., 2016), the kick velocity for a true age of 28 kyr is 830 km s^-1. in Figure 12. It shows that unless the birth period was only about 10 ms shorter than the spin period is today, the kick velocity should not exceed 2000 km s^-1. For a more realistic birth period of about 15 ms, we find that the kick velocity should be lower than 800 km s^-1. It is not probable that the true age is much larger than the upper limit set by Alsaberi et al. (2019), and thus the kick velocity much lower, based on the PWN and SNR morphology they present. Therefore, a Vela-like true braking index is unlikely.

The proper motion of the pulsar across the sky, extrapolated from the 700–800 km s^-1 kick velocity range from Alsaberi et al. (2019), does not induce a significant change in the period derivative that has been measured from timing, in part due to the large distance to the system. This change, known as the Shklovskii effect (Shklovskii, 1970), increases the observed period derivative by ($6\pm 1)\,\times\,10^{-20}$ s s^-1, which is 2 orders of magnitude smaller than the error on the period derivative resulting from our timing solution.

The spin-down luminosity of this pulsar is $\dot{E}_{\text{PSR}}=6.3\times 10^{36}$ erg s^-1. Alsaberi et al. (2019) measured the X-ray luminosity of the compact, hard X-ray emission region from the PWN and the pulsar in DEM S5 to be $L_{\text{PWN,X}}=1.5\times 10^{34}$ erg s^-1 at a distance of 60 kpc (Karachentsev et al., 2004). The ratio of $\frac{L_{\text{PWN,X}}}{\dot{E}_{\text{PSR}}}$ is thus about $2.4\times 10^{-3}$. This so-called ‘PWN X-ray efficiency’ is similar to that of other PWNe presented in Kargaltsev et al. (2008). The PWN X-ray power-law spectral index was measured in Alsaberi et al. (2019) to be 2$\pm$0.3. Kargaltsev et al. (2008) also relate the ‘PWN X-ray efficiency’ to the PWN X-ray power-law spectral index, and the values for the PWN of DEM S5 fit well within the scatter of the relationship they present.

We calculate the radio luminosity of the PWN of DEM S5, using the radio flux density and spectral index measurements from Alsaberi et al. (2019)’s Table 2, with Equation 4 from Frail & Scharringhausen (1997). This gives $L_{\text{PWN,radio}}=2.1\times 10^{34}$ erg s^-1, and a PWN radio efficiency of $3.3\times 10^{-3}$. This is similar to the Vela radio pulsar-PWN system (Frail & Scharringhausen, 1997). Rescaling the MeerKAT L-band discovery flux density of PSR J0040$-$7337 to 1400 MHz assuming a power-law radio spectral index of $-$1.60 (Jankowski et al., 2018), we find an approximate radio flux density of $S_{\text{PSR,1400\,MHz}}\simeq 14$ $\upmu$Jy. We can estimate the pulsar’s radio luminosity using Szary et al. 2014’s equation (2): $L_{\text{PSR,1400\,MHz}}\simeq 4\times 10^{29}$ erg s^-1 at 60 kpc. Thus the radio efficiency of the pulsar is $\eta_{\text{radio,PSR}}=\frac{L_{\text{PSR,1400\,MHz}}}{\dot{E}_{\text{PSR}}}\simeq 6\times 10^{-8}$. This fits within the large scatter of the $\eta_{\text{radio}}$-$\dot{E}$ relationship from Szary et al. (2014), among young pulsars with pulsed high-energy radiation and SNR association.

According to the expected X-ray efficiency of pulsars, where the ratio of pulsed X-ray luminosity to spin-down luminosity should be around $10^{-3}$ (Becker & Trümper, 1997), we can expect PSR J0040$-$7337 to have a pulsed X-ray luminosity of about $10^{33}$ erg s^-1 if its X-ray beam crosses our line of sight. This would be the case if the pulsar makes up for a few tenths of the total X-ray luminosity $L_{\text{PWN}}$ (which is of the order of $10^{34}$ erg s^-1). Therefore, our radio timing localisation, the low characteristic age and large spin-down luminosity of the pulsar all confirm that PSR J0040$-$7337 is the engine powering the DEM S5 PWN.

In paper I, PSR J0048$-$7317 was associated with a new PWN with no known parent SNR (Cotton et al., 2024). As for PSR J0040$-$7337, this was enabled by a seeKAT multi-beam localisation using the discovery observation. Furthermore, a preliminary estimate of the period derivative of the pulsar was obtained from two survey observations, which now matches the radio timing characteristic age of 38 kyr we derive here. We show the radio timing position¹⁴¹⁴14This supersedes the preliminary timing position with larger errors published in Cotton et al. (2024). of PSR J0048$-$7317 on the Cotton et al. (2024) radio continuum image of the new PWN in Figure 13. Just like PSR J0040$-$7337, the pulsar is situated at the head of the radio PWN near the possible bow-shock front, with a complex and structured tail extending south of the pulsar, presumably indicating the direction of motion. The pulsar has also a characteristic age and spin-down luminosity of the same order as PSR J0040$-$7337, consolidating this association.

Using the same calculations as for PSR J0040$-$7337, we find an approximate radio flux density of $S_{\text{PSR,1400\,MHz}}\simeq 49$ $\upmu$Jy, and a luminosity $L_{\text{PSR,1400\,MHz}}\simeq 1.3\times 10^{30}$ erg s^-1 at 60 kpc. The radio efficiency of the pulsar is $\eta_{\text{radio,PSR}}=\frac{L_{\text{PSR,1400\,MHz}}}{\dot{E}_{\text{PSR}}}\simeq 5\times 10^{-7}$. Again, this fits within the large scatter of the $\eta_{\text{radio}}$-$\dot{E}$ relationship from Szary et al. (2014), among young pulsars with pulsed high-energy radiation and SNR association, although the pulsar is in the high end of the radio luminosity distribution.

Using the radio timing position, we can constrain the X-ray flux of the pulsar-PWN system in an XMM-Newton 23 ks EPIC-pn observation (OBSID 0110000101). Assuming a power-law index of 1.7 and a Galactic foreground absorption of N${}_{\text{H}}=3\times 10^{20}$ cm^-2 the unabsorbed X-ray flux upper limit in the 0.2–12 keV band is $2.2\times 10^{-15}$ erg cm^-2 s^-1. We also take into account the sensitivity loss due to the position being off-axis in the observation. At a distance of 60 kpc, this yields an upper limit on the PWN X-ray luminosity of $L_{\text{PWN}}\lesssim 9.5\times 10^{32}$ erg s^-1. The ratio of $\frac{L_{\text{PWN}}}{\dot{E}_{\text{PSR}}}$ is thus lower than $4\times 10^{-4}$. This ‘PWN X-ray efficiency’ upper limit is on the low end of the distribution presented in Kargaltsev et al. (2008), but is still greater than e.g. the Vela PWN. If PSR J0048$-$7317 was to contribute 10 per cent of the PWN X-ray luminosity, we could expect its X-ray luminosity to be under $10^{32}$ erg s^-1, which is below the XMM-Newton limiting point source luminosity in the SMC (Haberl et al., 2012). Known X-ray rotation-powered pulsars with a spin-down luminosity of $\mathcal{O}(10^{36}$ erg s${}^{-1})$ have X-ray luminosities that range between $10^{30}$–$10^{33}$ erg s^-1(Kargaltsev et al., 2008; Shibata et al., 2016). Therefore, this non-detection of X-rays does not rule out emission from the system. We note that our upper limit is based on a choice of absorption value and may be increased for a higher absorption level (considering PSR J0048$-$7317 has the highest DM in the SMC).

As stated in Paper I, the localisation of PSR J0040$-$7335 is coincident with the SNR DEM S5, just north of the PWN of PSR J0040$-$7337. However, the DM of PSR J0040$-$7335 is nearly double that of PSR J0040$-$7337, so we had deemed it to be a chance background alignment. There is no multi-wavelength emission detected in the images from Alsaberi et al. (2019) at the new, more precise position provided by our radio timing solution for this young, glitching pulsar; and its characteristic age does not match that of the SNR.

In Paper I, we also established that PSR J0105$-$7208 is unlikely to be associated with DEM S128 because of the very high implied transverse velocity given the remnant age of 12 kyr (Leahy & Filipović, 2022). Here, we determine the characteristic age for the pulsar to be 206 kyr which, even given limitations on characteristic age interpretation, does not match the age of DEM S128.