# Timing six energetic rotation-powered X-ray pulsars, including the fast-spinning young PSR J0058-7218 and Big Glitcher PSR J0537-6910

Wynn C.G. Ho, Lucien Kuiper, Cristobal M. Espinoza, Sebastien Guillot, Paul S. Ray, D.A. Smith, Slavko Bogdanov, Danai Antonopoulou, Zaven Arzoumanian, Michal Bejger, et al.

► To cite this version:

Wynn C.G. Ho, Lucien Kuiper, Cristobal M. Espinoza, Sebastien Guillot, Paul S. Ray, et al. Timing six energetic rotation-powered X-ray pulsars, including the fast-spinning young PSR J0058-7218 and Big Glitcher PSR J0537-6910. Astrophys.J., 2022, 939 (1), pp.7. 10.3847/1538-4357/ac8743. hal-03676099

HAL Id: hal-03676099

https://hal.science/hal-03676099

Submitted on 11 Apr 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Timing Six Energetic Rotation-powered X-Ray Pulsars, Including the Fast-spinning Young PSR J0058-7218 and Big Glitcher PSR J0537-6910

Wynn C. G. Ho<fnref n="1" /> ID, Lucien Kuiper<fnref n="2" /> ID, Cristóbal M. Espinoza<fnref n="3" />,<fnref n="4" /> ID, Sebastien Guillot<fnref n="5" />,<fnref n="6" /> ID, Paul S. Ray<fnref n="7" /> ID, D. A. Smith<fnref n="8" /> ID, Slavko Bogdanov<fnref n="9" /> ID, Danai Antonopoulou<fnref n="10" />, Zaven Arzoumanian<fnref n="11" />, Michał Bejger<fnref n="12" />,<fnref n="13" /> ID, Teruaki Enoto<fnref n="14" /> ID, Paolo Esposito<fnref n="15" />,<fnref n="16" /> ID, Alice K. Harding<fnref n="17" /> ID, Brynmor Haskell<fnref n="13" />, Natalia Lewandowska<fnref n="18" /> ID, Chandreyee Maitra<fnref n="19" />, and Georgios Vasilopoulos<fnref n="20" /> ID

## Abstract

Measuring a pulsar’s rotational evolution is crucial to understanding the nature of the pulsar. Here, we provide updated timing models for the rotational evolution of six pulsars, five of which are rotation phase-connected using primarily NICER X-ray data. For the newly discovered fast energetic young pulsar, PSR J0058−7218, we increase the baseline of its timing model from 1.4 days to 8 months and not only measure more precisely its spin-down rate $\dot{\nu} = (−6.2324 \pm 0.0001) \times 10^{-11}$ Hz s⁻<fnref n="1" /> but also for the first time the second time derivative of its spin rate $\ddot{\nu} = (4.2 \pm 0.2) \times 10^{-21}$ Hz s⁻<fnref n="2" />. For the fastest and most energetic young pulsar, PSR J0537−6910 (with a 16 ms spin period), we detect four more glitches, for a total of 15 glitches over 4.5 yr of NICER monitoring, and show that its spin-down behavior continues to set this pulsar apart from all others, including a long-term braking index $n = −1.234 \pm 0.009$ and interglitch braking indices that asymptote to $\lesssim 7$ for long times after a glitch. For PSR J1101−6101, we measure a much more accurate spin-down rate that agrees with a previous value measured without phase connection. For PSR J1412+7922 (also known as Calvera), we extend the baseline of its timing model from our previous 1 yr model to 4.4 yr, and for PSR J1849−0001, we extend the baseline from 1.5 to 4.7 yr. We also present a long-term timing model of the energetic pulsar PSR J1813−1749, by fitting previous radio and X-ray spin frequencies from 2009–2019 and new ones measured here using 2018 NuSTAR and 2021 Chandra data.

*Unified Astronomy Thesaurus concepts: Neutron stars (1108); Pulsars (1306); Rotation powered pulsars (1408); X-ray sources (1822); Ephemerides (464)*

## 1. Introduction

Accurate measurements of the evolution of a pulsar’s rotation rate are vital for inferring properties of and classifying the pulsar. For example, the spin period $P(=1/\nu)$, where $\nu$ is the spin frequency) and spin period time derivative $\dot{P}$ enable estimates of a pulsar’s age (via the characteristic age $\tau_c \equiv P/2\dot{P}$) and magnetic field strength $B \approx 3.2 \times 10^{19}$ G ($P\dot{P})^{1/2}$], as well as give the rate of rotational energy loss $\dot{E} \approx 4.0 \times 10^{46}$ erg s⁻<fnref n="1" /> $\dot{P}/P^3$ (e.g.,

Shapiro & Teukolsky 1983). For neutron star classification, Figure 1 demonstrates that the ~3000 known pulsars reside in distinct regions in the $P-\dot{P}$ parameter space.

Regular monitoring observations of pulsars enable the detection of glitches, which are sudden changes in the spin frequency of a pulsar that mostly occur in pulsars with $\tau_c < 10^7$ yr (Espinoza et al. 2011; Yu et al. 2013; Liu et al. 2021; Lower et al. 2021; Basu et al. 2022), as well as timing irregularities (i.e., timing noise). Most glitches are thought to be due to unpinning of superfluid vortices in the neutron star (Anderson & Itoh 1975), and they can be used to infer properties of the superfluid interior (Link et al. 1999; Lyne et al. 2000; Andersson et al. 2012; Chamel 2013; Ho et al. 2015).

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA 19041, USA; wynnho@slac.stanford.edu</footnote>

<footnote n="2">²SRON-Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA, Leiden, The Netherlands</footnote>

<footnote n="3">³Departamento de Física, Universidad de Santiago de Chile (USACH), Av. Víctor Jara 3493, Estación Central, Chile</footnote>

<footnote n="4">⁴Center for Interdisciplinary Research in Astrophysics and Space Sciences (CIRAS), Universidad de Santiago de Chile, Santiago, Chile</footnote>

<footnote n="5">⁵IRAP, CNRS, 9 avenue du Colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France</footnote>

<footnote n="6">⁶Université de Toulouse, CNES, UPS-OMP, 31028 Toulouse, France</footnote>

<footnote n="7">⁷Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20735, USA</footnote>

<footnote n="8">⁸Laboratoire d'Astrophysique de Bordeaux, CNRS and Université Bordeaux, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="9">⁹Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="10">¹⁰Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Manchester, M13 9PL, UK</footnote>

<footnote n="11">¹¹X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="12">¹²INFN Sezione di Ferrara, Via Saragat 1, 44122 Ferrara, Italy</footnote>

<footnote n="13">¹³Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warsaw, Poland</footnote>

<footnote n="14">¹⁴Extreme Natural Phenomena RIKEN Hakubi Research Team, Cluster for Pioneering Research, RIKEN, Wako 351-0198, Japan</footnote>

<footnote n="15">¹⁵Scuola Universitaria Superiore IUSS Pavia, Palazzo del Broletto, piazza della Vittoria 15, 27100 Pavia, Italy</footnote>

<footnote n="16">¹⁶Istituto Nazionale di Astrofisica, IASF–Milano, via Alfonso Corti 12, 20133 Milano, Italy</footnote>

<footnote n="17">¹⁷Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA</footnote>

<footnote n="18">¹⁸Swarthmore College, 500 College Avenue, Swarthmore, PA 19081, USA</footnote>

<footnote n="19">¹⁹Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, 85748 Garching, Germany</footnote>

<footnote n="20">²⁰Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote>Received 2022 April 29; revised 2022 July 28; accepted 2022 July 28; published 2022 October 26</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Pulsar spin period *P* and spin period time derivative *Ṗ*. Circles denote pulsars whose values are taken from the ATNF Pulsar Catalogue (Manchester et al. 2005; version 1.66), and stars indicate pulsars considered in this work (see Table 1). Dashed lines indicate characteristic age $\tau_c \equiv P/2\dot{P}$, and dotted lines indicate magnetic field strength $B = 3.2 \times 10^{19} \text{ G} (P\dot{P})^{1/2}$.</figcaption></figure>

Pulsars are also potential sources of detectable gravitational waves that are continuously emitted for the lifetime of a pulsar and occur at gravitational-wave frequencies proportional to the pulsar's spin frequency. The most sensitive searches for continuous gravitational waves are those targeting known pulsars that have an accurate contemporaneous rotation phase-connected timing solution since such a model greatly reduces the parameter space of searches (Riles 2017; Sieniawska & Bejger 2019). Searches of data from the most recent 2019–2020 LIGO/Virgo observing run (O3) include five of the six pulsars studied here (see Table 1), with the exception being PSR J0058−7218 since its spin properties only became known in 2021. Constraints have been placed on the size of a gravitational wave-producing mountain in each pulsar (Abbott et al. 2021b, 2022a, 2022b) and on the size of an r-mode fluid oscillation in PSR J0537−6910 (Abbott et al. 2021a). A search was also made for gravitational waves produced by glitches of PSR J0537−6910 and PSR J1813−1749 (Abbott et al. 2022a). Maintaining X-ray monitoring of the pulsars studied here is necessary to enable stronger constraints, and even detections, in the more sensitive data that will be collected in future gravitational-wave observing runs, such as the next one (O4) scheduled to begin in 2023 March.

Here, we report on timing models for six pulsars (see Table 1), which are the result of monitoring data from Chandra (Weisskopf et al. 2002), NICER (Gendreau et al. 2016), and NuSTAR (Harrison et al. 2013). The spin pulsation of each pulsar is only detectable at X-ray energies, except the radio pulsation detection recently reported for PSR J1813−1749 (Camilo et al. 2021; see below). Our results include new long-term phase-connected timing models for PSR J0058−7218 and PSR J1101−6101 covering time baselines of 8 months and 1.7 yr, respectively. The baselines of the phase-connected timing models of PSR J0537−6910, PSR J1412+7922, and PSR J1849−0001 are extended up to 3 to 4 times those previously reported, and four

new glitches of PSR J0537−6910 are presented. For PSR J1813−1749, we fit spin frequency measurements made over 12 yr by Chandra, Green Bank, NICER, NuSTAR, and XMM-Newton, including two new measurements by Chandra in 2021, to obtain an updated spin-down rate, assuming a constant linear decline. We also report on the pulsation searches of Fermi Gamma-Ray Burst Monitor (GBM; Meegan et al. 2009) and Large Area Telescope (LAT; Atwood et al. 2009) data using our timing models for PSR J0058−7218, PSR J1101−6101, PSR J1412+7922, and PSR J1849−0001.

