Clock Error Impacts on X-Ray Pulsar Navigation Measurements

Kyle J. Houser; Demoz Gebre-Egziabher

doi:10.33012/navi.773

Abstract

X-ray pulsar navigation (XNAV) has the potential to provide autonomous navigation in deep space, but practical implementations are constrained, in part, by long signal integration times and single-pulsar observations. Under these conditions, navigation states (position, velocity, and clock error) are only partially observable, allowing the clock error to grow unchecked and potentially cause positioning, navigation, and timing algorithms to diverge. This work quantifies the impact of clock error on XNAV position and velocity estimates and identifies the oscillator performance required as a function of signal observation (integration) time. A maximum-likelihood estimator (MLE) is used to estimate phase and frequency from simulated pulsar observations corrupted by timing noise from oscillators of varying quality. The results demonstrate that typical temperature-controlled crystal oscillators are unsuitable for XNAV, whereas high-quality oven-controlled crystal oscillators or chip-scale atomic clocks provide sufficient stability to estimate and mitigate the clock error before clock noise dominates the MLE estimates.

Keywords

1 INTRODUCTION

With humanity preparing to return to the Moon and to send more spacecraft into deep space, navigation resource constraints are beginning to become a major factor in the planning of ongoing and future missions. However, the availability of the de facto system for deep space navigation, the Deep Space Network (DSN), may not be able to support this demand. Currently, the DSN is at times up to 40% over-subscribed, which imposes many operational constraints on missions trying to utilize this unique service (NASA Office of Inspector General, 2023). While upgrades and expansions to the DSN are currently being evaluated, another avenue of research being pursued in parallel aims to provide alternate navigation solutions to relieve strain on the DSN without having to physically modify it. In this regard, autonomous deep space navigation systems are being considered to fill gaps between DSN observations and to provide positioning, navigation, and timing (PNT) solutions.

This situation has motivated the consideration of signal-of-opportunity (SOP) systems for deep space PNT. As is the case with most SOP systems, the flexibility afforded by using signals that are already available comes at the cost of performance and capabilities (e.g., reduced accuracy). However, there are many real operational scenarios in which the reduced performance is an acceptable compromise. For example, one can generally accept this degraded performance in the cruise phases of missions through the deepest regions of space. In these regions, measurements have been historically unavailable, and thus, trajectory corrections are made infrequently. The ability to detect large trajectory errors early, an advantage that comes with an autonomous PNT system, allows the spacecraft to correct those errors earlier, which allows for lower fuel requirements.

X-ray navigation (XNAV) is an example of an autonomous PNT concept that has been proposed for deep space operations (Chen et al., 2020; Runnels & Gebre-Egziabher, 2017; Sheikh et al., 2006). In the XNAV concept, variable X-ray sources, specifically X-ray millisecond pulsars (MSPs), are used as SOPs. Pulsar signals are particularly well suited for PNT because they are periodic and unique from pulsar to pulsar, which makes them easily differentiable. MSPs have the added benefit of being particularly stable from a timing perspective. They have been extensively studied for scientific applications, such as the detection of gravitational waves (Agazie et al., 2023), resulting in the availability of detailed knowledge regarding their signal characteristics.

The phase and Doppler measurements of MSP signals can be exploited to form a position and velocity solution. By measuring the phase difference of a pulsar signal between a spacecraft and some reference point (such as the solar system barycenter [SSB] or the center of Earth), the range along the pulsar's line of sight between the spacecraft and reference point can be determined. This approach is similar to that of a differential global navigation satellite system (GNSS), where phase difference measurements provide observability of the baseline length between a user and a reference station. Unlike differential GNSS, however, the opportunistic use of these signals comes with many challenges. Most notably, owing to the weak nature of X-ray pulsar signals, special considerations are needed to be able to use these signals. One main challenge is the long integration times (or observation windows) required to form accurate measurements for PNT use.

Without contact with Earth or other external communications to provide a stable timing reference, these long integration times can become corrupted with measurement errors that arise from the onboard clocks. In principle, these errors may be corrected by using the XNAV system itself if one has simultaneous observations from four or more pulsars. For reasons that will be described later, this is not always possible. The impact of these clock errors on the XNAV phase and Doppler estimates has not been studied in depth in the XNAV literature. These timing errors vary drastically depending on the quality of the clock or oscillator used in the navigation system. For example, a standard crystal oscillator would add more timing error than a cesium clock. Understanding the impact of this difference in clock quality on forming phase and Doppler estimates is the focus of this paper.

1.1 Prior Work

