Pulsar data analysis

The tutorial is intended to provide you a basic introduction to the steps involved in the analysis of pulsar data, including searching for new pulsars/transients and timing a newly discovered pulsar. The filterbank data obtained from GMRT are converted to SIGPROC filterbank format using either the filterbank command from SIGPROC or the rficlean command from RFIClean. The tutorial will use data already converted to the SIGPROC filterbank format.

Introduction

You need to have PRESTO, SIGPROC, TEMPO2, TEMPO, and their dependencies installed on your machine. You also need the SIGPROC filterbank data to be available on your disk.

The header information of the data from a SIGPROC filterbank file can be inspected using the readfile command from PRESTO (an example using the header command from SIGPROC will be shown later).

readfile inputFile.fil

Data-inspection

A section of the raw data can be inspected by plotting the data as a function of time and frequency, e.g., using the waterfaller.py command from PRESTO. The command waterfaller.py several provisions that enable inspecting the data in several ways (e.g., before and after dedispersion, partial and full averaging over the bandwidth, partial averaging over time, etc.)

waterfaller.py inputFile.fil -d 0 --subdm 0 --show-ts -s 64 -T 0 -t 5

The single pulses from a reasonably strong pulsars might be visible after dedispersing the data at correct dispersion measure (DM).

waterfaller.py inputFile.fil -d <DM> --subdm <DM> --show-ts -s 64 -T 0 -t 5

Dedispersion and Folding

A small section of the raw data can be dedispersed and viewed using waterfaller.py as seen above. To dedisperse the whole data and also fold it over the rotation period of the pulsar, we can use prepdata from PRESTO.

prepfold -p 3.7452943 -dm 50.9 -topo -n 128 -nosearch  introTestd4.fil

prepfold is a powerful command, it can search for the most optimum parameters around the given fiducial parameters.

Searching for new pulsars

This part of the tutorial aims to demostrate the process of pulsar search using observations from GMRT. We will use band-4 (550-750 MHz) data of a test pulsar.

Introduction

Pulsars are fascinating extraterrestrial objects that produce pulse like emission with very precise periodicity. These objects are often detected in time domain radio observations. In this tutorial, we will learn how to search for these objects using the strengths of individual pulsed emission as well as their precise periodic behaviour. This tutorial will require PRESTO installed on your machine along with a filterbank file with pulsar in it.

Checking filterbank file

Filterbank file is a file-format for beam data from radio telescopes. This file has the metadata of the observation along with the raw beamdata recorded in the telescope. Use following syntax to check the metadata,

header <filterbank_file.fil>

Example of a filterbank header

Dedispersion

For an unknown pulsar, we don’t know the dispersion caused by the intersteller medium. We try a lot of trial DM values and perform search in each dedispersed time series. We will use prepsubband to dedisperse the filterbank file at a range of trial DM values,

prepsubband -lodm <smallest_trial_dm> -dmstep <seperation_between_dms> -numdms <number_of_trial_dms> -numout <number_of_output_samples> <input_fil_name> -o <output_base_name>

Please note that number of output samples should be an even number (in order to get a straight-forwand FFT in later stage). This will create a number of dedispersed time series with names output_base_name_dm.dat. These are the time series that will be used for single pulse search as well as periodicity search.

Exploring the time series

Now, since we have the dedispersed time series, we can try to visualize the pulses in the time series. Use exploredat to visualize the time series,

exploredat output_base_name.dat

This will generate an interactive plot,

Output plot of exploredat. One can zoom in, zoom out, and move around in this plot.

You can even try to estimate the period of the repeating signal from this plot.

Single pulse search

Though we have seen the single pulses from the pulsar directly in the time series, there is a formal way to search for single pulses in the time series.

single_pulse_search.py -t <significance_threshold> -m <maximum_width> <dedispersed_time_Series.dat>

This generates a text file with information of individual detected pulses for each time series. Now, we need to combine the results of each search in a single summary plot,

single_pulse_search.py -t <significance_threshold> <_dedispersed_time_Series.singlepulse>

Output of single pulse search. We get location of pulses, number of pulses, and single pulse S/Ns as a function of DM.

You can note down timestamp of some strong pulses from the text file and try to varify it’s presence in the time series using exploredat. The candidate plot of a single pulse is contains various diagonstic plots,

Diagnostic plot of a single pulse candidate.

False single pulse candidate.

Search for faint pulsars using periodicity search

There are pulsars that are weak and can not be seen in single pulses. To recover such puslars, we use their periodic property to get a good indication of their presence in the time series. Here, we will see one such method called FFT based periodicity search. The process is to take Fourier transform of the time series and search for strong signal in the power spectra. First we will see how FFT of a time series with repeating signal looks like. To compute the FFT of the time series,