An outline of the paper is as follows. Section 2 briefly summarizes the relevant information for each of our six pulsars. Section 3 describes the data and its processing. Section 4 presents our results, including the timing models for each of the six pulsars. Section 5 summarizes and discusses some implications of the work presented here.

## 2. Summary of Pulsars

PSR J0058−7218 is a newly identified fast-spinning young pulsar in the 14.7 kyr supernova remnant (SNR) IKT 16 in the Small Magellanic Cloud and is associated with a pulsar wind nebula (Owen et al. 2011; Maitra et al. 2015). Maitra et al. (2021) measured for the first time the timing properties of PSR J0058−7218 using XMM-Newton EPIC-pn in small window mode with a time resolution of 5.7 ms, including $\dot{v} = (-6.1 \pm 0.6) \times 10^{-11} \text{ Hz s}^{-1}$ (1$\sigma$ error) within the 118 ks exposure. The pulsar has a narrow single-peak pulse profile and a high pulsed fraction of $\approx 70\%$ in the 0.4–10 keV band. Its rapid spin-down rate means that PSR J0058−7218 has the fourth highest measured spin-down luminosity $\dot{E}$ among the 3000 known pulsars.

PSR J0537−6910 is the fastest-spinning young pulsar and is in the 1–5 kyr SNR N157B in the Large Magellanic Cloud (Wang & Gotthelf 1998; Chen et al. 2006). The pulsar has a narrow single-peak pulse profile and a pulsed fraction of $\sim 20\%$ in the 2–8 keV band (Marshall et al. 1998; Kuiper & Hermsen 2015). While its spin frequency decreases over the more than 17 yr of combined observations (1998–2011 with RXTE and 2017–2022 with NICER), a remarkable 60 glitches have been measured thus far, including 15 by NICER. This yields an average glitch rate of $\sim 3.5 \text{ yr}^{-1}$ and glitch magnitudes larger than those seen in most pulsars (Middleditch et al. 2006; Antonopoulou et al. 2018; Ferdman et al. 2018; Ho et al. 2020b; Abbott et al. 2021b); thus, PSR J0537−6910 is known as the Big Glitcher (Marshall et al. 2004). Its glitches are unusual in their predictability; in particular, there is a correlation seen between the size of its glitches $\Delta v$ and the time to its next glitch (Middleditch et al. 2006; Antonopoulou et al. 2018; Ferdman et al. 2018; Ho et al. 2020b). Its timing properties are also unusual, with a braking index $n \equiv \dot{v}\nu/\dot{v}^2 = -1.25 \pm 0.01$ (1$\sigma$ error) over the long term (17 yr duration) and a value that approaches $\lesssim 7$ over the short term ($\sim 100$ days) between glitches.

PSR J1101−6101 is associated with the hard X-ray source IGR J11014−6103 (Lighthouse nebula) and SNR MSH 11−61A (also known as G290.1−0.8). The remnant has an age of 10–30 kyr (García et al. 2012), which implies that PSR J1101−6101 is a fast-moving pulsar (Tomsick et al. 2012) with a velocity of 800–2400 km s$^{-1}$, which is consistent with its weak proper-motion upper limit (Pavan et al. 2016). The pulsar has a broad single-peak pulse profile, a pulsed fraction of $\approx 50\%$ in the 0.5–10 keV band, and a spin-down rate of $\dot{v} = (-2.17 \pm 0.13) \times 10^{-12} \text{ Hz s}^{-1}$.
<table><caption>**Table 1**

Properties of the Pulsars Considered in the Present Work</caption><thead><tr><th>Pulsar</th><th>P (ms)</th><th>P&#775; (10<sup>-14</sup> s s<sup>-1</sup>)</th><th>E&#775; (10<sup>37</sup> erg s<sup>-1</sup>)</th><th>B (10<sup>12</sup> G)</th><th>&tau;<sub>c</sub> (kyr)</th><th>SNR</th><th>SNR Age (kyr)</th><th>d (kpc)</th></tr></thead><tbody><tr><td>PSR J0058-7218</td><td>21.8</td><td>2.95</td><td>11</td><td>0.81</td><td>11.7</td><td>IKT 16</td><td>14.7 [1]</td><td>62 [2]</td></tr><tr><td>PSR J0537-6910</td><td>16.2</td><td>5.21</td><td>49</td><td>0.93</td><td>4.91</td><td>N157B</td><td>1-5 [3]</td><td>49.6 [4]</td></tr><tr><td>PSR J1101-6101</td><td>62.8</td><td>0.893</td><td>0.14</td><td>0.76</td><td>111</td><td>MSH 11-61A</td><td>10-30 [5]</td><td>7 &plusmn; 1 [6]</td></tr><tr><td>PSR J1412+7922</td><td>59.2</td><td>0.330</td><td>0.064</td><td>0.45</td><td>285</td><td>...</td><td>...</td><td>&le;3.3 [7]</td></tr><tr><td>PSR J1813-1749</td><td>44.7</td><td>12.7</td><td>5.7</td><td>2.4</td><td>5.58</td><td>G12.82-0.02</td><td>&lt;3 [8]</td><td>6-14 [9]</td></tr><tr><td>PSR J1849-0001</td><td>38.5</td><td>1.42</td><td>0.99</td><td>0.75</td><td>43.1</td><td>...</td><td>...</td><td>7 [10]</td></tr></tbody></table>

**Note.** Spin period *P* and spin period time derivative *P&#775*; spin-down luminosity *E&#775;* = 4.0 × 10⁴⁶ erg s⁻¹ *P&#775;/P³; magnetic field *B* = 3.2 × 10¹⁹ G (*P&#775;P*)¹/²; characteristic age *τc* ≡ *P*/(2*P&#775*); SNR association and age; and distance *d*.

**References.** [1]: Owen et al. (2011); [2]: Graczyk et al. (2020); [3]: Chen et al. (2006); [4]: Pietrzyński et al. (2019); [5]: García et al. (2012); [6]: Reynoso et al. (2006); [7]: Mereghetti et al. (2021); [8]: Brogan et al. (2005); [9]: Camilo et al. (2021); and [10]: Gotthelf et al. (2011).

(1σ error), which comes from two spin frequency measurements using XMM-Newton EPIC-pn in small window mode separated by about a year in 2013 and 2014 (Halpern et al. 2014). Kuiper & Hermsen (2015) measure a pulsed fraction of ~64% in the 2–10 keV band using the 2013 data.

PSR J1412+7922 (also known as Calvera) is a high-Galactic latitude pulsar, whose timing parameters suggest that the pulsar is relatively middle-aged and energetic (Zane et al. 2011; Halpern et al. 2013; Halpern & Gotthelf 2015). The pulsar could be a descendant of a member of the central compact object class of neutron stars (see De Luca 2017 for a review), whose SNR is no longer visible. Alternatively, PSR J1412+7922 may simply be a normal rotation-powered pulsar (Mereghetti et al. 2021). A major contributor to the uncertainty of its classification is its distance, with a lower limit of about 200 pc and a value of up to 3.8 kpc from spectral fitting (Halpern et al. 2013; Halpern & Gotthelf 2015; Mereghetti et al. 2021). Monitoring during the first year of NICER in 2017–2018 yielded a phase-connected timing model (Bogdanov et al. 2019), and subsequent work extended the timing model to over 3 yr using NICER data through 2021 February (Mereghetti et al. 2021).

PSR J1813–1749 is a highly energetic pulsar that produces a pulsar wind nebula and is associated with the gamma-ray/TeV source IGR J18135–1751/HESS J1813–178. The pulsar is located in the young (<3 kyr) SNR G12.82–0.02 (Brogan et al. 2005), and its proper motion from the center of the remnant implies an age of 1000–2200 yr (Dzib & Rodríguez 2021). PSR J1813–1749 has a broad single-peak X-ray pulse profile and a pulsed fraction of 50% in the 2–10 keV band (Gotthelf & Halpern 2009; Halpern et al. 2012; Kuiper & Hermsen 2015; see also Ho et al. 2020a). A phase-connected timing model spanning 37 days in 2019 was determined using NICER data (Ho et al. 2020a). More recently, Camilo et al. (2021) report radio observations at high frequencies using the Green Bank Telescope in which highly scattered radio pulses are detected for the first time for PSR J1813–1749. They also perform a linear fit of individual spin frequencies measured in 2009, 2011, and 2012 in X-ray and in 2012 and 2015 in radio and find a spin-down rate $\dot{\nu} = (-6.3364 \pm 0.0025) \times 10^{-11}$ Hz s⁻¹ (1σ error), which is within 3.2σ of the value determined in Ho et al. (2020a), which includes the 2019 NICER data but not the radio data. With confirmation of the pulsar in radio, a more reliable proper motion is measured using the Very Large Array (VLA; Dzib & Rodríguez 2021), and the apparent X-ray proper motion is likely due to brightness changes in the pulsar wind nebula, as pointed out in Ho et al. (2020a).

<table><caption>**Table 2**

Observation Log</caption><thead><tr><th>Telescope</th><th>Pulsar</th><th>Observation Date</th><th>Exposure (ks)</th></tr></thead><tbody><tr><td rowspan="2">Chandra</td><td rowspan="2">PSR J1813-1749</td><td>2021 Feb 10</td><td>20</td></tr><tr><td>2021 Jun 23</td><td>20</td></tr><tr><td rowspan="5">NICER</td><td>PSR J0058-7218</td><td>2021 Jun 1-2022 Jan 25</td><td>163</td></tr><tr><td>PSR J0537-6910</td><td>2017 Aug 17-2022 Feb 17</td><td>1374</td></tr><tr><td>PSR J1101-6101</td><td>2020 Apr 1-2021 Dec 16</td><td>370</td></tr><tr><td>PSR J1412+7922</td><td>2017 Sep 15-2022 Feb 8</td><td>1379</td></tr><tr><td>PSR J1849-0001</td><td>2018 Feb 13-2021 Nov 16</td><td>260</td></tr><tr><td rowspan="2">NuSTAR</td><td>PSR J1101-6101</td><td>2020 Nov 20</td><td>136</td></tr><tr><td>PSR J1813-1749</td><td>2018 Mar 25</td><td>26</td></tr></tbody></table>