Prior work has shown that maximum-likelihood estimation (MLE) paired with a digital phase-locked loop (DPLL) allows a dynamic observer to estimate and track the phase. It has been shown that for long integration times, this approach gives the theoretical Cramer–Rao performance bounds for phase and frequency estimates for a pulsar (Golshan & Sheikh, 2007). Additionally, it has been shown that the accumulated phase from a pulsar can be used to recursively determine the distance traveled by a spacecraft relative to its starting point to an accuracy of 50 km (3σ) (Runnels & Gebre-Egziabher, 2017). Chen et al. (2020) showed that by combining XNAV and radio pulsar data with trajectory models, a positioning solution with an accuracy of 8 km can be achieved during non-thrusting segments of the trajectory in the presence of clock errors similar to that of Global Positioning System (GPS) composite clocks. Brito et al. (2015) showed that accuracy on the order of 5 km can be achieved by using the radio signals from a pulsar after an hour of integration time. Even a small X-ray detector, compact enough to fit on a CubeSat-class vehicle, can provide usable navigation measurements. Runnels and Gebre-Egziabher (2021) showed that using a small detector to obtain time and angle-of-arrival measurements from X-ray pulsars and other bright X-ray sources can yield a six degree-of-freedom solution with position errors on the order of 1,000 km (which may be too large for many practical mission applications) and an attitude accuracy of approximately 5 arcseconds. Anderson et al. (2022) presented pulsar signal tracking with a DPLL and extended Kalman filter (EKF) to achieve a position error of 2–6 km depending on the pulsars being tracked. The first real-time on-orbit demonstration of XNAV was performed in 2017 as part of the Station Explorer for X-Ray Timing and Navigation (SEXTANT) mission integrated with the Neutron Star Interior Composition Explorer (NICER) observatory aboard the International Space Station, which demonstrated XNAV performance with an error of less than 10 km (Mitchell et al., 2018).

1.2 Problem Statement

The work in this paper aims to elucidate the impact of clock error on the estimates of phase and Doppler from X-ray pulsar measurements. XNAV measurements typically require long integration times over which clock noise can begin to have an increasingly large impact. In many prior XNAV simulations or studies that used real experimental data in post-processing, clock errors were either neglected or modeled as simple bias. For example, time synchronization by the DSN was used by Runnels and Gebre-Egziabher (2021) in the case of Chandra observatory data. GPS clock synchronization was used in the case of the NICER/SEXTANT work reported by Prigozhin et al. (2016).

Chen et al. (2020) explicitly addressed the estimation of clock errors in an XNAV system. However, the assumed clock errors were on the order of that of the composite clocks used in GPS satellites, which may be of unrealistically high accuracy for many satellite missions. Additionally, Chen et al. (2020) and many other studies assumed that all pulsars are simultaneously viewed as shown in Figure 1(a) and that, consequently, the measurements at any given epoch would be full rank (i.e., all position, velocity, and timing states would be estimated at the same time). In reality, the geometric diversity of pulsar locations in the sky and design constraints for the XNAV sensor would most likely require slewing to individual pulsars to make sequential observations as shown in Figure 1(b). This scenario leads to a diminished observability, meaning that long periods of time would pass before clock errors could be estimated and therefore corrected in the photon time-of-arrival measurements made by the XNAV system.

FIGURE 1

Concept of operations for XNAV

Given that XNAV systems may require large detector areas and/or long integration times, the concept of simultaneous ranging (a) may not always be feasible. Sequential ranging, as illustrated in (b), is more practical but introduces algorithmic challenges that are partially explored in this paper.

The aforementioned issues motivate the main question behind the current work: How long can a pulsar be viewed without estimating or correcting clock errors before timing noise dominates the measurement? To this end, the remainder of this paper is organized as follows: First, a description of the X-ray pulsar measurement and estimator is given. Following this, the clock error model is presented for generating timing noise in the pulsar measurements. Finally, a Monte Carlo simulation is performed to investigate the effects of clock errors on the navigation estimates, and its findings are discussed.

2 MODELS AND MEASUREMENTS FROM X-RAY PULSARS

The X-ray pulsar signals used in XNAV originate from extremely distant sources and are therefore incredibly weak when they are received in the solar system, often having count rates of less than one photon per second. Therefore, what the X-ray detector actually "sees" are discrete photon events instead of a continuous signal like one might imagine for the radio signals used in other PNT systems. This scenario is illustrated in Figure 2. The continuous signal is shown in blue, and the individual photons detected by the spacecraft as discrete energy pulses are shown in orange. If the signal from a particular pulsar were observed over a long period of time (i.e., with a long integration time) and a histogram of the number of received photons were formed as a function of phase, the result would approach the continuous blue line shown in the figure.

FIGURE 2

An example of how one period of the Crab pulsar would appear if viewed as a continuous signal shown as a continuous line

The units on the x-axis are in phase. The y-axis is normalized by the maximum value of the intensity. The individual marks symbolize what a spacecraft would observe in the solar system, i.e., discrete photons. Each mark denotes a single photon.

The actual measurements used by XNAV sensors do not correspond to energy pulses of the signal but rather the time at which the photon resulting in the energy pulse was received. From this, it is apparent why accurate timing is important when performing this measurement; clock errors directly corrupt the time-of-arrival measurements. Because the purpose of this section is to describe how time-of-arrival measurements are used to extract parameters relevant to PNT, for the moment, it will be assumed that a perfect clock is being used to sample the pulsar signal.

2.1 Pulsar Phase Modeling

The parameters estimated by an XNAV algorithm are the phase and frequency of the pulsar signal as observed by a spacecraft. These parameters can be translated into relevant PNT parameters such as range and velocity by assessing how the spacecraft dynamics influence the signal. To do this, a phase evolution function, which describes how the phase of a pulsar's signal evolves over time, is used. This phase is typically modeled using a polynomial function with the following form:

$ϕ (t) = ϕ_{0} + f (t - t_{0}) + \dots$ 1

where $ϕ (t)$ is the phase of the signal at time $t, ϕ_{0}$ is the initial phase offset at some initial reference time, t₀, and f is the rotation frequency of the pulsar. Higher-order terms, such as the spin-down frequency (frequency derivative), can be included in this function; however, for the purposes of this work, it will be sufficient to restrict the model to only the phase and frequency. It is important to note that this function can be used to model the phase observed at any point in space. Here, a reference point, such as the center of the Earth or the SSB, is used as the origin for the XNAV reference frame. Consequently, known values for $ϕ_{0}$ and f are available at some reference epoch t₀ for a stationary observer at the SSB.