realfft -fwd <dedispersed_time_series.dat>

This will create a file with .fft extension. To see this file in a plot,

explorefft <dedispersed_time_series.fft>

You will see signal at fundamental and subsequent harmonics.

Visualization of power-spectra of a time series. Strong signal can be seen at the fundamental and subsequent harmoics.

An algorithm called harmonic summing is used to add up signal power in the harmonics to increase the detection significance. This algorithm is implemented in accelsearch.

We use accelsearch to search for periodic signal, after cleaning the FFT of timeseries using rednoise

rednoise <dedispersed_time_series.fft>
aacelsearch -zmax 0 -numharm 16 <dedispersed_time_series_red.fft>

This will generate a text file containing the results of the search.

Result of accelsearch. First detection is the first harmonic of the pulsar signal.

Folding filterbank file

Once we know the correct period and DM of the pulsar, we can fold the filterbank file to generate characteristic plots of the pulsar. We use prepfold to fold a filterbank,

prepfold -p <period> -dm <DM> -nosearch -zerodm <filterbank_file.fil>

Result of prepfold. Profile of the pulsar along with subintegration vs phase, frequency vs phase, S/N vs DM, S/N vs period plots.

A false candidate. The pulse seen in the profile is caused by RFI in the time series.

There are lot of such plots for each trial DM in a blind search. All these plots are inspected either by eye or an machine learning based classifier to tell which one is a pulsar and which one is not. Once a pulsar is identified, one needs to check if it’s a new pulsar or an already known pulsar.

Timing of Pulsars

The tutorial is intended to provide you an introduction to the basic steps related to the timing analysis of pulsar using GMRT beam data. We will be using the band-3 (300 - 500 MHz) datasets for this purpose. We will only consider the case of an isolated pulsar for the exercise.

We will use PSRCHIVE and TEMPO2 Softwares.

Introduction

The GMRT beam data is in binary format and is converted to “filterbank” format (same binary file with a header). The “filterbank” is folded based on pulsar ephemeris, yielding three-dimensional (time, frequency, phase) folded pulsar data (“pfd”). For this tutorial, we will begin with “pfd” files.

You need to have PSRCHIVE and TEMPO2 installed on your machine. You also need the “pfd” files to be available on your disk.

Prediction of phases of pulsar using TEMPO2

Right Ascension (RA) and Declination (DEC) provide the source’s two-dimensional position in the sky, which is typically defined by the J2000 system. The pulsar’s spin frequency (F0) is defined as the reciprocal of its spining period. F1 is the first derivative of spin frequency with respect to time. Since the pulsar slows down over time, F1 has a negative value. The dispersion measure (DM) quantifies the line of sight column density of electrons. All of these parameters are listed in a file called the pulsar’s parameter file. Let’s start with a parameter file of pulsar J1646-2142, which contains the pulsar’s name, RA, DEC, F0, F1, and DM values.

cat J1646-2142.par

tempo2 -gr fake -f J1646-2142.par -nobsd 1 -ndobs 7 -randha y -ha 8 -start 59000 -end 60000 -rms 1e-3

This will generate time of arrivals (ToAs) for this pulsar ranging from MJD 59000 to 60000. ToAs are measured for weekly observations (i.e. one observation per seven days). A single ToA is generated for each observation. Each ToA will have a 1e-3 millisecond error (i.e. 1 microsecond). The observation hour angle is chosen at random over an 8-hour period.

At the end of this a file with name J1646-2142.simulate is created. The file J1646-2142.simulate contains the predicted ToAs of the pulsar.

head J1646-2142.simulate

Normally, the first column represents the file name; in this case, it is the same here because it was predicted using tempo2’s fake plugin. The frequency is represented in megahertz in the second column. The third and fourth columns contain the ToAs for the observation and the associated error. Note that we maintained the same error during prediction for all observations. The the last column contains the observatory site code.

The impact of parameter errors on timing data

First check the phases of the pulsar are derived accurately by running the command:

tempo2 -nofit -f J1646-2142.par  J1646-2142.simulate -gr plk

The command above displays the difference between the ToAs predicted by the parameter file J1646-2142.par and the ToAs listed in the timing file J1646-2142.simulate. The residuals must be randomly distributed around the zero line if the ToAs listed in the J1646-2142.simulate file matched the prediction from the J1646-2142.par file. The figure below demonstrates this.

Now let’s introduce a small error (say around an arc-sec) in postion (RA and DEC) in the parameter file to see the effect of wrong position on the timing data. After putting the offset in the parameter file run the above command again which will show that an error in position has yearly periodic variation in the residuals. This is demonstrated in figure below.

