© The Authors 2025

1. Introduction

International Atomic Time (TAI) is a timescale established and maintained by the Bureau International des Poids et Mesures (BIPM) using about 460 atomic clocks and several frequency standards from about 89 timing laboratories worldwide. TAI represents an ensemble timescale that is susceptible to accumulated errors capable of influencing its long-term stability. Pulsars emit steady periodic pulse signals outward, thereby establishing a timescale known as pulsar time (Petit 1996; Taylor 1991). They demonstrate remarkable long-term stability, with the highest stability surpassing 10⁻¹⁶ (Allan 1987), which renders them ideal for monitoring the enduring stability of TAI (Hartnett & Luiten 2010; Li et al. 2020). The stability of pulsar timing refers to the degree of random fluctuations in the rotation frequency of pulsars influenced by various sources of noise. It indicates the deviation of a pulsar’s rotation frequency or phase relative to the average frequency or phase over a specified time period, which serves as a crucial metric for assessing pulsar timing performance. Furthermore, pulsar timing can be employed to improve the solar system planetary ephemeris by detecting perturbations of the Earth (Foster & Backer 1990), while also offering essential evidence for the detection of gravitational waves (Weisberg & Taylor 1981; Jenet et al. 2005).

After recording the pulse times of arrival (TOAs) from a specific pulsar at a radio observatory, these pulse TOAs should be calibrated to accurately reflect the arrival times in an inertial reference frame (Edwards et al. 2006). Once the impact of both intrinsic and extrinsic noise sources are mitigated (including effects arising from the relative motion between the pulsar and observatory, instrumental errors, and signal propagation delays), the pulse TOAs can be regarded as an equivalent comparison between the pulsar’s clock and a reference atomic clock. Generally, public research on pulsar timing mainly relies on radio observations. In 2010, the International Pulsar Timing Array (IPTA) was formally established (Hobbs et al. 2010). The primary objective of IPTA is to detect and characterize low-frequency gravitational waves by timing approximately 100 millisecond pulsars (MSPs) globally and to establish a pulsar-based reference timescale. The current formal members of this organization include the Australian Parkes Pulsar Timing Array (PPTA) (Reardon et al. 2021), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (Alam et al. 2020), the European Pulsar Timing Array (EPTA) (Desvignes et al. 2016; EPTA Collaboration 2023; Chen et al. 2021), and the Indian Pulsar Timing Array (InPTA) (Paul et al. 2019; Tarafdar et al. 2022; Srivastava et al. 2023). Observing members consist of the MeerKAT Pulsar Timing Array (MPTA) in South Africa (Johnston et al. 2020; Miles et al. 2022; Spiewak et al. 2022) and the Chinese Pulsar Timing Array (CPTA) (Lee et al. 2012; Xu et al. 2023). The IPTA team has conducted long-term observations of over 60 MSPs (Manchester 2013), with some accumulating more than 20 years of data. The ensemble pulsar timescal (EPT) constructed using IPTA’s DR1 (Verbiest et al. 2016) or DR2 (Perera et al. 2019) data demonstrates long-term stability comparable to or slightly better than atomic time (Zharov et al. 2019; Rodin & Fedorova 2018). After removing linear and quadratic terms, it exhibits trends similar to TT (BIPM), which accurately reveals high-order systematic errors in TAI (Hobbs et al. 2012, 2020), and thereby demonstrates the viability of pulsar timing for monitoring atomic time (Zhao et al. 2023). Given that the precision of pulsar timing is largely dependent on the receiving antenna area, the completion of Five-hundred-meter Aperture Spherical radio Telescope (FAST) in recent years has provided an excellent opportunity for establishing high-precision ensemble pulsar timescale (EPT) (Zhou et al. 2021) simulated and evaluated ten-year timing datasets for several MSPs based on FAST observations, and achieved a stability of 1.39 × 10⁻¹⁶.

