The “Magnificent Seven” X-Ray Isolated Neutron Stars Revisited. I. Improved Timing Solutions and Pulse Profile Analysis

We present the first systematic X-ray pulse timing analysis of the six members of the so-called “Magnificent Seven” nearby thermally emitting isolated neutron stars (XINS) with detected pulsations. Using the extensive collection of archival XMM-Newton, Chandra, and NICER observations spanning over two decades, we obtain the first firm measurement of the spin-down rate for RX J2143.0+0654, while for the rest we improve upon previously published spin ephemerides and extend them by up to an additional decade. Five of the XINS follow steady spin-down with no indication of major anomalies in their long-term timing behavior; the notable exception is RX J0720.4−3125, for which, in addition to confirming the previously identified glitch, we detect a second spin derivative. The high-quality folded X-ray pulse profiles produced with the updated timing solutions exhibit diverse and complex morphologies, as well as striking energy dependence. These peculiarities cannot be readily explained by blackbody-like isotropic emission and simple hot-spot configurations, hinting at the presence of complex multitemperature surface heat distributions and highly anisotropic radiation patterns, such as may arise from a strongly magnetized atmospheric layer.

XINS exhibit a number of peculiarities in their observed properties that defy straightforward explanations.For instance, the confirmed optical and UV counterparts of all seven XINS (Kaplan et al. 2011) reveal that the measured photometry exceeds the flux predictions in the Rayleigh-Jeans tail extrapolated from the blackbody-like X-ray spectra.Furthermore, although XINS are typically much less luminous than younger thermally cooling NSs, their X-ray luminosities are too large to be explained by passive cooling alone; they also exceed their spin-down luminosities ( Ė ∝ Ṗ P −3 ), making them distinct from rotation-powered pulsars with comparable P and Ṗ values.This implies that the observed surface thermal radiation does not come at the expense of the rotational energy reservoir of the star and requires an additional source, most probably release of internal heat due to magnetic field decay, which may explain the unexpectedly strong surface magnetic fields of XINS.Regardless of the heating mechanism, the soft periodically modulated X-ray emission indicates the presence of a highly nonuniform heat distribution, which presumably traces the magnetic field topology at the surface and its vicinity.Thus, in principle, the observed surface X-ray emission harbors information about the magnetic field properties, including its surface distribution as well as its strength; the latter is also encoded in the beaming pattern of the surface atmospheric emission, which is strongly dependent on the magnetic field strength and field line orientation (see, e.g., Meszaros 1992;Pavlov et al. 1994;Ho et al. 2008;Perna et al. 2013;De Grandis et al. 2021).
At present, a self-consistent model of the multiwavelength spectral energy distribution for all XINSs is lacking and their surface composition and thermal and magnetic maps remain unconstrained, except for possibly RX J1856.5−3754(Ho et al. 2007;Ho 2007).To rectify this situation, we have commenced an in-depth analysis of the voluminous archival X-ray data set of the current sample of XINS.Mapping the surface X-ray emission as a function of time and energy for these sources by applying realistic neutron star models that consider all practically important physics has the potential to reveal the underlying physical processes responsible for their peculiar properties.A crucial prerequisite for such an undertaking is the availability of reliable long-term phase-coherent timing solutions which permit the folding of all existing event data with sufficient time resolution to obtain energy-resolved pulse profiles with the best possible photon statistics.
In this paper, the first in a series, we present a comprehensive and uniform analysis of the six2 XINS that exhibit pulsed emission to obtain greatly improved longterm timing solutions and to examine the properties of their pulsed emission more closely.To accomplish this, we make use of the wealth of X-ray Multi-Mirror Mission (XMM-Newton) and Chandra X-ray Observatory exposures collected over the past two decades, as well as more recent Neutron Star Interior Composition Explorer (NICER) observations.Here, we focus exclusively on timing and pulse profile analyses; discussions of spectral analyses, including those made by previous works, will be presented in the second paper in the series.The third paper will focus on realistic modeling of the surface thermal radiation.The present work is organized as follows.Section 2 summarizes the X-ray observations and reduction procedures used to extract the data products used in the analysis.In Section 3 we describe the timing procedure used in this work, and in Section 4 we lay out the resulting long-term timing analysis of the XINS, while Section 5 presents an examination of the properties of the pulsed emission as a function of energy and time.We offer a discussion in Section 6 and conclusions in Section 7.

OBSERVATIONS AND DATA ANALYSIS
The analysis presented herein makes use of archival Xray observations of the six bright XINS with confirmed X-ray pulsations.For the purposes of this work, we selected only observations with sufficient time resolution to permit a pulse profile to be constructed.Thus, in the case of RX J0420.0−5022we do not use XMM-Newton EPIC MOS1/2 full-frame (FF) mode observations due to the 2.6 s readout time, which is 75% of the neutron star spin period; likewise, we do not use Chandra ACIS-S data obtained in full imaging mode, which has a 3.2 s readout time.
The total filtered exposures for each XINS using each observatory/instrument combination used in this work are summarized in Table 1, while a detailed list of the individual observations including observation identifiers, UTC date and time, detector mode and filter, and the effective exposure time are presented in Tables 4-9 of the Appendix.The data reduction procedures employed for each telescope/instrument to arrive at the clean event lists used for timing are described below.