The goal is to determine the difference in phase and frequency between the signal observed at the SSB and a spacecraft in motion some distance away from the SSB. The measured initial phase offset of the spacecraft relative to the SSB, $ϕ_{S C}$ , depends on the spacecraft's range from the SSB. The frequency of the signal observed at the spacecraft, $f_{S C}$ , will be different from the frequency of the signal observed at the SSB, $f_{S S B}$ , owing to the Doppler effect caused by the spacecraft's range rate relative to the SSB. The phase and frequency of the signal observed at the spacecraft can be related to what is observed at the SSB by the following relationships:

$ϕ_{S C} = ϕ_{S S B} + \frac{r f_{S S B}}{c}$ 2

$f_{S C} = \sqrt{\frac{1 + \frac{\dot{r}}{c}}{1 - \frac{\dot{r}}{c}}} f_{S S B} \approx (1 + \frac{\dot{r}}{c}) f_{S S B}$ 3

where c is the speed of light and r and $\dot{r}$ are the range and range rate along the line of sight of the pulsar for the spacecraft, respectively. Note that in Equation (3), it is assumed that the spacecraft is traveling much slower than the speed of light. These equations can be inverted so that the range and range rate are written as functions of phase and frequency as follows:

$r = (ϕ_{S C} - ϕ_{S S B}) \frac{c}{f_{S S B}}$ 4

$\dot{r} = (\frac{f_{S C}}{f_{S S B}} - 1) c$ 5

Thus, at a given time t, the basic signal model for a particular pulsar (i.e., Equation (1)) provides the phase and frequency at the SSB, denoted as $ϕ_{S S B}$ and $f_{S S B}$ . At the spacecraft, an XNAV sensor measures the arrival time of photons, from which the corresponding phase and frequency at the spacecraft, $ϕ_{S C}$ and $f_{S C}$ , can be estimated using an algorithm described later in this paper. A key element of that algorithm is a mathematical model for the photon arrival time. These models, referred to here as photon count rate models, are discussed next.

2.2 Photon Count Rate Models

As mentioned above, the actual observation of the X-ray pulsar signal made by an XNAV sensor is not the amplitude of some continuous signal but rather the time of arrival of discrete X-ray photons. The number of photons observed over a given time interval can be described by a non-homogeneous Poisson process (NHPP). Rather than a standard Poisson process with a constant rate, the pulsar signal has a variable rate as it "pulses" over the duration of one signal period. Mathematically, the probability that N discrete photons arrive over some time interval corresponds to a Poisson process whose probability mass function is given by the following:

$\Pr [{\tilde{Λ}}_{k} = N ∣ {\bar{λ}}_{k}] = {\begin{array}{cc} e^{- {\bar{λ}}_{k}} \frac{{({\bar{λ}}_{k})}^{N}}{N!} & {N \in Z ∣ N \geq 0} \\ 0 & N < 0 \end{array}$ 6

where:

${\tilde{Λ}}_{k}$ is the number of events observed in a time window
${\bar{λ}}_{k}$ is the number of events expected to occur over the time window, calculated as follows:

${\bar{λ}}_{k} = \int_{t_{0}}^{t_{0} + T_{obs}} λ (t) d t$ 7

The time dependence on the rate function λ makes this process an NHPP rather than a homogeneous Poisson process. This rate function is defined using pulsar-specific parameters. These parameters characterize the pulsar's signal and noise count rates and are defined as follows:

$λ (t) = β + α h (ϕ (t))$ 8

where α and β are the average photon count rate for the pulsed signal and background, respectively, and $ϕ (t)$ is the phase evolution function. Background and source counts are dependent on many factors, such as the effective detector area, energy sensitivity, field of view, etc., and must be determined for a given detector design. The ratio $α / β$ serves as a proxy for the pulsar signal-to-noise ratio. This ratio is considered a proxy because β does not include sensor-specific noise sources or inefficiencies.

The function $h (ϕ)$ is called the pulse profile or light curve, an example of which is shown in Figure 2 as a solid line. This function describes the shape and relative intensity of the pulses from the pulsar signal. The pulse profile is defined such that it satisfies the following:

$\min_{ϕ \in (0, 1)} h (ϕ) = 0, \int_{0}^{1} h (ϕ) d ϕ = 1$

This rate function combined with the phase evolution function (Equation (1)) provides the information needed to accurately model the photon count statistics when viewing X-ray pulsars. Realizations of these processes can be generated in different ways described in standard simulation and modeling texts such as the work by Ross (1990). For example, thinning (accept/reject) methods or an inversion method can be used, as presented by Yu et al. (2020).

2.3 Maximum Likelihood Estimation

In this work, an MLE similar to the one described by Golshan and Sheikh (2007) is used to estimate the phase and frequency from the received pulsar signal. For completeness, some of the key points of MLE will be discussed below, but interested readers are encouraged to read the above work for the full derivation. The derivation of the MLE starts from the following fact: Whereas the number of photons arriving at the XNAV sensor in a given time window is described by an NHPP, the times at which those photons arrive are exponentially distributed. This distinction is important because the observation made by the XNAV sensor corresponds to photon arrival times $t_{k}$ for $k = 1, 2, \dots M$ , where M is some large number, rather than the number of photons N that arrive in a given interval. This exponential process is also non-homogeneous because of the variable rate function. The probability density function of the non-homogeneous exponential process (NHEP) is defined as follows:

$p (t_{k}; θ) = \exp [- \int_{t_{0}}^{t_{k}} λ (t; θ) d t] λ (t_{k}; θ)$ 9

where:

$t_{k}$ is the photon arrival time
$t_{0}$ is the starting observation time
θ is the vector of parameters to be estimated (phase and frequency)

Because λ depends on the phase ϕ, the above expression gives the probability that a photon was observed at time $t_{k}$ given the parameters that need to be estimated (phase and frequency). Note that it is assumed that each photon arrival time is an independent observation of this exponential process. Thus, the joint likelihood function for a set of M photon observations (denoted ${t_{k}}$ for $k \in [1 M]$ ) conditioned on θ is given by the following:

$p ({t_{k}}; θ) = \exp [- \int_{t_{0}}^{t_{0} + T_{obs}} λ (t; θ) d t] \prod_{k = 1}^{M} λ (t_{k}; θ)$ 10

where $T_{obs}$ is the total time of the observation. Taking the natural logarithm of this yields the following log-likelihood function (LLF):

$LLF ({t_{k}}; θ) = \sum_{k = 1}^{M} \log [λ (t_{k}; θ)] - \int_{t_{0}}^{t_{0} + T_{obs}} λ (t; θ) d t$ 11

This LLF is the cost function that the MLE aims to maximize. Golshan and Sheikh (2007) noted that the integral term exhibits little dependence on the parameters in θ and, thus, can be dropped from the LLF. Substituting the rate function along with the phase evolution function from above yields the following:

$LLF ({t_{k}}; ϕ, f) = \sum_{k = 1}^{M} \log [β + α h (ϕ + f (t_{k} - t_{0}))]$ 12

The MLE simply determines the arguments that maximize this function:

$(\hat{ϕ}, \hat{f}) = \underset{ϕ \in Φ, f \in Ω}{\arg max} \sum_{k = 1}^{M} \log [λ (t_{k}; ϕ, f)]$ 13

Equation (13) is a constrained optimization problem that can be solved in many ways. In this work, the above equation is solved numerically by minimizing the negative of the LLF and using a global solver method such as particleSwarm in MATLAB. Other approaches exist for solving this MLE, which may be more efficient to implement on real hardware (i.e., not a laptop or personal computer). For the purposes of this paper, however, using a numerical method as described above is deemed acceptable.

An important factor for numerically generating these estimates is bounding the search space for these parameters. Because te pulsar signal is periodic, a search space larger than one full cycle of the signal would yield multiple estimates that would be indistinguishable. This situation is similar to the integer ambiguity problem in carrier-phase differential GNSS. To this end, the search space in this paper is restricted to one full signal period because resolving these integer ambiguities is deemed to be outside the scope of this work.

There is more freedom in defining the frequency search space. The frequency will be shifted by Doppler effects as shown in Equation (3) and, therefore, is affected by velocity. Thus, a search space centered about the barycentric frequency, $f_{S S B}$ , can be used. The upper and lower bounds for this frequency search space can be set according to the minimum and maximum expected range rate for a given spacecraft.

2.4 Cramer-Rao Lower Bound Analysis

For the work in this paper, the Cramer-Rao lower bound (CRLB) is used as a lower bound on the variance of estimates generated by the MLE. It is well known that in the absence of timing errors, the MLE will approach the CRLB as the observation time increases. Accordingly, this paper proposes using the CRLB as a metric for assessing the impact of clock errors.

As long as the MLE does not diverge, the estimates will approach the CRLB with increasing observation time. However, in the presence of sufficiently large clock errors, the estimates will diverge away from the CRLB. Thus, for a given clock, the allowable integration time will be defined as the point at which clock noise begins to drive the MLE estimates away from the CRLB and causes divergence.

The CRLB is defined as follows (Kay, 1993):

$cov (\hat{θ}) \geq I^{- 1} (θ)$ 14

where I is the Fisher information matrix. For the multivariate case in this paper, each element in the Fisher information matrix is defined as follows:

$[I (θ)]_{i j} = - E [\frac{\partial^{2}}{\partial θ_{i} \partial θ_{j}} LLF (x; θ)]$ 15

The elements of the Fischer information matrix for the MLE in this paper are given by the following Snyder (1991):

$[I (θ)]_{i j} = \int_{t_{0}}^{t_{1}} \frac{\frac{\partial}{\partial θ_{i}} λ (t; θ) \frac{\partial}{\partial θ_{j}} λ (t; θ)}{λ (t; θ)} d t$ 16

A derivation for the CRLB when estimating phase and frequency has been given by Ashby and Golshan (2008). The final result for the CRLB at the end of this derivation is as follows:

$I^{- 1} (θ) = \frac{1}{L} [\begin{array}{cc} 4 T_{obs}^{- 1} & - 6 T_{obs}^{- 2} \\ - 6 T_{obs}^{- 2} & 12 T_{obs}^{- 3} \end{array}]$ 17

where the coefficient L is given by the following:

$L = \int_{t_{0}}^{t_{1}} \frac{{(α h^{'})}^{2}}{β + α h} d ϕ$

Because the integrand in L is always positive, the main driver for decreasing the lower bound of this MLE estimation error is the observation time. As observation time increases, the phase and frequency error variances are expected to decrease.

Once again, it should be noted that this bound corresponds to the error resulting from the stochastic nature of the NHPP/NHEP process. There has been no consideration regarding the impact of clock error on this estimate. Thus, this result corresponds to the optimal bound in the absence of any clock errors. In the following analyses on the impact of clock noise, the integration window length at which the estimator stops approaching and instead diverges away from the CRLB will be taken as the allowable integration window for a given clock. Another way to interpret this is that the point at which the MLE begins to diverge from the CRLB corresponds to the point at which clock errors have an impact greater than the inherent stochasticity of the photon arrival times.

3 CLOCK MODELS

In this work, a two-state, frequency random walk clock error model is used to simulate the timing errors. Specifically, the timing (phase) error, $δ_{t}$ , and frequency error, $δ_{f}$ , are modeled as a discrete random process driven by Gaussian noise:

${[\begin{matrix} δ_{t} \\ δ_{f} \end{matrix}]}_{k + 1} = \underset{A}{\underset{⏟}{[\begin{array}{cc} 1 & Δ t \\ 0 & 1 \end{array}]}} \underset{x_{k}}{\underset{⏟}{{[\begin{matrix} δ_{t} \\ δ_{f} \end{matrix}]}_{k}}} + \underset{w_{k}}{\underset{⏟}{{[\begin{matrix} w_{t} \\ w_{f} \end{matrix}]}_{k}}}$ 18

where $Δ t$ is the time between photon arrival times, $w_{t}$ is the phase noise, and $w_{f}$ is the average frequency noise over the time interval. The noise terms are distributed as $w_{k} \sim N (0, Q_{k})$ . The covariance $Q_{k}$ for this model is constructed by using Allan variance power law coefficients for a given oscillator. Van Dierendonck et al. (1984) derived the covariance matrix (with corrections from Brown and Hwang (2012)) of the process noise for a given timestep, $Q_{k}$ . The clock model here is based on the work of Van Dierendonck et al. (1984), with the only difference being that the flicker frequency noise is ignored, as suggested by Brown and Hwang (2012). The entries of this covariance matrix are given by the following:

$Q_{k} = [\begin{array}{ll} q_{11} & q_{12} \\ q_{21} & q_{22} \end{array}]$ 19

with:

$q_{11} = \frac{h_{0}}{2} Δ t + \frac{2}{3} π^{2} h_{- 2} Δ t^{3}$ 20a

$q_{12} = q_{21} = π^{2} h_{- 2} Δ t^{2}$ 20b

$q_{22} = 2 π^{2} h_{- 2} Δ t$ 20c

where h₀, and $h_{- 2}$ are the power law coefficients for white and random walk frequency noise, respectively.

To simulate noisy time measurements, $x_{0}$ (that is, $x_{k}$ at $k = 0$ in Equation (18)) is initialized to $x_{0} = {[\begin{array}{ll} 0 & 0 \end{array}]}^{T}$ . Then, a set of clean photon arrival time observations is generated and used to determine the times between photon arrivals, ${Δ t_{k}}$ . The value of $Δ t_{k}$ is used in both determining the covariance of the random noise and propagating the clock phase and frequency errors forward. The covariance, P, is also propagated in the same way as the states themselves. Using the same linear system in Equation (18) for this propagation yields the following:

$P_{k + 1} = A P_{k} A^{T} + Q_{k}$ 21

To confirm that this method generates errors with statistics consistent with the process noise matrix, clock error histories are generated for a given oscillator, and the resulting statistics are compared with the propagated covariance. The results of one such validation run are shown in Figure 3. The figure shows results from 1,000 trials, where each of the trials lasted 10⁴ s with a time step of $Δ t = 0.1 s$ . The process noise matrix was generated with the assumption that the clock is an oven-controlled crystal oscillator (OCXO) with the power law coefficients shown in Table 1.

FIGURE 3

Results of the simulation of an OCXO in both the time domain and in the Allan variance. (a) Results of 1,000 runs of timing (phase) and frequency noise generation. All runs are shown in gray, with an individual run highlighted in blue and 3σ bounds for this error shown in red, (b) Allan variance (AVAR) of timing noise (shown in blue). Both the white noise component and random walk noise component have distinct slopes. The power law fit of the Allan variance is shown by the red dashed line.

View this table:

Table 1 List of Oscillators Used in the Monte Carlo Experiment with Relevant Power Law Coefficients Describing Their Stability

For such a simulation, it is expected that $0.27 % (1, 000 \times 0.00027 \approx 3)$ of the runs may exceed the $\pm 3 σ$ bounds. It can be seen that only 2 of the 1,000 runs barely exceed the $\pm 3 σ$ bounds. This result suggests that the variance or standard deviation computed by Equation (21) properly describes the error being generated.

In the XNAV analysis described in the next section, this clock error model is used as follows: Given a time history of photon arrival times (referred to as the “clean” arrival time history in subsequent discussions), clock biases for each time step are determined using the above model and the clock characteristics given in Table 1. A time history of clock errors consistent with the photon arrival times is then generated. The clock error corresponding to each photon arrival time is subsequently added to the clean photon arrival times as follows:

${\tilde{t}}_{k} = t_{k} + δ_{t, k}$ 22

The impact of clock error depends on the interplay of signal flux, the quality of the clock, and the length of time during which the clock is relied upon for timing information. This latter aspect is a function of the $Δ t$ term in the clock error model given by Equation (18). For example, for a high-flux pulsar, where the average $Δ t$ is small, shorter-term frequency stability (i.e., low white frequency noise) may be important. In contrast, for low-flux pulsars, where the average $Δ t$ is high, low short-term instability might be tolerated in exchange for long-term stability (i.e., low random walk frequency noise).

In summary, the process of generating a history of photon arrival time measurements that accounts for clock errors is as follows:

Generate clean photon arrival times ( $t_{1}, t_{2}, \dots t_{N}$ ) consistent with the rate function for a given pulsar (Equation (8)).
Determine the set of times between each photon arrival, ${Δ t_{k} = t_{k} - t_{k - 1}}$ .
Generate multivariate white noise using each $Δ t_{k}$ with Equation (19) given a set of power law coefficients.
Propagate the timing errors according to Equation (18) from $t_{k - 1}$ to $t_{k}$ for $k = 1, 2, \dots M$ . Note that at $k = 1$ (the first time step), $t_{k - 1} = t_{0}$ and the clock error state is ${[\begin{array}{ll} 0 & 0 \end{array}]}^{T}$ .
Add each timing noise term $δ_{t, k}$ to the corresponding event time.

4 MONTE CARLO EXPERIMENT

To determine the effect of clock error on the pulsar measurements, a Monte Carlo simulation study was performed. For each run of the Monte Carlo simulation, a list of photon arrival times uncorrupted by clock error (i.e., a clean history of arrival times) was generated. The clock error model was then used to generate clock noise to be added to the clean photon arrival times in the process described in Section 3. These steps are applied for each clock model listed in Table 1. This process results in a set of five time histories; one clean set and four other sets corrupted by clock errors.

Each arrival time history is then run through the MLE independently and generate estimates for the initial phase and frequency of the pulsar signal observed by a spacecraft. This process is performed 100 times with different random observation noise realizations for a single observation duration (integration time) and clock model, and the variance of the estimates is recorded. Then, the observation duration is increased, and the process is repeated. Observation times ranging from 1,000 to $1 \times 10^{6} s$ are simulated. The history of the variance of each these runs is compared with the CRLB.

4.1 Clock Selection

Before discussing the results of the above-described simulations, a brief discussion of the clocks chosen for the simulation studies is warranted. Clock parameters were chosen to cover various qualities or broad classes of types of oscillators used in PNT. Parameters for a typical temperature-controlled crystal oscillator (TCXO), OCXO, and a cesium clock were taken from the work of Brown and Hwang (2012). Parameters for the Microchip SA65 chip-scale atomic clock (CSAC) were calculated using Allan deviation charts found in data sheets from the manufacturer. A list of these parameters is provided in Table 1.

The TCXO and OCXO were chosen because they are relatively common and inexpensive to integrate into a small satellite system; moreover, if they are stable enough to generate accurate estimates, they may prove beneficial for future XNAV systems. The cesium clock was chosen to represent the best available technology for space-borne timing applications. However, it should be noted that it would be difficult to practically integrate most space-qualified cesium standard clocks into a small CubeSat. The CSAC represents an interesting case. While it may not be able to provide the same stability as a traditional atomic clock, the CSAC provides an improvement in stability compared with the TCXO and OCXO while also being easy to integrate into a spacecraft system. The low size, weight, power, and cost (SWAP-C) of CSACs compared with that of a traditional cesium atomic clock presents an interesting option for integration into XNAV sensors.

4.2 Pulsar Selection

The pulsar used in this simulation study was carefully selected. Using a pulsar like the Crab pulsar would yield good estimates with shorter observation times, but the Crab pulsar is truly an outlier among MSPs. To test a more typical pulsar, PSR B1937+21 was used. The flux rate used for the pulsar in this study came from the data cataloged by the NICER instrument (Prigozhin et al., 2016). For this pulsar, the signal flux, α, is 0.029 photons per second, and the background flux, β, is 0.2 photons per second. The rotational period for this pulsar is approximately 1.56 ms or 641 Hz. These numbers are much more representative of pulsars that would be used for XNAV compared with the Crab pulsar, which, when viewed with the same instrument, has a signal and background flux of $α = 660.0$ and $β = 13860.2$ photons per second, respectively. Two cycles of the light curve, $h (ϕ)$ , for PSR B1937+21 are shown in Figure 4; these cycles were used for the Monte Carlo simulation study described above.

FIGURE 4

Two cycles of the light curve, $h (ϕ (t))$ , for PSR B1937+21

FIGURE 5

Results of the Monte Carlo experiment for all oscillators, showing the MLE standard deviation for each observation duration

The standard deviation for the estimates produced by the MLE using the perfect data is also included. The CRLB is shown as a dashed line.

4.3 Definition of MLE Search Space

Bounding the phase and frequency search space is an important step in solving the MLE optimization problem. Although the upper bound for the phase search space would be one period of the pulsar signal, it is inefficient and unnecessary to search the entire cycle. In contrast, restricting the search space to be small could artificially depress the estimate variance. Thus, the search space was bounded by using an inflated version of the CRLB itself. The CRLB for both phase and frequency was calculated for each observation time in the simulation. The search space for the problem was then defined to be the true value for phase and frequency plus or minus a padding equal to 10 times the CRLB.

5 RESULTS AND DISCUSSION

The results for each test case are presented below. The main result that is being sought here is when (or if) the standard deviation (or variance) of the estimation error for phase and frequency stops decreasing with observation window length and diverges away from the CRLB of the clean or uncorrupted arrival time history. The clean or uncorrupted set of time-of-arrival data shows that the MLE converges and rules out the presence of any effects caused by the signal itself. If the MLE using the corrupted data set diverges when the clean set does not, it can be concluded that timing noise from the clock has most likely caused this divergence. In these cases, sufficient error has accumulated from the clock noise that the estimate is no longer reliable.

The results for the least stable oscillator, the TCXO, are analyzed first. The TCXO's high degree of timing noise can be considered a worst-case scenario for clock errors. At 1,000 s, none of the estimates, not only the TCXO estimates, have converged to the CRLB because the estimator does not yet have enough information to reach this theoretical bound. However, increasing the observation time allows the other data sets to converge to and stay close to the CRLB. In the TCXO corrupted data set, it initially appears as though the errors start converging to the CRLB. However, this trend quickly reverses when the TCXO solution continues to diverge. The solution continues to diverge until it reaches the bounds of the MLE solver and remains at those bounds for the remainder of the observation times.

The OCXO provides an improvement over the TCXO in terms of stability. It can be seen that its performance is much better than that of the TCXO. As before, neither data set converges to the CRLB when 1,000 s of observation data is used. However, unlike the TCXO case, the OCXO data set continues to approach the CRLB with increasing observation time. This data set reaches the CRLB by approximately 2,000 s and remains close to the CRLB until approximately 20,000 s. At this point, the timing noise begins to add error to the estimate that cannot be mitigated by the additional data that come with increasingly longer observation windows. After this point, the standard deviation continues to increase until it too reaches the bounds of the MLE solver.

As would be expected, the CSAC performs better than the OCXO. The CSAC and OCXO initially show the same trend, with the errors converging to the CRLB after approximately 2,000 s. In the case of the CSAC, however, the errors remain converged to the CRLB until approximately 50,000 s, representing an improvement over the OCXO. Interestingly, over shorter durations of 1,000–2,000 s, the CSAC actually shows higher estimation variance before converging to the CRLB. This result is likely due to the higher white frequency noise on the CSAC compared with the OCXO used here.

Finally, the cesium clock shows the same behavior as the clean data (i.e., data in which the noise comes from the inherent stochasticity of photon arrival times). Once the estimate approaches the CRLB, it does not diverge as observed for the other clocks tested. The results are nearly indistinguishable from the clean data. Because the cesium clock is a high-stability oscillator, this result is not surprising, especially as these oscillators typically drift by only 10 ns over one day (Hutsell, 1994).

6 CONCLUSIONS

This paper presented an analysis of the impact of clock error on XNAV measurements and PNT estimates. Simulated XNAV observations consisting of X-ray photon arrival times were corrupted with timing noise of various levels based on common oscillator types. Low-quality oscillators are not suitable for XNAV systems, as they corrupt the photon arrival time data to the point at which that noise dominates the measurements. An OCXO may be good for short observation times in the range of 1,000–20,000 s, but longer observations would require a more stable clock, such as a CSAC. The cesium clock showed exceptional performance. Although cesium clocks are used in space for applications such as GPS satellites, it may not be feasible to integrate a cesium clock into some spacecraft with restricted SWAP-C constraints. However, as the results of this paper show, a CSAC or even a high-quality OCXO could be used for XNAV applications in which SWAP-C constraints are prohibitive.

The analysis above showed that timing errors must be accounted for in XNAV algorithms. Future work can include adding a representative clock model directly into the MLE to potentially mitigate the effect of some clock noise by estimating an average bias or drift in the clock. Other XNAV estimation methods, such as the correlation vector EKF reported by Runnels and Gebre-Egziabher (2017), are also good candidates for a similar analysis, as this filter uses the same photon arrival time measurements as the MLE used in this paper.

HOW TO CITE THIS ARTICLE:

Houser, K.J., & Gebre-Egziabher, D. (2026). Clock error impacts on X-ray pulsar navigation measurements. NAVIGATION, 73. https://doi.org/10.33012/navi.773

ACKNOWLEDGMENTS

The authors gratefully acknowledge support from the National Aeronautics and Space Administration (NASA) University Smallsat Technology Partnerships under cooperative agreement #80NSSC23M0233 and NASA/Minnesota Space Grant Consortium (MnSGC) for the work conducted in this paper. Additionally, the authors would like to acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing computing resources, which contributed to the simulation results shown in this paper. The authors would also like to acknowledge valuable input and advice from Dr. Lindsay Glesener (University of Minnesota - Twin Cities) as well as Dr. Paul Ray for help with processing NICER data. Although the authors gratefully acknowledge the support of the aforementioned organizations and individuals, the views and conclusions expressed in this paper are those of the authors alone and should not be interpret as representing the official policies, either expressed or implied, of the MnSGC, MSI, or NASA.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

↵
Agazie, G., Alam, M. F., Anumarlapudi, A., Archibald, A. M., Arzoumanian, Z., Baker, P. T., Blecha, L., Bonidie, V., Brazier, A., Brook, P. R., Burke-Spolaor, S., Bécsy, B., Chapman, C., Charisi, M., Chatterjee, S., Cohen, T., Cordes, J. M., Cornish, N. J., Crawford, F.,... Young, O. (2023). The NANOGrav 15 yr data set: Observations and timing of 68 millisecond pulsars. The Astrophysical Journal Letters, 951. https://doi.org/10.3847/2041-8213/acda9a
Google Scholar
↵
Anderson, K. D., Pines, D. J., & Sheikh, S. I. (2022). Investigation of X-ray pulsar signal phase tracking for spacecraft navigation. In AIAA Science and Technology Forum and Exposition, AIAA SciTech Forum 2022. https://doi.org/10.2514/6.2022-1589
Google Scholar
↵
Ashby, N., & Golshan, A. R. (2008, January). Minimum uncertainties in position and velocity determination using X-ray photons from millisecond pulsars. In Proceedings of the Institute of Navigation, National Technical Meeting (pp. 110–118). https://www.ion.org/publications/abstract.cfm?articleID=7668
Google Scholar
↵
Brito, D., Tavares, G., Fernandes, J., Noroozi, A., Verhoeven, C., & Lisboa, U. D. (2015, March). Radio pulsar receiver systems for space navigation. In 8th European Symposium on Aerothermodynamics for Space Vehicles.
Google Scholar
↵
Brown, R. G., & Hwang, P. Y. C. (2012). Introduction to random signals and applied Kalman filtering (4th ed.). Wiley.
Google Scholar
↵
Chen, P. T., Zhou, B., Speyer, J. L., Bayard, D. S., Majid, W. A., & Wood, L. J. (2020). Aspects of pulsar navigation for deep space mission applications. Journal of the Astronautical Sciences, 67, 704–739. https://doi.org/10.1007/s40295-019-00209-9
Google Scholar OpenURL
↵
Golshan, A. R., & Sheikh, S. I. (2007, April). On pulse phase estimation and tracking of variable celestial X-ray sources. In Proceedings of the 63rd Annual Meeting of the Institute of Navigation (2007) (pp. 413–422). https://www.ion.org/publications/abstract.cfm?articleID=7225
Google Scholar
↵
Hutsell, S. T. (1994, December). Fine tuning of GPS clock estimation in the MCS. In Proceedings of the 26th Annual Precise Time and Time Interval Systems and Applications Meeting (pp. 63–74).
Google Scholar
↵
Kay, S. M. (1993). Fundamentals of statistical signal processing, Volume I: Estimation theory. Pearson.
Google Scholar
↵
Mitchell, J. W., Winternitz, L. B., Hassouneh, M. A., Price, S. R., Semper, S. R., Yu, W. H., Ray, P. S., Wolff, M. T., Kerr, M., Wood, K. S., Arzoumanian, Z., Gendreau, K. C., Guillemot, L., Cognard, I., & Demorest, P. (2018). SEXTANT X-ray pulsar navigation demonstration: Initial on-orbit results. Advances in the Astronautical Sciences, 164. https://ntrs.nasa.gov/citations/20180001252
Google Scholar
↵
NASA Office of Inspector General. (2023). Audit of NASA’s Deep Space Network (tech. rep.). National Aeronautics and Space Administration.
Google Scholar
↵
Prigozhin, G., Gendreau, K., Doty, J. P., Foster, R., Remillard, R., Malonis, A., LaMarr, B., Vezie, M., Egan, M., Villasenor, J., Arzoumanian, Z., Baumgartner, W., Scholze, F., Laubis, C., Krumrey, M., & Huber, A. (2016). NICER instrument detector subsystem: Description and performance. Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, 9905, 99054X-1–99054X-7. https://doi.org/10.1117/12.2231718
Google Scholar
↵
Ross, S. M. (1990). A course in simulation (1st ed.). Academic Press.
Google Scholar
↵
Runnels, J. T., & Gebre-Egziabher, D. (2017). Recursive range estimation using astrophysical signals of opportunity. Journal of Guidance, Control, and Dynamics, 40, 2201–2213. https://doi.org/10.2514/1.G002650
Google Scholar OpenURL
↵
Runnels, J. T., & Gebre-Egziabher, D. (2021). Estimator for deep-space position and attitude using X-ray pulsars. IEEE Transactions on Aerospace and Electronic Systems, 57, 2149–2166. https://doi.org/10.1109/TAES.2021.3068432
Google Scholar OpenURL
↵
Sheikh, S. I., Pines, D. J., Ray, P. S., Wood, K. S., Lovellette, M. N., & Wolff, M. T. (2006). Spacecraft navigation using X-ray pulsars. Journal of Guidance, Control and Dynamics, 29, 49–63. https://doi.org/10.2514/1.13331
Google Scholar OpenURL CrossRef
↵
Snyder, D. L. (1991). Random point processes in time and space. Springer-Verlag.
Google Scholar
↵
Van Dierendonck, A. J., McGraw, J. B., & Brown, R. G. (1984, November). Relationship between Allan variances and Kalman filter parameters. In Proceedings of the 16th Annual Precise Time and Time Interval Systems and Applications Meeting (pp. 273–293). https://www.ion.org/publications/abstract.cfm?articleID=16168
Google Scholar
↵
Yu, W. H., Semper, S. R., Mitchell, J. W., Winternitz, L. B., Hassouneh, M. A., Price, S. R., Ray, P. S., Wood, K. S., Gendreau, K. C., & Arzoumanian, Z. (2020). NASA SEXTANT mission operations architecture. Acta Astronautica, 176, 531–541. https://doi.org/10.1016/j.actaastro.2020.06.040
Google Scholar OpenURL

PDF

Previous Next

Clock Error Impacts on X-Ray Pulsar Navigation Measurements

Abstract