Improved Pulsar Timing via Principle Component Mode Tracking
Hsiu-Hsien Lin Kiyoshi Masui Ue-Li Pen Jeffrey B. Peterson
astro-ph.IM astro-ph.HE


firstpage–lastpage
The [Y/Mg] clock works for evolved solar metallicity stars
Based on spectroscopic observations made with two telescopes: the Nordic Optical Telescope operated by NOTSA at the Observatorio del Roque de los Muchachos (La Palma, Spain) of the Instituto de Astrofísica de Canarias and the Keck I Telescope at the W.M. Keck Observatory (Mauna Kea, Hawaii, USA) operated by the California Institute of Technology, the University of California and the National Aeronautics and Space Administration.

    D. Slumstrup 1 F. Grundahl 1 K. Brogaard 1,2 A. O. Thygesen 3 P. E. Nissen 1 J. Jessen-Hansen 1 V. Van Eylen 4 
 M. G. Pedersen 5
    Received 3 July 2017 / Accepted 20 July2017 
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================








We present a principal component analysis method which tracks and compensates for short-timescale variability in pulsar profiles, with a goal of improving pulsar timing precision.
We couple this with a fast likelihood technique for determining pulse time of arrival, marginalizing over the principal component amplitudes. This allows accurate estimation of timing errors in the presence of pulsar variability.

We apply the algorithm to the slow pulsar PSR J2139+0040 using an archived set of untargeted raster-scan observations at arbitrary epochs across four years, obtaining an improved timing solution. The method permits accurate pulsar timing in data sets with short contiguous on-source observations, opening opportunities for commensality between pulsar timing and mapping surveys.  








pulsars: general – pulsars: individual: J2139+0040 – methods: data analysis






§ INTRODUCTION


While there are thought to be 100,000 pulsars in the Milky Way <cit.> only about 2600 have been cataloged <cit.>. A new generation of drift-scan interferometric telescopes, including CHIME <cit.>, HIRAX <cit.> and HERA <cit.>, will soon begin recording data with a primary goal of creating cosmological 21-cm intensity maps.  In addition, there are plans to use MeerKAT <cit.> to create intensity maps using a raster scan mode.  Such instruments could be used commensally to search for pulsars, and these instruments will have the collecting area, multi-beam capability, and on-sky integration time to substantially increase the pulsar discovery rate. 

The pulsar search will likely have to work within the observational parameters of the hydrogen surveys. For the drift-scan telescopes, the duration of a contiguous observation will be limited to the time sources take to drift through a beam, which for CHIME will be a few minutes. MeerKAT, with its raster scan observing mode will make shorter duration passes but will pass near a pulsar's coordinates more often. Integration time can be accumulated by combining multiple passes. Such a coherent pulsar search in time-disjointed data has been demonstrated in <cit.>. 

To further develop and test time-disjoint pulsar search algorithms, we used an existing hydrogen-mapping data set, the Green Bank Hydrogen Intensity Mapping (GBTIM) survey <cit.>. This is the data set in which fast radio burst FRB 110523 was discovered <cit.>. The set consist of raster scans, so it contains a time-disjointed set of near-passes of pulsar coordinates.

Our long-term goal is to search for new pulsars, but first we need to
test our algorithms and code on known examples, so we obtained a timing solution for PSR J2139+0040.
The pulsar was discovered by <cit.> at right ascension 21:39:42(16), declination +00:36(5), period 0.312470(3) s, and dispersion measure of 36(7) pc/cm^3. The pulsar is bright—we observe an average flux at 800 MHz of approximately 50 mJy—such that single pulses are readily detectable with the Green Bank Telescope. Our team began studying this pulsar because it lies within two degrees of FRB 110523 and could be used to constrain the Galactic component of the scintillation properties for this region of the sky.