XMM-Newton EPIC
In this work, we rely primarily on XMM-Newton European Photon Imaging Camera (EPIC) pn (Strüder et al. 2001) and MOS1/2 (Turner et al. 2001) data sets spanning over two decades.We do not consider any Reflection Grating Spectrometer (RGS) event data due to the substantially lower sensitivity compared to the imaging data and the greatly reduced signal-to-noise ratio of the dispersed source emission.
The X-ray events were reduced and extracted using the Science Analysis Software (SAS) package version xmmsas 20211130 0941-20.0.0 starting from the unprocessed (ODF) data products.After running the epproc and emproc pipelines for pn and MOS, respectively, the data were inspected for instances of strong background flaring which were excised by considering the count rate over the entire image in the 10-12 keV range for EPIC pn and >10 keV for EPIC MOS1/2 from light curves binned at 100 s intervals.Segments of the light curve with count rates exceeding 0.4 and 0.35 counts s −1 for pn and MOS1/2, respectively, were removed to produce a set of clean good time intervals.For the imaging observations, the source events were extracted from circular regions of radius 32 ′′ , which encircle ≈90% and ≈85% of the total point source energy at a See Tables 4-9 in the Appendix for a detailed list of all observations for each target.
b Only Chandra ACIS-S observations conducted in continuous clocking (CC3 GRADED) mode are considered.
c For RX J0420.0−5022,no usable Chandra or NICER data are available.
∼1.5 keV for the pn and MOS1/2, respectively, centered on the point source.For RX J1856.5−3754,there is a single 44 ks EPIC pn observation (ObsID 0201590101) obtained in Fast Timing mode, which foregoes one imaging dimension to permit fast readout (29 µs).The source events for this exposure were extracted from a rectangular region 13 pixels of width in the RAWX direction centered on row 37.The source event data were filtered using the recommended flag and pattern parameter values (PATTERN≤4 and #XMMEA EP for EPIC-pn and PATTERN≤12 and #XMMEA EM for EPIC MOS1/2) specified through the evselect tool in SAS.
We note that a portion of the XMM-Newton data sets for RX J0720.4−3125 , RX J1308.6+2127, and RX J1856.5−3754suffer from nonnegligible event pileup, as assessed by the epatplot task in SAS.However, although important for any detailed spectroscopic analyses, this effect does not significantly impact the pulse timing and cursory pulse profile analyses presented below.
The timestamps of the filtered source events were translated from the Terrestrial Time (TT) standard to Barycentric Dynamical Time (TDB) using the barycen tool in SAS with the DE405 JPL solar system ephemeris.

Chandra
The Chandra X-ray Observatory data set of the XINS used in this work consists of Advanced CCD Imaging Spectrometer (ACIS) exposures conducted in continuous clocking (CC) mode and High Resolution Camera (HRC) exposures originally intended for high-resolution dispersed spectroscopy with the Low Energy Transmission Grating (LETG) placed in the optical path.The data analysis was carried out using CIAO version 4.15 (Fruscione et al. 2006) and the accompanying calibration database CALDB 4.10.4.The barycentering of all Chandra ACIS and HRC source events was accomplished us-ing the axbary task in CIAO assuming the DE405 JPL solar system ephemeris.

ACIS-S Continuous Clocking Observations
A number of archival Chandra observations of XINS are available for which the back-illuminated ACIS-S3 chip was placed at the aim point.
To avoid photon pileup that would be caused by high source count rates, the detector was configured in Continuous Clocking (CC3 GRADED) mode, which affords a rapid (2.85 ms) readout time by sacrificing one dimension of spatial imaging.Following recommended procedures, the source events were extracted from a rectangular region of width 5 ′′ along the imaging direction.Due to the bright nature of the sources under consideration and the small point spread function of Chandra, the background contribution to the source count rate is virtually negligible.

HRC-S LETG Observations
For the Chandra HRC-S observations with the LETG, we extracted the zeroth-order (i.e., undispersed) pointsource emission using circular regions of radius 2 ′′ centered on the source position.In most instances, the first order dispersed spectrum is quite faint while also introducing higher background, and so is not used in the timing analysis.While the zeroth order provides no useful spectral information due to the poor intrinsic energy resolution of the HRC, because of the exceptionally soft spectra of all XINS we expect the vast majority of registered source events to have energies ≤1 keV.
The HRC event time tagging problem caused by a wiring error in the detector electronics 3 identified early in the Chandra mission results in a degraded time res-olution from the intrinsic 16 µs to ≈4 ms (determined mainly by the mean time between events over the entire detector).As this effective time resolution is still orders of magnitude smaller than the spin periods of the XINS, the HRC data are still suitable for the timing analyses that are the subject of this work.

NICER XTI
With the exception of RX J0420.0−5022,extensive NICER X-ray Timing Instrument (XTI; Gendreau et al. 2016) observations have been conducted of the sample of XINS starting as early as 2017 June.To obtain filtered event lists suitable for timing analyses, we used the NICERDAS software that is part of HEASoft version 6.31, following the procedures recommended by the NICER team.Specifically, the cutoff rigidity COR SAX parameter (as originally defined for the BeppoSAX mission), which serves as an indicator of ambient in-orbit particle background levels (see, e.g., Bogdanov et al. 2019), was restricted to values >1.5.In addition, the undershoot rate, which is a measure of the severity of contamination due to optical photons at the lowest energies, was limited to a median value of ≤200 during each good time interval.The rate of event overshoots, which are background contamination events caused by highenergy charged particles that deposit a large amount of charge on the detector, was restricted to values ≤1.5.The Potsdam planetary geomagnetic activity index4 k p , which is an estimate of fluctuation disturbances to the local magnetic fields, was limited to values ≤5.The filtered event arrival times were barycentered using the barycorr tool in HEASoft assuming the DE405 JPL solar system ephemeris.

PULSE TIMING PROCEDURE
To enable the timing and pulse profile analyses presented here, we first barycentered the XMM-Newton, Chandra, and NICER data sets assuming the best known position for each source and the JPL DE405 solar system ephemeris using the appropriate software tool for each observatory as described in Section 2. For RX J0720−3125, RX J1308.6+2127, and RX J1856.5−3754, which have measured proper motions, we accounted for the change in position over time by applying the barycentering for each observation using the proper-motion-corrected position at that epoch.In practice, due to the slow spins of the XINS, this position adjustment has a small effect on the corrected event timestamps and thus on the outcomes of the timing analyses.Moreover, setting the position and proper motion parameters free in the fit does not lead to an improvement in the timing solutions nor does it result in improved constraints on the astrometric parameters.
Initially, we attempted to fold the extracted source event data at the most recent published timing solutions for each XINS.However, since some of these ephemerides are over a decade old, they typically did not accurately predict the spin of the neutron star for the more recent data sets, as evidenced by phase drift and smearing of the accumulated profile.In light of this, we proceeded to carry out a formal timing analysis for the six XINS to obtain more up-to-date and precise timing solutions.We conducted this analysis with an unbinned event-based Bayesian likelihood evaluation approach as implemented in the event optimize code in the PINT pulsar timing software package (Luo et al. 2021), which has been thoroughly tested and is used extensively for timing of Fermi Large Area Telescope γ-ray pulsars (see, e.g., FERMI-LAT Collaboration et al. 2022;Smith et al. 2023).In this technique, the photon event data are not grouped into time of arrival (TOA) measurements but are instead evaluated collectively as an unbinned ensemble against a timing model (with likelihood evaluation based on the H-test; de Jager et al. 1989; de Jager & Büsching 2010) using a smooth pulse profile template constructed from one or more Gaussians (as necessary to reproduce the main features of the pulses).An initial profile template was constructed from the observation with the strongest pulse detection, which was most often the longest XMM-Newton exposure available.For parameter set sampling, the code uses the popular emcee affine-invariant Markov chain Monte Carlo (MCMC) ensemble sampler (Foreman-Mackey et al. 2013).For every inference run, we considered 200 sampler walkers and at least 1000 steps per chain to ensure convergence, with the first 100 steps used for burn-in.The photon based likelihood sampling technique is particularly beneficial for NICER XTI event data since it allows the inclusion of all short exposures (≲1 ks) that would otherwise be discarded in a conventional timing analysis as they do not provide sufficient photon statistics to produce a useful TOA measurement.In addition, for the fainter XINS such as RX J0420.0−5022 and RX J2143.0+0654, this approach is more advantageous due to the limited photon statistics that are obtained from the typical XMM and Chandra exposures available for these targets.The same holds true for RX J1856.5−3754,which, although the brightest in terms of total source count rate, is difficult to time due to the exceptionally low pulsed fraction.One drawback of this event-based likelihood procedure is that it is substantially more computationally intensive than a conventional TOA analysis in terms of processor and random access memory requirements, with computational time and memory usage scaling proportionally with the number of photons considered in the analysis.
As free parameters in the timing model we consider the spin frequency (ν), the first frequency derivative ( ν), and the arbitrary pulse phase offset relative to the ref- −61.9(4) b Astrometric parameters measured from optical or X-ray imaging observations and used for barycentering of the event data are from Haberl et al. (2004), Kaplan et al. (2002a), Walter et al. (2010), andSchwope et al. (2009).
c Reference epoch for measured and derived spin parameters quoted in the table.
d Single-trial H-test pulse detection significance of the barycentered and stacked event data sets folded at the highest-likelihood timing solution.
g Surface equatorial magnetic field strength B surf = 3.2 × 10 19 (− ν/ν 3 ) 1/2 G under the assumption of vacuum dipole magnetic braking for a neutron star with R = 10 km and M = 1.4 M ⊙ .
erence Modified Julian Date (MJD).For each target we also explored the possibility of a change in the spin-down rate by introducing a second frequency derivative (ν).In most cases, the 3σ credible intervals obtained for ν overlap with zero, suggesting that for the currently available data the addition of this parameter is not warranted.This is also indicated by the fact that including ν does not result in an improvement in the pulsed significance as determined by the H-test.As described later, for RX J0720.4−3125 the timing model required the introduction of additional glitch parameters as well as a second frequency derivative (ν).For all six XINS, the latest previously published ephemerides were used as the initial guess of the free parameters for the sampler.

TIMING RESULTS
The results of the timing analysis for five of the XINS are compiled in Table 2, and for RX J0702.4−3125 in Ta-ble 3. We quote measured parameters in frequency space (ν in Hz and ν in Hz s −1 ) but for convenience we also list the corresponding conversions to period space (with P in s and Ṗ in s s −1 ).For the event-based likelihood analysis results, we quote the customary 50 th percentile and the 16 th and 84 th percentiles credible intervals as the lower and upper uncertainties, respectively, of the posterior distribution.Figures 7 and 8 in the Appendix show the posterior distributions of the timing parameters in the form of corner plots.The reference epochs (in MJD) used for the quoted spin measurements correspond to the previously published values as listed in the ATNF pulsar catalog (Manchester et al. 2005), except for RX J0702.4−3125, for which the reference epoch of MJD 56,598 corresponds to the midpoint of the postglitch interval (see Section 4.2 for details).
The final timing solutions in the form of Tempo-style parameter (.par) files are provided as supplemental material to this paper; they contain the values of the fitted parameters corresponding to the highest likelihood sample found in the inference analysis for each XINS.In what follows, for each target we provide initial discovery and background information, a summary of previous work specific to pulse timing, and present the updated timing results.4.1.RX J0420.0−5022RX J0420.0−5022(1RXS J042003.1−502300) was identified as a neutron star candidate by Haberl et al. (1999).Follow-up ROSAT observations showed a remarkably soft thermal spectrum and yielded an improved position, both of which strengthened its classification as a neutron star.While initially a 22.7 s period was suggested, more sensitive observations with XMM-Newton found a period of 3.45 s instead, the shortest among the XINS (Haberl et al. 2004).The same observations also showed that RX J0420.0−5022 was the coolest INS, with blackbody kT ≃ 45 eV; it is also the faintest in terms of observed photon flux in the soft Xray band.The first coherent timing solution for this XINS was derived by Kaplan & van Kerkwijk (2011), who obtained a spin-down rate of ν = (−2.3± 0.2) × 10 −15 Hz s −1 based on the 2010-2011 XMM-Newton observing campaign and ν = (−2.314± 0.008) × 10 −15 Hz s −1 when observations dating back to 2002 were included.However, for the latter they noted a number of possible alternative timing solutions due to cycle count ambiguities, of which the quoted one was judged the best based on the χ 2 value of the TOA fit.The spin parameters imply B s = 1 × 10 13 G, τ c = 1.9 Myr, and by far the highest spin down luminosity among the XINS of Ė = 3 × 10 31 erg s −1 .
For the present timing analysis, we consider XMM-Newton data spanning the date range 2002 December 30 to 2019 May 22.The only available archival Chandra ACIS-S observations were obtained in full frame imaging mode, which has a readout time of 3.2 s, rendering them unusable for a timing analysis of this 3.45 s neutron star.No NICER XTI observations of RX J0420.0−5022 have been conducted to date, presumably due to its relatively faint nature.For the barycentering and timing analysis, we assume the source position obtained by Haberl et al. (2004) based on Chandra imaging observations.Despite the eight year gap in data, from our event based likelihood analysis we obtain a fairly precise spindown rate measurement of ν = (−2.4393± 0.0005) × 10 −15 Hz s −1 .This value is consistent with the 2010-2011 ν value from Kaplan & van Kerkwijk (2011) but not with their longer baseline (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) ν, suggesting that their chosen solution was not the correct one.This is expected since RX J0420.0−5022 has the shortest spin period among the XINS by a wide margin and is thus most susceptible to pulse cycle count ambiguities.3 for the timing model used and text for further details.
We note that a number of alternative, equally spaced solutions are also apparent in the ν − ν plane from our analysis, which likely arise due to the large gaps in data coverage from 2003-2010 and 2011-2019.Nevertheless, the quoted ν is strongly favored based on the strongest pulse detection significance determined from the H-test (12.09σ single trial).
4.2.RX J0720.4−3125RX J0720.4−3125(1RXS J072025.1−312554) is a XINS with a spin period of P = 8.39 s, discovered by Haberl et al. (1997) over the course of the ROSAT All-Sky Survey.After multiple attempts to constrain the spin-down rate (Kaplan et al. 2002b;Zane et al. 2002;Cropper et al. 2004), a fully phase-connected long-term timing solution was finally obtained by Kaplan & van Kerkwijk (2005a) using Chandra and XMM-Newton observations.A spin-down rate of Ṗ = 6.98 × 10 −14 s s −1 was measured, which implies a characteristic age of τ c = 2 Myr and a dipole magnetic field strength of B s = 2.4 × 10 13 G.Hohle et al. (2010, see also Haberl et al. 2006;van Kerkwijk et al. 2007;Hohle et al. 2009) used additional observations to extend the timing solution but identified residuals suggestive of long term variations in the spin behavior that could not be accounted for solely by a Ṗ .The interpretations offered for this behavior included  Ṗ (10 −14 s s −1 ) 7.151 +0.003 −0.003 1.9 Bs (10 13 G) 2.5 a The numbers quoted for the fitted parameters correspond to the 50 th percentile of the posterior distribution with the 16 th (−) and 84 th (+) percentiles given in parentheses.
b Astrometric parameters measured from optical imaging observations (Kaplan et al. 2007).
free precession of the neutron star or a glitch.However, a follow-up analysis (Hohle et al. 2012) concluded that the apparent timing irregularities are best explained by a glitch that occurred at some point during the gap in observations spanning the MJD interval 52, 866 ± 73.
van Kerkwijk et al. (2007) reported the first frequency derivative to be ν = (−1.04± 0.03) × 10 −15 Hz s −1 and the glitch permanent frequency increment to be ∆ν = (4.1 ± 1.2) × 10 −9 Hz, while the glitch permanent frequency first derivative increment was not well constrained with ∆ ν = (−4 ± 3) × 10 −17 Hz s −1 .RX J0720.4−3125 has been observed extensively over the past two decades (though with a sizable gap in coverage between 2012 and 2017), with a total on-source exposure second only to RX J1856.5−3754(see Table 1).XMM-Newton has targeted it in 21 separate observa-tions for which suitable on-source data is available, spanning from 2000 May 13 to 2019 October 16, for total effective net exposures of 476.4,547.4,and 524.1 ks for pn, MOS1, and MOS2, respectively.Archival Chandra observations consisting of 21 HRC-S LETG and 11 ACIS-S Continuous Clocking mode exposures that cover the period between 2000 February 1 and 2012 November 28 are also used in this analysis.Finally, there are 58 mostly brief NICER XTI exposures that add up to 104.6 ks of clean on-source time accumulated from 2017 June 30 to 2022 May 21.Collectively, the usable event data from the three observatories cover a period of 22 years, providing an additional decade of timing baseline coverage compared to the last published ephemeris.For the astrometric parameters, we adopt the coordinates and proper motion values measured by Kaplan et al. (2007) using optical observations.
As in prior published analyses of RX J0720.4−3125, to account for the anomalous long-term spin behavior we introduce three glitch parameters: the glitch spin phase increment ∆Φ, the permanent frequency increment ∆ν, and the permanent frequency first derivative increment ∆ ν.For illustrative purposes, Figure 1 shows the results of a TOA based timing analysis of RX J0720.4−3125 that considers a steady spin down model based on the most recent XMM-Newton and NICER data set that is propagated backwards in time.A clear discontinuity in the model is apparent around MJD 52886 (which we fix as the presumed epoch of the glitch), consistent with the glitch interpretation (see top panel of Fig. 1).Moreover, we find that even after accounting for the glitch, there are broad residuals that suggest the presence of a second frequency derivative ν, which we also include as a free parameter.
The updated timing solution for RX J0720.4−3125 using the event-based likelihood technique is presented in Table 3, while the posterior distributions of the free parameters are shown in Figure 8.The resulting values of ν = (−1.0156± 0.0004) × 10 −15 Hz s −1 , ∆ν = (5.02± 0.05)×10 −9 Hz, and ∆ ν = (−1.20±0.04)×10−17 Hz s −1 are in agreement with the values reported by van Kerkwijk et al. (2007) but with substantially improved precision, owing to the greatly extended time span of the current archival data set.We also obtain a fairly well constrained value of the second frequency derivative ν = 5.9 +0.9 −0.8 × 10 −27 Hz s −2 , making RX J0720.4−3125 the only XINS thus far with a significant ν.With this improved timing model we obtain a stacked pulse profile with a 204.73σ single-trial significance, by far the highest among the XINS.
RX J0806.4−4123 has the least amount of devoted on-source exposure among the XINS.Here, we mainly make use of XMM-Newton observations spanning the time period from 2000 November 8 to 2019 May 16 that yield net exposures of 98.6, 73.2, and 72.1 ks for pn, MOS1, and MOS2, respectively.Two Chandra observations of RX J0806.4−4123 have been conducted so far, one 18.1 ks ACIS-S Continuous Clocking mode exposure from 2010 March 21 and one 29.8ks HRC-I LETG exposure from 2018 September 1, which we include in the timing analysis.Finally, we consider NICER XTI data collected between 2019 March 11 and 2023 February 7 for a total clean on-source exposure of 104.6 ks.All observations are compiled in Table 6.For the barycentering of event data and timing analysis, the source sky coordinates measured by Haberl et al. (2004) from Chandra ACIS-S observations were considered.
The event-based likelihood timing analysis results in a spin frequency derivative ν = (−8.18± 0.01) × 10 −17 Hz s −1 , in full agreement with that found independently and using alternative methods by Posselt et al. (2024) but with substantially reduced uncertainties due to the wider time span of the data set .The inferred ν value is the smallest among the XINS, which implies a low spin-down luminosity of only Ė = 3 × 10 29 erg s −1 and a characteristic age of τ c = 17 Myr, an order of magnitude larger than the rest of the XINS.The combined XMM-Newton, Chandra, and NICER data yield a pulse detection significance of 29.11σ based on the Htest in the optimal 0.27-0.89keV energy range.

RX J1308.6+2127
RX J1308.6+2127(also cataloged as RBS 1223 and RXS J130848.6+212708) was identified as a strong neutron star candidate among the set of ROSAT bright sources by Schwope et al. (1999).A periodicity of 5.16 s was initially reported by Hambaryan et al. (2002), which later investigations identified to actually correspond to half of the true neutron star spin period of P = 10.31 s (Haberl 2004;Schwope et al. 2005).RX J1308.6+2127 is notable in that it is the only XINS to exhibit two clearly separated pulses per rotation, which suggest the presence of at least two distinct X-ray emitting regions on the stellar surface.The first and only published coherent timing solution was obtained by Kaplan & van Kerkwijk (2005b), who arrived at a frequency rate of change value of ν = (1.053± 0.003) × 10 −15 Hz s −1 .Hambaryan et al. (2011) have used XMM-Newton observations of RX J1308.6+2127 to place constraints on the emission geometry, surface chemical composition, magnetic field strength, and neutron star compactness.More recently, Borghese et al. (2017) claimed the detection of a narrow phase-dependent feature in the X-ray spectrum with an energy of ∼740 eV and an equivalent width of ∼15 eV.The purported feature was detected only in ∼1/5 of the phase cycle, and appears to be present for the entire time span covered by the observations examined (2001 December to 2007 June).The strong dependence on the pulsar rotation and the narrow width was speculated to be due to resonant cyclotron absorption/scattering in a confined high magnetic field structure close to the stellar surface.
For our updated timing analysis of RX J1308.6+2127we include a more recent XMM-Newton observation from 2019 December 16 (ObsID 0844140101), a single 87.2 HRC-S LETG observation (ObsID 4595) from 2005 April 12, and eight ACIS-S CC mode observations obtained between 2006 February 12 and 2006 August 3.To this we add a collection of NICER XTI exposures covering the period 2017 November 14 to 2018 June 26 with a combined clean effective exposure of 120.9 ks.In the timing analysis, the right ascension and declination for the source were assumed to correspond to the Chandra position reported in Kaplan et al. (2002a).
The median value we obtain for the rate of change of frequency for RX J1308.6+2127 is ν = (−1.05427± 0.00003) × 10 −15 Hz s −1 , in excellent agreement with the Kaplan & van Kerkwijk (2005b) timing solution but with a ≈100-fold improvement in precision, despite the decade-long gap in X-ray observations of this source (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017).By folding the entire collection of event data at the newly derived ephemeris we obtain a maximum pulse detection significance of 139.56σ as determined by the H-test in the 0.30-1.21keV band.4.5.RX J1856.5−3754RX J1856.5−3754(1RXS J185635.1−375433) is regarded as the canonical XINS and is the X-ray brightest.It was the first of the XINS to be identified as a strong neutron star candidate based on its soft spectrum and the lack of a bright optical counterpart (Walter et al. 1996).It is one of the nearest known NSs of any variety, with a parallax distance of only D = 123 +11 −15 pc (Walter et al. 2010).Despite these facts, it was the sixth XINS to have pulsations detected (Tiengo & Mereghetti 2007); although exceptionally bright in terms of photon flux, the thermal radiation features an extraordinarily low pulsed fraction of ∼1%, necessitating deep XMM-Newton observations to uncover the 7.055 s periodicity.This property can be plausibly interpreted as a thermally emitting region that is nearly aligned with the spin axis and/or a close alignment between the spin axis and the observer line of sight (Ho 2007).
Previously, van Kerkwijk & Kaplan (2008) conducted a timing analysis, which produced only an upper limit on the spin-down rate.A definitive measurement was obtained in a follow-up study by Kaplan & van Kerkwijk (2009b).Most recently, De Grandis et al. ( 2022) presented an updated timing solution for RX J1856.5−3754, with ν = −6.042(4)× 10 −16 Hz s −1 but using XMM-Newton EPIC-pn exposures only up to and including 2022 April 3 (ObsID 0810842001) and a portion of presently available NICER observations through 2019 June 18 (ObsID 2614010139).
There exist a large number of archival X-ray observations of RX J1856.5−3754 because it is observed with great regularity as an instrument calibration target.This has resulted in a bountiful archival data set with ≈2 Ms of XMM-Newton, ≈1.6 Ms of Chandra, and 186 ks of NICER XTI exposures accumulated over the course of the past two decades.Here, we make use of EPIC pn as well as MOS1/2 exposures and include two additional XMM-Newton observations (Ob-sIDs 0810842201 and 0810842301) from 2022 September 24 and 2023 April 1 (see Table 8 in Appendix) that were not used in De Grandis et al. (2022).We also analyze an extended set of NICER XTI observations up to 2022 October 10 that also include 33 calibration exposures acquired at off-axis angles5 that are usable for a timing analysis.Finally, we consider a number of Chandra ACIS-S Continous Clocking mode exposures as well as the zeroth order point source counts from the extensive set of Chandra HRC-S LETG observations, which were also not considered in De Grandis et al. (2022).For the astrometric parameters of RX J1856.5−3754,we use the measurements obtained by Walter et al. (2010).
The event-based likelihood timing analysis for RX J1856.5−3754results in ν = (−6.0373+0.0008 −0.0007 ) × 10 −16 Hz s −1 , which is fully consistent with the measurement from De Grandis et al. ( 2022) within uncertainties, although with significantly improved precision owing to the use of a larger data set that is also extended by an additional year.Despite the order of magnitude more exposure than the other XINS, resulting in ≈ 2 × 10 7 counts collected for this source, the pulse detection significance of 35.13σ from the H-test (in the optimal energy band 0.15-0.74keV) is substantially lower than those for RX J0720.4−3125 and RX J1308.6+2127 because of the exceptionally low pulse amplitude.4.6.RX J2143.0+0654RX J2143.0+0654(also known as RBS 1774 and 1RXS J214303.7+065419) was identified as a likely neutron star candidate by Zampieri et al. (2001).Its true nature was established by Zane et al. (2005), who uncovered a spin period of 9.437 s and a thermal, blackbody-like spectrum.Follow up X-ray timing observations by Kaplan & van Kerkwijk (2009b) yielded only a marginally significant spin-down rate of ν = (−4.6 ± 2.0) × 10 −16 Hz s −1 .
For the updated timing analysis we consider XMM-Newton observations covering the period 2004 May 31 to 2019 May 22, though with an 11 year gap in coverage between 2008 and 2019.The total filtered onsource exposures are 90.2,136.1 and 136.9 ks for pn, MOS1, and MOS2, respectively.The Chandra data set for RX J2143.0+0654consists of a single 12.1 ks ACIS-S CC mode observation from 2009 December 24 and four HRC-S LETG observations from 2018 August 18 to 2018 September 14 that total 139.7 ks.The 29 NICER XTI observations targeting this source were carried out between 2022 July 26 and 2023 August 28, yielding a substantial boost in event data with 272.3 ks of clean exposure.For the barycentering and timing analysis, we take the position published in Schwope et al. (2009) of the probable optical counterpart of RX J2143.0+0654.
With the greatly extended timing baseline obtained by including additional Chandra, XMM-Newton, and NICER observations spanning 2018 August to 2023 August, we obtain the first solid spin-down measurement ν = (−4.663+0.006 −0.004 ) × 10 −16 Hz s −1 or Ṗ = (4.145+0.005 −0.004 ) × 10 −14 s s −1 for RX J2143.0+0654,consistent with the initial constraint from Kaplan & van Kerkwijk (2009b).The posterior distribution is broad, asymmetric, and multipeaked (see Fig. 7 in the Appendix), which is presumably caused by the large gap in X-ray observations from 2009 to 2018.The current spindown rate measurement implies a surface dipole field of B s = 2 × 10 13 G, characteristic age of τ c = 3.6 Myr, and a spin-down luminosity of Ė = 2 × 10 30 erg s −1 , in line with the rest of the XINS sample.Folding the combined event list from all observatories results in an H-test pulse detection significance of 27.01σ in the 0.15-1.18keV band.

PULSE PROFILE ANALYSIS
With the updated timing solutions in hand, we can now construct high signal-to-noise pulse profiles and examine their morphology and energy dependence in greater detail than previously possible.Figure 2 shows a compilation of the normalized (counts per bin divided by the mean across all bins) broadband pulse profiles obtained by folding and stacking the event data from all telescopes considered in the timing analysis.We note that this stacking results in an average of the instrument response curves as a function of energy of the different telescopes, weighted based on the relative effective areas and exposure times.The energy range used for each XINS is the interval that was found to produce the maximum pulse detection significance: 0.16-1.15,0.15-1.21,0.27-0.89,0.30-1.21,0.15-0.74,and 0.15-1.18keV, for the six XINS in alphanumerical order, respectively.2 and 3. See text for the energy range used for each profile.A rotational cycle is repeated for clarity.
For some of the XINS, the high fidelity of the updated profiles reveal interesting morphology.RX J0720.4−3125,RX J0806.4−4123, and RX J1856.5−3754show a single, albeit clearly asymmetric, peak per rotation period.The pulse from RX J1856.5−3754 is characterized by a shallow leading edge and a steeper trailing edge.RX J0720.4−3125 has a clearly asymmetric "sawtooth" pulse, with a rapid rise and a more gradual decline, while RX J0806.4−4123exhibits an even more exaggerated triangular wave morphology, with a slower rise and a very fast decline of the pulse.RX J1308.6+2127 has two clearly separated peaks per rotation period, with comparable amplitude but with one noticably broader than the other.The current data set folded at the improved timing solution for J2143.0+0654shows the most complex broad band profile among the XINS, with a main broad peak flanked by two weaker and narrower pulses.In the case of RX J0420.0−5022,due to the paucity of pulsed photons the detailed structure of the pulse is not clearly evident.It appears to show a single pulse per rotation, but limited information can be extracted from the present data set.

Harmonic Decomposition
The complexity of the pulse shapes can be judged more quantitatively using the harmonic content of the periodic signal in the Fourier domain.In particular, the presence of significant harmonics in the power spectrum beyond the fundamental frequency provides an indication of how strongly the pulse shape deviates from a pure sinusoid.The total profiles for the six XINS with their Fourier components superposed are shown in Figure 3. Due to the limited data set, the profile RX J0420.0−5002 is adequately described by a single broad sine wave.RX J0720.4−3125 and RX J1856.5−3754 have a significant second harmonic (i.e., first overtone), which arise due to the pronounced skewness of the single pulse.RX J0806.4−4123requires the first four Fourier components to fully describe the highly asymmetric pulse shape.In RX J1308.6+2127, the second harmonic is predictably more prominent than the fundamental, as expected from a signal with two pulses per cycle, but additional higher harmonics are also present due to the differing widths of the two pulses.Despite the limited photon statistics for RX J2143.0+0654, its power spectrum is the richest in harmonic content due to the multiple features of the pulse profile.

Energy Dependence
In Figure 4 we show the profiles for the six XINS but split into different energy ranges.For RX J0720.4−3125 and RX J1308.6+2127, the data are divided into three, while for the rest into two adjoining bands.The en-ergy ranges were selected through a visual inspection of the folded data in a number of trial ranges; the bands selected were chosen because they exhibited the largest discrepancies between the profiles in the different energy ranges.The Chandra HRC events are excluded from these energy-resolved profiles as they provide no reliable spectral information.The wealth of pulsed source counts collected for RX J0720.4−3125,RX J1308.6+2127, and RX J1856.5−3754allows us to examine the energy dependence of the pulsed emission in even greater detail.For this purpose, in Figure 5   .Normalized pulse profiles of the six XINS split into multiple energy bands.For RX J0720.4−3125 and RX J1308.6+2127,profiles for three energy intervals are shown in orange, blue and purple for the soft, mid, and hard bands, respectively.For the rest, the binned profile data are split into two energy ranges with the softer and harder bands corresponding to the orange and blue curves, respectively.The superposed solid red, blue, and purple lines in the second rotational cycle show the best fit with a representation of the profile constructed from the empirical Fourier coefficients of the set of photon phases for the same energy range.
a function of both pulse phase and photon energy for these three XINS.To obtain the pulsed counts in every energy band considered we subtracted the DC (i.e., unpulsed) component, with the level of this DC component taken from the pulse phase bin with the minimum number of counts in each energy interval considered.
One of the outliers again is RX J0720.4−3125, which appears to undergo a dramatic transformation of the profile around 0.6 keV, with a narrower and offset pulse at higher energies (see Figure 5).This behavior is immediately apparent in Figure 4, where the X-ray profile is shown in three energy bands, 0.15-0.3,0.3-0.7,and 0.7-1.2keV.In Borghese et al. (2015), this peculiarity was interpreted as being due to an absorption feature based on a more limited data set (with only a single XMM-Newton observation used to draw that conclusion).An alternative interpretation for this feature is the presence of two distinct and widely separated X-ray emitting regions on the stellar surface with grossly different areas and temperatures, analogous to those of the Puppis A central compact object RX J0822−4300 (see Gotthelf & Halpern 2009) and the so-called "Three Musketeers" middle-aged rotation-powered pulsars (e.g., De Luca et al. 2005).
Similar behavior is observed in RX J1308.6+2127,where the relative strengths of the two distinct pulses, which have comparable amplitudes across the entire soft band, change substantially above ≈0.8keV.The broader pulse, in particular, also shifts ahead by ≈0.05 rotational cycles at higher energies, while the narrower pulse, which is slightly stronger in the softest band (0.15-0.5 keV), greatly diminishes in amplitude above ≈0.8keV.
The other XINS show more subtle but no less peculiar changes in their pulsed emission as a function of energy.For RX J0806.4−4123,J1856.5−3754, and perhaps RX J0420.0−5022,there is compelling (albeit presently marginal) evidence for a phase lag of the emission at higher energies as well as a possible narrowing of the pulse.In RX J2143.0+0654,there are hints that the relative strengths of the three pulses may change with energy and undergo phase shifts as a function of energy.However, in all these instances, the photon statistics of the pulsed emission are not sufficiently high for definitive conclusions to be drawn.

Hard Pulsed Emission
There have been recent reports of evidence for faint non-thermal components at higher energies in some XINS (Dessert et al. 2020;De Grandis et al. 2022).Hard X-ray radiation from neutron stars can be produced by a variety of mechanisms.Extended or compact, unresolved non-thermal radiation can be produced by a wind nebula or a bow shock driven by a particle outflow from the neutron star as it encounters the ambient medium (Gaensler 2005;Kargaltsev et al. 2015).A hard tail can also be caused by inverse Compton scattering of the surface thermal radiation by a population of intervening energetic particles; for low optical depth this results in a power-law tail (Nishimura et al. 1986).Finally, in rotation-powered neutron stars (especially young, energetic ones), pulsed emission with a power-law spectrum can arise from particle acceleration in the pulsar magnetosphere.The possible presence of a hard tail in the spectra of XINS has led to speculations that this high-energy excess could be due to exotic processes such as emission of axions from the core of the neutron star (Buschmann et al. 2021).
To test for one of these possibilities, namely pulsed non-thermal emission, we consider only events in the range 1.5-5 keV, where the contribution of the soft thermal radiation is expected to be negligible and the nonsource background to not be dominant.We only use XMM-Newton data since the NICER XTI background at higher energies is fairly high due to the lack of imaging capabilities and the Chandra exposures are substantially less sensitive.We find that in the 1.5-5 keV band none of the six XINS exhibit statistically significant pulsations, with H-test significances of 1.15σ, 2.47σ, 2.25σ, 2.49σ, 0.51σ, and 1.54σ, for RX J0420.0−5022,RX J0720.4−3125,RX J0806.4−4123,RX J1308.6+2127,RX J1856−3752, and RX J2143.0+0654,respectively.These non-detections translate into 3σ upper limits on  the pulsed fractions of the hard emission above 1.5 keV of 6.1%, 2.1%, 4.5%, 6.0%, 1.5%, and 5.9%, respectively.We defer an in-depth spectroscopic investigation of the apparent excess at high energies for all XINS to a subsequent paper.

DISCUSSION
The pulsed thermal radiation observed from the XINS presumably arises due to changing aspect of hot surface emitting regions as the star rotates.The peculiar features of the observed XINS profiles, their Fourier domain properties, and their variation with photon energy are not easily reproducible by one or more hot spots on the surface that radiate a simple blackbody spectrum with equal intensity in all directions.Therefore, any deviations of the profile from a sinusoidal shape and large changes as a function of energy would arise due to a combination of: i) complicated hot spot arrangements; ii) non-isotropic surface radiation, such as that produced by an atmospheric layer; and/or iii) phase and energy dependent attenuation, by processes such as cyclotron absorption (e.g., Suleimanov et al. 2012).We note that sharpened and/or distorted pulse shapes could also be produced by gravitational light bending for ultracompact neutron stars (R NS ≲ 8 km), although such small stellar radii are not favored by current constraints (Abbott et al. 2018;Miller et al. 2021;Riley et al. 2021).Moreover, for the comparatively slow spin periods under consideration here (3 − 11 s), asymmetries in the profile introduced by Relativistic Doppler boosting and aberration of the surface emission induced by the rotation of the star are insignificant (see, e.g., AlGendy & Morsink 2014).
The XINS dipole magnetic field strengths deduced from the spin-down rate are ∼10 13 G, which falls in a regime where the emergent radiation from a surface atmospheric layer is expected to be strongly affected by the field (Potekhin 2014).This manifests in a highly anisotropic radiation pattern, with a "pencil" beam in the direction of the magnetic field and a fan beam that peaks at intermediate angles (Meszaros 1992;Pavlov et al. 1994).The detailed beaming pattern has a strong energy dependence and is determined in part by the magnetic field strength and orientation, as well as on the chemical composition of the topmost emitting layer of the atmosphere (Mori & Ho 2007;Ho et al. 2008).The observed properties of the pulsed emission from the six XINS can, in principle, be explained by emission from such a magnetized surface layer (see, e.g., Ho et al. 2007;Ho 2007).RX J2143.0+0654, for instance, has a broad main peak flanked by two weaker and much narrower peaks, which is an unusual feature for thermal radiation, but could in principle be a consequence of such a pencil+fan angular dependence of the surface radiation.However, a formal analysis using a realistic neutron star model is required to verify this assertion.Such an investigation will be presented in a subsequent paper in this series.