PSR J1849–0001 is another young and energetic pulsar that produces a pulsar wind nebula and is associated with the gamma-ray/TeV source IGR J18490–0000/HESS J1849–000. The pulsar has a broad single-peak pulse profile and a high pulsed fraction of 77% in the 0.06–10 keV band (Kuiper & Hermsen 2015; see also Bogdanov et al. 2019). A phase-connected timing model spanning 20 days in 2010 was determined using RXTE data (Gotthelf et al. 2011; Kuiper & Hermsen 2015), and more recently a phase-connected timing model spanning about 1.5 yr in 2017–2018 was determined using Swift and NICER data (Bogdanov et al. 2019).

## 3. Data and Analysis Method

### 3.1. NICER Data

For five of the six pulsars considered here, we use and report new analyses of NICER data, with the exceptions being Chandra and NuSTAR data of PSR J1101–6101 and PSR J1813–1749 (see Sections 3.2 and 3.3). NICER data for these five pulsars are summarized in Table 2 and are processed following a similar procedure, which we outline, noting source-specific differences, below. More details are provided in Kuiper & Hermsen (2009) for PSR J0058–7218 and PSR J1101–6101, Ho et al. (2020b) for PSR J0537–6910, and Bogdanov et al. (2019) for PSR J1412+7922 and PSR J1849–0001.

We process and filter NICER data of each pulsar using HEASoft 6.22–6.29b (HEASARC 2014) and NICERDAS 2018-03-01_V003-2021-08-31_V008c. We exclude all events from “hot” detector 34, which gives elevated count rates in some circumstances, and portions of exposure accumulated during passages through the South Atlantic Anomaly. While NICER is sensitive to 0.25–12 keV photons, we make an energy cut and
extract only events within a specific energy range optimized for pulsation searches. These are 0.8–5 keV for PSR J0058−7218, 1–7 keV for PSR J0537−6910, 1.5–10 keV for PSR J1101−6101, 0.37–1.97 keV for PSR J1412+7922, and 1.89–6 keV for PSR J1849−0001. For PSR J0058−7218 and PSR J1101−6101, we use cleaned standard event data, which are subsequently screened to exclude events that were collected during high background levels as determined from light curves at 12–15 keV; in instances when the exposure time of cleaned event data is too short to derive an accurate time-of-arrival (TOA) measurement (see below), unfiltered event data with background filtering are used instead. For PSR J0537−6910, PSR J1412+7922, and PSR J1849−0001, we ignore time intervals of enhanced background affecting all detectors by constructing a light curve binned at 16 s and removing intervals strongly contaminated by background flaring when the count rate exceeds a threshold value. The thresholds are $10 \text{ cs}^{-1}$ for PSR J0537−6910, $4.5 \text{ cs}^{-1}$ for PSR J1412+7922, and $5 \text{ cs}^{-1}$ for PSR J1849−0001, and they are the same as previously used in Ho et al. (2020b) for PSR J0537−6910 and in Bogdanov et al. (2019) for PSR J1412+7922 and PSR J1849−0001. Using these filtering criteria, we obtain clean data for pulse timing analysis.

We combine sets of individual ObsIDs into merged observations, with each merged observation yielding a single TOA measurement. ObsIDs are combined such that there is sufficient exposure to confidently detect the spin frequency of each pulsar, with typical total exposures of 20–30 ks for PSR J0058−7218, 4–9 ks for PSR J0537−6910, 30 ks for PSR J1101−6101, 6–12 ks for PSR J1412+7922, and 5–9 ks for PSR J1849−0001. Merged ObsIDs are those acquired usually within a 3–4 day span and on rare occasions within 6–7 days. Before performing a pulsation search, we use barycorr to transform between Terrestrial Time (TT), used for event time stamps, and Barycentric Dynamical Time (TDB) and to account for the effects of satellite motion with respect to the barycenter. In all timing analyses performed here, unless otherwise noted (in particular, Sections 4.4 and 4.6), source positions (and proper motions, if measured previously) are held fixed at the values given in the corresponding tables below, along with our adopted solar system ephemeris.

For PSR J0058−7218 and PSR J1101−6101, the search for pulsations, generation of TOAs, and determination of timing models are conducted by following procedures described in Kuiper & Hermsen (2009). For PSR J0537−6910, PSR J1412+7922, and PSR J1849−0001, we conduct acceleration searches using PRESTO (Ransom et al. 2002), with searches using a time bin and including a number of harmonics that are specific to each pulsar. In particular, these are 0.5 ms and 8 harmonics for PSR J0537−6910, 3 ms and 4 harmonics for PSR J1412+7922, and 1 ms and 4 harmonics for PSR J1849−0001. Pulsations at the spin frequency of each pulsar are usually the strongest detected. Data are folded at the candidate pulse frequency using prepfold, and a refined frequency is determined. On occasion, further iterations are performed to obtain a more robust measurement. Finally, we produce a pulse profile template by fitting a set of NICER pulse profiles with a Gaussian shape; this template is then used to determine the TOA of each merged observation following the unbinned maximum likelihood technique described in Ray et al. (2011). We use TEMPO2 (Hobbs et al. 2006) to fit the TOAs with a timing model and to measure the glitch parameters in the case of PSR J0537−6910.

## 3.2. NuSTAR Data

NuSTAR observed PSR J1101−6101 on 2020 November 20 (ObsID 30601029002) for 136 ks and PSR J1813−1749 on 2018 March 25 (ObsID 30364003002) for 26 ks (see Table 2). For PSR J1101−6101, we use data processed with NuSTAR-DAS v2.0.0 and CALDB v20200811, and we barycenter cleaned event data extracted from a circular region of 50" radius. For PSR J1813−1749, we process the data following the standard procedure with NuSTARDAS v2.1.1 and CALDB v20210210 and use cleaned event data. For event selection, we use an extraction circle of 64" (26 pixels) radius and in the 3–50 keV energy range. Barycentric correction is done using barycorr. We use PRESTO and prepfold to perform a pulsation search and to determine the spin frequency of PSR J1813−1749.

## 3.3. Chandra Data

Chandra observed PSR J1813−1749 using the ACIS-S detector in continuous clocking (CC) mode on 2021 February 10 (ObsID 23545) and 2021 June 23 (ObsID 23546) for 20 ks on each date (see Table 2). We reprocess these data following the standard procedure with chandra_repro of the Chandra Interactive Analysis of Observations (CIAO) package version 4.14 and Calibration Database 4.9.6 (Fruscione et al. 2006). We extract events from the one-dimensional CC data sets along a 3" (3.1 pixels) length centered on the pulsar position (see Section 4.5) and in the 2.5–8 keV energy range. We transform the selected event time stamps from TT to TDB using the axbary tool and the pulsar position calculated at the epoch of each observation with the proper motion reported in Dzib & Rodríguez (2021). As with the analyses of NICER data for the other pulsars, we use PRESTO and prepfold to perform a pulsation search and to determine the spin frequency of PSR J1813−1749.

# 4. Results

## 4.1. PSR J0058−7218

NICER observations of PSR J0058−7218 began on 2021 June 1, and we are able to obtain a phase-connected timing model using data through 2022 January 25. Figure 2 shows timing residuals of the 15 TOAs used to obtain our best-fit timing model, which is given in Table 3. While the initial timing model only has a precision of 10% in $\dot{v}$ (Maitra et al. 2021), the new longer-time span data yield a precision of 0.002% in $\dot{v}$, as well as a $\ddot{v}$ and thus braking index $n = 50 \pm 2$.

Figure 3 shows the pulse profile at 0.8–5 keV from the combined NICER observations. The pulse profile can be fit by a Gaussian with an FWHM of 0.14 in phase or 3 ms. Pulsations are only weakly detected using NICER at 0.3–0.8 keV and are not detected above background at 5–10 keV due to contamination from the underlying SNR.

No clear spin-up glitches have been detected in NICER data so far, although we can estimate a glitch wait time (albeit applicable for pulsars with $|\dot{v}| < 3 \times 10^{-11} \text{ Hz s}^{-1}$; Fuentes et al. 2017) of $1/(420 \text{ Hz}^{-1}|\dot{v}|) \approx 15$ months, which is longer than our current time span of observations. On the other hand, when we extrapolate the spin frequency using the NICER timing model back to the epoch of the XMM-Newton observation on 2020 March 15 (MJD 58,924), the result is 57.2 µHz higher than that reported in Maitra et al. (2021), with
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Timing residuals of PSR J0058−7218 from a best fit of NICER pulse TOAs using the timing model given in Table 3. Errors are $1\sigma$ uncertainty.</figcaption></figure>

<table><caption>**Table 3**
Timing Parameters of PSR J0058−7218</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R.A. α (J2000)</td><td>00h58m16s85</td></tr><tr><td>Decl. δ (J2000)</td><td>−72°18′05″60</td></tr><tr><td>Solar system ephemeris</td><td>DE405</td></tr><tr><td>Range of dates (MJD)</td><td>59,366−59,604</td></tr><tr><td>Epoch t<sub>0</sub> (MJD TDB)</td><td>59,408</td></tr><tr><td>Frequency ν (Hz)</td><td>45.940434278(6)</td></tr><tr><td>Frequency 1st derivative ṽ (Hz s<sup>−1</sup>)</td><td>−6.2324(1) × 10<sup>−11</sup></td></tr><tr><td>Frequency 2nd derivative ṽ (Hz s<sup>−2</sup>)</td><td>4.2(2) × 10<sup>−21</sup></td></tr><tr><td>rms residual (μs)</td><td>352</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>9.94/11</td></tr><tr><td>Number of TOAs</td><td>15</td></tr></tbody></table>

**Note.** Numbers in parentheses are the $1\sigma$ uncertainty in the last digit. The position is from a Chandra ACIS-S image (MJD 56,332), with a 90% confidence level uncertainty of 0.6 (Maitra et al. 2015).

an uncertainty of 0.4 µHz from the XMM-Newton measurement and 0.2 µHz from our $\dot{v}$ uncertainty. This difference in spin frequency could be due to one or more spin-up glitches occurring in the 15 months between the times of the XMM-Newton and NICER observations. In fact, the large braking index ($n=50$) suggests recovery from a recent large glitch, similar to the behavior seen in PSR J0537−6910. Detection of future large glitches in PSR J0058−7218 would validate this possibility.

## 4.2. PSR J0537−6910

PSR J0537−6910 has been observed by NICER since early in the start of the mission in mid-2017. Data on PSR J0537−6910 from 2017 August 17 to 2020 April 25 and their timing results, including eight glitches during this period, are reported in Ho et al. (2020b). Data from 2020 May 12 to October 29 and their results, including three more glitches during this period,

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Pulse profile of PSR J0058−7218 from 153 ks of NICER data at 0.8−5 keV. Errors are $1\sigma$ uncertainty. Two rotation cycles are shown for clarity.</figcaption></figure>