In contrast to radio observations, X-ray observations conducted in outer space avoid atmospheric dispersion effects, thereby enhancing the precision of pulse TOA measurements and making X-ray data more suitable for constructing comprehensive pulsar timing. Due to the extremely faint signals emitted by most MSPs in the X-ray band and the inherent limitations of X-ray observations, such as the difficulty in increasing the detector area, several pulsars known for their exceptional stability in the radio band (e.g., J0437-4715) exhibit significantly reduced signal-to-noise ratios (S/N) in the X-ray band. Additionally, many MSPs display thermal X-ray pulsations with considerably broader peaks compared to those observed in the radio band (Deneva et al. 2019; Bogdanov 2017), which further impacts the precision of pulse TOA measurements. However, two notable exceptions, PSR B1937+21 and PSR B1821-24, exhibit remarkably narrow nonthermal X-ray pulsations, making them highly suitable candidates for EPT (Kawai et al. 1998). These two pulsars have undergone extensive long-term observations by several prominent radio telescopes, which has led to significant successes in evaluating their stability (Backer et al. 1982; Kaspi et al. 1994; Lyne 1987; Cognard et al. 1996; Hazboun et al. 2022). Nonetheless, elevating their stability and constructing EPT in the X-ray spectrum still remain challenging due to weak X-ray emissions and limited detector areas. The situation has seen significant improvement with the emergence of the Neutron Star Interior Composition Explorer (NICER), due to its robust X-ray observation capabilities and a remarkable time resolution of less than 100 ns (Gendreau et al. 2016; Prigozhin et al. 2016; Gendreau & Arzoumanian 2017). Deneva et al. 2019 assessed the one-year stabilities of PSR B1937+21 and PSR B1821-24 using NICER observations, and obtained stabilities of 3.1 × 10⁻¹⁴ and 1.5 × 10⁻¹², respectively. However, longer-term stability investigations have not been conducted due to the limited duration of the data. Sun et al. 2023 determined that the five-year NICER stability of PSR B1937+21 is 1.7 × 10⁻¹⁴, which demonstrates that the pulsar’s rotation exhibits exceptional stability over an extended observation period.

To address the current gap in the research on constructing EPT for X-ray data, this paper introduces a comprehensive X-ray pulsar timing processing framework and establishes EPT over a 7.3-year period for PSR B1937+21 and PSR B1821-24. The key contributions of this study can be summarized as follows. Firstly, this study processed and obtained the 7.3 years NICER timing residuals for the two pulsars, PSR B1937+21 and PSR B1821-24. Secondly, an X-ray pulsar timing processing framework was introduced, which includes a comprehensive noise analysis component that detected significant evidence of red noise present in both pulsars using Bayesian factors, and employed noise model fitting to better model ensemble timing residuals, demonstrating a substantial improvement over the classical cubic spline interpolation method. Lastly, an enhanced Wiener filtering algorithm was proposed for constructing ensemble pulsar timing, specifically designed to mitigate the high white noise levels inherent in X-ray observations.

The remainder of this paper is structured as follows. Section 2 details the processing procedures for NICER data. Section 3 describes the framework and algorithms developed for X-ray pulsar timing. Section 4 presents an analysis of the experimental results. Finally, Section 5 provides the conclusions and discussion.

2. Observations and data processing

Most MSPs typically exhibit low flux in the X-ray band, which results in a low signal-to-noise ratio of the pulse profile and subsequently affects the accuracy of pulse TOAs. Thus, the implementation of certain preprocessing techniques becomes essential to augment the S/N of pulsar profiles, which include the screenings of the energy and flux rate ranges, and S/N (Sun et al. 2023). In terms of optimal energy range screening, we employed a step size of 0.05 keV and calculated the chi-square values of pulsar profiles across various energy segments. The optimal energy ranges determined for the two pulsars are 1.1–6.55 keV for PSR B1937+21 and 0.8–5.5 keV for PSR B1821-24. For optimal flux rate range screening, we employed an eight-second observation interval with a flux rate threshold of 2 cnts/s. Segments with a photon flux rate that surpasses this threshold were excluded, since it indicates prevalent noise. Additionally, we excluded individual observations with notably low S/Ns or inadequate photon counts. The data preprocessing results for PSR B1937+21 and PSR B1821-24 are presented in Table 1.