If you fit for the RA/DEC for which the error is introduced, you will again get the randomly sampled residuals around the zero line in the post-fit residuals.

Now, reset the parameter file to its default values and try introducing a small error in F0. The residuals for the F0 error show a straight line pattern with no periodicity, as shown below.

Again, if you fit for the F0, you will get the randomly sampled residuals around the zero line in the post-fit residuals.

Try the same thing again with error in F1, and this time you’ll notice a parabolic long-term variation like the one in the image below.

Understanding the pattern in the residuals caused by parameter error will assist us in solving a newly discovered pulsar. Once accurate (phase-coherent) solutions for a pulsar are obtained, one can establish its clock-like behaviour and predict its phases at any point in time.

Time to solve a Newly Discovered Pulsar

This time we will start with the real GMRT data set for a newly discovered pulsar J1244-4708 which is discovered in the GMRT High Resolution Southern Sky Survey.

This pulsar was localised through imaging of a single observational epoch which provides us with its RA and DEC which are accurate to a few arc-seconds. The F0 and DM values of the pulsar were estimated using data from an epoch observation in which the pulsar was strongly detected. We used the reciprocal of the period value (in seconds) given in the waterfall plot of that observation to calculate F0. Similarly, from the same plot, we noted the DM and epoch of the observation. Using these values, we created a parameter file for this pulsar (in this case, J1244-4708.par) that contains an initial guess of its parameters.

cat J1244-4708.par

In this exercise, we will use GMRT’s incoherent array (IA) beam data. This pulsar has been observed 13 times over a baseline of 9 months. The GMRT IA beam data is first converted to “filterbank” format, then dedispersed and folded to “pfd” format. We’ll start with the “pfd” files right away.

Time of arrival estimation

The “pfd.ps” files contain waterfall plots for all “pfd” files. The first step is to go through all of the profiles and find the cleanest one with the sufficiently bright integrated pulse. We will use that profile as our reference or template profile, assuming it is the true profile of the pulsar. Use the following command to view all the profiles.

evince *ps

After selecting the template profile, the next step is to estimate how many ToAs should be extracted from each observation. We will derive a single ToA for full observation for weak detections, two ToAs for mildly detectect profiles, and four ToAs for bright detections. This approach will be adequate to find the solution in this specific case. To determine the number of ToAs for individual observations, check all the profiles again using the command above.

Let’s say you’ve selected k ToAs for a particular observation A.pfd and B.pfd is your chosen epoch for the template. Then run the following command to get the ToAs for observation A.pfd.

get_TOAs.py -n k -t B.bestprof A.pfd >> J1244-4708.tim

This command will find the ToAs for observation A.pfd by cross-correlating it with template B.bestprof. The ToAs will be saved in the J1244-4708.tim timing file.

The value of k varies depending on the detection strength of the pulsar for each observation. So, for each observation, use the same command, but change the file names A.pfd and the corresponding k value while keeping the template same (i.e., B.bestprof). Once all observations’ ToAs have been generated, run the command below.

tempo2 -nofit -f J1244-4708.par J1244-4708.tim -gr plk

This command displays the difference (or residuals) between the observed ToAs (in the J1244-4708.tim file) and the predicted ToAs (predicted using J1244-4708.par file). Because the ToAs for a newly discovered pulsar are not phase-connected, we will not see any general systematic pattern in the residuals. This is illustrated in the figure below.

Now, take a sample of densely sampled points and start fitting from F0. Once you’ve identified a pattern in the ToAs, try fitting RA and DEC. Finally, you can try fitting F1, but keep in mind that the value of F1 for 9 months of data will be highly unreliable. Once phase coherent solutions are obtained, ToAs will exhibit near-random behaviour around the zero line, as shown in the figure below.

Examine the F0, F1, RA, and DEC values. Since RA and DEC are obtained through imaging, the post-fit RA and DEC should be within a few arc-seconds of the initial guess (i.e., pre-fit values). The F1 value after fitting should be negative and small (compared to error in F0). Below are the fitted values for this exercise.

Finally, on the graphical interface, select “new par” to create the new parameter file. Now that the pulsar’s phase-coherent solution has been found, you will be able to predict the ToAs in future observations. The prediction accuracy improves as the number of observations and the timing baseline increase.

Acknowledgements and Credits

The Searching for new pulsars and Timing of pulsars sections of the tutorial were entirely contributed by Shubham Singh and Shyam Sunder, respectively. The data sets used in these parts were kindly provided by Dr. Jayanta Roy and Dr. Bhaswati Bhattacharyya. Yogesh Maan is responsible for the overall pulsar tutorial coordination as well as for making available the data used in the Introduction section.