Fermi Gamma-ray Space Telescope

Pulsar Tools Anatomy

This article describes data processing procedures by the pulsar analysis tools distributed as a part of a Science Tools package. The versions of the tools described here is the ones included in Science Tools version v9r8p2: gtbary version v6r2p4, gtephem version v8r2, gtpulsardb version v8r2, gtpphase version v8r3, gtophase version v8r3, gtpsearch version v10r5, gtpspec version v10r5, gtptest version v10r5. A different version of those tools may or may not behave the same way as described in this article. All of the above tools report their version numbers on the first line of the tools' textual output. If you find the tool behave differently than described in this article, please check the version number of the tool in its textual output against the ones listed above.

Chapter 1. Data processing procedures of pulsar analysis tools

In this chapter, data processing procedures are described step by step in great detail for gtephem, gtpphase, gtophase, gtpsearch, gtpspec, and gtptest. Those tools are designed to play central roles in users' pulsar analysis. In this chapter, data processing procedures by those tools are described in detail step by step.

gtephem

  1. The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
  2. The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
  3. The tool computes a reference time of computation, or a moment in time that reftime parameter points to as in a time format specified by timeformat parameter in a time system specified by timesys parameter.
  4. The tool chooses the best spin ephemeris for the reference time as described in Appendix C. If strict parameter is set to true, ephemeris validity is required. Otherwise, ephemeris validity is not required.
  5. The tool displays the chosen spin ephemeris.
  6. The tool chooses the best orbital ephemeris for the reference time as described in Appendix C.
  7. The tool displays the chosen orbital ephemeris.
  8. With the chosen spin ephemeris, the tool computes and displays the following quantities at the reference time, as described in Appendix A: the pulsar's Right Ascension and Declination, a pulse phase, a pulse frequency, the first time-derivative of pulse frequency, and the second time-derivative of pulse frequency.
  9. The tool searches for an ephemeris remark that is effective between the reference time and the reference epoch of the chosen ephemeris. If found, the tool produces a warning message listing ephemeris remarks that are found.

gtpphase

  1. The tool opens input event file(s) as described in Step 1 of Chapter 3.
  2. The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
  3. The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
  4. The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
  5. The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
  6. The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
  7. The tool reads pphasefield parameter from the parameter interface of the tool, and use it for the name of a FITS column which the tool writes pulse phase values.
  8. For all given event files, the tool creates a new FITS column with the name if it doesn't exist.
  9. The tool reads pphaseoffset parameter from the parameter interface of the tool, and keeps the value during the process.
  10. For each event in the given event file(s), the tool takes the following steps to assign a pulse phase to the event.
    1. The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
    2. The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
    3. The tool chooses the best spin ephemeris for the corrected arrival time as described in Appendix C requiring ephemeris validity.
    4. With the chosen spin ephemeris, the tool computes a pulse phase at the corrected time as described in Appendix A, adds pphaseoffset parameter value to it, then truncates the integer part of the value so that the result is non-negative and less than one (1).
    5. The tool writes the truncated phase value to a FITS column specified by pphasefield parameter.

gtophase

  1. The tool opens input event file(s) as described in Step 1 of Chapter 3.
  2. The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
  3. The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
  4. The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
  5. The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
  6. The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
  7. The tool reads ophasefield parameter from the parameter interface of the tool, and use it for the name of a FITS column which the tool writes orbital phase values.
  8. For all given event files, the tool creates a new FITS column with the name if it doesn't exist.
  9. The tool reads ophaseoffset parameter from the parameter interface of the tool, and keeps the value during the process.
  10. For each event in the given event file(s), the tool takes the following steps to assign an orbital phase to the event.
    1. The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
    2. The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
    3. The tool chooses the best orbital ephemeris for the corrected arrival time as described in Appendix C.
    4. With the chosen orbital ephemeris, the tool computes an orbital phase at the corrected time as described in Appendix A, adds ophaseoffset parameter value to it, then truncates the integer part of the value so that the result is non-negative and less than one (1).
    5. The tool writes the truncated phase value to a FITS column specified by pphasefield parameter.