We performed time corrections using the “barycorr” command within the HEASoft software. The Einstein delay was modeled using 127 sine periodic terms, with a precision of up to 100 ns (Fairhead & Bretagnon 1990). Due to the relatively lower precision of pulse TOAs in X-ray timing data, an enhanced parameter fitting framework was adopted to iteratively refine timing parameters by incorporating TOAs across the entire observation span (Sun et al. 2023). The timing residuals for the two pulsars, obtained by fitting the timing parameters, are illustrated in Figure 1. The residuals for PSR B1937+21 show a smaller amplitude of red noise, with a root mean square (RMS) of 4.17 μs, while PSR B1821-24 exhibits significant fluctuations and strong red noise, with a RMS of 25.07 μs. The fitted timing parameters are detailed in Table 2.

Fig. 1. Timing residuals after fitting the timing model: (a) PSR B1937+21, (b) PSR B1821-24.

To compare the stability between NICER and radio data, we analyzed NANOGrav data for PSR B1937+21 and PPTA data for PSR B1821-24, as depicted in Figure 2. For PSR B1937+21, the timing residuals in the radio data are generally within ±5 μs, although certain data points exhibit larger error bars exceeding 15 μs. The timing residuals for PSR B1821-24 fall within the range of ±30 μs. In the radio timing residuals of PSR B1821-24, previous studies have demonstrated the presence of discernible red noise (Vivekanand 2019; The NANOGrav Collaboration 2019).

Fig. 2. Radio timing residuals for PSR B1937+21 and PSR B1821-24. (a) PSR B1937+21, (b) PSR B1821-24.

3. Ensemble pulsar timing analysis

Pulsar timing residuals are affected by various sources of stochastic noise, with many of these noise components being independent from each other. To mitigate the influence of independent noise in timing data, the ensemble averaging technique can be implemented, for which the most classic algorithm is the weighted average method (Kaspi 1995; Petit et al. 1993). By conducting this method, the constructed EPT achieved a stability of 2 × 10⁻¹⁴ over a period of 2.5 years, outperforming the stability calculated from individual pulsars (Petit 1996). Afterward, wavelet filtering and Wiener filtering algorithms were successively applied to EPT. Following the evaluation of seven-year data from the Arecibo Observatory for PSR B1855+09 and PSR B1937+21, the Wiener filtering algorithm exhibits excellent performance, achieving a stability of 1.7 × 10⁻¹⁵ (Zhong & Yang 2007). Furthermore, a joint algorithm combining wavelet decomposition with Wiener filtering was proposed and achieved a level of stability surpassing that of other algorithms (Zhong & Yang 2007). In addition, a bispectral filtering algorithm was introduced for constructing ensemble pulsar timing. The results show that the bispectral filtering algorithm can effectively reduce the influence of observation noise and improve the stability of ensemble pulsar time compared to the classic weighting method (Zhou et al. 2021). By utilizing 8.2-year data obtained from the above two pulsars, Rodin improved the stability to (0.8 ± 1.9) × 10⁻¹⁵ (Rodin 2008).

3.1. Ensemble pulsar timing framework

In the construction of ensemble pulsar timing, it is essential to have simultaneous observational data from different pulsars. Radio observations typically employ cubic spline interpolation to connect these data points. Given the high precision of radio TOAs, this method introduces minimal error. However, X-ray observations, which are subject to significant white noise and large gaps between observations, present a challenge. In such cases, cubic spline interpolation can result in substantial errors due to the large gaps between the data points. Moreover, this method generates only one set of deterministic timing residuals for stability analysis, which renders it insufficient for accurately representing the pulsar’s stability. To address this, we developed a framework for generating simulated residuals that incorporates red noise estimation, as illustrated in Figure 3. The framework consists of modeling both red noise and white noise, as well as extracting residual trends (such as clock offsets and gravitational wave signals). After removing the red noise, we rely on the randomness of the white noise and fixed trend components to generate simulated residuals, and we conduct multiple Monte Carlo experiments to achieve more accurate statistical results.

Fig. 3. Framework for generating ensemble pulsar timing residuals.