are reported in Abbott et al. (2021b). Here, we report on data from 2020 November 10 to 2022 February 17. During this period, four new glitches were detected. The timing model and glitch parameters are given in Tables 4 and 5, respectively. We use the previous naming convention where each segment is separated by a glitch and is labeled by glitch number, with segment 1 occurring after glitch 1, which is the first NICER-detected glitch. Note that glitch 8 was first reported in Ho et al. (2020b), but revised timing parameters for segment 8 and glitch 8 were given in Abbott et al. (2021b) after the accumulation of more data for the segment; similarly, here we revise the parameters first given in Abbott et al. (2021b) for segment 11 and glitch 11. We find the epoch of glitch 15 to be MJD 59,522, even though a glitch epoch between MJD 59,529 and 59,556 produces similar fit results; we favor the earlier epoch since periodicity analyses of an observation on MJD 59,529 suggest that the glitch had already occurred by this date.

Figure 4 shows $\dot{v}$ for each segment, as well as interglitch $\dot{v}$ values measured using RXTE from Antonopoulou et al. (2018). A simple linear fit to only the NICER set of $\dot{v}$ gives $\dot{v} = (−8.2 \pm 0.5) \times 10^{-22}$ Hz s$^{-2}$ and a long-term braking index $n = −1.27 \pm 0.08$, while a fit to the entire NICER and RXTE set of $\dot{v}$ gives $\dot{v} = (−7.92 \pm 0.06) \times 10^{-22}$ Hz s$^{-2}$ and $n = −1.234 \pm 0.009$; all are in agreement with values found in Ho et al. (2020b) and values found using only RXTE data in Antonopoulou et al. (2018) and Ferdman et al. (2018). Thus, there is not strong evidence for a change in the long-term braking index of PSR J0537−6910.

As discussed in Ho et al. (2020b; see also Middleditch et al. 2006; Andersson et al. 2018; Antonopoulou et al. 2018; Ferdman et al. 2018), the short-term spin-down behavior (i.e., the behavior in between glitches) is much different than the long-term behavior. In particular, the interglitch braking index $n_{ig}$ is non-negative and much greater than the canonical value of 3 that is the result of conventional spin-down by electromagnetic dipole radiation at constant magnetic field and moment of inertia. This is illustrated in Figure 5, which
<table><caption>**Table 4**
Timing Parameters of PSR J0537−6910</caption><thead><tr><th>Segment</th><th>Epoch (MJD)</th><th>Start (MJD)</th><th>End (MJD)</th><th>ν (Hz)</th><th>(10<sup>-10</sup> Hz s<sup>-1</sup>)</th><th>(10<sup>-20</sup> Hz s<sup>-2</sup>)</th><th>n<sub>ig</sub></th><th>rms (μs)</th><th>χ<sup>2</sup>/ dof</th><th>TOAs</th></tr></thead><tbody><tr><td>11</td><td>59,195</td><td>59,107.7</td><td>59,283.4</td><td>61.904295359(1)</td><td>-1.997588(2)</td><td>0.59(2)</td><td>9.2(3)</td><td>179.9</td><td>8.5</td><td>31</td></tr><tr><td>12</td><td>59,318</td><td>59,286.9</td><td>59,349.9</td><td>61.902180422(4)</td><td>-1.99789(2)</td><td>0.7(4)</td><td>11(6)</td><td>82.4</td><td>2.6</td><td>10</td></tr><tr><td>13</td><td>59,399</td><td>59,352.6</td><td>59,446.6</td><td>61.900794415(3)</td><td>-1.997802(9)</td><td>1.6(1)</td><td>25(2)</td><td>132.4</td><td>4.6</td><td>17</td></tr><tr><td>14</td><td>59,487</td><td>59,461.5</td><td>59,518.4</td><td>61.899292039(3)</td><td>-1.99828(3)</td><td>[1]<sup>a</sup></td><td>...</td><td>72.7</td><td>3.0</td><td>6</td></tr><tr><td>15</td><td>59,578</td><td>59,529.6</td><td>59,626.6</td><td>61.897742999(2)</td><td>-1.998105(6)</td><td>1.04(8)</td><td>16(1)</td><td>110.8</td><td>3.2</td><td>20</td></tr></tbody></table>

**Note.** Columns are interglitch segment number, timing model epoch, segment start and end dates, spin frequency and its first two time derivatives, interglitch braking index, timing model residual, goodness-of-fit measure, and number of TOAs. Numbers in parentheses are the $1\sigma$ uncertainty in the last digit. Segments 1–7 are in Ho et al. (2020b), and segments 8–10 are in Abbott et al. (2021b). The position of R.A. = 05h37m47s416, Decl. = −69°10′19″88 (J2000) is from a Chandra ACIS-I image (MJD 51,442), with a $1\sigma$ uncertainty of ~0.6 (Townsley et al. 2006). The solar system ephemeris used is DE421.

a $\nu$ is fixed at $10^{-20}$ Hz s$^{-2}$.

<table><caption>**Table 5**
Glitch Parameters of PSR J0537−6910</caption><thead><tr><th>Glitch</th><th>Glitch Epoch (MJD)</th><th>ΔΦ (cycle)</th><th>Δν (μHz)</th><th>(10<sup>-13</sup> Hz s<sup>-1</sup>)</th><th>(10<sup>-20</sup> Hz s<sup>-2</sup>)</th></tr></thead><tbody><tr><td>11</td><td>59,103(5)</td><td>0.5(6)</td><td>33.9(4)</td><td>-1(1)</td><td>-3(2)</td></tr><tr><td>12</td><td>59,285(2)</td><td>-0.26(1)</td><td>7.872(8)</td><td>-0.94(3)</td><td>...</td></tr><tr><td>13</td><td>59,351(2)</td><td>0.51(2)</td><td>12.27(3)</td><td>-0.8(2)</td><td>0.9(4)</td></tr><tr><td>14</td><td>59,454(8)</td><td>0.31(1)</td><td>16.60(1)</td><td>-1.71(4)</td><td>...</td></tr><tr><td>15</td><td>59,522(8)</td><td>-0.30(2)</td><td>22.08(1)</td><td>-0.61(3)</td><td>...</td></tr></tbody></table>

**Note.** Columns are glitch number, glitch epoch, change in rotation phase, and changes in spin frequency and its first two time derivatives at each glitch. Numbers in parentheses are the $1\sigma$ uncertainty in the last digit. Glitches 1–7 are in Ho et al. (2020b), and glitches 8–10 are in Abbott et al. (2021b).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolution of the spin frequency time derivative $\dot{\nu}$ of PSR J0537−6910. $\dot{\nu}$ are measured by fitting a timing model to TOAs in each interglitch segment (see Table 4 and Ho et al. 2020b; Abbott et al. 2021b for NICER and Table 1 of Antonopoulou et al. 2018 for RXTE). Errors are the $1\sigma$ uncertainty. Dashed line shows a linear fit of NICER and RXTE data with the best-fit $\dot{\nu} = -7.92 \times 10^{-22}$ Hz s$^{-2}$.</figcaption></figure>

shows $n_{ig}$ measured using NICER and RXTE, with the latter values taken from Antonopoulou et al. (2018), and the time since the last glitch is the epoch of the segment minus the

epoch of the corresponding glitch (e.g., time since glitch 11 = 59195–59103 = 92 days). It is clear that large values of $n_{ig}$ are measured for short times after a glitch and that small $n_{ig}$ are measured after long post-glitch times. In other words, it appears there is recovery back to a rotational behavior that is characterized by a braking index $\lesssim 7$ after disruption by a glitch. Fits to an exponential decay yield decay timescales of 19–44 days, with a longer timescale for a lower asymptotic braking index (Ho et al. 2020b). Braking indices of 5 and 7 are expected for spin-down by gravitational-wave quadrupole and r-mode emission, respectively (see Ho et al. 2020b; Abbott et al. 2021b, and references therein).

Figure 6 shows the glitch parameters $\Delta\nu$ and $\Delta\dot{\nu}$ for the 15 glitches measured so far using NICER (see Table 5 and Ho et al. 2020b; Abbott et al. 2021b). Note that the alternating sizes of $\Delta\nu$ and the pairing of $\Delta\dot{\nu}$ seen in the first eight glitches do not continue in more recent glitches.

Glitch activity can be characterized by the parameter $A_g \equiv \sum_i (\Delta\nu/\nu)_i / T_{\text{obs}}$, where the summation is over each glitch $i$ and $T_{\text{obs}}$ is the time over which the pulsar is monitored (McKenna & Lyne 1990). For glitches detected using NICER and $T_{\text{obs}} = 4.5$ yr, we find $\sum_i \Delta\nu_i = (254.6 \pm 0.6)\mu\text{Hz}$ and $A_g = (9.15 \pm 0.02) \times 10^{-7} \text{yr}^{-1}$. Figure 7 plots the cumulative fractional glitch magnitude $\Delta\nu/\nu$ over the RXTE and NICER eras. Combining RXTE and NICER glitches produces an activity parameter $A_g = (8.918 \pm 0.009) \times 10^{-7} \text{yr}^{-1}$.

The glitches of PSR J0537−6910 are unique in that the time to the next glitch is correlated with the size of the preceding glitch (see also Middleditch et al. 2006; Antonopoulou et al. 2018; Ferdman et al. 2018; Ho et al. 2020b). This is illustrated in Figure 8. The correlation can be fit by time to the next
<figure><img src="image_8.png" /></figure>

**Figure 5.** Interglitch braking index $n_{ig}$ of PSR J0537−6910 calculated from the spin parameters of each segment between glitches as a function of time since the last glitch. Large and small circles denote NICER and RXTE values, respectively (from here and Antonopoulou et al. 2018; Ho et al. 2020b; Abbott et al. 2021b). Errors in $n_{ig}$ are $1\sigma$. Horizontal dotted lines indicate braking indices $n=3$, 5, and 7, which are expected for pulsar spin-down by electromagnetic dipole radiation, gravitational wave-emitting mountain, and gravitational wave-emitting r-mode oscillation, respectively.

<figure><img src="image_9.png" /></figure>

**Figure 6.** Glitch $\Delta\nu$ (top) and $\Delta\dot{\nu}$ (bottom) as functions of time (see Table 5 and Ho et al. 2020b; Abbott et al. 2021b). Errors in $\Delta\nu$ are the $1\sigma$ uncertainty.