gtpsearch

  1. The tool opens input event file(s) as described in Step 1 of Chapter 3.
  2. The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
  3. The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
  4. The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
  5. The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
  6. The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
  7. The tool chooses the best spin ephemeris for the time origin computed in Step 5-2 of Chapter 3.
  8. With the chosen spin ephemeris, the tool computes a pulse frequency at the time origin.
  9. The tool computes the Fourier resolution of the given event data in the following procedure.
    1. The tool computes the start time and the stop time of the given event data as described in Appendix B.
    2. The tool applies arrival time corrections to the start time and the stop time as described in Appendix D. If it is not successful, then the tool produces an error message and terminates the process.
    3. The tool computes a time difference between the corrected start time and the corrected stop time in the time system that is determined to be used for time series analysis in Step 5-2 of Chapter 3.
    4. The tool computes an inverse of the time difference computed in the previous step, and takes it as the Fourier resolution of the given event data.
  10. The tool determines trial frequencies in the following procedure.
    1. The tool reads numtrials parameter from the parameter interface of the tool, and sets it to the number of trial frequencies.
    2. The tool uses the pulse frequency at the time origin computed above to determine the center of trial frequency range. If the number of trial frequencies is odd, the computed frequency is set to the center of trial frequency range. If the number of trial frequencies is even, the computed frequency is set to the first frequency after the center of trial frequencies. For example, if the number of trial frequencies is 11, the computed frequency is set to the 6th trial frequency. If the number of trial frequencies is 12, the computed frequency is set to the 7th trial frequency.
    3. The tool reads scanstep parameter from the parameter interface of the tool, and multiplies it to the Fourier resolution computed in the previous step, and uses it as a frequency step size, or a difference in frequency between neighboring trial frequencies.
    4. The tool computes and sets all trial frequencies based on the central frequency and the frequency step size, both of which are computed above.
  11. The tool evaluates algorithm parameter, chooses a statistical test to perform on the given data, and initializes the chosen statistical test as shown in the following table.
    algorithm value Statistical test Initialization
    CHI2 χ2 test The tool reads numphase parameter from the parameter interface of the tool, and uses it for the number of phase bins of a folded light curve.
    RAYLEIGH Rayleigh test Same as Zn2 test for n = 1.
    Z2N Zn2 test The tool reads numharm parameter from the parameter interface of the tool, and uses it for the number of harmonics to sum up (i.e., the value of the suffix "n" of Zn2 test).
    H H test The tool reads maxharm parameter from the parameter interface of the tool, and uses it for the maximum number of harmonics to evaluate in searching for the optimal number of harmonics to sum up.

    Those statistical tests are well-known and their definitions and descriptions are found in the literature. For example, Leahy et al. (1983) described χ2 test in detail, and Vaughan et al. (1994) pointed out importance of noise-signal interactions in χ2 test. Also, Buccheri et al. (1983) explained Zn2 test as a part of their data analysis, and De Jager et al. (1989) defined and discussed H test.

    References
    Buccheri et al. 1983, A&A 128, 245
    De Jager et al. 1989, A&A 221, 180
    Leahy et al. 1983, ApJ 266, 160
    Vaughan et al. 1994, ApJ 435, 362

  12. For each event in the given event file(s), the tool takes the following steps to accumulate contribution(s) of the event to the test statistic for the chosen statistical test at all the trial frequencies.
    1. The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
    2. The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
    3. The tool computes an elapsed time at the corrected arrival time since the time origin computed in Step 5-2 of Chapter 3. The computation is performed in a time system that is determined to be used for time series analysis in Step 5-2 of Chapter 3.
    4. For each of all the trial frequencies, the tool computes a contribution of this event to the test statistic of the chosen statistical test. First, a pulse phase is computed for each trial frequency. Then, if χ2 test was chosen in the previous step, one event is added to an appropriate phase bin of a folded light curve for the trial frequency. If Zn2 test, Rayleigh test, or H test was chosen, the tool uses the computed pulse phase to compute the sine and the cosine components of Zn2 test for all harmonics that are needed to compute the chosen test statistic, and add them to the sum of those components over events that have already been processed.
  13. The tool computes values of the test statistic for the chosen statistical test for all trial frequencies.
  14. The tool identifies the trial frequency that gives the maximum statistic, and evaluates a chance probability that a value of the test statistic of the chosen statistical test exceeds the maximum value found under the given test condition, under the assumption that the no pulsations exist in the given time series data. The procedure to compute it can be further broken down as follows.
    1. The tool computes a single-trial chance probability, p, which is a probability that a value of the test statistic of the chosen statistical test exceeds the maximum statistic identified in the previous step after one trial. If no pulsations exist in the given time series data, the test statistic follows a probability distribution shown in the table below, where numphase and numharm are values of the named parameters taken from the parameter interface of the tool.
      algorithm value The test statistic follows...
      CHI2 χ2 distribution with (numphase - 1) degrees of freedom
      RAYLEIGH χ2 distribution with 2 degrees of freedom
      Z2N χ2 distribution with numharm×2 degrees of freedom
      H special probability distribution for H test, defined by De Jager et al. (1989)

      For a technical reason, the tool computes a range of probability that the true value of the single-trial chance probability lies. The algorithm to compute a probability range for a given statistic value is described in a separate document titled "Numerical integration of the chi-squared probability density function" (PDF file, 148 kB) for your reference.

      References
      De Jager et al. 1989, A&A 221, 180

    2. The tool computes the number of independent trials, N, which is the number of trial frequencies that are considered statistically independent of each other. Letting fmin be the smallest trial frequency, fmax the largest, and fFourier the Fourier resolution of the given event data, the tool first computes the maximum number of independent trial frequencies that can fit in the range of trial frequency, Nmax , by:

      Nmax = ceil((fmax - fmin) / fFourier)

      where ceil(x) is the ceiling function of x whose value is the smallest integer that is greater than or equal to x. Then, the tool compares it with the actual number of trial frequencies, Ntrial (which is specified by numtrials parameter in the parameter interface of the tool), and determines N to be the smaller of them, or:

      N = min(Nmax , Ntrial)

    3. The tool computes a multi-trial chance probability, Q, which is defined as:

      Q = 1 - (1 - p)N

      whose computation procedure is described in a separate document titled "Numerical computation of false alarm probability for multiple trials" (PDF file, 148 kB) for your reference. Note that Q is computed for each of the lower bound of p and the upper bound, and the tool takes them as the lower and the upper bounds of Q, respectively.
  15. The tool evaluates outfile parameter, and creates an output FITS file unless the value of the parameter is NONE. The value of outfile parameter is used for the name of the output file, and the output file contains the test result, including values of the test statistic at all trial frequencies, a multi-trial chance probability computed in the previous step, and other information related to the statistical test.
  16. The tool reports the same test result as in an output FITS file for human reading.
  17. The tool evaluates plot parameter and title parameter to create a plot showing the result.
    1. If the value of plot parameter is a logical true, the tool creates a graphical output that plots values of the test statistic against trial frequencies. Otherwise, the tool does nothing in this step.
    2. If the value of title parameter is "DEFAULT", the tool uses a pre-defined title for a title of the plot. Otherwise, the value of title parameter is used for the plot title.

gtpspec

  1. The tool opens input event file(s) as described in Step 1 of Chapter 3.
  2. The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
  3. The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
  4. The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
  5. The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
  6. The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
  7. The tool computes the start time and the stop time of the given event data as described in Appendix B.
  8. The tool applies arrival time corrections to the start time and the stop time as described in Appendix D. If it is not successful, then the tool produces an error message and terminates the process.
  9. The tool reads binwidth parameter from the parameter interface of the tool to define time bins for DFFT (Discrete Fast Fourier Transform) algorithm. A time bin is a time interval that the tool counts the number of photon arrival times that fall in. In this step, the tool defines time bins with the length specified by binwidth parameter, starting at the corrected start time.
  10. The tool reads numbins parameter from the parameter interface of the tool, and defines time segments for DFFT computations. A time segment is a group of successive time bins that DFFT is performed at once. In this step, the tool defines time segments such that the number of time bins in one time segment is equal to the value of numbins parameter, starting at the corrected start time.
  11. For each event in the given event file(s), the tool takes the following steps to determine a time segment and a time bin that the event belongs to.
    1. The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
    2. The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
    3. The tool computes an elapsed time at the corrected arrival time since the corrected start time. The computation is performed in a time system that is determined to be used for time series analysis in Step 5-2 of Chapter 3. Then:
      • if the corrected photon arrival time is at or later than the corrected start time and at or earlier than the corrected stop time, the tool divides the computed elapsed time by the length of a time segment (= binwidth×numbins), and takes the ratio as the segment number and the reminder the time bin number.
      • Otherwise, the event is discarded, i.e., the tool stops processing this event and moves on to the next event.
  12. For each time segment, the tool takes the following steps to compute a power spectral density using DFFT algorithm.
    1. The tool counts the number of events that belongs to each time bin in a time segment being processed to create a light curve, or an array whose element stores the number of events in a time bin.
    2. The tool computes the average number of events per time bin, and subtract the result from the element of the light curve computed in the previous step.
    3. The tool performs DFFT on the time segment.
    4. At every frequency of the computed Fourier transform of the light curve, the tool computes a Fourier power, or a squared sum of the real part and the imaginary part of the Fourier transform.
  13. The tool adds the Fourier power spectra for all the time segments computed in the previous step.
  14. The tool reads lowfcut parameter, and uses it for the lower boundary of frequency for a pulsation search. The tool identifies the lowest frequency bin in the computed Fourier transform whose frequency is not less than the lowfcut parameter value. In all the procedure described below, the tool ignores all frequency bins below the lowest frequency bin.
  15. The tool identifies the trial frequency that gives the maximum statistic, and evaluates a chance probability that a summed Fourier power exceeds the maximum value found under the given test condition. The procedure to compute it can be further broken down as follows.
    1. The tool computes a single-trial chance probability, p, which is a probability that a summed Fourier power exceeds the maximum Fourier power identified in the previous step after one trial. If no pulsations exist in the given time series data, a summed Fourier power follows the χ2 distribution with 2×M degrees of freedom, where M is the number of time segments that contains at least one event. For a technical reason, the tool computes a range of probability that the true value of the single-trial chance probability lies. The algorithm to compute a probability range for a given statistic value is described in a separate document titled "Numerical integration of the chi-squared probability density function" (PDF file, 148 kB) for your reference.
    2. The tool counts the number of frequency bins that are searched for the maximum Fourier power, N, starting at the lowest frequency bin identified in the earlier step and ending at the last frequency bin which corresponds to the Nyquist frequency.
    3. The tool computes a multi-trial chance probability, Q, which is defined as:

      Q = 1 - (1 - p)N

      whose computation procedure is described in a separate document titled "Numerical computation of false alarm probability for multiple trials" (PDF file, 148 kB) for your reference. Note that Q is computed for each of the lower bound of p and the upper bound, and the tool takes them as the lower and the upper bounds of Q, respectively.
  16. The tool evaluates outfile parameter, and creates an output FITS file unless the value of the parameter is NONE. The value of outfile parameter is used for the name of the output file, and the output file contains the test result, including the summed Fourier powers at and above the lowest frequency determined in the earlier step, a multi-trial chance probability computed in the previous step, and other information related to the pulsation search.
  17. The tool reports the same test result as in an output FITS file for human reading.
  18. The tool evaluates plot parameter and title parameter to create a plot showing the result.
    • If the value of plot parameter is a logical true, the tool creates a graphical output that plots the summed Fourier powers against pulse frequencies. Otherwise, the tool does nothing in this step.
    • If the value of title parameter is "DEFAULT", the tool uses a pre-defined title for a title of the plot. Otherwise, the value of title parameter is used for the plot title.

