New Timing Results of MSPs from NICER Observations

Millisecond pulsars (MSPs) are known for their long-term stability. Using six years of observations from the Neutron Star Interior Composition Explorer (NICER), we have conducted an in-depth analysis of the X-ray timing results for six MSPs: PSRs B1937+21, B1821$-$24, J0437$-$4715, J0030+0451, J0218+4232, and J2124$-$3358. The timing stability parameter $\sigma_z$ has been calculated, revealing remarkable timing precision on the order of $10^{-14}$ for PSRs B1937+21 and J0437$-$4715, and $10^{-13}$ for PSRs B1821$-$24, J0218+4232, and J0030+0451 over a timescale of 1000 days. These findings underscore the feasibility of autonomous in-orbit timekeeping using X-ray observations of MSPs. In addition, the consistency of long-term spin-down noise in the X-ray and radio bands has been investigated by comparison with IPTA radio data.


Introduction
Pulsars, particularly MSPs, exhibit long-term stability in their pulsation periods, which rivals that of current atomic timescales [1][2][3][4][5][6].Consequently, high-precision pulsar timing can constrain the mass-radius relationship of neutron stars (NSs) and elucidate their emission mechanisms [7,8].Moreover, it can facilitate the detection of nanohertz gravitational waves from supermassive black hole binaries well in advance of system coalescence [9,10], and enable precise measurement of solar system planetary masses [11].Furthermore, given the separate emission regions for radio and X-ray radiation, investigation of the consistency of spin-down noise in the X-ray and radio bands will underscore the understanding the origins of pulsar emissions, the characteristics of timing noise, and the evolutionary processes involved [12][13][14].
Pulsars also play an important role in positioning, navigation, and timing (PNT).As the number of spacecraft launched into Earth's orbit, the vicinities of the Moon and Mars, and beyond into deep space continues to grow, the demand for PNT capabilities in space is becoming ever more critical.Current navigation and timekeeping for spacecraft in low-Earth orbits predominantly depend on Global Navigation Satellite Systems (GNSSs) and ground-based tracking, whereas deep-space missions rely primarily on radio technologies and space-Earth time comparisons.Pulsar-based navigation, a vital auxiliary solution, has been proposed and successfully tested in orbit [15][16][17][18][19][20].In addition, pulsars will offer the potential for timekeeping services as a 'cosmic clock', compensating for biases in satellitemounted atomic clocks or oscillators due to varying space environments, long-term drift, inadequate relativity corrections, or infrequent synchronization with Earth's high-stability clocks [21].Pulsars stand as prime candidates to provide an external absolute time standard, particularly for space missions in diverse orbits.A pulsar-based timescale, utilizing a pulsar timing array, has been proposed ,employing high-precision times of arrival (TOAs) detected by ground-based radio telescopes [6].X-ray timing offers unique advantages over radio due to its negligible dispersion measure (DM) effects, miniaturized detectors, and independence from ground-based systems, despite sharing similar pulsar parameters and intrinsic noise profiles with radio observations.Frequency steering of spaceborne clocks has been demonstrated based on XPNAV-1 observations [22].However, the accuracy of X-ray TOA detections is still a bottleneck, and a higher-precision demonstration still needs to be conducted.
The Neutron Star Interior Composition Explorer (NICER), with its unparalleled capability to detect millisecond pulsars, enables detailed analyses of the time-stability of X-ray pulses, assessing their utility for on-orbit timekeeping.The six X-ray MSPs, including PSRs B1937+21, B1821−24, J0437−4715, J0030+0451, J0218+4232, and J2124−3358, demonstrate high X-ray fluxes, small rotation periods, and minimal period derivatives, coupled with high stability (although PSR B1937+21 and B1821−24 show low-frequency timing irregularities dominate the timing residuals [3,23]), rendering them excellent candidates for timekeeping purposes.Extensive observational data of NICER for these six MSPs span over six years.Here, the new timing results of the MSPs with more than six years of observations from NICER have been obtained and demonstrate the feasibility of using pulsars for autonomous space navigation.