$\text{glitch} = 49.2 \text{ d} (\Delta\nu/10 \mu\text{Hz}) + (25 \pm 18)\text{d}$, in agreement with that found in Ho et al. (2020b). This correlation enables the prediction of when glitches will occur in PSR J0537−6910. In particular, glitch 16 should occur on 2022 March 18, although there is a large uncertainty ($\pm 26$ d), due in part to the uncertain epoch of glitch 15.

<figure><img src="image_10.png" /></figure>

**Figure 7.** Fractional glitch magnitude $\Delta\nu/\nu$ of PSR J0537−6910 shown as a cumulative sum over each previous glitch. RXTE values are from Table 2 of Antonopoulou et al. (2018). Dashed line indicates a line with a slope of $8.918 \times 10^{-7} \text{ yr}^{-1}$, which is the glitch activity $A_g \equiv \sum_i (\Delta\nu/\nu)_i / T_{\text{obs}}$, where $T_{\text{obs}}$ is the time over which the pulsar is monitored. NICER values are offset by $\Delta\nu/\nu = 16.7 \times 10^{-6}$, which is the extrapolated value of the RXTE-only glitch activity at the epoch of NICER segment 0 at MJD 58,020.

<figure><img src="image_11.png" /></figure>

**Figure 8.** Correlation between time interval to the next glitch $\Delta T$ and size of glitch $\Delta\nu$ of PSR J0537−6910. Large and small circles denote NICER and RXTE values, respectively (from here and Antonopoulou et al. 2018; Ho et al. 2020b; Abbott et al. 2021b). Errors in $\Delta\nu$ are $1\sigma$. The vertical dotted line indicates the size of NICER glitch 15, which is the most recent glitch (on 2021 November 4) and for which time to next glitch is not known yet. Dashed line shows the linear fit result $\Delta T = 49.2 \text{ d} (\Delta\nu/10 \mu\text{Hz}) + 25 \text{ d}$.

### 4.3. PSR J1101−6101

NICER observations of PSR J1101−6101 began on 2020 April 1, and we are able to obtain a phase-connected timing model using data through 2021 December 16. Figure 9 shows timing residuals
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Timing residuals of PSR J1101-6101 from a best fit of NICER pulse TOAs with the timing model given in Table 6. Errors are the $1\sigma$ uncertainty.</figcaption></figure>

<table><caption>**Table 6**
Timing Parameters of PSR J1101-6101</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R.A. α (J2000)</td><td>11h01m44s915</td></tr><tr><td>Decl. δ (J2000)</td><td>-61°01′38″66</td></tr><tr><td>Solar system ephemeris</td><td>DE405</td></tr><tr><td>Range of dates (MJD)</td><td>58,940–59,564</td></tr><tr><td>Epoch t<sub>0</sub> (MJD TDB)</td><td>59,171</td></tr><tr><td>Frequency ν (Hz)</td><td>15.9230402557(5)</td></tr><tr><td>Frequency 1st derivative ṷ (Hz s<sup>-1</sup>)</td><td>-2.26504(4) × 10<sup>-12</sup></td></tr><tr><td>rms residual (ms)</td><td>1.60</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>7.80/9</td></tr><tr><td>Number of TOAs</td><td>12</td></tr></tbody></table>

**Note.** Number in parentheses are the $1\sigma$ uncertainty in the last digit. The position is from a Chandra ACIS-I image (MJD 56,211), with a 90% confidence level uncertainty of 0.64 (Pavan et al. 2014).

of the 12 TOAs used to obtain our best-fit timing model, which is given in Table 6, and Figure 10 shows the 1.5–10 keV pulse profile from the combined NICER observations. Our measured spin-down rate of ṷ = (-2.26504 ± 0.00004) × 10⁻¹² Hz s⁻¹ is consistent with and significantly improves upon the precision of the previously measured incoherent timing model value of ṷ = (-2.17 ± 0.13) × 10⁻¹² Hz s⁻¹ from Halpern et al. (2014). The addition of ṷ to the timing model yields a fit improvement of only Δχ²=0.5 and an unconstrained ṷ = (-1 ± 18) × 10⁻²⁴ Hz s⁻².

The time span of our NICER timing model overlaps with a NuSTAR observation taken on 2020 November 20. We use the NICER timing model to extract the pulsed emission from the NuSTAR data. The resulting pulse profiles are shown in Figure 10. It is clear that the two sets of pulse profiles bear strong resemblance in shape and have consistent alignment. More detailed analysis of the NuSTAR data can be found in Klingler et al. (2022).

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Pulse profile of PSR J1101-6101 from 370 ks of NICER data at 1.5–10 keV (top), while lower panels show pulse profiles from 136 ks of NuSTAR data at 1.6–10, 10–20, 20–60, and 1.6–60 keV. Errors are $1\sigma$ uncertainty. Two rotation cycles are shown for clarity.</figcaption></figure>

### 4.4. PSR J1412+7922

NICER observations of PSR J1412+7922 began on 2017 September 15, and Bogdanov et al. (2019) report a 1 yr phase-connected timing model using NICER data through 2018 October 3. Mereghetti et al. (2021) report a 3.4 yr phase-connected timing model using NICER data through 2021 February 26. Here, we extend the timing model to over 4.4 yr using NICER data through 2022 February 8. Our analysis procedure yields 138 TOAs that are barycentered, but not corrected for proper motion, with respect to the pulsar’s position as measured by Halpern & Gotthelf (2015). While we are able to successfully obtain a phase-connected timing model, the resulting timing residuals over the 4.4 yr span of data display a systematic wave-like behavior on a timescale of several months, even with the addition of ṷ and ṷ̇.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Timing residuals of PSR J1412+7922 from a best fit of NICER pulse TOAs with the timing model given in Table 7. Errors are $1\sigma$ uncertainty.</figcaption></figure>

<table><caption>**Table 7**
Timing Parameters of PSR J1412+7922</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R.A. α (J2000)</td><td>14<sup>h</sup>12<sup>m</sup>56<sup>s</sup>05(3)</td></tr><tr><td>Decl. δ (J2000)</td><td>+79°22′03″68(7)</td></tr><tr><td>Solar system ephemeris</td><td>DE405</td></tr><tr><td>Range of dates (MJD)</td><td>58,014.2–59,616.9</td></tr><tr><td>Epoch t<sub>0</sub> (MJD TDB)</td><td>58,750</td></tr><tr><td>Frequency ν (Hz)</td><td>16.8921082712(1)</td></tr><tr><td>Frequency 1st derivative ṷ (Hz s<sup>-1</sup>)</td><td>-9.40547(4) × 10<sup>-13</sup></td></tr><tr><td>Frequency 2nd derivative ṷ̇ (Hz s<sup>-2</sup>)</td><td>-2.83(3) × 10<sup>-23</sup></td></tr><tr><td>rms residual (ms)</td><td>1.449</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>303.7/132</td></tr><tr><td>Number of TOAs</td><td>138</td></tr></tbody></table>

**Note.** Numbers in parentheses are the $1\sigma$ uncertainty in the last digit. The position epoch is the same as the timing model epoch. No proper motion is assumed.

terms in the timing model (see, e.g., Figure 1 of Mereghetti et al. 2021). Therefore, we consider a timing model that includes the pulsar’s position as a fit parameter. To construct such a model, we calculate 138 spacecraft topocentric, not barycentric, TOAs from the same data set, and we fit these TOAs using PINT (Luo et al. 2021). Figure 11 shows timing residuals of the TOAs used to obtain our best-fit timing model, which is given in Table 7. While the timing model of Bogdanov et al. (2019) only required a spin-down ṷ term to achieve a timing residual of 1.36 ms, our four times longer time baseline requires ṷ̇ for a comparable timing residual of 1.45 ms. The addition of proper motion as a fit parameter to the timing model only yields a fit improvement of $\Delta\chi^2 = 13$.