gtptest

  1. The tool opens input event file(s) as described in Step 1 of Chapter 3.
  2. The tool checks whether all the given event files contain a FITS column with the name specified by pphasefield parameter. If it is missing in one of the event files, the tool produces an error message and terminates the process.
  3. The tool prepares and initializes all the statistical tests that are currently implemented: χ2 test, Rayleigh test, Zn2 test, and H test.
    • The tool reads numphase parameter from the parameter interface of the tool, and the parameter value is used for the number of phase bins of a folded light curve in χ2 test.
    • The tool reads numharm parameter from the parameter interface of the tool, and the parameter value is used for the number of harmonics to sum up in Zn2 test (i.e., the value of the suffix "n" of Zn2 test).
    • The tool reads maxharm parameter parameter from the parameter interface of the tool, and the parameter value is used for the maximum number of harmonics to evaluate in searching for the optimal number of harmonics to sum up in H test.
  4. For each event in the given event file(s), the tool takes the following steps to accumulate contribution(s) of the event to the test statistics for all the statistical tests.
    1. The tool reads a FITS column specified by pphasefield parameter.
    2. The tool computes a contribution of this event to each test statistic of all the statistical test. For χ2, one event is added to an appropriate phase bin of a folded light curve. For Zn2 test, Rayleigh test, and H test, the tool uses the pulse phase to compute the sine and the cosine components of Zn2 test for all harmonics that are needed to compute the test statistics, and add them to the sum of those components over events that have already been processed.
  5. The tool computes values of all the test statistics.
  6. The tool evaluates a chance probability that a value of each test statistic exceeds the value computed from the given event file(s), under the assumption that the no pulsations exist in the given time series data. Specifically speaking, null hypotheses for all the implemented statistical test are following.
    • For χ2 test, the test statistic is assumed a sample from the χ2 distribution with (numphase - 1) degrees of freedom, where numphase is the value of numphase parameter taken from the parameter interface of the tool.
    • For Rayleigh test, the test statistic is assumed a sample from the χ2 distribution with 2 degrees of freedom, where numharm is the value of numphase parameter taken from the parameter interface of the tool.
    • For Zn2 test, the test statistic is assumed a sample from the χ2 distribution with (numharm × 2) degrees of freedom, where numharm is the value of numphase parameter taken from the parameter interface of the tool.
    • For H test, the test statistic is assumed a sample from the special probability distribution for H test, defined by De Jager et al. (1989).
    For a technical reason, the tool computes a range of probability that the true value of the single-trial chance probability lies. The algorithm to compute a probability range for a given statistic value is described in a separate document titled "Numerical integration of the chi-squared probability density function" (PDF file, 148 kB) for your reference.

    References
    De Jager et al. 1989, A&A 221, 180

  7. The tool evaluates outfile parameter, and creates an output FITS file unless the value of the parameter is NONE. The value of outfile parameter is used for the name of the output file, and the output file contains the test result, including the folded light curve created for χ2 test, values of the test statistics, and other information related to the statistical tests.
  8. The tool reports the same test result as in an output FITS file for human reading.
  9. The tool evaluates plot parameter and title parameter to create a plot showing the result.
    • If the value of plot parameter is a logical true, the tool creates a graphical output that plots the folded light curve created for χ2 test. Otherwise, the tool does nothing in this step.
    • If the value of title parameter is "DEFAULT", the tool uses a pre-defined title for a title of the plot. Otherwise, the value of title parameter is used for the plot title.

Chapter 2. Other supporting tools for pulsar analysis

There are two more pulsar-related tools in a Science Tools package: gtbary and gtpulsardb. Those tools use the same components for fundamental computations such as arrival time corrections, time system conversions, and pulsar ephemerides loading, as the ones used by the tools described in the previous chapter. But, due to their distinct functionality, their procedural steps are different. In this chapter, data processing procedures by the two supporting tools are described in detail step by step.

gtbary

  1. The tool evaluates tcorrect parameter to determine a type of arrival time corrections to perform. If tcorrect parameter is set to BARY, then the tool performs barycentric corrections. If tcorrect parameter is set to GEO, then the tool performs geocentric corrections.
  2. The tool opens an input event file specified by evfile parameter. It is not allowed to give more than one event file at a time.
  3. The tool copies the input event file to a temporary output file. The name of the temporary output file is the value of outfile parameter with ".tmp" attached at the end. For example, if outfile parameter is set to "gtbary_output.fits", then the temporary output file name is "gtbary_output.fits.tmp".
  4. The tool reads solareph parameter from the parameter interface of the tool, and determines a solar system ephemeris to use in barycentric corrections.
  5. The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
  6. The tool reads ra and dec parameters from the parameter interface of the tool, and uses them as the pulsar's Right Ascension and Declination for which geocentric or barycentric corrections are performed. Note that the values of ra and dec parameters are interpreted in the same celestial reference frame as the one used by the solar system ephemeris chosen in the earlier step. To be specific, FK5 is used by JPL DE200, and ICRS by JPL DE405.
  7. The tool updates values of the following header keywords of all extensions in the temporary output file.
    Header keyword The value is set to...
    TIMESYS TDB if tcorrect = BARY, and TT if tcorrect = GEO.
    TIMEREF SOLARSYSTEM if tcorrect = BARY, and GEOCENTRIC if tcorrect = GEO.
    RA_NOM Right Ascension used for geocentric or barycentric corrections.
    DEC_NOM Declination used for geocentric or barycentric corrections.
    RADECSYS the name of celestial reference frame (FK5 or ICRS).
    PLEPHEM solar system ephemeris used for geocentric or barycentric corrections (JPL-DE200 or JPL-DE405).
    TIMEZERO 0. (zero)
    CREATOR gtbary v6r2p4
    DATE file creation date in YYYY-MM-DDThh:mm:ss format.
    TIERRELA the same value as in the input event file, and 1.e-9 if it is missing in the input event file.
  8. The tool reads scfile, sctable, and angtol parameters from the parameter interface of the tool, and stores them in a memory space for later use.
  9. The tool reads the following time expressions in the input event file, applies geocentric or barycentric corrections to them, and writes the result to the temporal output file.
    • TSTART header keyword value
    • TSTOP header keyword value
    • DATE-OBS header keyword value
    • DATE-END header keyword value
    • FITS column specified by timefield parameter of all extensions other than GTI extensions
    • START column of GTI extensions
    • STOP column of GTI extensions
    Barycentric times are computed as described in Appendix D. Geocentric times are computed in the same way as barycentric times, except that a geocentric time, tgeo defined by the following formula, is computed in place of tbary in Appendix D.

    tgeo = tobs + rsat · npsr / c

    where rsat is the vector pointing from the center of the Earth to the spacecraft. See Appendix D for definitions of the other quantities in the above formula.
  10. The tool renames the temporal output file to a file name specified by outfile parameter.

