Timing noise

Timing noise is defined as a continuous, noise like flucutation in rotation rate. It is well known that young pulsars show large "amounts" of timing noise, whereas millisecond pulsars show relatively little. We wish to explore statistical methods for analysing pulsar timing noise and to apply these methods to both young pulsars and the sample of millisecond pulsars in the PPTA sample.

The work described below was undertaken by Bill Coles and George Hobbs.

The research is now being continued by Sheila Kanani - a summer student at the ATNF. Details of her work are available here

Analysis methods

Problem definition: For a given data-set of pulsar timing residuals we would like to 1) search for periodic variations (such as might indicate the presence of an undetected planet) and 2) estimate the spectral exponent of the red noise process that characterises the low frequency behaviour of the residuals.

Technique: Normally we would use a standard power spectrum analysis which is optimal for detection of periodic signals and can easily be used to estimate the exponent of a power-law spectrum. However, our time series are not equally sampled, so classical spectrum analysis cannot be applied immediately. Also, the time series often have very steep spectra so leakage in spectral estimation is very important.

Normally when analysing time series with very steep spectra one would use a "prewhitening" technique and perhaps a low sidelobe spectral window (e.g. the Hanning window). It is even more important to use some form of prewhitening when the data are unequally spaved because this inevitably causes high sidelobes in the spectral window. There are many ways to deal with the problem, but there is no univeral solution because it depends so much on the nature of the actual sampling and the characteristics of the actual data. There appear to be nearly as many algorithms as there are practitioners of the art!

In our application the samples are clustered as tight groups of a few daily observations spaced by a few weeks. The spectra appear to be very steep at low frequencyes, with spectral exponents of -6 or steeper, but they flatten at a frequency of the order of a cycle/year to an exponent of -2 or so. We have developed an effective procedure that is easily understood and modified that doesn't appear to introduce significant artifacts in the estimated spectra. The idea is to separate the steep and flatter spectrum components so they can be analysed separately.

We begin by interpolating the raw data to an equally spaced grid using a piecewise cubic hermite interpolating polynomial (PCHIP in Matlab). This has some advantages for our purposes, over a cubic spline. In particular it has no overshoots and less oscillation if the data are not smooth. We then lowpass filter this interpolated series using a simple block smoother of adjustable length and we discard the transients at either end of the series. The smoothing block is typically of the order of a year, so the smoothed series is heavily bandlimited. As we intend to subtract the smoothed component from the raw data, we need a linear-phase filter. Any symmetric FIR filter would satisfy this constraint. We chose a simple block smoother becauase it appears to be adequate for our purpose and it has the minimum transient length, so we lose the minimum amount of the time series possible for a given bandwidth. This preserves maximum frequency resolution and maximum sensitivity to periodic components. Thus we obtain the low-frequency of "red" component.

We then interpolate the red component back to the orginal sample spacing using the same PCHIP interpolator. We could use a bandlimited interpolation here, but any interpolator is adequate because the red component is much narrower in bandwidth than the original data. We subtract the red component from the original data to obtain the high-frequency or white component at the original sample spacing.

We spectrum analyse the equally spaced red component first prewhitening by a second difference, then forming the periodogram classically, finally postdarkening by the appropriate sin^4(pi f/fny) function. We use two windows, a rectangular and a Hann window in each case. Comparing these two windows gives us warning when leakage is important. The periodigram samples are well-behaved statistically. They all have chi-squared distributions. Those from the rectangular window are independent, and those from the Hann window are partially dependent as that window is twice the width of the rectangular window.

We spectrum analyse the unequally spaced which component using the Lomb-Scargle scheme, modified in that we do not normalise by the variance. Effectively this is the same as a least squares fit of equally spaced sinusoidal components. The statistics of these spectral estimates are also chisquared and they are roughly independent. We do not oversample in frequency, or sample past the Nyquist frequency.

To estimate the significance of any "spikes" in the data such as would indicate periodic components, we fit smooth functions to the estimated spectra and plot the 95% confidence limits around these smooth functions. For the red component we fit Pr(f) = (Af^-exp1 + Bf^-exp2)*sinc^2(fNS)/ Where A and B are found using linear least squares. Typically exp1 = -7 and exp2 = -2, but we adjust these exponents roughly to best fit the data. The white component roughly fits this model, i.e. Pw(f) = (A f^-exp1 + Bf^-exp2)*(1-sinc(fNS))^2, but this model is not adequate to give good confidence limits. We obtain a better smoothed spectrum simply by smoothing the actual spectra estimates by a rectangular block of length 25.

An example of this process is shown below. This is a real data set from B1929+10 with an artificial planet of period 500 d. The planet was adjusted to be the smallest we could detect. The time series is shown first. The upper panel has raw data, interpolated data and red component. The lower panel has the white component with actual samples overlaid. An expanded plot is shown beneath the first one.