Observation
NICER, an X-ray observatory affixed to a movable arm on the exterior of the International Space Station (ISS), has been operational since June 2017 [24].Designed specifically to investigate the X-ray emissions of NSs, NICER's primary objective is the high-precision timing of X-ray-emitting pulsars to constrain NSs' mass-radius relationships and to study their high-energy emission mechanisms.The exposures of PSRs B1937+21, B1821−24, J0437−4715, J0030+0451, J0218+4232, and J2124−3358 were around several mega-seconds over the period from June 2017 to September 2023, as listed in Table 1, which could be utilized to supply more accurate timing results.
The events from the six pulsars have been successfully filtered with the above criteria.For PSR B1937+21, a total of 543,133 photons meet our selection criteria across 1.65 Ms of net exposure time, yielding an average count rate of 0.33 counts/s.In the case of PSR B1821−24, we have identified 407,839 photons over 0.67 Ms of clean exposure, with an average count rate of 0.61 counts/s.For PSR J0218+4232, the data comprise 819,464 photons across 1.85 Ms of net exposure, averaging 0.44 counts/s.PSR J0437−4715 yielded 4,681,958 photons over 2.50 Ms of clean exposure, resulting in an average count rate of 1.87 counts/s.For PSR J0030+0451, we obtained 2,679,284 photons over 3.70 Ms of clean exposure, with an average count rate of 0.72 counts/s.Lastly, PSR J2124−3358 contributed 1,018,017 photons over 2.09 Ms of clean exposure, averaging 0.49 counts/s.It is important to note that these count rates include both source and background contributions.

ToA Calculation
To obtain ToAs of each pulsar, the arrival time of each photon at the local observatory is corrected to the solar system barycenter (SSB) and then folded with the pulsar's ephemeris.The initial pulsar ephemerides are sourced from the International Pulsar Timing Array (IPTA) Data Release 2 (DR2) [26] and refined with NICER data [27].Owing to decades of monitoring, the positions, proper motions, parallaxes, and binary orbits of these pulsars have been determined with high precision by radio telescopes.Consequently, we fix these parameters in our analysis.The rotational period and its derivative serve as initial inputs and are updated with the X-ray observations.Our timing analysis employs the DE436 JPL ephemeris, Barycentric Dynamical Time (TDB) units, and the TT(BIPM2022) clock realization.Additionally, we have incorporated the GNU Multiple Precision Arithmetic Library (GMP) (https://gmplib.org/(accessed on 13, December, 2023)) to facilitate highprecision data recording at the 1 ns level.
For each pulsar, we have constructed a high-precision "standard" pulse profile using six years of data.Due to the substantial variation in photon counts across different observation IDs (ObsIDs)-with some containing fewer than 100 photons-we have segmented the entire dataset based on count rates.Each segment comprises approximately 50,000 to 200,000 photons within a 30-day span.We obtain pulsed profiles for each segment and calculate ToAs through cross-correlation analysis with the standard profile.The ToA error is calculated by performing Gaussian sampling across the pulse profile.The pulsar period parameters and timing residuals are then refined using Tempo2 [28,29].