gtpulsardb

  1. The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
  2. The tool evaluates psrdbfile parameter, and creates a list of file names from the parameter value. If the value of psrdbfile parameter starts with an at-mark ("@"), the part that follows the at-mark is interpreted as the name of a text file which contains a list of pulsar ephemerides database file names, and the contents of the text file is set to the file name list. Otherwise, the value of psrdbfile parameter is interpreted as a pulsar ephemerides database file name, and set to the first and only element of the file name list.
  3. The tool loads the entire contents of all the pulsar ephemerides database files in the file name list, in the order of their appearance in the file name list.
  4. The tool evaluates filter parameter, and filter the loaded pulsar ephemerides as follows. Then, the tool removes entries in OBSERVERS extension and ALTERNATIVE_NAMES extension that become unnecessary as a result of the filtering.
    • If filter = NAME, the tool reads psrname parameter from the parameter interface of the tool, and filters the loaded pulsar ephemeris data by pulsar name. Spin ephemerides, orbital ephemerides, and ephemeris remarks for the pulsar specified by psrname parameter are sub-selected, and those for other pulsars are removed.
    • If filter = TIME, the tool reads tstart and tstop parameters from the parameter interface of the tool, and filters the loaded pulsar ephemeris data by time. Spin ephemerides whose validity window overlaps with the time interval specified by tstart and tstop parameters are sub-selected, and other spin ephemerides are removed. See Appendix C for an explanation of validity of a spin ephemeris. Also, ephemeris remarks that are effective in (a part of) the time interval are sub-selected, and other ephemeris remarks are removed.
    • If filter = SOLAREPH, the tool reads solareph parameter from the parameter interface of the tool, and filters the loaded pulsar ephemeris data by solar system ephemeris name. Spin ephemerides and orbital ephemerides that are derived using the solar system ephemeris specified by solareph parameter are sub-selected, and those using a different solar system ephemeris are removed.
    • If filter = NONE, the tool does not filter the loaded pulsar ephemeris data.
  5. The tool reads author parameter from the parameter interface of the tool, in order to use its value in the next step for AUTHOR header keyword of the output pulsar ephemerides database.
  6. The tool creates an output pulsar ephemerides database file with a file name specified by outfile parameter, and writes all the remaining contents of the loaded pulsar ephemerides database into the output file.

Chapter 3. Details of common data processing procedures

Despite of differences in functionality and capability of the tools, most of them share common procedures such as pulsar ephemeris loading and photon arrival time corrections. This chapter describes common procedures of data processing by most of the pulsar analysis tools. Not all the steps are performed by all those tools, and not all the tools behave in the same way as others in all steps. Most of the tools, however, share procedures to prepare for their data processing, which this chapter describes in greater detail.

Each section of this chapter starts with a list of the tools that actually perform procedures described in the section. The tool names are listed after "Performed by:" tag at the beginning of each section. Also, it is summarized in the following tabulated matrix, where "Yes" means the step is performed by the tool, and "No" mean it is not. The matrix also outlines which part of the procedures each of the tools performs for your reference.

gtbary gtephem gtpulsardb gtpphase gtophase gtpsearch gtpspec gtptest  
No No No Yes Yes Yes Yes Yes Step 1. Open event files
Yes Yes Yes Yes Yes Yes Yes No Step 2. Determine leap second file name
No No No Yes Yes Yes Yes No Step 3. Interpret user's request on arrival time corrections
No Yes No Yes Yes Yes Yes No Step 4. Load pulsar ephemerides
No No No Yes Yes Yes Yes No Step 5. Initialize arrival time corrections
No No No Yes Yes Yes Yes No Step 6. Summarize ephemeris status

Step 1. Open event files

Performed by: gtpphase, gtophase, gtpsearch, gtpspec, gtptest

The tool evaluates evfile parameter, and creates a list of file names from the parameter value. If the value of evfile parameter starts with an at-mark ("@"), the part that follows the at-mark is interpreted as the name of a text file which contains a list of event file names, and the contents of the text file is set to the file name list. Otherwise, the value of evfile parameter is interpreted as an event file name, and set to the first and only element of the file name list.

Each file in the event file name list is then opened for an event extension and a Good Time Interval (GTI) extension in the order of appearance in the list. The name of an event extension to open is specified by evtable parameter, while a GTI extension name is fixed to "GTI".

When opening those extensions, the tool reads several FITS header keywords to determine a data format of each file. For example, the tool may read TELESCOP and INSTRUME header keywords to identify a mission name and an instrument type, as well as TIMEFRAM and TIMESYS header keywords to judge which arrival time corrections have already been applied to various kinds of times recorded in the file. If the data format is not supported by the tool, it produces an error message and terminates the process.

Step 2. Determine leap second file name

Performed by: gtbary, gtephem, gtpulsardb, gtpphase, gtophase, gtpsearch, gtpspec

The tool reads leapsecfile parameter from the parameter interface of the tool, and takes its value as the name of a leap second file to use in the later stage of the process, unless the value of leapsecfile parameter is DEFAULT. If leapsecfile parameter is set to DEFAULT, the tool looks for a default leap second file, although it is recommended to explicitly set leapsecfile parameter to the name of a leap second file of your choice if your analysis is sensitive to leap seconds, for the following reasons.

  • An insertion or a removal of a leap second is unpredictable by nature, and a delivery of the latest leap second file does not necessarily synchronize with a software release. As a result, a default leap second file can be obsolete for your data set, and may not contain all the leap seconds necessary for your analysis.
  • A default leap second file is a file named leapsec.fits under a directory specified by TIMING_DIR environment variable. However, using the environment variable to control a search for a leap second file could be misleading, because the value of the environment variable can be modified by a wrapper script, for instance, before the flow of execution reaches to the point where the tool reads it.

The tool needs a leap second file when UTC time system is involved in computation. Examples of such computations include, but are not limited to, interpreting a time expression in UTC and computing a time difference between two moments in time, one of which is expressed in UTC. The tool likely performs a UTC-involved computation when timesys or usersys is set to UTC, or a pulsar ephemeris is given as defined in UTC.

Step 3. Interpret user's request on arrival time corrections

Performed by: gtpphase, gtophase, gtpsearch, gtpspec

The tool evaluates tcorrect parameter to interpret a user's request on which arrival time corrections should be applied to photon arrival times prior to computations of physical or statistical quantities from time series data, with an exception of gtophase which assumes that a user requests a pre-determined combination of arrival time corrections.

There are three levels of request for each type of arrival time corrections: suppressed, allowed, or required. If an arrival time correction is requested to be suppressed, the tool will not apply the arrival time correction on photon arrival times. If an arrival time correction is requested to be allowed, the tool searches for ephemeris data for the arrival time correction after loading all given pulsar ephemerides in the later stage of the process. If a set of ephemeris data is found for the arrival time correction, then the arrival time correction will be applied. Otherwise, it will not. If an arrival time correction is requested to be required, the tool searches for ephemeris data for the arrival time correction after loading all given pulsar ephemerides in the later stage of the process. If a set of ephemeris data is found for the arrival time correction, then the arrival time correction will be applied. Otherwise, the tool produces an error message and terminates the process.

The following table summarizes tool's interpretations of user's request on arrival time corrections for all possible values of tcorrect parameter. Due to differences in functionality and purposes of each tool, some tools have a slightly different way to interpret a user's request than others. Also in the table is gtophase, which does not have tcorrect parameter, showing that gtophase behaves as if a user requested arrival time corrections as indicated in the table.