The framework consists of two main components: noise modeling and residual simulation. Understanding timing noise is crucial for studying pulsars and advancing pulsar timing applications (Arzoumanian et al. 2018). The first step is to determine the presence of red noise in timing residuals by calculating the Bayesian factor. White noise has a uniform power spectral density, whereas red noise is concentrated at lower frequencies. To accurately identify variations related to the reference clock, it is essential to extract and mitigate the influence of red noise. If red noise is detected, it is first modeled using 100 harmonic components. The optimal phase combination is then determined by minimizing the chi-square value, which enables the reconstruction of the red noise time-domain signal. Furthermore, to overcome two major limitations of traditional methods, we construct both a white noise model and a residual low-frequency signal model, which are then used to generate simulated residuals for Monte Carlo analysis. Typically, red noise can be represented using a power-law model

$$\begin{array}{c}\hfill P_{\text{RN}}=\frac{A_{\text{RN}}^2}{12{\pi }^2}{\left(\frac{f_T}{f_{\text{yr}}}\right)}^{-\gamma }{\text{yr}}^3.\end{array}$$(1)

Here, P_RN represents the power spectral density (PSD) of red noise, A_RN denotes the amplitude of red noise, γ is the spectral index, f_T is taken as (1/T, n/T), T signifies the total length of observational data, n is typically set to 30, yr represents the number of seconds in a year, and f_yr = 1/yr.

Achromatic red noise in pulsar data is a stochastic signal, making Gaussian process regression (Rasmussen & Williams 2006) a widely used modeling approach. To characterize this noise, a normal kernel process is applied in the Fourier domain, and follows a power-law model. This method ensures an accurate statistical representation of the noise and its temporal evolution. Bayesian factors are utilized to assess the relative goodness of fit between the red noise model and the white noise model (Hazboun et al. 2022). When the Bayesian factor exceeds three, we obtain strong evidence for the presence of red noise in the timing data.

Specifically, we construct the noise model using the Enterprise (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE)¹ library and employ Dynesty to calculate the evidence for the red-noise likelihoods. Enterprise is a pulsar timing analysis code that is designed for noise analysis, gravitational-wave searches, and timing model analysis (Ellis et al. 2019). For parameter estimation, we utilize the standard Markov Chain Monte Carlo (MCMC) sampler, PTMCMCSAMPLER.

The power-law models primarily describe the spectral index and amplitude of red noise, and effectively characterize its power spectrum. However, they lack critical phase information, which is essential for accurately reconstructing the red noise in the timing data. This limitation makes it difficult to capture the full complexity of the red noise using power-law models alone.

To address this, we use harmonic fitting, which provides a more comprehensive approach. The amplitudes and frequencies of the harmonics can be directly mapped to the PSD of the red noise, which allows for a more precise representation of its characteristics. Furthermore, by minimizing the chi-square value, we can accurately determine the optimal phase combination of the harmonics, which enables a more precise reconstruction of the red noise time-domain signal. This method not only captures the full complexity of red noise but also enhances the accuracy of timing residual analysis, particularly in identifying clock-related variations.

3.2. Mirror-symmetric Wiener filtering method based on Levinson-Durbin iteration

The current ensemble pulsar timing algorithms are primarily suitable for radio timing data with strong red noise and weak white noise characteristics. However, X-ray timing data exhibits strong white noise. To address this, we used an enhanced time-domain Wiener filtering algorithm based on Levinson-Durbin iteration, which effectively reduces white noise by estimating AR model parameters directly from the observed signal (Rasmussen & Williams 2006). By accurately modeling the signal’s autocorrelation structure, the Wiener filter can differentiate between the signal of interest and noise components, which are assumed to be uncorrelated and exhibit a flat power spectral density (white noise). This approach is widely applied in fields such as audio processing, telecommunications, and image enhancement due to its ability to tailor its response to the specific statistical properties of the observed signal, leading to significant noise reduction and improved signal quality.

We set r = s + v, where r represents the timing residuals, s represents the signal to be estimated, and v represents the white noise. Then the signal to be estimated is denoted as $\widehat{\mathit{s}}={\mathit{w}}^{\ast \text{T}}\mathit{x}$, where ${\mathit{w}}^{\ast }={\widehat{\mathit{R}}}_{\mathit{rr}}^{-1}{\widehat{\mathit{r}}}_{\mathit{rs}}$. Here, ${\widehat{\mathit{R}}}_{\mathit{rr}}$ is the autocorrelation matrix of the signal r, and autocorrelation sequence ${\widehat{r}}_{\mathit{rr}}\left(k\right)$ can be described as