Timing Results
The frequency and its first-order derivative in each pulsar ephemeris have been updated by refitting all ToAs, while the other parameters remain fixed as in IPTA.The higher-quality pulse profiles of the six pulsars are shown in Figure 1.With the six years of observational data, we have been able to obtain a highly accurate pulse profile.PSRs B1937+21, B1821−24, and J0218+4232 consist of a main pulse and an interpulse approximately 180 degrees apart in phase.For PSRs J0437−4715, J0030+0451, and J2124−3358, the profiles exhibit an approximate sinusoidal shape, indicative of radiation primarily from one or more hot spots on the surface of the neutron star.Additionally, the asymmetric pulse profiles and the distinct "hump" observed in PSRs J0437−4715 and J2124−3358 are attributed to a secondary hot spot [30].
The timing residuals obtained with the updated pulsar ephemerides may reflect various factors, including pulsar timing noise, the gravitational wave background, clock errors, solar system ephemeris and pulsar ephemeris errors, and other uncalibrated instrumental effects.The dispersion measure (DM) and its variation have a negligible impact for X-ray observation as mentioned above.The gravitational wave background is also considered negligible [10].The event times have been calibrated and corrected with a precision of about 40 ns and are referenced to UTC with an accuracy better than 100 ns [31].The influence of the solar system ephemeris is expected to reach about 100 ns over typical periods of 10 years [32].Furthermore, inaccuracies in the IPTA radio pulsar ephemerides arise from covariances between certain parameters and variations in the DM, which could potentially introduce unaccounted-for errors in our timing measurement.Consequently, the timing residuals of the six pulsars, containing all the factors mentioned above, are shown in Figure 2, ranging from a few microseconds to tens of microseconds.For PSR B1937+21, we observe the lowest X-ray band residuals, with the RMS within 1 µs due to the smallest spin period and the sharpest (and narrowest) profile, as illustrated in Figure 2. Therefore, PSR B1937+21 holds significant potential for applications in X-ray timekeeping.Low-frequency timing irregularities have been observed over a period of more than three years [3].However, if these timing irregularities can be accurately modeled and removed, it will raise the possibility of producing an essentially stable X-ray clock [12].
The second derivative of the spin frequency, ν, which may have intrinsic or kinematic origins in MSPs, is an important consideration for high-precision pulsar timing.Therefore, we have refitted the ephemerides of the six pulsars, including the second derivatives of frequency, using TEMPO2.A statistically significant non-zero second derivative of the frequency (with a significance greater than 3σ) is detected only in PSRs B1937+21 and B1821−24.As shown in Figure 2, the cubic structures in the timing residuals are suppressed effectively with ν included, especially for PSR B1821−24.The root mean square of the timing residuals of PSR B1937+21 decreases from 0.9 µs to 0.6 µs while it decreases from 16.7 µs to 6.7 µs for PSR B1821−24.

Evaluation of Rotational Stability
Clock stability is typically described by the Allan variance, σ 2 y , which is designed to quantify instability in clocks operating at a constant frequency.However, pulsars are known to have a frequency derivative due to the gradual slowdown of their rotation.According to the timing analysis, we can remove the effects of pulsar frequency, frequency derivative, proper motion, parallax, and other systematic influences.The remaining timing residuals include the measurement error, uncalibrated instrumental error, and intrinsic timing noise, the latter providing insight into the rotational stability of the pulsar and determining its potential as a precise timekeeping clock.To evaluate the rotational stability, another parameter, σ 2 z , has been introduced to evaluate the rotational stability by calculating the variations in the timing residuals [2].
For an interval of length τ starting at time t 0 , we can fit a cubic polynomial to the timing residuals in that interval, where X(t) minimizes the sum of [(x i − X(t i ))/σ i ] 2 over all TOA t i with residuals t i and uncertainties σ i .Then, where angle brackets denote the weighted average over the third-order coefficients of the best-fit polynomials of all non-overlapping intervals of length τ within the set of timing residuals.
The timing stability, σ z , is plotted as a function of the data span, τ, for the six millisecond pulsars and is presented in Figures 3-8.For a direct comparative analysis, the results are juxtaposed with observations from the radio band.Furthermore, in the scenario where the residuals are solely influenced by white noise, we expect σ z to scale with τ as σ z ∝ τ −1.5 , as described by [2].This relationship is depicted by the dashed lines included in each of the figures for reference.Despite some deviation at lower time spans, it demonstrate comparable stability in both the radio and X-ray bands for PSRs B1937+21, B1821−24, and J0218+4232, showing the consistence of the timing between the X-ray and radio bands.However, for PSRs J0030+0451, J0437−4715, and J2124−3358, the stability indices in the X-ray band are lower than those in the radio band, with the largest discrepancies exceeding a factor of 100 [25,36].As shown in Figure 2, the time residuals for the latter three pulsars are larger due to the broad and sinusoidal-like profiles, which suggests that the timing stability is likely still limited by white noise.However, other possible explanations corresponding to the X-ray observations will be considered.Consequently, a varying number of detectors from NICER were utilized in the analysis to evaluate the timing stability in relation to different detector areas.As illustrated in Figure 9, an increase in detector area leads to enhanced timing stability, thereby reinforcing the predominance of white noise.Consequently, for PSRs J0030+0451, J0437−4715, and J2124−3358, a larger detection area, more observation time, and reduced background noise are required to suppress white noise and achieve results comparable to those in the radio band.