Pulsar tool tcorrect value Barycentric correction Binary demodulation P-dot cancellation
gtpphase NONE Suppressed Suppressed Suppressed
AUTO Allowed Allowed Suppressed
BARY Required Suppressed Suppressed
BIN Required Required Suppressed
ALL Required Required Suppressed
gtophase     Required Required Suppressed
gtpsearch
gtpspec		 NONE Suppressed Suppressed Suppressed
AUTO Allowed Allowed Allowed
BARY Required Suppressed Suppressed
BIN Required Required Suppressed
PDOT Required Suppressed Required
ALL Required Required Required

Note that not all possible combinations of arrival time corrections can be requested, because some combinations are physically inappropriate and should not be allowed. For example, because the barycentric correction must be applied before binary demodulation, a user cannot request to have the barycentric correction suppressed and the binary demodulation required.

Step 4. Load pulsar ephemerides

Performed by: gtephem, gtpphase, gtophase, gtpsearch, gtpspec

The tool reads pulsar ephemerides and related remarks, and stores them in a memory space for later use. There are two types of pulsar ephemeris for the tool to load: a spin ephemeris and an orbital ephemeris. Naively speaking, a spin ephemeris is a set of parameters that describes temporal changes in pulsar's spin frequency, while an orbital ephemeris is a set of parameters that characterize propagation delays in a binary system that a pulsar of interest belongs to. Not all the pulsar analysis tools listed above loads the both types of pulsar ephemeris. An ephemeris remark can be any statement by a provider of pulsar ephemerides who wants to call user's attention, such as a report of a glitch or excessive timing noise during a certain period of time. For chatter = 4 or larger, if any information (spin ephemerides, orbital ephemerides, or ephemeris remarks) is loaded from pulsar ephemerides database file(s), the tool compiles a brief summary of the creation history of the loaded pulsar ephemerides database file(s), and reports the summary for human reading.

Step 4-1. Load spin ephemerides

Performed by: gtephem, gtpphase, gtpsearch

The tool evaluates ephstyle parameter, and loads spin ephemerides as described in the following table. Because of the functionality and the capability of the tool, gtephem does not have ephstyle parameter, and loads spin ephemerides as if ephstyle were set to DB. For each spin ephemeris loaded, the tool chooses an appropriate numerical model of temporal changes in pulsar's spin frequency as described in Appendix A.

Pulsar tool ephstyle value Spin ephemerides are loaded...
gtpphase
gtpsearch		 DB from SPIN_PARAMETERS extension(s) of pulsar ephemerides database file(s) specified by psrdbfile parameter. Spin ephemerides are sub-selected by pulsar name given by psrname and, if requested, by solar system ephemeris name given by solareph parameter, before being stored in a memory space for later use. Sub-selection of spin ephemerides by solar system ephemeris name will be performed only when matchsolareph parameter is set to PSRDB or ALL.
FREQ from parameters ephepoch, timeformat, timesys, ra, dec, phi0, f0, f1, and f2. For a tool without phi0 parameter in the parameter interface, phi0 parameter is assumed to be zero (0.0).
PER from parameters ephepoch, timeformat, timesys, ra, dec, phi0, p0, p1, and p2. For a tool without phi0 parameter in the parameter interface, phi0 parameter is assumed to be zero (0.0).
gtephem     in the same way as gtpphase and gtpsearch do for ephstyle = DB.
Step 4-2. Load orbital ephemerides

Performed by: gtephem, gtpphase, gtophase, gtpsearch, gtpspec

Unless binary demodulation is requested to be suppressed (see Step 3 for types of request for arrival time corrections), orbital ephemerides are loaded from ORBITAL_PARAMETERS extension(s) of pulsar ephemerides database file(s) specified by psrdbfile parameter. Orbital ephemerides are sub-selected by pulsar name given by psrname and, if requested, by solar system ephemeris given by solareph parameter, before being stored in a memory space for later use. Sub-selection of orbital ephemerides by solar system ephemeris name will be performed only when matchsolareph parameter is set to PSRDB or ALL. For each orbital ephemeris loaded, the tool chooses an appropriate numerical model of pulsar's binary orbit as described in Appendix A.

Step 4-3. Load ephemeris remarks

Performed by: gtephem, gtpphase, gtophase, gtpsearch, gtpspec

The tool evaluates reportephstatus parameter, and if the value of reportephstatus parameter is a logical true, then the tool loads ephemeris remarks from REMARKS extension of pulsar ephemerides database file(s) specified by psrdbfile parameter.

Step 5. Initialize arrival time corrections

Performed by: gtpphase, gtophase, gtpsearch, gtpspec

The tool determines which arrival time corrections to apply to photon arrival times on the fly, defines a time axis for time series analysis, and initializes arrival time corrections that are determined to be applied. If chatter parameter is set to 3 or larger, the tool summarizes its determinations on arrival time corrections, and reports the summary for human reading.

Step 5-1. Determine arrival time corrections to apply

Performed by: gtpphase, gtophase, gtpsearch, gtpspec

The tool determines which arrival time corrections to apply to photon arrival times by a combination of a user's request and availability of ephemeris data for requested arrival time corrections. The following table summarizes tool's determination on whether each type of arrival time correction should be applied.

Arrival time correction Pulsar tool User's request The arrival time correction is determined...
binary demodulation gtpphase
gtophase
gtpsearch
gtpspec		 suppressed not to be applied.
allowed to be applied if and only if at least one orbital ephemeris is loaded.
required to be applied if at least one orbital ephemeris is loaded. If none is loaded, the tool produces an error message and terminates the process.
p-dot cancellation gtpsearch suppressed not to be applied.
allowed to be applied if and only if at least one spin ephemeris is loaded.
required to be applied if at least one spin ephemeris is loaded. If none is loaded, the tool produces an error message and terminates the process.
gtpspec suppressed not to be applied.
allowed
required		 to be applied.
barycentric correction gtpphase
gtophase
gtpsearch
gtpspec		 suppressed not to be applied. If binary demodulation or p-dot cancellation are determined to be applied, the tool produces an error message and terminates the process.
allowed
required		 to be applied.
Step 5-2. Define time axis for time series analysis

Performed by: gtpsearch, gtpspec

The tool determines a time origin and a time system for time series analysis and stores in a memory space for later use. There are multiple ways to specify a time origin, and a user can choose a desired one by setting an appropriate value to timeorigin parameter.

timeorigin value A time origin is set to...
START the start time of a given data set. See Appendix B for definition of the start time of a data set.
STOP the stop time of a given data set. See Appendix B for definition of the stop time of a data set.
MIDDLE the mid-point of the start time and the stop time of a given data set. See Appendix B for definition of a start and a stop time of a data set. The mid-point of the two moments in time is computed in a time system specified by TIMESYS header keyword of the given event file, or the first event file in the event file list given through @-mark convention if more than one event files are given.
USER a moment in time that usertime parameter points to as in a time format specified by userformat parameter in a time system specified by usersys parameter.

A time system for time series analysis is determined depending on arrival time corrections (barycentric correction, binary demodulation, and p-dot cancellation) that are requested. If at least one of them is requested, TDB system is chosen for time series analysis. If no arrival time corrections are requested, the tool reads a TIMESYS header keyword value from all the given event files, and checks whether they are all identical to each other. If they are identical to each other, then the time system is chosen for time series analysis. Otherwise, the tool produces an error message and terminates the process.

Step 5-3. Initialize p-dot cancellation

Performed by: gtpsearch, gtpspec

If p-dot cancellation is determined to be applied in the earlier section of this step, (see above for the determination procedure), the tool determines values of the following parameters for p-dot cancellation and a time system in which p-dot cancellation computations are performed, and stores them in a memory space for later use. If p-dot cancellation is determined not to be applied, the tool does nothing in this section.