$$\begin{array}{c}\hfill {\widehat{r}}_{\mathit{rr}}\left(k\right)=\frac{1}{N-k}{\displaystyle \sum _{n=0}^{N-k-1}}r\left(k\right)r_1\left(k\right),k=0,\dots ,M-1,\end{array}$$(2)

and can be equivalently represented as

$$\begin{array}{c}\hfill {\widehat{r}}_{\mathit{rr}}\left(k\right)=\frac{1}{N-k}{r}^{\text{T}}\left(k\right)r_1\left(k\right),k=0,\dots ,M-1,\end{array}$$(3)

where N is the number of residuals, r(k) = [r(0), ..., r(N−k−1)]^T, r₁(k) = [r(k), ..., r(N−1)]^T, and M represents the order of the filter. In the formula, ${\widehat{\mathit{r}}}_{\mathit{rs}}$ is the cross-correlation matrix between r and s, which can be equivalently represented as the auto-correlation matrix ${\widehat{\mathit{r}}}_{\mathit{ss}}$ of s, that is,

$$\begin{array}{c}\hfill {\widehat{r}}_{\mathit{ss}}\left(k\right)={\widehat{r}}_{\mathit{rr}}\left(k\right)-r_{\mathit{vv}}\left(k\right),\end{array}$$(4)

where, r_vv(k) = σ_v²δ(k), k = 0, ..., M − 1, δ(k) is the Kronecker function. When constructing the time-domain Wiener filter, two issues need to be considered. The first is that filtering at the tap locations is not possible. We can address this by applying a second filter using the method of mirror symmetry. The formula is r_i = r_n − i + 1, where r_i is the i-th timing residuals, i ranges from 1 to n, and n is its length. The second issue is that the variance of white noise needs to be known. We construct the covariance matrix of the timing residuals and apply Cholesky decomposition to whiten them. The corresponding formula is

$$\begin{array}{c}\hfill \mathit{C}=\mathit{L}{\mathit{L}}^{\text{T}},\end{array}$$(5)

where C is the covariance matrix of the timing residuals and L is the lower triangular matrix. The whitened timing residuals can be expressed as

$$\begin{array}{c}\hfill {\mathit{R}}_w=\mathit{r}{\mathit{L}}^{-1},\end{array}$$(6)

where R_w is the whitened timing residuals, while r refers to the initial timing residuals. By using this approach, the covariance matrix of the whitened data becomes an identity matrix, thereby removing the correlation among data points.

4. Experimental results

4.1. Timing red noise analysis

Bayesian evidence calculations were carried out with the Enterprise and Dynesty nested sampling Python package. This approach was applied separately to red-noise models and models without red noise for comparison. The search range for A_RN was set between 10⁻¹⁴ and 10⁻¹¹, while γ varied from 0 to 7. The results are summarized in Table 3. Subsequently, we conducted 10⁵ Monte Carlo simulations using PTMCMCSAMPLER to obtain the posterior distributions of the red noise parameters, as illustrated in Figure 4.

Fig. 4. Results of two-dimensional search for Bayesian factors for PSR B1937+21 and PSR B1821-24.

We find strong evidence for the presence of red noise in both PSR B1937+21 and PSR B1821-24 using 7.3 years of NICER data, with Bayesian factors of 19.182 ± 1.213 and 257.617 ± 1.188, respectively. For PSR B1937+21, the median value of log₁₀A is −12.439 with an uncertainty range of [ − 0.082, +0.081], while the spectral index γ has a median of 0.884 with an uncertainty range of [ − 0.433, +0.508]. Similarly, for PSR B1821-24, the median log₁₀A is −11.527 with an uncertainty range of [ − 0.055, +0.059], and the median γ is 2.121 with an uncertainty range of [ − 0.301, +0.333]. After obtaining evidence of the presence of red noise, extracting it from the timing data becomes a critical step in generating ensemble residuals. To achieve this, we employed the previously introduced harmonic fitting method for red noise extraction, as shown in Figure 5.

Fig. 5. Red noise extracted from PSR B1937+21 and PSR B1821-24 using harmonic fitting method.