Summary and Conclusions
X-ray pulsar-based timekeeping presents a promising application for space missions, with the long-term stability of millisecond X-ray pulsars (MSPs) being a critical parameter.Utilizing six years of observations from NICER, the timing stability, σ z , has been assessed for six candidate pulsars: PSRs B1937+21, B1821−24, J0437−4715, J0030+0451, J0218+4232, and J2124−3358.Despite the presence of red noise and/or white noise, long timing stability of 10 −14 for PSRs B1937+21 and J0437−4715, and 10 −13 for PSRs B1821−24, J0218+4232, and J0030+0451 has been achieved over a span of 1000 days, supporting the feasibility of timekeeping.It is important to note that a thorough investigation of the red noise and efforts to mitigate or minimize its impact on timing residuals will further enhance the utility of this approach.Moreover, by comparison with IPTA radio data, the consistency of long-term spin-down noise in the X-ray and radio bands has been investigated.It shows a similar upturn for PSRs B1937+21, B1821-24, and J0218+4232, which could provide meaningful information on the emission and spin-down mechanisms of neutron stars.

Data Availability Statement:
All the data used in this paper can be found at NASA's HEASARC website, https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl.

Figure 3 .
Figure3.Measure of timing stability, σ z , vs. the data span over which a third-order polynomial is fitted to the timing residuals of PSR B1937+21.Black circle points show σ z for NICER data, and red square points show the radio results from[25].A dashed line shows the slope for the case where timing precision is limited by white noise only.

Figure 6 .
Figure 6.Same as Figure 3 but for PSR J0030+0451.The red square points show the radio results with IPTA DR2 dataset.

Figure 7 .
Figure 7. Same as Figure 3 but for PSR J0437−4715.The red square points show the radio results with IPTA DR2 dataset.

Figure 8 .
Figure 8. Same as Figure 3 but for PSR J2124−3358.The red square points show the radio results with IPTA DR2 dataset.

Figure 9 .
Figure 9. Same as Figure3but for PSR J0437−4715.A varying number of detectors from NICER were employed in the analysis: black squares represent data from a single detector, black triangles correspond to data from ten detectors, and black circles indicate data from all but three of the hottest detectors, totaling fifty-three.

Author Contributions:
Conceptualization, S.Z.(Shijie Zheng), D.H., and M.G.; methodology, S.Z.(Shijie Zheng), H.X., K.L., J.Y., and L.Z.; software, H.W.; validation, Y.L. and Y.Y.; formal analysis, X.M.; writing-original draft preparation, S.Z.(Shijie Zheng) and D.H.; writing-review and editing, M.G., Y.C., and S.Z.(Shuangnan Zhang).All authors have read and agreed to the published version of the manuscript.Funding: This work is supported by the National Key R&D Program of China (2022YFF0711404, 2021YFA0718500), and the National Natural Science Foundation of China (Grant No. 12333007).Heng Xu is supported by the Project funded by China Postdoctoral Science Foundation No. 2023M743518, Kejia Lee and Heng Xu are supported by National SKA Program of China 2020SKA0120100.

Table 2 .
The timing parameters of the six MSPs.