r1: ratio of the first time-derivative of spin frequency over the spin frequency at the time origin
r2: ratio of the second time-derivative of spin frequency over the spin frequency at the time origin

Due to differences in functionality and capability of the tools, each tool has its own way to determine them.

  • gtpsearch chooses one spin ephemeris for the time origin from the loaded spin ephemerides as described in Appendix C. Validity of a chosen spin ephemeris is not required. See Appendix C for an explanation of validity of a spin ephemeris. Using the chosen spin ephemeris, the tool computes spin frequency and its first and second time-derivatives at the time origin (see Appendix A for how to compute those quantities). The tool then computes a ratio of the first derivative over the spin frequency and another ratio of the second derivative over the spin frequency, and takes them as r1 and r2. Finally, the tool identifies an appropriate time system in which ephemeris computations with the chosen spin ephemeris are performed (see Appendix A for how to determine a time system for ephemeris computations), and takes it as a time system to be used for p-dot cancellation.

  • gtpspec evaluates ephstyle parameter. If ephstyle = FREQ, the tool reads f1f0ratio and f2f0ratio parameters from the parameter interface of the tool, and takes them as r1 and r2, respectively. In other words, r1 and r2 are set as:

    r1 = f1f0ratio
    r2 = f2f0ratio

    where f1f0ratio and f2f0ratio are parameter values taken from the parameter interface of the tool. If ephstyle = PER, the tool reads p1p0ratio and p2p0ratio parameters from the parameter interface of the tool, and computes r1 and r2 as:

    r1 = - p1p0ratio
    r2 = 2 × (p1p0ratio)2 - p2p0ratio

    where p1p0ratio and p2p0ratio are parameter values taken from the parameter interface of the tool. In both cases, the tool performs p-dot cancellation computations in a time system for time series analysis determined in Step 5-2.
Step 5-4. Initialize barycentric correction

Performed by: gtpphase, gtophase, gtpsearch, gtpspec

The tool reads scfile, sctable, and angtol parameters from the parameter interface of the tool, and stores them in a memory space for later use. See Appendix D for how the tool uses those parameter values. The tool also evaluates matchsolareph parameter, and if matchsolareph = EVENT or ALL, then the tool checks all event files in which event times are already barycentric, in order to enforce consistent use of a single solar system ephemeris in barycentric corrections throughout all the given event file(s). For other values of matchsolareph parameter, the tool does not perform the checking. For each event file whose event times are already barycentric, the tool performs the checking as follows.

  • The tool identifies the name of solar system ephemeris with which barycentric times in the event file were computed. Typically, it is recorded as a PLEPHEM header keyword value in the event file.
  • The tool compares it with the value of solareph parameter. If they match, the tool continues processing the event file. Otherwise, the tool produces an error message and terminates the process.

In addition, gtophase and gtpspec reads ra and dec parameters from the parameter interface of the tool, and stores them in a memory space for later use as a direction in the sky for which barycentric corrections will be computed.

Step 5-5. Initialize binary demodulation

Performed by: none

The tool does nothing to initialize binary demodulation.

Step 6. Summarize ephemeris status

Performed by: gtpphase, gtophase, gtpsearch, gtpspec

The tool examines the ephemeris information that are loaded in Step 4 (i.e., spin ephemerides, orbital ephemerides, and ephemeris remarks), and produces a warning message if any noteworthy ephemeris status is identified. A warning message produced in this step is only for human reading, and does not affect data processing by the tool at all. It is designed merely to notify a user availability of ephemeris data throughout anticipated data processing (i.e., continuity of ephemeris coverage without a gap), or to call user's attention to a certain special condition such as a glitch.

The loaded ephemeris information is examined within a time interval derived from a time coverage of given event data. First, the tool computes the start time of a given data set. See Appendix B for definition of the start time of a data set. Then, the start time is applied arrival time corrections that are determined to be applied to photon arrival times in Step 5-1. See Appendix D for how the tool applies arrival time corrections. If it is successful, then the tool starts ephemeris examination at the corrected start time. If it is unsuccessful (which can happen due to lack of ephemeris data for the moment in time, for example), the tool starts ephemeris examination at the uncorrected start time. Similarly, the tool determines a stopping time of ephemeris examination based on the stop time of a given data set.

Note that unsuccessful arrival time corrections in this step do not produce an error message or a warning. That is because unsuccessful arrival time corrections on the start time and the stop time do not necessarily mean arrival time corrections will be unsuccessful in the later process of data.

Also note a warning message from the ephemeris examination performed in this step may not be complete. If arrival time corrections are unsuccessful in computing the start time or the stop time of the examination time interval, ephemeris gap(s) or ephemeris remark(s) may be missed near the uncorrected time(s), for example.

The tool then examines the loaded ephemeris information within the time interval determined above, and reports the result of examination. Not all the tools examine all kinds of ephemeris status. Depending on the tool's functionality, each tool has a set of ephemeris status to examine, as described below.

  • gtpphase examines and reports ephemeris gaps within the time interval and ephemeris remarks that overlaps with the time interval.
  • gtophase, gtpsearch, and gtpspec examines and reports only ephemeris remarks that overlaps with the time interval.

For each ephemeris status reported in a warning message, the tool also reports a time interval during which the ephemeris status persists. Based on a warning message produced in this step, a user may choose to further sub-select events by event time, in order to avoid potential problem(s) in user's analysis caused by the reported ephemeris status.

Appendix

Appendix A. Numerical models of pulsar's timing properties

Pulsar's timing properties reflect two types of pulsar's physical motions: a rotational motion and an orbital motion. In a conventional time series analysis, a rotational motion of a pulsar is characterized by a numerical model that describes temporal changes in pulsar's spin frequency. There are a few models for temporal changes in spin frequency. Each spin-frequency model is parametrized by a set of parameters in its own way, and a spin frequency and its time derivatives at an arbitrary time can be computed with the parameter set. Also, it is common in a conventional time series analysis to characterize a pulsar's orbital motion by a numerical model that describes a propagation time for a photon that is emitted at a pulsar and travels in a binary system that the pulsar belongs to. Various physical and empirical models have been developed for an orbital motion of a pulsar in a binary system. Each binary model is parametrized by a set of parameters in its own way, and a photon propagation time in a binary system can be computed with the parameter set.

This section lists the spin-frequency models and the binary orbital models that are implemented in the pulsar analysis tools, and explains how the pulsar analysis tools compute relevant quantities for time series analysis with given spin and orbital parameters.

Ephemeris computations with spin parameters

A spin-frequency model can be represented by a function of time that computes the number of revolutions that the pulsar will turn until a given point of time since a reference time. When the tool loads spin ephemerides, the tool identifies an appropriate model of spin-frequency history for each spin ephemeris, and chooses a function of time to compute a pulsar position, a pulse phase, a pulse frequency, and time derivatives of pulse frequency whenever necessary in data processing.