The Curious Case of RX J0720.4−3125
The presence of a significant second spin derivative in RX J0720.4−3125allows us to compute a pulsar braking index as n = ν ν/ ν2 .In the idealized model of a rotation-powered neutron star as a magnetic dipole spinning down in a vacuum, pulsar electrodynamics predicts a value of n = 3 (Gunn & Ostriker 1969).More realistic treatments of pulsar magnetospheres that incorporate wind outflows, magnetic field obliquity, and field evolution can produce values less than 3 (e.g., Blandford & Romani 1988;Melatos 1997;Harding et al. 1999;Ho 2015).However, the braking index for RX J0720.4−3125 based on the inferred ν, ν, and ν is an unphysical value of n = 682.This implies that the measured ν is not due to magnetic braking and could instead be due to timing noise (see, e.g., Parthasarathy et al. 2020, and references therein).An additional indicator of timing noise evident in Figure 1 is the significant scatter of the TOAs around the best fit model even after the inclusion of ν.
While previous investigations have found that most of the XINS show no evidence of long-term variability down to levels of a few percent (e.g., Kaplan & van Kerkwijk 2009a;Sartore et al. 2012;De Grandis et al. 2022), RX J0720.4−3125 has been reported to undergo spectral and profile changes on timescales of years (de Vries et al. 2004;Hohle et al. 2012;Borghese et al. 2015).Figure 6 shows the XMM-Newton profile in the 0.15-1.2keV range split into three MJD ranges: the first covering the interval prior to the assumed glitch epoch (MJD 51,575), the second after the glitch and up to the midpoint of the post glitch data span (MJD 56,598), and the third from this midpoint onward.It is apparent that the three profiles are not identical, with systematic differences in pulse amplitude, shape, and peak phase (bottom three panels Fig. 6), which cannot be explained merely by statistical fluctuations.The difference in peak phase is consistent with the level of scatter observed in the timing residuals shown in Figure 1.Thus, the likely explanation for the apparent ν and the scatter in the timing residuals is this gradual change in the pulsed emission over time, though the underlying process responsible for this atypical behavior remains to be understood.A thorough analysis of the spectral and profile variations in RX J0720.4−3125 that may shed more light onto this mystery is deferred to a subsequent work.
Based on an analysis of the set of XMM-Newton EPIC pn Full Frame mode observations, Hambaryan et al. (2017) proposed that the true spin period of RX J0720.4−3125 may actually be twice that of the commonly reported value.To verify this claim, we examined the power spectra of the deepest (65.1 ks) XMM-Newton observation from 16 October 2019 (ObsID 0852980201, obtained in Small Window mode) and the collection of NICER XTI exposures.We find no evidence for significant signal power around 0.0596 Hz, even when dividing the event data in the energy bands considered in Hambaryan et al. (2017) (see in particular their Fig. 5).We also carried out a timing analysis of the NICER XTI data alone restricted to events above 0.7 keV, where there is an apparent difference between the two pulses per rotation for the presumed 0.0596 Hz spin frequency.The NICER XTI data was chosen as it provides a uniform, high quality data set far away from the glitch and its span is such that a ν is not important.The highest likelihood solution found (ν = 0.0596 Hz, ν = −5.78× 10 −16 Hz s −1 ) results in an H-test significance that is lower than the best solution for the nominal spin at 0.119 Hz, implying that there is no significant contribution to the signal power from the pulsar at a frequency of 0.0596 Hz.This is confirmed by harmonic decomposition of the folded data, which reveals that the Fourier component of the putative fundamental frequency 0.0596 Hz is not required to accurately repro-duce the observed pulses.Moreover, at twice the nominal spin period, the two pulses per rotation are statistically indistinguishable and separated by exactly 0.5 in spin phase in all energy ranges considered.For reference, in the case of RX J1308.6+2127(see Section 4.4), the two pulses are not half a rotational cycle apart and their peak-to-peak separation changes depending on the choice of energy range (see Fig. 4).Based on these findings, we conclude that the commonly reported spin of P = 8.39 s is in fact the true rotation period of the RX J0720.4−3125neutron star.