It is notable that the position we measure using timing data is different from that measured by Halpern & Gotthelf (2015) using Chandra HRC-I imaging data from 2007 and 2014. In particular, their position from the 2014 data (MJD 56,749), which is much longer and closer to the time of our data than the 2007 data, is $(\alpha, \delta) = (14^{h}12^{m}55^{s}815 \pm 0^{s}011, +79^{\circ}22'03''697 \pm 0''030)$,

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Spin frequency $\nu$ of PSR J1813−1749 (top) and the difference between the best-fit linear model and data (bottom) as functions of time. Measurements of $\nu$ are made using XMM-Newton (triangles), Chandra (squares), Green Bank Telescope in radio (stars), NuSTAR (cross), and NICER (circle). Dashed line shows a linear fit of all $\nu$ measurements with best-fit $\dot{\nu} = -6.3442 \times 10^{-11}$ Hz s$^{-1}$.</figcaption></figure>

and thus our position differs by $\Delta\alpha = +0^{\circ}23 \pm 0^{\circ}04$ and $\Delta\delta = -0^{\circ}02 \pm 0^{\circ}10$. The implied proper motion of PSR J1412+7922 between its 2014 and 2019 positions is $\mu_{\alpha} \cos \delta = +120 \pm 20$ mas yr$^{-1}$ and $\mu_{\delta} = -3 \pm 20$ mas yr$^{-1}$, which is also at odds with $\mu_{\alpha} \cos \delta = -40 \pm 30$ mas yr$^{-1}$ and $\mu_{\delta} = -56 \pm 21$ mas yr$^{-1}$ as measured by Halpern & Gotthelf (2015) from changes of the 2007 and 2014 positions. Differences in the two techniques used, i.e., imaging versus timing, could contribute to the different measured positions of PSR J1412+7922, including the impact of timing noise on the timing model fits. A future Chandra observation may resolve these discrepancies.

### 4.5. PSR J1813−1749

Our analysis of PSR J1813−1749 only considers the spin frequencies measured using various radio and X-ray telescopes over a 12 yr time span. In particular, we fit the evolution of $\nu$ with a linear decline in order to determine $\dot{\nu}$. This is done because the data for PSR J1813−1749 are sparse and a phase-connected timing analysis like that done in Ho et al. (2020a) requires significant NICER observing time.

We obtain three new measurements of the spin frequency of PSR J1813−1749. In particular, we find $\nu = 22.3537223 \pm 0.0000008$ Hz on MJD 58,202.39 from the 2018 NuSTAR observation, and $\nu = 22.3479859 \pm 0.0000039$ Hz on MJD 59, 255.45 and $\nu = 22.3471999 \pm 0.0000039$ Hz on MJD 59, 388.14 from the two 2021 Chandra CC observations. Meanwhile, the previous spin frequency measurements are from 2009 and 2011 using XMM-Newton; from 2012 using Chandra (Halpern et al. 2012); the 2012 and 2015 measurements at radio wavelengths (Camilo et al. 2021); and a 2019 measurement using NICER (Ho et al. 2020a). These are all shown in Figure 12. Fitting a simple linear decline in spin frequency yields a best-fit spin-down rate $\dot{\nu} = (-6.34425 \pm 0.00072) \times 10^{-11}$ Hz s$^{-1}$,
<table><caption>**Table 8**
Incoherent Timing Parameters of PSR J1813−1749</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R.A., α (J2000)</td><td>18<sup>h</sup>13<sup>m</sup>35<sup>s</sup> 173</td></tr><tr><td>Decl., δ (J2000)</td><td>−17°49′57″ 75</td></tr><tr><td>Solar system ephemeris</td><td>DE405</td></tr><tr><td>Range of dates (MJD)</td><td>54,918.14−59,388.14</td></tr><tr><td>Epoch t<sub>0</sub> (MJD TDB)</td><td>58,681.04</td></tr><tr><td>Frequency ν (Hz)</td><td>22.35108384(2)</td></tr><tr><td>Frequency 1st derivative ˙ν (Hz s<sup>−1</sup>)</td><td>−6.3442(7) × 10<sup>−11</sup></td></tr><tr><td>Proper motion µ<sub>α</sub> cos δ (mas yr<sup>−1</sup>)</td><td>−5.0</td></tr><tr><td>Proper motion µ<sub>δ</sub> (mas yr<sup>−1</sup>)</td><td>−13.2</td></tr></tbody></table>

**Note.** Numbers in parentheses are the $1\sigma$ uncertainty in the last digit. The position is from VLA data (MJD 58,119) with uncertainties of ~0.009 and ~0.011, and the proper motion has uncertainties of 3.7 and 6.7 mas yr$^{-1}$ in $\mu_\alpha \cos \delta$ and $\mu_\delta$, respectively (Dzib & Rodríguez 2021).

which is in agreement with that determined in Ho et al. (2020a) because the fit is driven by the two most precise measurements by XMM-Newton in 2009 and NICER in 2019.

Table 8 presents the resulting 12 yr timing model, albeit one that is incoherent in contrast to the 37 d phase-connected timing model presented in Ho et al. (2020a). Residuals from the timing model suggest that PSR J1813−1749 undergoes glitches with sizes as large as a few tens of µHz (see the bottom panel of Figure 12), and the spin-down rate suggests a glitch wait time of ~14 months (Fuentes et al. 2017; see Section 4.1). Such large glitches are uncommon but can occur multiple times in a single pulsar over a period of a few years, such as in the Vela pulsar and PSR J0537−6910 (see, e.g., Figure 6). Analysis of the 2019 NICER data of PSR J1813−1749 indicates a spin-down rate of $-6.428 \times 10^{-11}$ Hz s$^{-1}$ during the 37 day observing period and a possible glitch with $\Delta\nu \approx 3$ µHz near the end of the observation (see Ho et al. 2020a for more discussion). Note that the spin frequency difference of 786 µHz between the two 2021 Chandra observations separated by 133 days implies a spin-down rate of $-6.84 \times 10^{-11}$ Hz s$^{-1}$.

### 4.6. PSR J1849−0001

NICER observations of PSR J1849−0001 began on 2018 February 13, and Bogdanov et al. (2019) report a 1.5 yr phase-connected timing model using NICER data through 2018 September 29 as well as a Swift observation on 2017 March 19. Here, we extend the timing model to nearly 4.7 yr using NICER data through 2021 November 16. Figure 13 shows timing residuals of the 41 TOAs used to obtain our best-fit timing model, which is given in Table 9. While the timing model of Bogdanov et al. (2019) only required a spin-down $\dot{\nu}$ term to achieve a timing residual of 602 µs, our three times longer time baseline requires up to $\ddot{\nu}$ with a comparable timing residual of 569 µs.

With our time baseline of several years, we attempt to measure a change of the position of PSR J1849−0001, such as could be caused by the pulsar's proper motion, using spacecraft topocentric TOAs. In this case, the best-fit timing model produces an rms residual of 513 µs, an F-test probability of 0.03, and position differences of $\Delta\alpha = -0.008 \pm 0.003$ and $\Delta\delta = +0.016 \pm 0.011$ compared to the timing model given in Table 9. Since the fit improvement is marginal and the resulting position and other timing parameters, e.g., $\nu$ and $\dot{\nu}$, are all within $1\sigma$ uncertainties of

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Timing residuals of PSR J1849−0001 from a best fit of NICER (circles) and one Swift (cross) pulse TOAs, with the timing model given in Table 9. Errors are $1\sigma$ uncertainty.</figcaption></figure>

<table><caption>**Table 9**
Timing Parameters of PSR J1849−0001</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R.A., α (J2000)</td><td>18<sup>h</sup>49<sup>m</sup>01<sup>s</sup> 632</td></tr><tr><td>Decl., δ (J2000)</td><td>−00°01′17″ 45</td></tr><tr><td>Solar system ephemeris</td><td>DE421</td></tr><tr><td>Range of dates (MJD)</td><td>57,832.1−59,533.9</td></tr><tr><td>Epoch t<sub>0</sub> (MJD TDB)</td><td>58,682</td></tr><tr><td>Frequency ν (Hz)</td><td>25.9586535424(3)</td></tr><tr><td>Frequency 1st derivative ˙ν (Hz s<sup>−1</sup>)</td><td>−9.53597(1) × 10<sup>−12</sup></td></tr><tr><td>Frequency 2nd derivative ˙ν (Hz s<sup>−2</sup>)</td><td>8.3(1) × 10<sup>−23</sup></td></tr><tr><td>Frequency 3rd derivative ˙ν (Hz s<sup>−3</sup>)</td><td>−2.25(3) × 10<sup>−30</sup></td></tr><tr><td>Frequency 4th derivative ˙ν (Hz s<sup>−4</sup>)</td><td>6.6(5) × 10<sup>−38</sup></td></tr><tr><td>rms residual (µs)</td><td>569</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>90.7/35</td></tr><tr><td>Number of TOAs</td><td>41</td></tr></tbody></table>

**Note.** Numbers in parentheses are the $1\sigma$ uncertainty in the last digit. The position is from a Chandra HRC-S image (MJD 55,885), with a 90% confidence level uncertainty of $0.016$ (Kuiper & Hermsen 2015).

the values shown in Table 9, we do not report this alternative timing model here.

### 5. Discussion

Using long-term monitoring observations made by NICER over the past several years, we have calculated rotation phase-connected timing models for five pulsars that are only known to be visible at X-ray energies, as well as determined the long-term spin-down rate of the highly energetic PSR J1813−1749 by making use of recent Chandra and NuSTAR observations. These timing models have time spans that greatly exceed those of previous models, thus providing more reliable characterizations of the spin properties of these rapidly rotating, mostly young pulsars. Continued monitoring of these pulsars is needed for searches at other energies (e.g., gamma-ray energies using
Fermi; see below) and crucially for gravitational-wave searches of more sensitive data to be obtained in upcoming observing runs (see Section 1). The high-cadence monitoring of PSR J0537−6910 also enables detection of its glitches, which provides a unique probe of superfluidity.

## 5.1. Gamma-Ray Searches

The timing models presented here enable us to search for pulsed gamma-ray emission from PSR J0058−7218, PSR J1101−6101, PSR J1412+7922, and PSR J1849−0001 in Fermi LAT data. We neglect PSR J0537−6910 and PSR J1813−1749 because the complex timing behaviors of these pulsars and strong gamma-ray emission from nearby pulsar wind nebula and SNRs make detection of a gamma-ray pulsar component very difficult (see, e.g., Fermi LAT Collaboration et al. 2015). We conduct two types of searches here. The first follows the methodology described in Smith et al. (2019), and the second follows that described in Kuiper et al. (2018).

For the first set of searches, we use our timing models to gamma-ray phase-fold each pulsar six times, i.e., using three values of the gamma-ray photon-weighting parameter $\mu_w$ described by Bruel (2019) and used in Smith et al. (2019) and either LAT data during only the epoch range of each timing model or all 12.6 yr of currently available LAT data. We find no statistically significant ($>1\sigma$) deviation from a uniform phase distribution, and we place limits on phase-integrated flux above 100 MeV of $7.2 \times 10^{-12} \text{ erg s}^{-1} \text{ cm}^{-2}$ for PSR J1101−6101, $1.8 \times 10^{-12} \text{ erg s}^{-1} \text{ cm}^{-2}$ for PSR J1412+7922, and $1.2 \times 10^{-11} \text{ erg s}^{-1} \text{ cm}^{-2}$ for PSR J1849−0001. These limits are based on pulsar sensitivity estimates made using the second Fermi pulsar catalog (see in particular Section 8.2 of Abdo et al. 2013), but updated for the 4FGL-DR3 12 yr data set (T. Burnett, private communication). For PSR J0058−7218, we measure a total flux above 100 MeV of $(1.6 \pm 0.3) \times 10^{-12} \text{ erg s}^{-1} \text{ cm}^{-2}$, which is likely dominated by emission from other nearby sources such as the star-forming region NGC 346 (see also Maitra et al. 2021). For PSR J1101−6101 and PSR J1849−0001, detections of their pulsations could be hindered by several LAT sources that lie within the 1° LAT angular resolution above 500 MeV of each pulsar. The nondetection of PSR J1412+7922, despite its high $\dot{E}$ relative to the LAT flux sensitivity at the pulsar's high Galactic latitude and low gamma-ray background, could be due to an unfavorable viewing geometry, as illustrated by Johnston et al. (2020) and the angles inferred from Mereghetti et al. (2021). Importantly, the methodology of Smith et al. (2019) may not be suitable for discovering pulsed emission from soft gamma-ray pulsars.

For the second set of searches using the methodology of Kuiper et al. (2018), we barycenter LAT timing data of PSR J0058−7218 collected during the 8 month period of the model. We apply event selections similar to those for the LAT timing analysis of PSR J1846−0258 presented in Kuiper et al. (2018). We then fold in phase the data using the timing parameters specified in Table 3. For energies above ~500 MeV, we obtain a pulse phase distribution with a narrow single-peaked shape that is coincident in phase with the prominent pulse in the X-ray band. However, the overall $Z_n^2$-based significance varies between 2.2–2.8$\sigma$ for $n$ between 3 and 8. Below ~500 MeV, the distribution is statistically flat. Further LAT exposure and NICER monitoring are needed to determine whether this possible

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Pulse profiles of PSR J1101−6101 from Fermi LAT data at 0.3–3 GeV (solid histogram) and from NICER data at 1.5–10 keV (dashed curve with arbitrary normalization; see also Figure 10). Errors are $1\sigma$ uncertainty. Two rotation cycles are shown for clarity.</figcaption></figure>

detection of pulsed high-energy gamma-rays from PSR J0058−7218 is real or a statistical fluctuation.

For PSR J1101−6101, we use a similar strategy to search LAT data collected during the 1.7 yr period of the ephemeris given in Table 6. We obtain pulse phase distributions for logarithmically binned energy bands of 30–300 MeV, 0.3–3 GeV, and 3–30 GeV. The 30–300 MeV and 3–30 GeV bands do not show any evidence for a significant pulsed signal. In the intermediate band, we find a potential signal at a significance of about $3.2\sigma$ by adopting a $Z_3^2$-test with the bulk of the enhanced emission (see Figure 14) aligned with the X-ray pulsations detected by NICER and NuSTAR (see Figure 10). Furthermore, a restricted frequency search in the 0.3–3 GeV band, near the NICER ephemeris prediction and keeping the frequency derivative fixed at the value given in Table 6, shows only one prominent, but weak, maximum coinciding with the prediction of the NICER ephemeris. Further Fermi LAT exposure and X-ray timing in the future are required to corroborate this tentative detection of pulsed gamma-ray emission from PSR J1101−6101.

For PSR J1412+7922 and PSR J1849−0001, we again do not detect pulsed emission using Fermi LAT data. On the other hand, we detect pulsations of PSR J1849−0001 up to ~150 keV using the Fermi GBM NaI detectors, with a $> 4\sigma$ significance for the 100–150 keV band (L. Kuiper et al., in preparation). Note that even though we detect none (or potentially two) of the six sources in pulsed gamma-rays, this does not necessarily mean that each is gamma-ray-quiet, given the difficulties described that hamper possible detections of their pulsations.

## 5.2. Superfluidity and Glitch Predictability from PSR J0537−6910

In the conventional two-component model for glitches (Anderson & Itoh 1975; Alpar et al. 1984), glitch activity is related to the ratio of superfluid moment of inertia $I_{sf}$ to total stellar moment of inertia $I$, such that $I_{sf}/I \ge 2\tau_c A_g$ (Link et al. 1999),
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Schematic of a glitch (see also Ray et al. 2019). The pulsar spin frequency $\nu$ (solid line) is observed to decrease at a rate $\dot{\nu}$, due to energy loss from electromagnetic radiation, while a superfluid component within the star rotating at $\nu_{sf}$ is pinned and does not spin down with the rest of the star ($\dot{\nu}_{sf} = 0$; dashed line). When the superfluid component unpins and couples to the nonsuperfluid component after time $T_g$, the superfluid transfers angular momentum ($\propto \Delta\nu_{sf}$) to the rest of the star, and a spin-up glitch $\Delta\nu$ is observed. The pulsar then continues to spin down when the superfluid becomes pinned again.</figcaption></figure>

where $A_g \equiv \sum_i (\Delta\nu/\nu)_i / T_{obs}$; for simplicity, we here neglect entrainment, which can be taken into account by a factor of order unity on the right hand side (Andersson et al. 2012; Chamel 2013). An equivalent measure of glitch activity to $A_g$ is $\dot{v}_g \equiv \sum_i \Delta\nu / T_{obs} \approx \nu A_g$ (Lyne et al. 2000) since each glitch induces only a small change in $\nu$, such that $I_{sf}/I \ge \dot{v}_g / |\dot{\nu}|$. Many works find that $2\tau A_g$ and $\dot{v}_g / |\dot{\nu}|$ can be as large as $\sim 0.01$ and use this to constrain superfluid properties and even to infer the mass of isolated pulsars (Link et al. 1999; Andersson et al. 2012; Chamel 2013; Ho et al. 2015; Fuentes et al. 2017, 2019; Montoli et al. 2020). Note that the mass inferred from glitch activity depends on the nuclear equation of state, but the uncertainty in mass is only weakly dependent; the mass uncertainty in the case of PSR J0537–6910 is dominated by the uncertainty in its age and, hence, the temperature of its superfluid/nonsuperfluid regions (see Ho et al. 2015, for details).

One might expect that pulsars with long observing times $T_{obs}$ that are seen to glitch many times ($N_g \gg 1$, where $N_g$ is the number of observed glitches) would show a correlation between the size of their glitches $\Delta\nu$ and the time between glitches $T_g$,$^{21}$ which could be due to surpassing a critical threshold or

building up or depleting an angular momentum reservoir. Such a correlation can be derived by considering the simple spin-down behavior of normal and superfluid components, as illustrated in Figure 15. A spin rate lag between these two components builds up during the time between glitches $T_g$, due to a difference in spin-down rates of the components, with $\dot{\nu}$ being the spin-down rate of the normal component and $\dot{v}_{sf} (= 0)$ being that of the superfluid. At the glitch,

$$ \Delta\nu + |\Delta\nu_{sf}| = \nu_{sf} - \nu = |\dot{\nu}|T_g, \quad (1) $$

while angular momentum conservation implies

$$ |\Delta\nu_{sf}| = \frac{I}{I_{sf}} \Delta\nu, \quad (2) $$

where $I_{sf} \ll I$ is assumed. Combining Equations (1) and (2) yields

$$ \Delta\nu = \frac{I_{sf}}{I} |\dot{\nu}| T_g, \quad (3) $$

which can also be related to glitch activity by taking $T_{obs} \approx N_g T_g$ and $\sum \Delta\nu \sim N_g \Delta\nu$, such that $\dot{v}_g / |\dot{\nu}| = I_{sf} / I$. Equation (3) can be used to constrain the superfluid moment of inertia, although in this case individual glitches are used rather than an ensemble average that the glitch activity parameters imply. Figure 16 illustrates this constraint for $\dot{v}_{ig}$ equal to $\dot{\nu}$ measured between glitches and for three values of $I_{sf} / I$ (although recall that we neglect entrainment, which would result in a scaling factor of order unity in the constraint; see above), as well as for glitch data for PSR J0537–6910 and some other pulsars. Note that another superfluid constraint can be obtained by assuming that the spin-down torque on the entire star at a glitch is equal to the torque on the normal (nonsuperfluid) component before or after the glitch, which results in

$$ \left| \frac{\Delta \dot{\nu}}{\dot{\nu}} \right| = \frac{I_{sf}}{I} \quad (4) $$

(see also Alpar et al. 1981).

Another correlation can be obtained simply by assuming that a glitch-induced change in spin-down rate $\Delta\dot{\nu}$ fully recovers with time linearly at a rate $\dot{v}_{ig}$, where $\dot{v}_{ig}$ is $\dot{\nu}$ measured between glitches (see also Akbal et al. 2017). This then implies

$$ |\Delta\dot{\nu}| = \dot{v}_{ig} T_g. \quad (5) $$

Note that if the spin-down rate does not fully recover, then the above can be replaced by $k|\Delta\dot{\nu}| = \dot{v}_{ig} T_g$, where $k < 1$. Equation (5) and glitch data from Lower et al. (2021) and for PSR J0537–6910 are plotted in Figure 17. Lower et al. (2021) show that the glitches of their 16 pulsars follow the correlation, albeit with scatter. We see here that the addition of glitches of PSR J0537–6910 extends the apparent correlation to much higher values of $|\Delta\dot{\nu}|/T_g$ and $\dot{v}_{ig}$. Also note that we do not observe a strong correlation between $\dot{v}_{ig}$ and $\delta t^2$, where $\delta t$ is the time since previous glitch, as expected from Haskell et al. (2020). However, we cannot make a firm conclusion since $\dot{v}_{ig}$ is likely strongly affected by glitch recovery effects in the first $\sim 80$ days after a glitch (see Figure 5). All the above and the issues discussed in Section 4.2 illustrate the importance of continuing X-ray timing observations and gravitational-wave

<footnote>$^{21}$ It is often assumed that the time between glitches $T_g$, the time since the previous glitch, and the (wait) time to the next glitch are all equivalent to each other. However, these times are not necessarily the same in glitch models (see, e.g., Carlin & Melatos 2021, and references therein) and do not appear to be equivalent in observations. For a few pulsars, there is a correlation between observed glitch sizes and the times to the next glitch (see, e.g., Figure 8 for PSR J0537–6910). But there is little support for a correlation between observed glitch sizes and the times to previous glitch in the case of PSR J0537–6910 (Middleditch et al. 2006; Antonopoulou et al. 2018; Ferdman et al. 2018) and in other glitching pulsars (Melatos et al. 2018; Fuentes et al. 2019; Lower et al. 2021). Nevertheless, for the simple scenario outlined here, we assume that these times are equivalent; in particular, we use the time to the next glitch as $T_g$ for plotting glitch data in this section.</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Glitch size $\Delta\nu$, interglitch spin-down rate $\dot{\nu}_{ig}$, and time to next glitch $T_g$ for PSR J0537−6910 from NICER (large circles) and RXTE (small circles), with the latter from Antonopoulou et al. (2018), and for 16 other pulsars (crosses; errors not shown since they are smaller than the symbols in most cases) from Lower et al. (2021). Dotted lines are $I_{st}/I = 0.001, 0.01, 0.1$ [see Equation (3)], where $\dot{\nu}_{ig}$ is $\dot{\nu}$ between glitches and entrainment is neglected.</figcaption></figure>

searches of PSR J0537−6910 for its impact on understanding glitches and revealing properties of dense nuclear matter.

We thank S. Dzib for providing the radio position of PSR J1813−1749 and A. Królak, S. Mastrogiovanni, M. Pitkin, K. Riles, and G. Woan for support. W.C.G.H. acknowledges support through grants 80NSSC21K0091 and 80NSSC21K1907 from NASA and Chandra award SAO GO1-22061X. Chandra grants are issued by the Chandra X-ray Center (CXC), which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. C.M.E. acknowledges support from ANID FONDECYT 1211964. S.G. acknowledges the support of the Centre National d'Etudes Spatiales (CNES). S.B. acknowledges support from NASA grant 80NSSC20K0275. D.A. acknowledges support from an EPSRC/ STFC fellowship (EP/T017325/1). M.B. is partially supported by National Science Center (NCN) Poland grants 2016/22/E/ST9/ 00037 and 2017/26/M/ST9/00978. B.H. acknowledges support from NCN Poland via SONATA BIS grant 2015/18/E/ST9/ 00577.

This research made use of data obtained from the Chandra Data Archive and the Chandra Source Catalog, and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and Sherpa. This work is supported by NASA through the NICER mission and the Astrophysics Explorers Program and uses data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. This work made extensive use of the NASA Astrophysics Data System (ADS) Bibliographic Services and arXiv.

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Interglitch $\ddot{\nu}_{ig}$, glitch size $\Delta\dot{\nu}$, and time to next glitch $T_g$ for PSR J0537−6910 from NICER (large circles) and RXTE (small circles), with the latter from Antonopoulou et al. (2018), and for 16 other pulsars (crosses, i.e., all data shown with $\ddot{\nu}_{ig} < 3 \times 10^{-21}$ Hz s$^{-2}$) from Lower et al. (2021). Dotted line is $\dot{\nu}_{ig} = |\Delta\dot{\nu}|/T_g$ [see Equation (5)].</figcaption></figure>

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include NASA and the Department of Energy (DOE) in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the High Energy Accelerator Research Organization (KEK), and the Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work is performed in part under DOE Contract DE-AC02-76SF00515.

Facilities: Chandra, Fermi, NICER, NuSTAR.

**Software:** CIAO (https://cxc.harvard.edu/ciao/), HEAsoft (https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/), PINT (https://github.com/nanograv/pint), PRESTO (https://github.com/scottransom/presto), TEMPO2 (Hobbs et al. 2006).

## ORCID iDs

Wynn C. G. Ho https://orcid.org/0000-0002-6089-6836
Lucien Kuiper https://orcid.org/0000-0002-7889-6586
Cristóbal M. Espinoza https://orcid.org/0000-0003-2481-2348
Sebastien Guillot https://orcid.org/0000-0002-6449-106X
Paul S. Ray https://orcid.org/0000-0002-5297-5278
D. A. Smith https://orcid.org/0000-0002-7833-0275
Slavko Bogdanov https://orcid.org/0000-0002-9870-2742
Michał Bejger https://orcid.org/0000-0002-4991-8213
Teruaki Eno https://orcid.org/0000-0003-1244-3100
Paolo Esposito https://orcid.org/0000-0003-4849-5092

Alice K. Harding https://orcid.org/0000-0001-6119-859X

Natalia Lewandowska https://orcid.org/0000-0003-0771-6581

Georgios Vasilopoulos https://orcid.org/0000-0003-3902-3915

## References

Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, ApJ, 922, 71

Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, ApJL, 913, L27

Abbott, R., Abbott, T. D., Acernese, F., et al. 2022a, ApJ, 932, 133

Abbott, R., Abe, H., Acernese, F., et al. 2022b, ApJ, 935, 1

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

Akbal, O., Alpar, M. A., Buchner, S., & Pines, D. 2017, MNRAS, 469, 4183

Alpar, M. A., Anderson, P. W., Pines, D., & Shaham, J. 1981, ApJL, 249, L29

Alpar, M. A., Pines, D., Anderson, P. W., & Shaham, J. 1984, ApJ, 276, 325

Anderson, P. W., & Itoh, N. 1975, Natur, 256, 25

Andersson, N., Antonopoulou, D., Espinoza, C. M., Haskell, B., & Ho, W. C. G. 2018, ApJ, 864, 137

Andersson, N., Glampedakis, K., Ho, W. C. G., & Espinoza, C. M. 2012, PhRvL, 109, 241103

Antonopoulou, D., Espinoza, C. M., Kuiper, L., & Andersson, N. 2018, MNRAS, 473, 1644

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Basu, A., Shaw, B., Antonopoulou, D., et al. 2022, MNRAS, 510, 4049

Bogdanov, S., Ho, W. C. G., Enoto, T., et al. 2019, ApJ, 877, 69

Brogan, C. L., Gaensler, B. M., Gelfand, J. D., et al. 2005, ApJL, 629, L105

Bruel, P. 2019, A&A, 622, A108

Camilo, F., Ransom, S. M., Halpern, J. P., & Roshi, D. A. 2021, ApJ, 917, 67

Carlin, J. B., & Melatos, A. 2021, ApJ, 917, 1

Chamel, N. 2013, PhRvL, 110, 011101

Chen, Y., Wang, Q. D., Gotthelf, E. V., et al. 2006, ApJ, 651, 237

De Luca, A. 2017, JPhCS, 932, 012006

Dzib, S. A., & Rodríguez, L. F. 2021, ApJ, 923, 228

Espinoza, C. M., Lyne, A. G., Stappers, B. W., & Kramer, M. 2011, MNRAS, 414, 1679

Ferdman, R. D., Archibald, R. F., Gourgouliatos, K. N., & Kaspi, V. M. 2018, ApJ, 852, 123

Fermi LAT Collaboration, Ackermann, M., Albert, A., et al. 2015, Sci, 350, 801

Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, Proc. SPIE, 6270, 62701V

Fuentes, J. R., Espinoza, C. M., & Reisenegger, A. 2019, A&A, 630, A115

Fuentes, J. R., Espinoza, C. M., Reisenegger, A., et al. 2017, A&A, 608, A131

García, F., Combi, J. A., Albacete-Colombo, J. F., et al. 2012, A&A, 546, A91

Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, Proc. SPIE, 9905, 99051H

Gotthelf, E. V., & Halpern, J. P. 2009, ApJL, 700, L158

Gotthelf, E. V., Halpern, J. P., Terrier, R., & Mattana, F. 2011, ApJL, 729, L16

Graczyk, D., Pietrzyński, G., Thompson, I. B., et al. 2020, ApJ, 904, 13

Halpern, J. P., Bogdanov, S., & Gotthelf, E. V. 2013, ApJ, 778, 120

Halpern, J. P., & Gotthelf, E. V. 2015, ApJ, 812, 61

Halpern, J. P., Gotthelf, E. V., & Camilo, F. 2012, ApJL, 753, L14

Halpern, J. P., Tomsick, J. A., Gotthelf, E. V., et al. 2014, ApJL, 795, L27

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103

Haskell, B., Antonopoulou, D., & Barenghi, C. 2020, MNRAS, 499, 161

HEASARC 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

Ho, W. C. G., Espinoza, C. M., Antonopoulou, D., & Andersson, N. 2015, SciA, 1, e1500578

Ho, W. C. G., Espinoza, C. M., Arzoumanian, Z., et al. 2020b, MNRAS, 498, 4605

Ho, W. C. G., Guillot, S., Saz Parkinson, P. M., et al. 2020a, MNRAS, 498, 4396

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

Johnston, S., Smith, D. A., Karastergiou, A., & Kramer, M. 2020, MNRAS, 497, 1957

Klingler, et al. 2022, ApJ, submitted

Kuiper, L., & Hermsen, W. 2009, A&A, 501, 1031

Kuiper, L., & Hermsen, W. 2015, MNRAS, 449, 3827

Kuiper, L., Hermsen, W., & Dekker, A. 2018, MNRAS, 475, 1238

Link, B., Epstein, R. I., & Lattimer, J. M. 1999, PhRvL, 83, 3362

Liu, H.-Y., Zhou, S.-Q., Zhang, Y.-Q., Feng, Z.-W., & Zhou, X. 2021, RAA, 21, 154

Lower, M. E., Johnston, S., Dunn, L., et al. 2021, MNRAS, 508, 3251

Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45

Lyne, A. G., Shemar, S. L., & Smith, F. G. 2000, MNRAS, 315, 534

Maitra, C., Ballet, J., Filipović, M. D., et al. 2015, A&A, 584, A41

Maitra, C., Esposito, P., Tiengo, A., et al. 2021, MNRAS, 507, L1

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993

Marshall, F. E., Gotthelf, E. V., Middleditch, J., Wang, Q. D., & Zhang, W. 2004, ApJ, 603, 682

Marshall, F. E., Gotthelf, E. V., Zhang, W., Middleditch, J., & Wang, Q. D. 1998, ApJL, 499, L179

McKenna, J., & Lyne, A. G. 1990, Natur, 343, 349

Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791

Melatos, A., Howitt, G., & Fulgenzi, W. 2018, ApJ, 863, 196

Mereghetti, S., Rigoselli, M., Taverna, R., et al. 2021, ApJ, 922, 253

Middleditch, J., Marshall, F. E., Wang, Q. D., Gotthelf, E. V., & Zhang, W. 2006, ApJ, 652, 1531

Montoli, A., Antonelli, M., & Pizzochero, P. M. 2020, MNRAS, 492, 4837

Owen, R. A., Filipović, M. D., Ballet, J., et al. 2011, A&A, 530, A132

Pavan, L., Bordas, P., Pühlhofer, G., et al. 2014, A&A, 562, A122

Pavan, L., Pühlhofer, G., Bordas, P., et al. 2016, A&A, 591, A91

Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Natur, 567, 200

Ransom, S. M., Eikenberry, S. S., & Middleditch, J. 2002, AJ, 124, 1788

Ray, P. S., Guillot, S., Ho, W. C. G., et al. 2019, ApJ, 879, 130

Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17

Reynoso, E. M., Johnston, S., Green, A. J., & Koribalski, B. S. 2006, MNRAS, 369, 416

Riles, K. 2017, MPLA, 32, 1730035

Shapiro, S. L., & Teukolsky, S. A. 1983, Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects (New York: Wiley)

Sieniawska, M., & Bejger, M. 2019, Univ, 5, 217

Smith, D. A., Bruel, P., Cognard, I., et al. 2019, ApJ, 871, 78

Tomsick, J. A., Bodaghee, A., Rodriguez, J., et al. 2012, ApJL, 750, L39

Townsley, L. K., Broos, P. S., Feigelson, E. D., Garmire, G. P., & Getman, K. V. 2006, AJ, 131, 2164

Wang, Q. D., & Gotthelf, E. V. 1998, ApJ, 494, 623

Weisskopf, M. C., Brinkman, B., Canizares, C., et al. 2002, PASP, 114, 1

Yu, M., Manchester, R. N., Hobbs, G., et al. 2013, MNRAS, 429, 688

Zane, S., Haberl, F., Israel, G. L., et al. 2011, MNRAS, 410, 2428