The harmonic fitting method used to extract red noise contains some high-frequency components, which leads to less smooth results. To address this, we applied a low-pass filter to remove signals with frequencies greater than 0.5 years, which resulted in a smoother reconstruction that successfully captured the fluctuations of the timing residuals across the entire observation period. However, slight discrepancies remained at both ends of the time series. These discrepancies are due to edge effects caused by the finite duration of the observation period, which limits the accuracy of the red noise reconstruction at the start and end of the data. The red noise amplitude for PSR B1937 is ±3 μs, while for PSR B1821 it is ±20 μs, differing by an order of magnitude.

We acknowledge that the red noise estimation for B1937, as illustrated in Figure 5, remains marginal relative to the dominant white noise component. Our analysis benefits from 3.8 additional years of data. Specifically, our dataset spans 7.3 years, in contrast to the 3.5 years used in Hazboun et al. (2022). Nevertheless, we recognize that the conclusions presented by Hazboun et al. (2022) remain fundamentally valid. Their results, which are derived from simulations informed by real observational data, constitute a credible and important benchmark. Although our analysis provides a reasonable estimate of the red noise in B1937, further observational data is required to more precisely characterize and confirm its presence.

After removing the extracted red noise, we fitted the remaining low-frequency signals containing clock errors and calculated the white noise amplitudes, with standard deviations of σ_B1937 = 4.17 μs and σ_B1821 = 25.07 μs, to simulate the generation of the residuals.

4.2. Single pulsar timing stability analysis

Presently, two notable issues plague the calculation methods for pulsar stability. Firstly, the reliability of stability calculations using the σ_z method diminishes at the last data point, which corresponds to the longest observation time. This issue arises because only one set of determined timing residuals is utilized at this juncture, thereby introducing uncertainty. Secondly, the calculation of stability error bars is fundamentally limited, as it relies solely on the distribution of observed timing residuals, which lacks theoretical and practical significance. Consequently, incorporating the uncertainty of timing residuals into simulations and repeating the calculations enabled us to derive statistically validated stability metrics and error bars. This approach emphasizes that more meaningful stability assessments can be derived from larger sample statistics. Figure 6 depicts the stability results derived from 500 Monte Carlo simulations using 7.3 years of NICER data from two pulsars. In this context, σ_z stability is derived from the average of 500 Monte Carlo experiments, with the error bars obtained from the standard deviation. The stability values at different timescales are summarized in Tables 4 and 5.

Fig. 6. Pulsar stability of PSR B1937+21 and PSR B1821-24. (a) Stability analysis using NICER and NANOGrav data for PSR B1937+21. (b) Stability analysis of PSR B1821-24 using NICER data and PPTA data.

Pulsar stability serves as a metric for gauging the variance in a pulsar’s rotational frequency fluctuations. Ideally, the stability derived from observations in both the X-ray and radio bands should be consistent when disregarding any errors. However, constrained by factors such as limited detector size and X-ray flux, the measurement uncertainties in the X-ray band notably exceed those in the radio band. Consequently, in the short term, X-ray stability is inferior to that derived from radio observations. Nonetheless, as observation time increases, the impact of observational uncertainties on stability gradually diminishes, with red noise emerging as the primary factor that influences stability calculations. Consequently, over extended observation periods, the stability derived from X-ray and radio observations closely aligns, with the long-term stability in X-rays potentially even surpassing that of radio data, owing to the sensitivity of the σ_z method for timing red noise. From this perspective, utilizing X-ray timing observations for stability calculations holds significant implications for the long-term assessment of stability. In 2019, Deneva analyzed the one-year NICER data stability of PSR B1937+21 and PSR B1821-24. Due to the limited amount of observational data, the author fitted the stability curve and obtained fitted results of one-year stability for these two pulsars, with values of approximately 1.5 × 10⁻¹³ and 6.5 × 10⁻¹³, respectively, which are close to the results presented in this study (Deneva et al. 2019).

4.3. Ensemble pulsar timing stability analysis