CONCLUSIONS
We have presented a comprehensive timing reanalysis of all available X-ray data of the six nearby XINS with confirmed pulsations.In specific cases, the resulting timing solutions obtained provide orders of magnitude improvements in precision of the spin frequency and spin-down rate measurements compared to previous published results, while also extending the baseline to span over two decades.For RX J2143.0+0654, which previously had only a marginal measurement on its spin down, we obtained the first firm measurements of ν.Five of the XINS are well behaved in the sense that their long-term timing behavior is consistent with steady spin-down, with no clear evidence for measurable timing noise, higher order spin frequency derivatives, or glitches.The exception is RX J0720.4−3125, which is already known to have undergone a glitch around MJD 52,886, and exhibits broad residuals that can be accommodated by a second spin derivative (ν).The erratic spin behavior is possibly related to the previously identified long-term variations in the X-ray spectrum and pulse profile of this source, though the underlying physics responsible for these irregularities remains unknown.
The final timing solutions resulting from the analysis presented herein are included as supplementary products to this paper in the form of Tempo-style pulsar parameter (.par) files.The updated spin ephemerides allowed us to produce folded X-ray pulse profiles that reveal previously unseen details.The pulsed emission from the XINS shows surprising complexity, with pronounced asymmetries in the otherwise broad pulses expected from thermally emitting neutron stars.Moreover, the profiles show peculiar energy dependence, with significant shape changes, phase shifts in the pulse peak and, in some instances, enhancements in pulsed fraction at higher energies.These features are difficult to reconcile with a simple model of an isotropic blackbodylike radiation from small hot spots on the stellar surface.Thus, it may be necessary to invoke more sophisticated models that include the detailed beaming patterns of magnetic atmospheres and more complicated multitemperature hot regions on the neutron star surface.
In a follow up paper, we will present a detailed investigation of the spectra (both phase averaged and phase resolved) and their variation with time to gain additional insight into the surface temperature distribution, chemical composition, and magnetization of the stellar surface.Ultimately, we intend to employ detailed modeling methods of the extensive X-ray data sets to fully map the surface temperature distribution on the stellar surface for each XINS.
In the meantime, continued monitoring of the XINS sample is needed to maintain phase coherence of these timing solutions and to monitor for glitches and other timing anomalies.For RX J0420.0−5022,RX J0806.4−4123, and RX J2143.0+0654, in particular, additional data is necessary to improve the photon statistics of their pulse profiles, which would lead to tighter constraints on their surface emission properties.When coupled with the prospect of additional XINS discoveries with eROSITA (see, e.g, Kurpas et al. 2023Kurpas et al. , 2024)), this would provide new opportunities to glean insight into the properties of these puzzling objects and lead to an improved understanding of the Galactic neutron star population and uncertain aspects of neutron star evolution.
We thank Deven Bhakta for implementing parallel processing in the PINT event optimize code, and Paul Ray and Scott Ransom for helpful discussions.(Hunter 2007), PINT (Luo et al. 2021), SAS (Gabriel et al. 2004), Tempo2 (Hobbs et al. 2006)

