Pulse Arrival Time Computation Algorithms

The technique of pulsar timing uses the predictability of radio pulsar signals as a clock tick that can provide detailed information about both the pulsar itself and the system in which it resides. In most cases, timing can yield parameters like the dispersion measure along the line of sight to the pulsar, provide a basic estimate of the spindown age, and so on, probing the evolution and distribution of stars and plasma in the galaxy. In rare circumstances, timing can been used to test general relativity through the detection of effects like orbital period decay in the relativistic binary system B1913+16. In this case, the rate of binary period change (as measured through pulse timing) was seen to match that expected due to loss of energy by gravitational radiation.

Pulsar timing involves a complicated analysis chain, from sophisticated instrumentation at the front end of a telescope, all the way through to detection and processing using high-speed computers and custom-made software. Each step in the process must be carefully examined and tested in order to ensure that the results of any experiment are valid. Part of the work done at Swinburne involves testing the software routines provided in the PSRCHIVE scheme to ensure their validity under a wide range of circumstances.

The focus of this section is our investigation into methods for computing the arrival time of a pulse, given an archive with an accurate time stamp. This process has two aspects, the first is how to choose an appropriate standard template profile against which to time the observation, and the second is how to use the template to extract an arrival time from the data.

Template Creation

Analytic Templates

• Noiseless by definition
• Flexible to changing profile structure
• Cannot fully recover fine detail

Integrated Templates

• Contain additional or copied noise (self standarding)
• Static objects... not always suitable for evolving profiles
• Match the "true" profile well in a small parameter space

Examples of analytic standards include the method described by Kramer et al. (1999, ApJ 520, p324-334) that involves decomposition of the profile into a sum of gaussian components, some of whose parameters can be allowed to vary. While static, integrated templates are likely to be the best possible representation of the characteristic pulse profile, this is only true over small bandwidths, especially if the profile of the pulsar evolves rapidly with frequency (or any other parameter).

Once a standard template profile has been built, it must somehow be compared with the observed profile in order to determine the offset from the time stamp recorded in the data file (which may represent the first phase bin of the profile for example). This usually involves some form of cross-correlation, the maximum point of which will represent the relative offset in pulse phase space. The complexity is introduced when we realise that it is possible, through the application of clever algorithms, to time a pulsar to significantly greater precision than the width of one binning unit. This can be done using one of two methods:

Arrival Time Computation

• Mathematical interpolation of the cross-correlation function
• Fitting for a phase gradient in the Fourier domain
• (Both methods can be extended to use full polarimetric information)

The cross correlation function will measure the relative shift to within one phase bin of precision, but it is possible to extrapolate across the gaps in the discrete CCF using a simple mathematical function which can then be used to estimate the true shift. It is also possible to move into the Fourier domain and use that fact that two profiles that differ only by a relative shift will have fourier transforms that differ only by a linear phase gradient. It is possible to fit for this phase gradient and so determine the original time domain shift. The strengths and weaknesses of the two methods are still unclear and we are in the process of developing simulation code to examine the algorithms in more detail. It should be noted that the simplest way to compute an arrival time is to use the Stokes I (total intensity) profile. However, in the case where more polarimetric information is available, the methods above can be extended into additional dimensions that may provide an even more precise shift measurement.

The PSRCHIVE scheme currently supports both schemes. The interpolation method makes use of a parabolic fit to the point of maximum correlation and it's two nearest neighbours. The apex of this parabola is taken to be the new point of maximum correlation. The error associated with this fit is somewhat difficult to determine analytically as it involves a component due to the intrinsic width of the parabola and also to the goodness of fit. The figure below illustrates this method:

Testing arrival time routines requires a data set whose parameters are precisely known. Relevant parameters include signal to noise ratio, profile shape, width and the distribution of shifts from some model. The only way to satisfy all these requirements is to simulate a data set by building artificial profiles. This is the job of the SyntheticProfile class in the PSRCHIVE scheme. This class can create Gaussian profiles with injected random noise that precisely match a given phase location, width and height relative to the RMS of the baseline. An ensemble of these synthetic profiles can be created and timed to see if our arrival time algorithms can re-produce the parameters we begin with.

The following plots show the first results from our simulations of a set of Gaussian profiles with identical structure and shape, as a function of increasing signal to noise. The first plot makes use of the interpolation algorithm and the second uses the Fourier shift algorithm.

Each simulation step increases the signal to noise of the profile ensemble by unity. Only the lower left-hand plot integrates with each step, the others simply illustrate the workings of the program. Each ensemble is given a set of pre-defined shifts, either an analytic ramp (as shown), a random distribution, or a single value. The PSRCHIVE timing algorithms are asked to report on the shift of each profile in the ensemble, using a standard profile constructed to have the same shape, but much higher signal to noise. The sucess or failure of the different algorithms can be measured in a variety of ways that are still being perfected. Preliminary results indicate that interpolation-based methods are better when the signal to noise of the observed or standard profiles is low, but Fourier techniques are superior when the profiles are not noise dominated.

408 MHz map of the galactic plane, courtesy of NASA