Measuring clock jumps using pulsar timing

TL;DR

This paper presents a Bayesian statistical method to detect and measure instant clock jumps using pulsar timing data, achieving high precision comparable to satellite-based time transfer services.

Contribution

The paper introduces an analytic marginalization technique for pulsar timing systematics within a Bayesian framework, enabling efficient and precise clock jump detection.

Findings

Achieved 80-ns precision in clock jump measurement

Demonstrated method's applicability with real pulsar data

Potential to aid in developing autonomous pulsar time scales

Abstract

In this paper, we investigate the statistical signal-processing algorithm to measure the instant local clock jump from the timing data of multiple pulsars. Our algorithm is based on the framework of Bayesian statistics. In order to make the Bayesian algorithm applicable with limited computational resources, we dedicated our efforts to the analytic marginalization of irrelevant parameters. We found that the widely used parameter for pulsar timing systematics, the `Efac' parameter, can be analytically marginalized. This reduces the Gaussian likelihood to a function very similar to the Student's $t$-distribution. Our iterative method to solve the maximum likelihood estimator is also explained in the paper. Using pulsar timing data from the Yunnan Kunming 40m radio telescope, we demonstrate the application of the method, where 80-ns level precision for the clock jump can be achieved. Such a precision is comparable to that of current commercial time transferring service using satellites. We expect that the current method could help developing the autonomous pulsar time scale.

Equations

$$R_{i,j}=s_j+n_{i,j},$$

$$f(R_{i,j})\propto \prod _{i=1}^N\prod _{j=1}^{M_i}\frac{1}{{\eta }_i{\sigma }_{i,j}}{e}^{-\frac{1}{2}{\left(\frac{R_{i,j}-s_{i,j}-{\mu }_i}{{\eta }_i{\sigma }_{i,j}}\right)}^2}.$$

$$s_{i,j}=\left\{\begin{array}[]{l}s_{0},\textrm{for }T_{i,j}>t_{0}\,,\\ 0,\textrm{ for alternatives}\,,\end{array}\right.$$

$$f(R_{i,j})\propto \prod _{i=1}^N\prod _{j=1}^{\mathrm{before}}\frac{1}{{\eta }_i{\sigma }_{i,j}}{e}^{-\frac{1}{2}{\left(\frac{R_{i,j}-{\mu }_i}{{\eta }_i{\sigma }_{i,j}}\right)}^2}\prod _{i=1}^N\prod _{j=\mathrm{after}}^{M_i}\frac{1}{{\eta }_i{\sigma }_{i,j}}{e}^{-\frac{1}{2}{\left(\frac{R_{i,j}-s_0-{\mu }_i}{{\eta }_i{\sigma }_{i,j}}\right)}^2}.$$

$$s_0,t_0,{\mu }_i,{\eta }_i=\mathrm{argmax}[f(R_{i,j})].$$

$$f^{\prime }(R_{i,j})\equiv \int f(R_{i,j})p(\eta )\mathrm{d}\eta ,$$

$$f^{\prime }(R_{i,j})\propto \prod _{i=1}^N{\left[\sum _{j=1}^{\mathrm{before}}{\left(\frac{R_{i,j}-{\mu }_i}{{\sigma }_{i,j}}\right)}^2+\sum _{j=\mathrm{after}}^{M_i}{\left(\frac{R_{i,j}-{\mu }_i-s_0}{{\sigma }_{i,j}}\right)}^2\right]}^{-M_i/2}.$$

$$\widehat{{\mu }_i}=\frac{{\sum }_{j=1}^{\mathrm{before}}{R}_{i,j}{\sigma }_{i,j}^{-2}}{{\sum }_{j=1}^{\mathrm{before}}{\sigma }_{i,j}^{-2}}.$$

$$\sum _{i=1}^N\frac{{\sum }_{j=\mathrm{after}}^{M_i}\frac{R_{i,j}-{\mu }_i-\widehat{s_0}}{{\sigma }_{i,j}^2}}{\frac{1}{M_i}\left[{\sum }_{j=1}^{\mathrm{before}}{\left(\frac{R_{i,j}-{\mu }_i}{{\sigma }_{i,j}}\right)}^2+{\sum }_{j=\mathrm{after}}^{M_i}{\left(\frac{R_{i,j}-{\mu }_i-\widehat{s_0}}{{\sigma }_{i,j}}\right)}^2\right]}=0.$$

$${\sigma }_{i,j}^{\prime 2}={\sigma }_{i,j}^2\frac{1}{M_i}\left[\sum _{j=1}^{\mathrm{before}}{\left(\frac{R_{i,j}-{\mu }_i}{{\sigma }_{i,j}}\right)}^2+\sum _{j=\mathrm{after}}^{M_i}{\left(\frac{R_{i,j}-{\mu }_i-\widehat{s_0}}{{\sigma }_{i,j}}\right)}^2\right],$$

$${\widehat{s_0}}^{\prime }=\frac{{\sum }_{i=1}^N{\sum }_{j=\mathrm{after}}^{M_i}\left(R_{i,j}-{\mu }_i\right){\sigma }_{i,j}^{\prime -2}}{{\sum }_{i=1}^N{\sum }_{j=\mathrm{after}}^{M_i}{\sigma }_{i,j}^{\prime -2}}.$$

$$\mathrm{Var}[s_0]\equiv ⟨{\widehat{s_0}}^2⟩-⟨\widehat{s_0}{⟩}^2\simeq A+B,$$

$$A$$

$$B$$

Full text

\PACS

97.60.Gb, 05.45.Tp, 95.55.Sh

1 Introduction

Pulsars are stable celestial clocks, particularly, the long-term stability of millisecond pulsars (MSPs) is comparable to that of the current atomic time-scales [2, 1]. On a short time-scale, the noise in pulsar timing signals becomes signal-to-noise ratio limited, and it is widely accepted that the atomic clock is more stable [1]. At first glance, pulsar timing would not help much, if one focuses on the short term timing stability. Indeed, most of the previous studies aimed at developing a pulsar time-scale with long-term stability [6, 5, 3, 7, 4], while the short-term stability is not extensively investigated.

The traditional local time-frequency sources based on the \Authorfootnoteatomic clock to generate the frequency standard, GPS technology have two major components: (i) a local and (ii) a GPS receiver as a servo to regularize the long-term timing stability [8]. The atomic frequency standard easily reaches a daily Allan variance of $10^{-13}$, while the GPS system performs rather well in time broadcasting, and the precision of 10 ns or better can be achieved using the commercial systems [9]. Clearly, it is unreasonable to expect that pulsar timing will surpass such high stability on short-timescales.

However, the above mentioned stability of the atomic frequency standard is based on the assumption that the atomic frequency standard is well behaved. In practice, due to environmental variations, frequency standard tuning, or other factors (e.g. accidentally changing clock configuration, replacing cables, or interference in the GPS receiver), the local time standard sometimes shows unexpected jitter. In most of the cases, the time signal jumps for a few micro- to milli-seconds. This is a rather common phenomenon, and appears nearly at all astronomical observatories (for an example, see [10]). Clearly, the very first step to construct a standalone pulsar timescale aiming at real scenario applications would be the identification of those signal jitters and related corrections.

In this paper, we propose to use the timing data of multiple pulsars to diagnose and correct such local time jitter. This is very similar to a pulsar timing array (PTA) experiment, which primarily aims at detecting gravitational waves at the nano-Hz band [11]. We introduce the algorithm to compute the local clock jumps from pulsar timing data in Section 2. In Section 3, we demonstrate the algorithm using pulsar timing data from the Kunming 40m (KM40m) radio telescope. The conclusions and discussions are made in Section 4.

2 Waveform estimation for clock jumps

In the common pulsar timing practice, one measures the TOA at the telescope site. The TOA is then converted to the TOA as seen by an observer at the pulsar’s co-rotating reference frame. During the conversion, one accounts for effects such as Earth motion in the Solar system, electromagnetic wave dispersion, pulsar motion with respect to the Solar-system barycenter. The orbital dynamics of the pulsar also need to be corrected, if it is in a binary system. The next step is to compute the pulsar timing residuals. Here, an integer phase is subtracted from the pulsar-frame TOA, and the difference between the TOA and the modeled integer phase is the timing residual. All the unmodeled effects, such as clock jitter, reside in these timing residuals. We refer the interested readers [12] and references therein for an extensive description of pulsar timing techniques.

The clock jitter or jump introduces a common mode in all pulsars’ timing data, i.e. all the TOAs will be late for a given value at the same epoch, if the clock ticking leaps forward. Such a common mode is very different from other possible noise signals, e.g. pulsar intrinsic noise will not be correlated among pulsars. The clock-induced common mode is identical for all pulsars. We need to identify this identical signal and estimate the waveform in order to compute the clock jump.

Our algorithm to extract the clock jump is based on the maximum likelihood estimator (MLE). For the waveform estimation purpose, the MLE was proven to be asymptotically optimal [13]. Such a property guarantees that the error of waveform estimation approaches the best possible values when the signal-to-noise ratio becomes large.

Clock jitter happens on timescales usually less than a few minutes. We can use pulsar timing data around the jitter epoch to determine the clock jitter amplitude. Neglecting the red noise contribution for short timescales111We can neglect the red noise, if its power is much lower than measurement error. For data span of less than one year, we can neglect red noise of most PTA MSPs such as that measured by [7] for long-term datasets., the timing residuals contain mainly two parts, the clock jitter and measurement error, i.e.

[TABLE]

where we denote the timing residual of the $i$-th pulsar at the $j$-th epoch as $R_{i,j}$, $s_{i,j}$ is the waveform of the clock jitter and $n_{i,j}$ is the measurement error, following the same convention for the indices as for the timing residual ($R$). The noise ($n_{i,j}$) is a zero mean Gaussian random variable. In this paper, we assume that there are extra systematics in determining the errorbar of each TOA. We summarize the systematics using the ‘Efac’ parameters similar to previous studies [14]. The likelihood for the timing residual thus becomes

[TABLE]

Here, the function $f$ is the likelihood function, i.e. the probability density for the given data ($R_{i,j}$). $N$ and $M_{i}$ are the total number of pulsars and epochs for the $i$-th pulsar. $\mu_{i}$ is the local mean of the residual, $\sigma_{i,j}$ is the errorbar of each data point. The ‘Efac’ parameter $\eta_{i}$ is defined for each pulsar to include the minimum practical modeling for the measurement of systematics.

Focusing on the clock jitter, the signal will be identical for all pulsars. We model the clock jitter using a step function of the form

[TABLE]

where $s_{0}$ is the amplitude of the clock jitter, which allows for both positive and negative values corresponding to the two possible ways of clock jitter, i.e. delay or advance. The time epoch of the data points is described by $T_{i,j}$.

We can now split the likelihood into a multiplication of two independence parts, the one before the clock jitter and the one after it, i.e.

[TABLE]

By inserting Equation 3 into Equation 4, one can see that the likelihood is controlled by $2+2N$ parameters. There are two parameters for the clock jitter, the event’s amplitude and epoch ($s_{0}$ and $t_{0}$), $N$ pulsar ‘Efac’ parameters ($\eta_{i}$) and finally $N$ local mean values ($\mu_{i}$). The MLEs for all the parameters are derived according to

[TABLE]

The statistical inference seems to be straightforward, e.g. one can perform Monte-Carlo Bayesian inference with the likelihood function (Equation 4) to measure all $2+2N$ parameters and related errors.

However, for the current problem, such a full Monte-Carlo Bayesian method is not ideal. It has a very high computational cost due to the potentially large number of parameters. On the other hand, the Bayesian method would perform the inference on all parameters, for some of which we are not interested in their values, e.g. $\eta_{i}$ and $\mu_{i}$ are parameters for individual pulsars. Furthermore, we only need to measure the epoch of jitter to the precision of observation cadence due to the sparseness of pulsar timing observations. The requirement for jitter-epoch estimation is rather rough. As we will show in the following paragraphs, it is possible to compute the MLE in a very computationally efficient way, if we perform a one-dimensional search for the jitter epoch ($t_{0}$). That is, we can compute the MLE and error for the clock jump amplitude using a direct method, at each trial jitter epoch.

Since ‘Efac’ parameters are not required, we can marginalize the them by integrating the likelihood function ($f$) under Bayesian prior, and the reduced likelihood ($f^{\prime}$) is then defined as

[TABLE]

where function $p(\eta)$ is the prior of $\eta$. The least informative prior will be $p(\eta)\propto 1/\eta$, i.e. when the $\eta$ probability distributes uniformly in logarithmic space [15]. Inserting Equation 4 into Equation 6, one can show that 222We have used integral of $\int_{0}^{\infty}\eta^{-n-1}\exp\left[-a/(2\eta^{2})\right]\,d\eta=2^{n/2-1}a^{-n/2}\Gamma(n/2)$ [16], where $\Gamma$ is the Gamma-fucntion.

[TABLE]

The above result is very similar to the Student’s $t$-distribution with $M-1$ degrees of freedom, which describes the mean of a normally distributed population with finite sample size and unknown standard deviation. Note that for any time epoch between two available data points, the data defining ‘before’ and ‘after’ are the same, and as such for all time epochs between these points the likelihood value will be the same. In this way, the best measurement precision for the epoch of clock jump is the length between the two nearest data points containing the jump. This becomes more apparent in case where the jump is associated with data gaps longer than the average data cadence, as we will also see in Section 3.

The MLE for $\mu_{i}$ can be readily computed from $\partial f^{\prime}/\partial\mu_{i}=0$, which leads to

[TABLE]

Here the hat symbol $\hat{X}$ stands for the estimator for any given parameter, $X$. The condition for the MLE for $s_{0}$ is $\partial f^{\prime}/\partial s_{0}=0$, which gives an implicit expression

[TABLE]

We can use a two-step iterative method to solve the MLE for $s_{0}$ from Equation 9, as follows: step i. start with the initial value of $s_{0}=0$ and compute the re-scaled errorbar using

[TABLE]

step ii. update value of $s_{0}$ using

[TABLE]

Here, Equation 10 acts to estimate the errorbars of each pulsar using the total standard deviation after removing the jump. Equation 11 uses the re-scaled errorbar to compute the error-weighted average of the jump. The experiment shows that 3 to 5 iterations are usually enough to get the $\sigma_{i,j}^{\prime}$ stable to within 15 digits with the our updating scheme. The above results (Equation 8 and 9) agree with the ad-hoc $\chi^{2}=1$ approach, which is a rather common method to deal with the systematics in errorbar [17].

We now proceed to compute the error of the MLEs. From Equation 11, one can show that the variance of the estimator $\hat{s_{0}}$ has two parts

[TABLE]

where the part $A$ and $B$ are

[TABLE]

Part $A$ comes from the error in estimating the local $\mu_{i}$ using data before the jump, and part $B$ originates from the error of the data after the jump. The error of $\hat{s_{0}}$ is $\delta\hat{s_{0}}=\sqrt{{\rm Var}[s_{0}]}$. From Equation 12, one can see that the final error in the estimator agrees with the optimal weighting averaging. It is clear that we have used all the available information in the data when estimating the clock jitter amplitude.

The above MLE, $\hat{s_{0}}$ and $\hat{\mu_{i}}$, are the estimators for the clock jump amplitude and residual leveling given a pre-determined clock jump epoch. In practical situations, we will not know the clock jump epoch beforehand, and we need to search for it. As shown in the next section, we can simply calculate $\hat{s_{0}}$ on a dense grid of clock-jump epoch trials. The most likely jump is the largest among the trials.

While in this discussion we focus on the clock jump signal, in practice one also needs to account the pulsar timing-parameters fitting. One can implement the timing-model fitting by replacing $\mu_{i}$ in Equation 7 with the linearised timing waveform. In this paper, we take a simpler iterative approach, where we repeat the step of fitting for the clock jump and for timing parameters several times, until the result converges. Note that in each iteration, we subtract the fitted jump from the site arrival times (TOAs at the observatory) when fitting the timing parameters. Such an iterative method also helps to derive the clock jump referring to the telescope site.

3 Example using pulsar timing data from the Yunnan observatory

We first validate our algorithm using simulated data. We simulate timing data of two pulsars (PSRs J1713+0747 and J0437-4715, which are the most accurately timed pulsars in our campaign) with similar length and cadence as the real data, which we are going to discuss right after, in this section. The data is simulated with precision of 5 $\mu$s. We inject a clock jump at MJD 58380 with an amplitude of $5$ $\mu$s. The data and recovered clock jump value are shown in Figure 1. As one can see, the current algorithm recovers the clock jump amplitude well.

We now turn to the real data. The Yunnan Astronomical Observatory (YNAO) of the Chinese Academy of Science operates the Kunming 40-meter radio telescope (KM40m). The telescope was built in 2006 for the Chinese lunar-probe mission. It is located in the south-west of China (${\rm N}25^{\circ}01^{\prime}38^{\prime\prime}$, ${\rm E}102^{\circ}47^{\prime}45^{\prime\prime}$), approximately 15 kilometers away from the city of Kunming. The total collecting area of KM40m is 1250 m2. Our timing data were collected using a room-temperature receiver (70 K system temperature), originally designed for communication purposes. The center frequency is 2.5 GHz and the bandwidth is 300 MHz.

We have been regularly timing five millisecond pulsars since 2017. The observing schedule is much more intense compared to a typical pulsar timing program (e.g. see [18]), as we observe the same set of pulsars almost daily. Part of our timing efforts aim at helping to identify the potential issues in the frequency-clock signal chain of telescope systems. Following the standard timing pipeline, we measured the TOAs using the software PSRCHIVE [19], and computed the timing residual using TEMPO2 [20, 14].

On the 25th of August 2018, the clock signal distributor was broken during an intense lightning storm. Switching to the back-up signal path led to a jump in the pulsar timing signal, as shown in Figure 2.

The search of the clock jump is shown in Figure 3, where we compute the clock-jump value on a grid of jump epoch. The maximum clock jump value appears in the span between MJD 58350 and MJD 58365. After correcting the jump at MJD 58350, no significant clock jump can be further detected. The measured jump is $12.71\pm 0.08$ $\mu$s. We have determined the clock jump with a precision of 80 ns, which is 8 times lower than the original pulsar timing precision (600 ns for PSR 0437-4715). The 80-ns precision mainly comes from the statistical algorithm presented in this paper.

4 Conclusion and Discussion

In this work, we have investigated a computationally efficient method to measure the observatory’s local clock jumps using pulsar timing data as the first step to build the autonomous pulsar time-scale. In order to minimize errors in determining the clock jump, we used Bayesian statistics, with which the unknown systematics can be marginalized with an analytic calculation. The estimator we adopted is the MLE, which is thought to be asymptotically optimal. We found that the Gaussian likelihood after marginalization over the noise scales leads to a distribution similar to the Student’s $t$-distribution. Using the method, we measured the clock jump at Yunnan observatory using the pulsar timing data, where 80-ns precision could be achieved using the observed data. The current precision is worse than that of widely used GPS time delivery service. Using a better receiver (e.g. 20 K noise temperature) with wider observation bandwidth (e.g. 900MHz) will bring the precision down to an 8-ns level, which will be better than that of the current GPS time delivery service. In addition, the precision can also be improved if some larger facilities can be used, for example, the Five-hundred-meter Aperture Spherical radio Telescope (FAST). FAST has just complete its initial commissioning [21] and first studies on pulsars have also been made by using it [22, 23, 24, 25, 26]. Moreover, in general, the clock precision scales as $\sim N_{\rm psr}^{1/2}$, with $N_{\rm psr}$ the number of pulsar observed. If more telescope time is available, we can time more pulsars, which would allow substantial improvement on the precision of pulsar time-scales. In this way, pulsar timing may be used as an independent method to check the local clock stability, if the GPS signal is not available or if the autonomy is the priority.

Due to the intrinsic phase confusion of the pulsar timing technique, the clock-jump value ($s$) seen in the timing residuals is the remainder of the division of the real clock jump ($s_{\rm real}$) by pulsar period $p$, i.e. $s_{\rm real}=s\,({\rm mod\,}p)$. In this way, one will not recover a clock jump that is larger than the rotational period of the pulsar in the array with the longest one. This difficultly can be resolved by observing a long period pulsar in the array. The long-period pulsar can help to derive an initial guess of the jump value. For our analysis, we used PSR B0329+54, which has a rotational period of 0.7 s. The MSP timing data can be preconditioned to the coherent state, and the final clock jump without phase confusion can be measured with the MSP ensemble.

The measured clock jump in the timing residuals refers to the barycentric dynamical time. However, in practice we need the clock jump with respect to the local terrestrial time (most of the times, Coordinated Universal Time). In this paper, we perform the reference frame correction using an iterative approach, where we can measure the clock jump in the barycentric frame, correct the local TOA with the measured jump, reform the timing residual (with timing-parameter fitting), and repeat the iteration. Once the process converges, i.e. when no further clock correction is needed, one determines the clock jump in the terrestrial frame. The iterative method thus also enables us to fit for the pulsar timing parameters.

In the current paper, we focus on the clock jump, so the timing-parameter modeling is not explicitly shown. One can include the timing parameter modeling in the Equation 7 by replacing the central value $\mu_{i}$ with the linearised timing model [14]. Here, the iterative procedure automatically takes care of the timing parameter fitting.

Pulsar timing residuals often contain red noise, i.e. time-correlated noise. In the current paper, we consider white noise only, which reduces the computational cost significantly. There are two major reasons that legitimize our decision to take such a simplified approach. Firstly, the clock jump is highly localized in time, therefore we do not need an extended dataset to measure its value. The contribution of red noise in such a short time scale is small and usually the white noise dominates. Secondly, we have marginalized the ‘Efac’ parameters, which take care of the noise systematics.

In order to get a robust error estimation, we include the ‘EFac’ in our model. We neglected the jitter noise modeling because the radiometer noise (error bars) is much larger than the jitter level. For example, the 1-hr weighted precision of J0437$-$4715 is 600 ns at YN40m, which is 6-20 times larger than the jitter noise measured [27]. Similarly, the current accuracy for the clock jump is 80 ns based on 150-day observation. The jitter noise contribution is further averaged out by another order of magnitude with the current data processing algorithm.

*This work was supported by XDB23010200, 11690024, NSFC U15311243, 2015CB857101 and funding from TianShanChuangXinTuanDui and the Max-Planck Partner Group. *