When observational data contain gaps, cubic spline interpolation shows significant fluctuations and distortions, and thereby fails to accurately reflect the true trends and variations. To address these issues, we proposed a new framework for generating ensemble residuals that reliably produces stable timing residuals. We compared the results of cubic spline interpolation with those obtained from the proposed framework for 7.3 years of NICER data for PSR B1937+21, as depicted in Figure 7. We observe that the simulated residuals generated by our framework exhibit a distribution range consistent with the observed residuals (±10 μs), whereas the traditional method (the cubic spline interpolation method) introduces significant errors in observation gaps, with residuals reaching up to 150 μs. Additionally, we proposed an enhanced Wiener filtering algorithm tailored for X-ray data, and the resulting timing residuals are shown in Figure 8. The standard deviation values of the two pulsars decreased significantly after applying the enhanced Wiener filtering, with PSR B1937 reducing from 3.68 μs to 2.31 μs and PSR B1821 reducing from 21.61 μs to 17.00 μs.

Fig. 7. Comparison of simulated timing residuals from the traditional cubic spline method and the proposed framework using 7.3 years NICER observations of PSR B1937+21.

Fig. 8. Timing residuals of PSR B1937+21 before (red) and after (blue) applying the enhanced Wiener filtering algorithm.

We applied the proposed framework to conduct 500 Monte Carlo simulations, and generated simulated timing residuals based on 7.3 years of NICER data for the two pulsars. We then used the long-term timing residuals of the two pulsars to compare the ensemble pulsar stability obtained through the classical weighted method, Wiener filtering method, wavelet decomposition method, and enhanced Wiener filtering method, as illustrated in Figure 9. The stability of the four ensemble pulsar timing algorithms over various observation periods is presented in Table 6. The enhanced Wiener filtering algorithm exhibits its effectiveness in mitigating white noise from X-ray pulsar timing data, thereby elevating stability. Over a 7.3-year timescale, the stability achieved with the enhanced Wiener filtering method exceeds that of the conventional methods by an order of magnitude. These results suggest that the method used significantly enhances the timing stability of X-ray pulsars.

Fig. 9. Ensemble pulsar stability calculated by using the four algorithms.

5. Conclusions and discussion

In this study, we establish X-ray pulsar timing scales using 7.3 years of NICER observation data from PSR B1937+21 and PSR B1821-24, with stability values of 5.21 × 10⁻¹⁴ and 1.14 × 10⁻¹², respectively. By proposing a framework to generate simulated timing residuals, we provided a robust and precise evaluation of the long-term stability for X-ray ensemble pulsar timing. Addressing the prominent presence of white noise in X-ray pulsar timing residuals, we introduced a mirror-symmetric Wiener filtering method based on the Levinson-Durbin iteration technique. This approach comprehensively surpasses all other methods, and achieves the best stability over a 7.3-year period, with a precision of 4.19 × 10⁻¹⁵. Indeed, stability estimators already account for white noise by incorporating more samples over longer lags, thus averaging down the white noise. Using the enhanced Wiener method to reduce white noise may result in a “double averaging” issue. This practice might result in stability comparisons that lack practical significance when contrasted with other comprehensive methods that do not employ such averaging. However, in practical applications, particularly when using pulsar timing for atomic clock calibration, the optimal time-domain Wiener filter can play a significant role. It can enhance calibration accuracy by effectively reducing noise components that are not related to the true signal, thus providing a cleaner and more precise signal for calibration purposes.

The X-ray data exhibit closer consistency with intrinsic noise because they are immune to the propagation effects of the interstellar medium that are contingent on radio-frequency dependence. Hazboun et al. detected red noise in PSR B1821+24 using 3.5 years of NICER observation data. They also predicted that with an observation time of five years and a TOA error of 5 μs, the red noise in pulsar timing can be observed from PSR B1937+21 (Hazboun et al. 2022). In this study, we provided strong evidence for detecting red noise in both PSR B1937+21 and PSR B1821+24 using 7.3 years of NICER observation data. We adopt the approach of fitting harmonic components and successfully reconstruct the red noise time-domain signals of the two pulsars.