APPENDIX
A. LOGS OF X-RAY OBSERVATIONS In Tables 4-9 we list the XMM-Newton, Chandra, and NICER observations used for the timing analysis presented in Section 3 for each target.Each table provides the numeric observation identifier (ObsID), the calendar date and UTC start time, the instrument, imaging mode/filter used (where applicable), and the net exposure time after applying the data filtering procedures described in Section 2. For a number of XMM-Newton observations, one or more of the detectors have no scientifically useful data due to the detector not being active or the filter wheel being closed, and are thus not listed.

B. TIMING SOLUTION POSTERIOR DISTRIBUTIONS
In Figure 7 and 8 we present the one and two-dimensional posterior distributions of ν and ν from the timing inference analyses in the form of corner plots.For RX J0720.4−3125, the corner plot also shows ν and the glitch parameters listed in Table 3.In each figure, the square marks the highest likelihood solution found.In the timing solution parameter (.par) files provided as supplementary material to this paper, the values for the fitted parameters correspond to these highest likelihood points.b For the Chandra HRC-S LETG observation, only the zeroth order (undispersed) point source emission was extracted for the timing analysis.and Small Window (SW).
b For the Chandra HRC LETG observation, only the zeroth order (undispersed) point source emission was extracted for the timing analysis.b NICER observations marked with an asterisk were acquired with a 1. ′ 47 pointing offset from RX J1856.5−3754 for instrument calibration purposes.b For the Chandra HRC-S LETG observation, only the zeroth order (undispersed) point source emission was extracted for the timing analysis.(bottom) from the event-based likelihood analysis described in Section 3. The contours correspond to the 0.5, 1, 1.5, and 2σ credible interval levels.The blue square marks the sample with the maximum likelihood that was used to fold the event data and produce the pulse profiles shown throughout the paper.The median value and credible intervals are reported in Table 2.  (ν), the glitch spin phase increment ∆Φ, the glitch permanent frequency increment ∆ν, the glitch permanent frequency first derivative increment ∆ ν, and the arbitrary pulse phase offset from the event-based likelihood analysis described in Section 3. The blue square marks the sample with the maximum likelihood that was used to fold the event data and produce the pulse profiles shown in Figure 2. The contours correspond to the 0.5, 1, 1.5, and 2σ credible interval levels The median value and credible intervals are reported in Table 3.