Letting α(t) and δ(t) be the pulsar's Right Ascension and Declination, respectively, N(t) the number of pulsar's revolutions at time t, φ(t) a pulse phase at time t, f(t) a pulse frequency at time t, and f(n)(t) the n-th time derivative of pulse frequency at time t, φ(t) is computed as the fractional part of N(t) independent of spin-frequency model. By convention, φ(t) is computed to be greater than or equal to zero (0.0) and less than unity (1.0). Other quantities are computed as described below, depending on a spin-frequency model that the tool chooses for a given spin ephemeris.

  • If a spin ephemeris is loaded from a SPIN_PARAMETERS extension of a pulsar ephemerides database file with header keyword EPHSTYLE = "FREQ", the quantities are computed in TDB time system as:

    α(t) = RA
    δ(t) = DEC
    N(t) = φ0 + F0 × (t - t0) + 1/2 × F1 × (t - t0)2 + 1/6 × F2 × (t - t0)3
    f(t) = F0 + F1 × (t - t0) + 1/2 × F2 × (t - t0)2
    f(1)(t) = F1 + F2 × (t - t0)
    f(2)(t) = F2
    f(n)(t) = 0 for n ≥ 3

    where RA, DEC, F0, F1, and F2 are field values of FITS columns of such names, φ0 a pulse phase constant which is further explained below, t0 a moment in time that the sum of EPOCH_INT and EPOCH_FRAC points to as an MJD number in TDB time system. Pulse phase constant φ0 is the fractional part of N0, which is then computed as

    N0 = - F0 × (ta - t0) - 1/2 × F1 × (ta - t0)2 - 1/6 × F2 × (ta - t0)3,

    where ta is a moment in time that the sum of TOABARY_INT and TOABARY_FRAC points to as an MJD number in TDB time system. Combining the definitions above, one obtains φ(ta) = 0.0 (zero). In other words, a pulse phase at a time of pulse arrival is always zero (0.0), allowing multiple sets of spin parameters to be used to assign pulse phase values to photon arrival times in a consistent manner, to the extent of reliability of ephemeris data used.
  • If a spin ephemeris is loaded from the parameter interface of the tool with ephstyle = FREQ, the quantities are computed in a time system specified by timesys parameter as:

    α(t) = ra
    δ(t) = dec
    N(t) = phi0 + f0 × (t - t0) + 1/2 × f1 × (t - t0)2 + 1/6 × f2 × (t - t0)3
    f(t) = f0 + f1 × (t - t0) + 1/2 × f2 × (t - t0)2
    f(1)(t) = f1 + f2 × (t - t0)
    f(2)(t) = f2
    f(n)(t) = 0 for n ≥ 3

    where ra, dec, phi0, f0, f1, and f2 are parameter values taken from the parameter interface of the tool, t0 a moment in time that ephepoch parameter points to as in a time format specified by timeformat parameter in a time system specified by timesys parameter.
  • If a spin ephemeris is loaded from the parameter interface of the tool with ephstyle = PER, the quantities are computed in a time system specified by timesys parameter as:

    α(t) = ra
    δ(t) = dec
    N(t) = phi0 + I(t - t0, p0, p1, p2)
    f(t) = 1 / q0
    f(1)(t) = - q1 / q02
    f(2)(t) = - q2 / q02 + 2 q12 / q03
    f(n)(t) = Dn(q0, q1, p2)

    where ra, dec, phi0, p0, p1, and p2 are parameter values taken from the parameter interface of the tool, t0 a moment in time that ephepoch parameter points to as in a time format specified by timeformat parameter in a time system specified by timesys parameter, I(t, t0, p0, p1, p2) the definite integral of f(t) over time range from t0 to t, q0 and q1 are a pulse period and its first time derivative at time t as defined by:

    q0 = p0 + p1 × (t - t0) + 1/2 × p2 × (t - t0)2
    q1 = p1 + p2 × (t - t0)

    and finally Dn(q0, q1, p2) the n-th time derivative of f(t) expressed as a polynomial of q0, q1, and p2. The explicit forms of the functions I(t, t0, p0, p1, p2) and Dn(q0, q1, p2) are derived and described in a separate document titled "Computing pulse phase and pulse frequency using pulse period ephemeris" (PDF file, 92 kB) for your reference.
Ephemeris computations with orbital parameters

An orbital model of a binary system can be represented by a function of time that computes a time delay in a photon arrival time caused by an orbital motion of a pulsar in a binary system. When the tool loads orbital ephemerides, the tool identifies an appropriate binary orbital model for each orbital ephemeris, and chooses a function of time to compute a photon propagation time in a binary system. Based on the chosen function, the tool computes a photon arrival time at which a photon of interest would have arrived at the spacecraft if the photon left at the center of gravity of the binary system, instead of the pulsar. This computation, called binary demodulation, is performed as a part of photon arrival time corrections, if so requested. Binary demodulation is an iterative process to search for a corrected photon arrival time (see Appendix D for details), and each orbital model defines the maximum tolerance in time difference as convergence criteria in the iterative search.

  • If an orbital ephemeris is loaded from a ORBITAL_PARAMETERS extension of a pulsar ephemerides database file with header keyword EPHSTYLE = "DD", time delays are computed in TDB time system based on the Damour-Deruelle model (Damour and Deruelle 1985 and 1986) with no aberration (i.e., A = B = 0). The explicit formula of time delay was originally derived by Damour and Deruelle (1986), and can be found as Equation 7 by Taylor and Weisberg (1989), for example, where a brief description of the model and an excellent summary of other binary models are also given. The maximum tolerance in time difference for binary-demodulation is set to 1 nano-second (1 × 10-9 seconds).

    References
    Damour, T., and Deruelle, N. 1985, Ann. Inst. H. Poincare (Physique Theorique), 43, 107.
    Damour, T., and Deruelle, N. 1986, Ann. Inst. H. Poincare (Physique Theorique), 44, 263.
    Taylor, J. H., and Weisberg, J. M. 1989, ApJ 345, 454.

Appendix B. Time boundaries of event data

Time boundaries of a give set of data often plays an important roll in time series analysis. In the pulsar analysis tools, the time boundaries of given event data are defined as described below.

Start time
a moment in time that a START column value of the first FITS row in GTI extension of the given event file points to. If more than one event files are given, the tool computes the start time for each event file, and takes the earliest one among them as the start time for the entire data set.
Stop time
a moment in time that a STOP column value of the last FITS row in GTI extension of the given event file points to. If more than one event files are given, the tool computes the stop time for each event file, and takes the latest one among them as the stop time for the entire data set.

Appendix C. Best ephemeris for a given moment in time

In several occasions in data processing by the pulsar analysis tools, the tool needs to choose a spin ephemeris or an orbital ephemeris that is best used to compute a physical quantity (such as a spin frequency or a photon propagation delay in a binary system) at a given moment in time. Which ephemeris is best among the loaded ephemerides depends on a quantity to compute and a purpose of the computation, as well as a moment in time for which an ephemeris is chosen. This section describes the ephemeris choosing procedures that are implemented in the pulsar analysis tools.

Best spin ephemeris when validity required

Each spin ephemeris comes with a time interval during which temporal changes of spin frequency is well reproduced by the ephemeris. The time interval is called a validity window of the spin ephemeris, and the spin ephemeris is considered valid during the time interval. When validity of a chosen ephemeris is required, the tool chooses a spin ephemeris from the loaded spin ephemerides in the following procedure.

  1. The tool searches for a spin ephemeris which is valid at the given moment in time.
  2. If only one spin ephemeris is found, the tool chooses it as the best spin ephemeris.
  3. If more than one spin ephemerides are found, the tool takes the following steps to compare them until only one spin ephemeris is left. Then, the tool chooses it for the best spin ephemeris.
    1. Ephemeris with a validity window that starts latest
    2. Ephemeris with the longest validity window
    3. Ephemeris that appears latest in the loaded spin ephemerides
  4. If none is found, the tool reports an error.
Best spin ephemeris when validity not required