Despite the sporadic scan-mode observations in the data set, we show that a timing solution for PSR J2139+0040 can be obtained, and we present improved timing parameters. We find timing residuals are reduced by use of a novel principal component analysis (PCA) technique to fit the time-variable pulse waveforms. The PCA technique automatically compensates for the rapid mode switching of the pulse shape, allowing a more precise timing solution. We describe this technique in detail since it may be more widely useful in pulsar timing.




§ DATA AND PROCESSING


Here we describe our data and observations as well as the processing of the data into folded pulse profiles.



 §.§ Observations

The data set we used came from an 21-cm intensity mapping survey, in which we raster-scanned four different Wiggle-Z (WZ) fields (1hr, 11hr, 15hr, and 22hr)<cit.> recording spectra from 700-900 MHz with 4096 spectral frequency channels, using integration intervals of 1.024 ms. The full survey is comprised of 680 hours of observations between 2011 and 2015. PSR J2139+0040 is located in the 22hr field, which was observed as part of GBT projects 10B-336 and 14B-339. Scans of this field have the beam crossing close to the pulsar often, but for short periods of time, with the beam typically moving across the sky at a rate of several degrees per minute. The set consists of sporadically-spaced observing sessions of several hours. This provides a sample of pulsar data with a wide range of crossing angles and with baseline drifts in the data due to variable ground spill and sky brightness structure. The cadence of the observing epochs is variable. It is therefore challenging to extract a precise pulsar timing solution from the set, providing a test of our processing techniques.

Our WiggleZ data are stored in blocks of 2048 time samples (PSRFITS records <cit.>). From our bulk data, we select all such blocks where at the block midpoint the telescope boresight is within 0.15^∘ from the published pulsar position. This corresponds roughly to a half width at half maximum of the telescope beam at 800 MHz. The resulting data set includes 1975 seconds of integration time, which we term the WZ data set.

In addition to the raster scan data described above, we obtained a single pointed observation of PSR J2139+0040 with a duration of one hour on MJD 57178. The frequency range of the pointed data is 720-920 MHz with 2048 spectral frequency channels, and an integration interval 8.192e-5 s, grouped into records of 4096 samples. We dub this the pointed data set.