Figure 1 .
Figure1.Timing residuals expressed in seconds for RX J0720.4−3125 from a TOA analysis in Tempo2 using XMM-Newton (black), Chandra (red), and NICER (blue) observations.The top panel shows the best fit timing solution of all observations after MJD 52866 (marked by the vertical dashed line) without a glitch to illustrate the sudden change in spin behavior.The bottom panel shows the best fit with a glitch included with an assumed epoch of MJD 52886.See Table3for the timing model used and text for further details.

Figure 2 .
Figure2.Normalized pulse profiles of the six XINS obtained by folding and stacking all XMM-Newton EPIC pn and MOS1/2, Chandra HRC-S and ACIS, and NICER (where available) X-ray event data at the timing solutions reported in Tables2 and 3. See text for the energy range used for each profile.A rotational cycle is repeated for clarity.

Figure 3 .
Figure3.The folded profiles from Figure2but with the most significant Fourier components superposed.Error bars are omited to improve clarity and a single rotational cycle is shown.The solid red lines show the best fit with a model of the profile constructed from the empirical Fourier coefficients of the set of photon phases.The dashed lines show the sinusoids corresponding to the first harmonically related components of the fit (blue, orange, purple, and cyan for the first, second, third, and fourth harmonic, respectively).The bottom panel for each pulsar shows the residuals from the fit.
we present two-dimensional histograms of pulsed counts as