Due to certain constraints, such as low radiation flux and limited detection area, the precision of X-ray TOAs is typically one to two orders of magnitude lower than that of radio TOAs, which results in a distinct manifestation of white noise within the X-ray timing residuals. Nevertheless, by utilizing the extraction technique of red noise, it is feasible to extract red noise out of the dominant white noise component. The magnitude of the extracted red noise is consistent with what existed in radio observations (since the X-ray and radio data do not correspond in time, we can only compare the magnitude of the red noise components). Furthermore, we believe that the red noise extracted from long-term X-ray timing residuals is unaffected by interstellar medium noise (Hazboun et al. 2022). When conducting X-ray and radio observations within the same time period, we can compare the trends of red noise in both X-ray and radio spectral bands to confirm whether there is any remaining low-frequency noise in the radio observations (caused by atmospheric dispersion and interstellar medium). This cross-comparison facilitates the understanding and potential correction of any such type of additional noise sources in the radio observations.

It is evident that the long-term stability of pulsar timing may be comparable to or even better than atomic time, whereas the short-term stability of ensemble pulsar timing is significantly lower due to the larger measurement errors in pulsar TOAs. Based on the results obtained in this study for these two pulsars, the one-month stability of ensemble pulsar timing approximates 5.5 × 10⁻¹² order of magnitude. In contrast, some commercial Cesium (Cs) clocks can achieve a one-month stability of 1 × 10⁻¹⁴ order of magnitude, and high-performance optical clocks can even reach a one-month stability of 1 × 10⁻¹⁶ order of magnitude (Hartnett & Luiten 2010). Integrating pulsar timing with existing high-precision atomic time poses tremendous challenges. However, it has been demonstrated that the observations from multiple pulsars can not only be used to detect significant errors in TT (BIPM17) (Hobbs et al. 2012), but to monitor and calibrate atomic clock jumps (Hartnett & Luiten 2010). Although pulsar-based timescales are unlikely to make significant contributions to the stability of atomic timescales in the coming ten years, they will be serving as a valuable independent means for verifying atomic timescales (Hobbs et al. 2012).

Identifying pulsars that exhibit high stability in both the radio and X-ray bands is highly valuable and will serve as one of our future research directions. Future space X-ray observation projects are expected to see improvements in both detection sensitivity and observational coverage, which will contribute to obtaining higher-quality datasets. Projects such as eXTP (Zhang et al. 2016) and STROBE-X (Wilson-Hodge et al. 2017), with their advanced detection capabilities and high temporal and spectral resolution, are expected to greatly enhance X-ray pulsar timing precision. With the deployment of these high-performance X-ray detectors, the precision of X-ray pulse TOA measurements will continue to improve. By designing more precise pulsar timing filtering algorithms, there is huge potential to further advance the applications of pulsar timing in fields such as atomic time error correction and solar system ephemeris refinement.

Acknowledgments

This work was supported by the National Natural Science Foundation of China No. 62103313. This research utilized the data obtained from the NICER Data Archive and software resources made available by the High Energy Astrophysics Science Archive Research Center (HEASARC).

References

All Tables

All Figures

	Fig. 1. Timing residuals after fitting the timing model: (a) PSR B1937+21, (b) PSR B1821-24.
In the text

	Fig. 2. Radio timing residuals for PSR B1937+21 and PSR B1821-24. (a) PSR B1937+21, (b) PSR B1821-24.
In the text

	Fig. 3. Framework for generating ensemble pulsar timing residuals.
In the text

	Fig. 4. Results of two-dimensional search for Bayesian factors for PSR B1937+21 and PSR B1821-24.
In the text

	Fig. 5. Red noise extracted from PSR B1937+21 and PSR B1821-24 using harmonic fitting method.
In the text

	Fig. 6. Pulsar stability of PSR B1937+21 and PSR B1821-24. (a) Stability analysis using NICER and NANOGrav data for PSR B1937+21. (b) Stability analysis of PSR B1821-24 using NICER data and PPTA data.
In the text

	Fig. 7. Comparison of simulated timing residuals from the traditional cubic spline method and the proposed framework using 7.3 years NICER observations of PSR B1937+21.
In the text

	Fig. 8. Timing residuals of PSR B1937+21 before (red) and after (blue) applying the enhanced Wiener filtering algorithm.
In the text

	Fig. 9. Ensemble pulsar stability calculated by using the four algorithms.
In the text