Using the PRESTO[<http://www.cv.nrao.edu/ sransom/presto/>] software package on the pointed data set yields DM = 31.7262 pc cm^-3, which we use for all subsequent analysis except Section <ref> where we perform a full fit for the DM.



 §.§ Preprocessing



Data pre-processing produces an initial calibration of the data, removes flux from  a noise injection source, mitigates radio-frequency interference (RFI), long-time-scale noise, and other sky signals such as point sources moving through the beam. 
To remove the system bandpass response function, we divided the recorded intensity by its time mean, converting to units of the system temperature.


To mitigate RFI and long-time-scale signals we apply a stack of filters to the data consisting of: 1. deleting frequencies where the time variance is anomalously high, 2. subtracting the time mean from the data, 3. subtracting the time-linear component from the data, 4. deleting time samples where the frequency mean is outlying, 5. subtracting the frequency mean from the data, and finally repeating step 1. We use 5σ for all thresholds.

The above filters result in a zero-mean artifact in the folded profiles where the pulse profile is pushed negative just outside the main pulse. This is visible in the middle panel of Figure <ref>. However, this artifact is a result of linear operations on the data and as such does not bias our subsequent timing analysis.



 §.§ Barycentric Time and Folding


We choose to time the pulsar in the Barycentric coordinate system because we are working toward searching such data sets for weaker pulsars. This step is not critical to the mode-tracking technique presented in Section <ref>, but we describe it for completeness. We converted time stamps from Universal Time to Barycentric Time, τ(t, α, δ), defined as the time a signal from a far-off source at right ascension α and declination δ arriving at the observatory (here GBT) at time t, would arrive at the solar system barycenter.



We perform this conversion using the TEMPO[<http://tempo.sourceforge.net/>] software package as invoked by the bary command in PRESTO. For the right ascension and declination we initially use the previously published location of the pulsar, which we later refine. As described in Section <ref>, refining this position requires the partial derivatives of the Barycentric times with respect to right ascension and declination, which we calculate by finite difference.

To form pulse profiles, we group our sporadically spaced WZ data into groups spanning five minutes or less. These are typically very sparsely sampled with most groups containing only a few seconds of on-target data. The five-minute duration was selected to reduce the set to a manageable size of fixed cadence, with reasonable signal-to-noise for each folded waveform, and to ensure little worry of phase drift within a group.
The pointed data, in contrast, is contiguously sampled and we divide it into 640 groups of duration 5.369 s, which we use to characterize the time variability, and later stack into 40 pulse profiles of duration 86 s used for timing.

We fold the time stream on the pulsar period 0.312470 s <cit.>.
We use 200 phase bins for the WZ data, and 800 phase bins for pointed data.

We dedisperse the folded data, and then average it over frequency yielding pulse profiles as shown in Figure <ref>. This yields 289 such profiles although a number of these are noise dominated and show no evidence of the pulsar. We discard these, resulting in 232 pulse profiles in the WZ data set. Combined with the 40 profiles from the pointed data set, there are a total of 272 pulse profiles used for the timing analysis. Figure <ref> below shows how these profiles are distributed in time.







§ ANALYSIS AND RESULTS


Here we describe the analysis of the folded pulse profiles. Our analysis accounts for the pulsar's short-time-scale variability, estimates the pulsar phase, and fits a timing solution.



 §.§ Principal component analysis of pulse shapes


Examining the set of folded pulse profiles, we find substantial variation in shape. In particular, some profiles have a stronger first peak, while others have a stronger second peak as shown in Figure <ref>.
This is the dominant mode of time variability for the profiles, with the preference for one peak over the other appearing to be correlated over timescales of several minutes.
This is thus a form of mode switching, which is common for a substantial fraction of pulsars on a range of time scales <cit.>.
To find the best-fit arrival time for each folded waveform, one could simply fit each waveform to an average profile, but waveform shape variations would then cause timing errors, since the fit would favor earlier arrival times when the first peak is stronger. To address this, we use a principal component analysis of the pulse profiles (PCA), which automatically tracks the pulsar mode. This PCA Mode Tracking compensates for mode-switching waveform-shape variation.



We use the 640 pulse profiles from the pointed data (each containing 5.369 s of integration) to construct the PCA. Using Fourier techniques to achieve sub-bin alignment, we align these profiles according to a preliminary, linear timing solution valid for the pointed data. We create a matrix of 640 pulse profiles called d_E i where E is the epoch, an index of the mid-point time
of each folding interval, and i is the 800-point phase bin within the waveform. We carry out a Singular Value Decomposition:

    d_Ei =∑_m U_EmS_mmV_mi,
 
where for each mode m: U_Em is the eigenfunction in epoch E, V_mi the eigenfunction in phase i, and S_mm is the singular value. We show V_mi for the first ten m modes in Figure <ref>. 



Inspecting the phase eigenfunctions, the first mode V_0 i appears similar to the average waveform. We tested this equivalence by averaging the 640 pulse profiles into single template. We compared this average template with the mode V_0 i of our PCA technique, and found the two consistent within the noise.

The second mode V_1i allows the two peaks in the waveform to depart in amplitude from the average. Essentially, U_E1 tracks the variation with epoch of the waveform mode, changing sign when the mode switches. 


Since we intend to use the V modes as templates for fitting pulsar phases, a concern is that noise in these modes will bias the phase measurements. To reduce this noise for the modes above the first, we set them to zero outside of the main pulse (all phase bins except the central 160 of 800) since we see no evidence for pulsar flux outside this region. For the primary mode, the filtering artifact mentioned in Section <ref> results in a smooth structure outside this region and so we spline the profile outside the central 160 phase bins.




 §.§ Pulsar time of arrival estimation from direct integration of the likelihood

The natural next step would be to simultaneously fit mode amplitudes and a pulse phase to each of our individual pulse profiles, using V_mi as templates. However, we find that for this high-dimensional parameter space, the likelihood often has multiple maxima, and so the traditional least-squares fitting method fails. The multi-modal nature of the parameter space is due to the fact that for some profiles, there are multiple combinations of the V_mi modes that, with different phases, may adequately describe the profiles such as those shown in Figure <ref>. To deal with this, we employ a new pulse profile fitting technique that fully samples the parameter space.

We denote an individual measured folded pulsar profile d_i where the index i runs over the phase bins (we suppress the index E for epoch in this section, since the analysis is performed independently for each profile). We model the measured profile as

    d_i = ∑_m A_m V_mi(ϕ) + n_i,

where V_mi (ϕ) are the profile template modes (with m running over mode number), A_m is the mode amplitude, ϕ is the finely adjusted pulsar phase, and n_i is the noise contribution. The dependence of V_mi on ϕ describes the rotation of the templates to match the data. For notational brevity, we now switch to vector notation, where the above equation becomes:

    d =  A^T  V(ϕ) +  n.

We assume Gaussian noise that is uniform and uncorrelated: ⟨ n_i n_j ⟩ = δ_ijσ^2.
The template modes have already been measured so the parameters of the model are ϕ and the amplitudes, A.

Of these parameters we are chiefly interested in ϕ. All the information about this parameter is in p(ϕ| d), the posterior probability distribution of ϕ given the measurements d. For instance, the measurement mean and variance of the phase are the first and second moments of this distribution. This can be calculated from the posterior of all the parameters marginalized over the mode amplitudes:

    p(ϕ| d) = ∫ d^N A  p(ϕ, A| d),

where N is the length of the vector A.


For flat priors on the initial parameters, Bayes' Theorem states that the posterior distribution for the parameters is proportional to the likelihood function, p(ϕ, A| d) ∝ p( d | ϕ, A). Since the data is Gaussian-distributed, the latter is given by:

    p( d | ϕ, A)    ∝exp[
            -1/2χ^2( d, ϕ, A)]
        
    χ^2( d, ϕ, A)    = 
        1/σ^2∑_i [
            d_i - ∑_m  A_m V_mi(ϕ)
        ]^2.


Combining the above, we have:

    p(ϕ| d) ∝∫ d^N A exp{-1/2σ^2∑_i [ d_i - ∑_m A_m  V_mi(ϕ) ]^2 }

This is a multidimensional integral over the vector space of A, which would be prohibitively expensive to do numerically. A key insight is that since the model is linear in A, the likelihood is Gaussian not only in the data, but in A as well (but not in ϕ). This permits the integral to be done analytically.
This same insight was used in <cit.> where a similar integral appears over the frequency dependant template amplitudes in wide-band pulsar timing. However, while in wide-band timing this is a computational convenience, here being able to perform this integral is essential due to the multi-modal nature of the likelihood space. Similarly, the integral appears in <cit.> when marginalizing over epoch dependant profile amplitudes in profile domain timing analysis.


When performing the integral, there is a factor that depends on the expression

    ∑_i V_mi(ϕ) V_ni(ϕ).

This expression is independent of ϕ because the dependence of V_mi on ϕ is just a shift of the template modes, which disappears after summing over phase bins. Thus, there is no need to carry out this part of the calculation. Note that this is only true if the noise is uniform.
The final expression is

    p(ϕ| d) 
           ∝exp[
            -1/2σ^2∑_m ∑_i V_mi(ϕ)  d_i ∑_j V_mj(ϕ) d_j,
        ]

or equivalently

    p(ϕ| d)
           ∝exp[
            -1/2σ^2 
            ( V  d)^T
             V  d
        ].


which is proportional to

    p(ϕ| d)
           ∝exp{-1/2χ^2[ d, ϕ,Â( d, ϕ)]},

where Â( d, ϕ) is the linear-best-fit value for the template amplitudes at fixed ϕ. Linear regression gives

    Â( d, ϕ) = [ V(ϕ)  V^T(ϕ)]^-1 V(ϕ)  d.




So succinctly, the procedure for estimating the time of arrival is:

    
  * At fixed ϕ perform a linear fit for the mode amplitudes
    A, which we found in Section <ref>.
    
  * Evaluate χ^2 for these parameters and take e^-χ^2/2 as the likelihood. 
    
  * Repeat for a range of ϕ values, covering the region of low χ^2 that dominates the likelihood. (We repeat this step for a range of -20 to +20 phase bins from the fixed ϕ, with steps of 0.02 phase bins.)
    
  * Take the zeroth, first, and second moments to get the normalization, phase estimate, and variance, respectively. In Section <ref>, we use the phase estimate and variance to find the timing solution.


While our procedure has been described in the profile space, our actual analysis is performed in the Fourier domain, as in <cit.>. We use only the first through 90th harmonics to limit contamination from noise in the template modes and observed extraneous noise near 300 Hz. We use N=6 template modes in our fits. An example of a profile fit is shown in Figure <ref>.

Finally, properly estimating the error on ϕ requires an estimate of the noise power σ^2. We use the value of σ^2 that results in a reduced chi-squared of unity (χ^2/ν, where ν is the number of degrees of freedom) for the best fit (maximum likelihood) parameters.



 §.§ Timing solution



We proceed to adjust the parameters of the pulsar timing model:

    ϕ_E =
    ϕ_0
    + Δ P/P^2(τ_E - τ_0)
    + Ṗ/2 P^2(τ_E - τ_0)^2
    + 1/P[-∂τ_E/∂αΔα
          - ∂τ_E/∂δΔδ]

where, ϕ_E is the barycentric phase of the pulse profile indexed by E, τ_E is its barycentric time, τ_0 is the reference epoch, P is the period of the pulsar used for folding, ∂τ_E/∂α is the derivative of the barycentric time with respect to source right ascension, and ∂τ_E/∂δ is derivative of the barycentric time with respect to source declination. There are also five free parameters in the timing model, including the period derivative Ṗ, period correction Δ P, initial phase offset ϕ_0, right ascension correction Δα, and declination correction Δδ. We proceed to adjust these parameters such that the above equation fits our phase measurements using standard weighted least squares and extract best fit parameters with uncertainties.

One challenge to obtaining a timing solution with this data set is the large gap in our data. We have phase measurements for epochs over several months in 2011 and 2015 but none in between. Compared to more uniformly sampled data, this distribution of measurements weakly constrains the total number of pulsar rotations in the gap period. Indeed we find multiple χ^2 minima corresponding to changing the parameter Δ P by integer multiples of ∼ 8× 10^-10 s, which changes the number of rotations in the gap. However, the lowest of these minima is δχ^2=44.5 smaller than the next lowest. Thus this minimum is strongly preferred over the others, and our timing parameters are well constrained.

Inspecting the timing residuals, we find one phase measurement that is an extreme, 9σ, outlier. Inspecting the full likelihood curve for the profile phase fit (as described in Section <ref>), we find the likelihood to be very non-Guassian, making such an outlier roughly 0.1% likely. This is not terribly improbable given that we have 272 such data points, and is far more likely than a 9σ outlier for Gaussian data. Since our least-squares fit to the timing solution implicitly assumes Gaussian-distributed phase errors, we exclude this outlying data point in the timing solution fit.

The timing residuals are shown with the fitting results in Figure <ref>.
The fit has a reduced chi-squared of χ^2/ν=1.26 for 266 degrees of freedom, which is reasonable given the non-Gaussian nature our phase measurements which we have represented with 1-σ error bars. We inflate the uncertainties on the timing parameters by this factor.



The parameters for the final timing solution for PSR J2139+0040 are given in Table <ref>.  We find new right ascension and declination, period, period derivative, reference time of arrival (derived from ϕ_0), and also provide the epoch used for the timing model. We use the period and period derivative to estimate the magnetic field strength (B), the characteristic age (t_c), and the spin-down luminosity (Ė)<cit.>. 





To check for evidence of proper motion of the pulsar's sky position, we fix all parameters other than the position, then separately fit the 2011 and 2015 data. We find the two output positions are consistent with each other and the solution using the full data set, and as such we find no evidence for proper motion.

We used the TEMPO pulsar timing software to verify our timing solution. For this we converted our pairs of epoch-phase measurements to barycentric times of arrival then inverted PRESTO's bary command to convert to topocentric times of arrival. Feeding these and our timing solution parameters into TEMPO we find agreement in the value of χ^2. Letting TEMPO refit our timing-model parameters we find no significant shifts. 




 §.§ Dispersion Measure


We used the pointed data to determine DM. We align the 40 folded profiles in phase based on our timing solution, then stacked them into a single profile. 
We fit this profile for the dispersion using the first V mode as a template according to: 

    d_if   = A(f/800  MHz)^βV_0i(ϕ_f)
    ϕ_f    ≡ϕ_900 +  DM/P4148.808  s MHz^2/ pc / cm^3[1/f^2-1/(900  MHz)^2]

Here, d_if is the frequency-dependant profile data, which are fit with parameters for the amplitude A, power-law slope β, and frequency dependant phase ϕ_f for each frequency f. The frequency dependant phase ϕ_f contains an offset ϕ_900 and the dispersion phase delay.

 We subsequently calculate the dispersion-measure-based distance using Cordes-Lazio NE2001 Galactic Free Electron Density Model <cit.> website[<https://www.nrl.navy.mil/rsd/RORF/ne2001>] using the new right ascension, declination, and DM.
 The results of this analysis are given in Table <ref>




 §.§ Timing improvement available from PCA Mode Tracking



To assess the impact of PCA mode tracking we compared the timing residuals using long and short integrations. We carried out these parallel analyses on sets of 16 individual 5.4 s profiles from the pointed data. First, we stacked these profiles into a single pulse profile with a total integration time of 86 s and fit this profile with the technique described in Section <ref>. Second, we fit the 16 profiles individually, then combined their resulting posterior distributions. We did this for the 40 sets of 16 profiles covering the full hour of pointed data. The first technique is similar to the conventional technique of taking long averages to smooth out pulse-to-pulse variability. The second technique uses fine grained data allowing the PCA algorithm to track and compensate for the variability.

We find that using PCA mode tracking and fitting fine-grained 5.4 second averages results in a ∼20% smaller uncertainty on average in the pulse time of arrival (calculated from the second moment of the posterior as described in Section <ref>). While the improvement is substantial for this pulsar, it may not be generic. Pulsars are highly individual: they vary substantially in the degree and time scales mode switching, so we anticipate the improvement available from PCA mode tracking will also vary strongly from case to case.




§ DISCUSSION AND CONCLUSIONS


Pulsar timing observations often make use of profiles averaged for minutes. These long averages smooth out the effect of rapid mode switching. Since it is in general sub-optimal to average over a time-variable signal, we suggest that PCA Mode Tracking could be applied broadly, by intentionally using short averages that display rapid mode switching—which the PCA then follows and compensates for.  
In particular, it may be useful to test this technique using data from millisecond pulsars, since improvement of timing precision of these pulsars allows more precise tests of general relativity, along with tighter constraints on long-wavelength gravitational wave backgrounds. 

We have obtained a substantially improved timing solution for PSR J2139+0040, using untargeted observations in conjunction with a ∼1 hour duration targeted observation. Our data are well-described by a simple five-parameter model,  including its sky location, reference phase, period, and period derivative. The measured slow-down rate of Ṗ = (7.64±0.13)×10^-18(s s^-1) is well below average, with ∼ 99% of slow pulsars having a higher rate (<cit.>) as shown in Figure <ref>.  This indicates that PSR J2139+0040 has an abnormally weak magnetic field.



For many of our pulse profiles, we find a highly non-Gaussian likelihood with multiple maxima for the phase. However, by analytically integrating over the template amplitudes, we are able to fully sample the marginalized posterior for phase and calculate the distribution's mean and standard deviation. The subsequent fit to a timing model achieves a reduced chi-squared of 1.26, which is reasonable given the non-Gaussian nature of our phase errors. A more complete treatment would be to use the phase measurement posteriors directly when fitting a timing solution, rather than reducing them to their first and second moments. This would optimally extract the timing information from the phase measurements but would substantially complicate fitting a timing solution, since least-squares would no longer be applicable.


Our PCA mode tracking method addresses issues similar to those addressed via  "profile-domain pulsar timing" <cit.>, particularly, what these authors refer to as “low-frequency stochasticity” and “phase-correlated stochasticity”. The profile-domain techniques uses pre-defined uncorrelated shapelet components to contribute to the individual profiles.
In contrast, the PCA mode tracking method is non-parametric. There is no need to guess the shaplets, the PCA finds them automatically. Compared to the profile-domain strategy we anticipate that fewer components will be needed, yielding less degeneracy with the pulse phase and likely reducing timing error.  The PCA mode tracking technique derives its mode templates from the data, and so for weak pulsars these waveforms may have noise that is absent when using the pre-defined templates of the profile domain technique. 

A PCA technique was employed to characterize pulsar variability in <cit.>, where, rather than fitting for the time of arrival and component amplitudes in high cadence data, the authors used the empirical correlation between the timing residuals and component amplitudes to retroactively correct for systematic timing error. 

Our results demonstrate the feasibility of timing pulsars commensally with mapping surveys. In particular, upcoming hydrogen intensity mapping surveys using MeerKAT <cit.> and Phase 1 of the Square Kilometre Array <cit.> will map large fractions of the southern sky with thousands of hours of telescope time. The telescopes will operate as a collection of single dishes and will thus need to scan rapidly to overcome 1/f noise. There is the potential to obtain pulsar timing solutions for free (in terms of telescope time) using the data from these surveys in much the same way as we have done here. The challenge is in storing the data at the rapid cadence required for pulsar studies, but, in single dish mode, data volumes are modest compared to interferometric mode. Mapping surveys spend a small fraction of their time pointing at known pulsars, however the mapping survey gets this data at many epochs, for every pulsar in its survey field, providing a large volume of timing data.

One of our goals in this work is to enable pulsar searches commensal with intensity mapping experiments. Having obtained a precise timing solution for a known pulsar in the GBTIM data, we have demonstrated that accumulation of pulsar data over five years in few-pulse snippets can be accomplished, although this may be more difficult for weaker or more erratic pulsars. 

This analysis was possible because the GBTIM data used short 1 ms integrations. The intensity mapping data from CHIME <cit.> and HIRAX <cit.> will use longer integrations for their intensity mapping data, but these instrument will have additional transient-search backend hardware allowing high cadence analysis of the data streams. New search algorithms have recently been proposed by <cit.> to reduce the cost of such long-term pulsar waveform assembly by several orders of magnitude. If these techniques could be employed at upcoming transit survey instruments such as CHIME and HIRAX the pay off could be a substantially increased rate of pulsar discovery. 



§ ACKNOWLEDGEMENTS


We thank Kendrick Smith, Ingrid Stairs, Maura McLaughlin and Alexander Roman for valuable discussions.

K. W. M. is supported by the Canadian Institute for Theoretical Astrophysics National Fellows program.

U.-L. P. acknowledges support from the Natural Sciences and Engineering Research Council of Canada.
Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada  and by the Province of Ontario through the Ministry of Research, Innovation and Science.

J. B. P. acknowledges support from NSF Award 1211777. 

Computations were performed on the GPC supercomputer at the SciNet HPC Consortium.



mnras