Figure 4
Figure4.Normalized pulse profiles of the six XINS split into multiple energy bands.For RX J0720.4−3125 and RX J1308.6+2127,profiles for three energy intervals are shown in orange, blue and purple for the soft, mid, and hard bands, respectively.For the rest, the binned profile data are split into two energy ranges with the softer and harder bands corresponding to the orange and blue curves, respectively.The superposed solid red, blue, and purple lines in the second rotational cycle show the best fit with a representation of the profile constructed from the empirical Fourier coefficients of the set of photon phases for the same energy range.

Figure 5 .
Figure 5. Two-dimensional histograms of pulsed counts versus rotational phase and photon energy for RX J0720.4−3125,RX J1308.6+2127, and RX J1856−3754 (from top to bottom, respectively).The color bar shows the number of pulsed counts in each pixel with counts increasing from darker to lighter colors.The horizontal striping is an artifact produced by instrumental edge features.Two pulse phase cycles are shown for clarity.

Figure 6 .
Figure6.Top panel: XMM-Newton EPIC pn+MOS1/2 pulse profiles of RX J0720.4−3125 from three MJD ranges: one including all data prior to the glitch (orange) and the other two using events from the first (blue) and second (purple) halves of the post-glitch data set.The solid red, blue, and purple lines show the model based on the empirical Fourier coefficients of the phase folded event data.Bottom three panels: Differences in normalized counts between the profiles from the second and first (upper panel), third and first (middle panel), and third and second (lower panel) MJD intervals from the top panel.