When validity of a chosen ephemeris is not required, the tool chooses a spin ephemeris from the loaded spin ephemerides in the following procedure.

  1. The tool tries to find the best spin ephemeris which is valid at the given moment in time in the procedure described above.
  2. If found, the tool chooses it as the best spin ephemeris.
  3. If none is found (i.e., the tool reports an error in finding the best spin ephemeris which is valid at the given moment in time), the tool searches for a spin ephemeris that is closest to the given moment in time in the following procedure.
    1. For each of the loaded spin ephemeris, the tool computes in TDB system the time difference between the start time of its validity window and the given moment in time, and the time difference between the stop time of its validity window and the given moment in time.
    2. The tool compares the lengths of the two time differences, and takes the shorter one as the time distance between the spin ephemeris and the given moment in time.
    3. The tool searches for a spin ephemeris with the shortest time distance to the given moment in time, and takes it as the closest spin ephemeris to the given moment in time. If more than one spin ephemerides are tied, the tool takes the one that appears latest in the loaded spin ephemerides.
  4. The tool chooses the closest spin ephemeris to the given moment in time for the best spin ephemeris.
Best orbital ephemeris

Unlike a spin ephemeris, an orbital ephemeris does not come with a validity window, based on which the best spin ephemeris can be defined for an arbitrary moment in time. Instead, the tool defines a reference epoch for each orbital model, and uses it to choose the best orbital ephemeris among the loaded orbital ephemerides as described below. Unless otherwise noted, the tool uses a time of periastron passage as a reference epoch of an orbital ephemeris, since it is commonly seen in orbital parameters of most binary models.

  1. For each of the loaded orbital ephemeris, the tool computes in TDB system the time difference between the reference time of the orbital ephemeris and the given moment in time, and takes it as the time distance between the orbital ephemeris and the given moment in time.
  2. The tool searches for an orbital ephemeris with the shortest time distance to the given moment in time, and takes it as the closest orbital ephemeris to the given moment in time. If more than one orbital ephemerides are tied, the tool takes the one that appears latest in the loaded orbital ephemerides.
  3. The tool chooses the closest orbital ephemeris to the given moment in time for the best orbital ephemeris.

Appendix D. Arrival time corrections

Three arrival time corrections are implemented for the pulsar analysis tools: the barycentric correction, binary demodulation, and p-dot cancellation. When requested, the tool applies those corrections to a given moment in the order listed above. This section describes each of the arrival time corrections in detail.

Barycentric correction

For an event time (or any other representation of time) in an event file, the tool computes a barycentric time, which is a moment in time when the photon would have arrived at the solar system barycenter, as described below.

  1. The tool computes the pulsar position (Right Ascension and Declination) in the sky for a given moment in time. gtpphase and gtpsearch chooses the best spin ephemeris for the given moment in time (see Appendix C for ephemeris choosing procedures), and computes the pulsar position at the given moment in time with the spin ephemeris (see Appendix C for ephemeris computation procedures). gtophase, gtpspec, and gtbary read ra and dec parameters in the parameter interface of the tool, and use them as the sky position for which barycentric corrections are applied.
  2. If it is indicated that the event time is already a barycentric time (typically by TIMEREF header keyword of the event file being SOLARSYSTEM), the tool checks consistency in pulsar position, by comparing the pulsar position for which the barycentric time was computed with the one for which the barycentric corrections would have been applied to the event time if it were not barycentric. If a separation angle between the two pulsar positions is less than or equal to the value of angtol parameter, then the tool determines they coincide with each other. Otherwise, the tool produces an error message and terminates the process.
  3. The tool reads spacecraft data from a table (usually a FITS extension) specified by sctable parameter in a spacecraft file specified by scfile parameter, and computes the spacecraft position in reference to the center of the Earth at the given moment in time.
  4. The tool obtains or computes the following quantities at the given moment in time.
    • robs : the vector pointing from the solar system barycenter to the spacecraft
    • npsr : a unit vector pointing to the pulsar
    • Tsys : the time difference between TT and TDB systems (TDB - TT).
    • rsat : the vector pointing from the geocenter (the center of the Earth) to the spacecraft
    • vearth : the velocity of the Earth with respect to the solar system barycenter
    • θ : the pulsar-Sun-Earth angle, i.e., the angle formed between a line connecting the pulsar and the center of the Sun and a line connecting the Earth and the center of the Sun. Note it is the center of the Sun to use here, not the solar system barycenter.
  5. The tool computes a barycentric time for the given moment in time. Let tobs be the given moment in time, tbary a barycentric time for it, then the tool computes it as:

    tbary = tobs + Δtg + ΔtE - ΔtS

    where Δtg is a geometric time delay in the solar system (i.e., a light travel time without gravitational effects), ΔtE is the Einstein corrections, ΔtS is the Shapiro delay due to the gravitational field of the solar system. Each term on the right-hand side is computed as:

    Δtg = robs · npsr / c
    ΔtE = Tsys + rsat · vearth / c2
    ΔtS = - (2GMsun / c3) log(1 + cosθ)

    where c is the speed of light, G the gravitational constant, and Msun the solar mass.
Binary demodulation

For a barycentric time, the tool computes a binary-demodulated time, which is a moment in time when the photon would have arrived at the solar system barycenter if it were emitted at the center of gravity of the binary system that the pulsar belongs to. The tool performs binary demodulation on the fly if it is determined to be applied in Step 5-1 of Chapter 3. Binary demodulation is an arrival time correction for an orbital delay, i.e., a time delay caused by the fact that the photon was emitted and traveled within a binary system, with respect to the case where the photon were emitted at the center of the gravity of the binary system. Note that computations of orbital delays are dependent on a binary orbital model, and may or may not include gravitational effects depending on a binary orbital model.

As described below, the procedure is iterative, but typically a convergence is achieved in a few iterations.

  1. The tool computes initial guesses on a binary-demodulated time and an orbital delay for the photon.
    1. The tool takes the given barycentric time as an initial guess on a binary-demodulated time.
    2. The tool chooses the best orbital ephemeris for the given barycentric time.
    3. With the chosen orbital ephemeris, the tool computes an orbital delay for a photon emitted at the given barycentric time, and takes it as an initial guess on an orbital delay.
  2. The tool iterates the following procedures until it converges (the convergence criteria is also described below). In the first iteration, the tool takes the initial guess of a binary-demodulated time as a candidate time, and the initial guess of an orbital delay as a tentative orbital delay.
    1. The tool subtracts the tentative orbital delay from the candidate time, and takes the result as a new candidate time.
    2. The tool chooses the best orbital ephemeris for the new candidate time.
    3. With the chosen orbital ephemeris, the tool computes an orbital delay for a photon emitted at the new candidate time, and takes it as a new tentative delay.
    4. The tool adds the new tentative delay to the new candidate time, and compares the result with the given barycentric time. If the time difference is less than or equal to the maximum tolerance in time difference for the chosen binary orbital model (see Appendix A for details), it is deemed converged. Otherwise, take another iteration with the new candidate time and the new tentative orbital delay.
  3. The tool takes the candidate time at the convergence as a binary-demodulated time.
P-dot cancellation

The tool modifies a barycentric time or a binary-demodulated time to cancel or reduce the effect of frequency derivatives by the given amount, in order to make potential pulsations in the time series data appear to have a constant frequency over time.

Letting tunc be a photon arrival time that is to be corrected for p-dot cancellation, tpdot a corrected time, then the tool computes it in a time system determined in Step 5-3 of Chapter 3 as:

tpdot = tunc + 1/2 × r1 × (tunc - torigin)2 + 1/6 × r2 × (tunc - torigin)3

where r1 and r2 are the p-dot cancellation parameters as defined in Step 5-3 of Chapter 3, and torigin is the time origin that is computed in Step 5-2 of Chapter 3.