Figure 7 .
Figure 7. Frequency (ν), frequency derivative ( ν), and pulse phase posterior distributions inferred for RX J0420.0−5022(top left), RX J0806.4−4123(top right), RX J1308.6+2127(middle left), RX J1856.5−3754(middle right), and RX J2143.0+0654(bottom) from the event-based likelihood analysis described in Section 3. The contours correspond to the 0.5, 1, 1.5, and 2σ credible interval levels.The blue square marks the sample with the maximum likelihood that was used to fold the event data and produce the pulse profiles shown throughout the paper.The median value and credible intervals are reported in Table2.

Figure 8 .
Figure8.Posterior distributions inferred for RX J0720.4−3125 for frequency (ν), frequency first ( ν) and second derivative (ν), the glitch spin phase increment ∆Φ, the glitch permanent frequency increment ∆ν, the glitch permanent frequency first derivative increment ∆ ν, and the arbitrary pulse phase offset from the event-based likelihood analysis described in Section 3. The blue square marks the sample with the maximum likelihood that was used to fold the event data and produce the pulse profiles shown in Figure2.The contours correspond to the 0.5, 1, 1.5, and 2σ credible interval levels The median value and credible intervals are reported in Table3.

Table 1 .
Total filtered X-Ray Exposures (in Kiloseconds) Used in This Work
′′ 5 µα cos δ (mas yr −1 ) The numbers quoted for the fitted parameters correspond to the 50 th percentile value of the posterior distribution and the −/+ values correspond to the 16 th (−) and 84 th (+) percentiles of the posterior credible intervals.Numbers quoted in parentheses give the uncertainty in the last digit of the nominal value.
This work was funded in part by NASA Astrophysics Data Analysis Program (ADAP) grant No. 80NSSC20K0884.We acknowledge use of computing resources from Columbia University's Shared Research Computing Facility project, which is supported by NIH Research Fa-cility Improvement Grant 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) Contract C090171, both awarded April 15, 2010.The work presented is based in part on observations with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member states and NASA, and on observations from the NICER mission, which is funded by NASA.The XMM-Newton Science Analysis Software (SAS) is developed and maintained by the Science Operations Centre at the European Space Astronomy Centre and the Survey Science Centre at the University of Leicester.This research has made use of data obtained from the Chandra Data Archive, contained in the Chandra Data Collection (CDC) 224 doi:10.25574/cdc.224, and software provided by the Chandra X-ray Center (CXC) in the application package CIAO, as well as data products 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.We acknowledge extensive use of NASA's Astrophysics Data System (ADS) bibliographic services and the ArXiv.

Table 4 .
Log of X-ray observations of RX J0420.0−5022 a Listed exposure times are values after filtering for instances of strong background flaring.EPIC MOS1/2 observations acquired in Full Window mode are not used due to the coarse 2.6 s time resolution.The imaging mode abbreviations correspond to Full Window (FW) and Small Window (SW).

Table 7 continued
Table 7 (continued) Listed exposure times are values after filtering for instances of strong background flaring.The imaging mode abbreviations correspond to Full Window (FW)

Table 8 continued
Table 8 (continued) Listed exposure times are values after filtering for instances of strong background flaring.The imaging mode abbreviations correspond to Full Window (FW), Large Window (LW), Small Window (SW), and Fast Timing (FT).

Table 9 .
Log of X-ray observations of RX J2143.0+0654Listed exposure times are values after filtering for instances of strong background flaring.The imaging mode abbreviation SW is for Small Window.