Improved Pulsar Timing via Principle Component Mode Tracking
Hsiu-Hsien Lin, Kiyoshi Masui, Ue-Li Pen, Jeffrey B. Peterson
MMNRAS , 1–9 (2017) Preprint 14 November 2018 Compiled using MNRAS L A TEX style file v3.0
Improved Pulsar Timing via Principle Component ModeTracking
Hsiu-Hsien Lin, (cid:63) Kiyoshi Masui, Ue-Li Pen, , , , Jeffrey B. Peterson McWilliams Center for Cosmology, Carnegie Mellon University, Department of Physics, 5000 Forbes Ave, Pittsburgh, PA, 15213,USA Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Rd, Vancouver, BC, V6T 1Z1, Canada Canadian Institute of Theoretical Astrophysics, 60 St George St, Toronto, ON M5S 3H8, Canada Canadian Institute for Advanced Research, CIFAR Program in Cosmology and Gravity, Toronto, ON, M5G 1Z8 Dunlap Institute for Astronomy & Astrophysics, University of Toronto, 50 St George St, Toronto, ON, M5S 3H4, Canada Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
14 November 2018
ABSTRACT
We present a principal component analysis method which tracks and compensatesfor short-timescale variability in pulsar profiles, with a goal of improving pulsar timingprecision. We couple this with a fast likelihood technique for determining pulse time ofarrival, marginalizing over the principal component amplitudes. This allows accurateestimation of timing errors in the presence of pulsar variability.We apply the algorithm to the slow pulsar PSR J2139+0040 using an archived setof untargeted raster-scan observations at arbitrary epochs across four years, obtainingan improved timing solution. The method permits accurate pulsar timing in data setswith short contiguous on-source observations, opening opportunities for commensalitybetween pulsar timing and mapping surveys.
Key words: pulsars: general – pulsars: individual: J2139+0040 – methods: dataanalysis
While there are thought to be 100,000 pulsars in the MilkyWay (Yusifov & K¨u¸c¨uk 2004) only about 2600 have been cat-aloged (Manchester et al. 2005). A new generation of drift-scan interferometric telescopes, including CHIME (Banduraet al. 2014), HIRAX (Newburgh et al. 2016) and HERA(DeBoer et al. 2017), will soon begin recording data with aprimary goal of creating cosmological 21-cm intensity maps.In addition, there are plans to use MeerKAT (Fonseca et al.2017) to create intensity maps using a raster scan mode.Such instruments could be used commensally to search forpulsars, and these instruments will have the collecting area,multi-beam capability, and on-sky integration time to sub-stantially increase the pulsar discovery rate.The pulsar search will likely have to work within theobservational parameters of the hydrogen surveys. For thedrift-scan telescopes, the duration of a contiguous observa-tion will be limited to the time sources take to drift througha beam, which for CHIME will be a few minutes. MeerKAT,with its raster scan observing mode will make shorter dura- (cid:63)
E-mail: [email protected] tion passes but will pass near a pulsar’s coordinates more of-ten. Integration time can be accumulated by combining mul-tiple passes. Such a coherent pulsar search in time-disjointeddata has been demonstrated in Anderson (1993, Chapter 4).To further develop and test time-disjoint pulsar searchalgorithms, we used an existing hydrogen-mapping dataset, the Green Bank Hydrogen Intensity Mapping (GB-TIM) survey (Chang et al. 2010; Masui et al. 2013; Switzeret al. 2013). This is the data set in which fast radio burstFRB 110523 was discovered (Masui et al. 2015). The setconsist of raster scans, so it contains a time-disjointed set ofnear-passes of pulsar coordinates.Our long-term goal is to search for new pulsars, but firstwe need to test our algorithms and code on known exam-ples, so we obtained a timing solution for PSR J2139+0040.The pulsar was discovered by Camilo et al. (1996) atright ascension 21:39:42(16), declination +00:36(5), period0.312470(3) s, and dispersion measure of 36(7) pc/cm . Thepulsar is bright—we observe an average flux at 800 MHz ofapproximately 50 mJy—such that single pulses are readilydetectable with the Green Bank Telescope. Our team be-gan studying this pulsar because it lies within two degreesof FRB 110523 and could be used to constrain the Galactic © a r X i v : . [ a s t r o - ph . I M ] D ec H.-H. Lin et al. component of the scintillation properties for this region ofthe sky.Despite the sporadic scan-mode observations in the dataset, we show that a timing solution for PSR J2139+0040 canbe obtained, and we present improved timing parameters.We find timing residuals are reduced by use of a novel prin-cipal component analysis (PCA) technique to fit the time-variable pulse waveforms. The PCA technique automaticallycompensates for the rapid mode switching of the pulse shape,allowing a more precise timing solution. We describe thistechnique in detail since it may be more widely useful inpulsar timing.
Here we describe our data and observations as well as theprocessing of the data into folded pulse profiles.
The data set we used came from an 21-cm intensity mappingsurvey, in which we raster-scanned four different Wiggle-Z (WZ) fields (1hr, 11hr, 15hr, and 22hr)(Glazebrooket al. 2007) recording spectra from 700-900 MHz with 4096spectral frequency channels, using integration intervals of1.024 ms. The full survey is comprised of 680 hours of ob-servations between 2011 and 2015. PSR J2139+0040 is lo-cated in the 22hr field, which was observed as part of GBTprojects 10B-336 and 14B-339. Scans of this field have thebeam crossing close to the pulsar often, but for short pe-riods of time, with the beam typically moving across thesky at a rate of several degrees per minute. The set consistsof sporadically-spaced observing sessions of several hours.This provides a sample of pulsar data with a wide rangeof crossing angles and with baseline drifts in the data dueto variable ground spill and sky brightness structure. Thecadence of the observing epochs is variable. It is thereforechallenging to extract a precise pulsar timing solution fromthe set, providing a test of our processing techniques.Our WiggleZ data are stored in blocks of 2048 time sam-ples (
PSRFITS records (Hotan et al. 2004)). From our bulkdata, we select all such blocks where at the block midpointthe telescope boresight is within 0.15 ◦ from the publishedpulsar position. This corresponds roughly to a half width athalf maximum of the telescope beam at 800 MHz. The re-sulting data set includes 1975 seconds of integration time,which we term the WZ data set.In addition to the raster scan data described above, weobtained a single pointed observation of PSR J2139+0040with a duration of one hour on MJD 57178. The frequencyrange of the pointed data is 720-920 MHz with 2048 spec-tral frequency channels, and an integration interval 8.192e-5 s, grouped into records of 4096 samples. We dub this thepointed data set.Using the PRESTO software package on the pointed dataset yields DM = 31.7262 pc cm − , which we use for all sub-sequent analysis except Section 3.4 where we perform a fullfit for the DM. Data pre-processing produces an initial calibration of thedata, removes flux from a noise injection source, mitigatesradio-frequency interference (RFI), long-time-scale noise,and other sky signals such as point sources moving throughthe beam. To remove the system bandpass response func-tion, 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 astack of filters to the data consisting of: 1. deleting frequen-cies where the time variance is anomalously high, 2. sub-tracting the time mean from the data, 3. subtracting thetime-linear component from the data, 4. deleting time sam-ples where the frequency mean is outlying, 5. subtracting thefrequency mean from the data, and finally repeating step 1.We use σ for all thresholds.The above filters result in a zero-mean artifact in thefolded profiles where the pulse profile is pushed negative justoutside the main pulse. This is visible in the middle panelof Figure 1. However, this artifact is a result of linear oper-ations on the data and as such does not bias our subsequenttiming analysis. We choose to time the pulsar in the Barycentric coordinatesystem because we are working toward searching such datasets for weaker pulsars. This step is not critical to the mode-tracking technique presented in Section 3.1, but we describeit for completeness. We converted time stamps from Univer-sal Time to Barycentric Time, τ ( t , α, δ ) , defined as the timea signal from a far-off source at right ascension α and decli-nation δ arriving at the observatory (here GBT) at time t ,would arrive at the solar system barycenter.We perform this conversion using the TEMPO softwarepackage as invoked by the bary command in PRESTO . For theright ascension and declination we initially use the previ-ously published location of the pulsar, which we later refine.As described in Section 3.3, refining this position requiresthe partial derivatives of the Barycentric times with respectto right ascension and declination, which we calculate byfinite difference.To form pulse profiles, we group our sporadically spacedWZ data into groups spanning five minutes or less. These aretypically very sparsely sampled with most groups containingonly a few seconds of on-target data. The five-minute dura-tion was selected to reduce the set to a manageable size offixed cadence, with reasonable signal-to-noise for each foldedwaveform, and to ensure little worry of phase drift within agroup. The pointed data, in contrast, is contiguously sam-pled and we divide it into 640 groups of duration 5.369 s,which we use to characterize the time variability, and laterstack into 40 pulse profiles of duration 86 s used for timing.We fold the time stream on the pulsar period 0.312470 s(Camilo et al. 1996). 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 overfrequency yielding pulse profiles as shown in Figure 1. This http://tempo.sourceforge.net/ MNRAS , 1–9 (2017) mproved Pulsar Timing Figure 1. Typical Pulse profile for PSR J2139+0040.
Up-per panel: Dedispersed dynamic spectra. Middle panel: Blue dotsare the pulse profile averaged across frequency, and red line is thebest-fit model. Barycentric Arrival Times are used in these plots.Bottom panel: Fit residuals. yields 289 such profiles although a number of these are noisedominated and show no evidence of the pulsar. We discardthese, 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 timinganalysis. Figure 5 below shows how these profiles are dis-tributed in time.
Here we describe the analysis of the folded pulse profiles. Ouranalysis accounts for the pulsar’s short-time-scale variability,estimates the pulsar phase, and fits a timing solution.
Examining the set of folded pulse profiles, we find substan-tial variation in shape. In particular, some profiles have astronger first peak, while others have a stronger second peakas shown in Figure 2. This is the dominant mode of timevariability for the profiles, with the preference for one peak
Figure 2. Pulse shape variations.
The two waveforms showfolded data from the WZ set, focusing on the phase range -0.2 to0.2. The pulsar shows mode-switching behavior: The ratio of thefluence between the first and second peak varies substantially. Ifthese two waveforms were naively fit to an average profile, therewould be a substantial timing error because of this profile varia-tion. over the other appearing to be correlated over timescalesof several minutes. This is thus a form of mode switching,which is common for a substantial fraction of pulsars ona range of time scales (Lorimer & Kramer 2004). To findthe best-fit arrival time for each folded waveform, one couldsimply fit each waveform to an average profile, but wave-form shape variations would then cause timing errors, sincethe fit would favor earlier arrival times when the first peakis stronger. To address this, we use a principal componentanalysis of the pulse profiles (PCA), which automaticallytracks the pulsar mode. This PCA Mode Tracking compen-sates 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 thePCA. Using Fourier techniques to achieve sub-bin align-ment, we align these profiles according to a preliminary, lin-ear timing solution valid for the pointed data. We create amatrix of 640 pulse profiles called d Ei 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 carryout a Singular Value Decomposition: d Ei = (cid:213) m U Em S mm V mi , (1)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 singularvalue. We show V mi for the first ten m modes in Figure 3.Inspecting the phase eigenfunctions, the first mode V i appears similar to the average waveform. We tested thisequivalence by averaging the 640 pulse profiles into singletemplate. We compared this average template with the mode V i of our PCA technique, and found the two consistentwithin the noise.The second mode V i allows the two peaks in the wave-form to depart in amplitude from the average. Essentially, U E tracks the variation with epoch of the waveform mode,changing sign when the mode switches. MNRAS000
The two waveforms showfolded data from the WZ set, focusing on the phase range -0.2 to0.2. The pulsar shows mode-switching behavior: The ratio of thefluence between the first and second peak varies substantially. Ifthese two waveforms were naively fit to an average profile, therewould be a substantial timing error because of this profile varia-tion. over the other appearing to be correlated over timescalesof several minutes. This is thus a form of mode switching,which is common for a substantial fraction of pulsars ona range of time scales (Lorimer & Kramer 2004). To findthe best-fit arrival time for each folded waveform, one couldsimply fit each waveform to an average profile, but wave-form shape variations would then cause timing errors, sincethe fit would favor earlier arrival times when the first peakis stronger. To address this, we use a principal componentanalysis of the pulse profiles (PCA), which automaticallytracks the pulsar mode. This PCA Mode Tracking compen-sates 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 thePCA. Using Fourier techniques to achieve sub-bin align-ment, we align these profiles according to a preliminary, lin-ear timing solution valid for the pointed data. We create amatrix of 640 pulse profiles called d Ei 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 carryout a Singular Value Decomposition: d Ei = (cid:213) m U Em S mm V mi , (1)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 singularvalue. We show V mi for the first ten m modes in Figure 3.Inspecting the phase eigenfunctions, the first mode V i appears similar to the average waveform. We tested thisequivalence by averaging the 640 pulse profiles into singletemplate. We compared this average template with the mode V i of our PCA technique, and found the two consistentwithin the noise.The second mode V i allows the two peaks in the wave-form to depart in amplitude from the average. Essentially, U E tracks the variation with epoch of the waveform mode,changing sign when the mode switches. MNRAS000 , 1–9 (2017)
H.-H. Lin et al.
Figure 3.
The 1st to the 10th modes of V are in the sequenceof red, green, blue, yellow, cyan, black and dark to light greys,respectively. We set the offset of the 1st to the 10th modes bysteps of -0.3. Since we intend to use the V modes as templates forfitting pulsar phases, a concern is that noise in these modeswill bias the phase measurements. To reduce this noise forthe modes above the first, we set them to zero outside ofthe 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 inSection 2.2 results in a smooth structure outside this regionand so we spline the profile outside the central 160 phasebins. The natural next step would be to simultaneously fit modeamplitudes and a pulse phase to each of our individual pulseprofiles, using V mi as templates. However, we find that forthis high-dimensional parameter space, the likelihood oftenhas multiple maxima, and so the traditional least-squaresfitting method fails. The multi-modal nature of the param-eter space is due to the fact that for some profiles, there aremultiple combinations of the V mi modes that, with differentphases, may adequately describe the profiles such as thoseshown in Figure 2. To deal with this, we employ a new pulseprofile fitting technique that fully samples the parameterspace.We denote an individual measured folded pulsar profile d i where the index i runs over the phase bins (we suppressthe index E for epoch in this section, since the analysis isperformed independently for each profile). We model themeasured profile as d i = (cid:213) m A m V mi ( φ ) + n i , (2)where V mi ( φ ) are the profile template modes (with m run-ning over mode number), A m is the mode amplitude, φ isthe finely adjusted pulsar phase, and n i is the noise contri-bution. The dependence of V mi on φ describes the rotationof the templates to match the data. For notational brevity, we now switch to vector notation, where the above equationbecomes: d = A T V ( φ ) + n . (3)We assume Gaussian noise that is uniform and uncorrelated: (cid:104) n i n j (cid:105) = δ ij σ . The template modes have already been mea-sured so the parameters of the model are φ and the ampli-tudes, A .Of these parameters we are chiefly interested in φ . Allthe information about this parameter is in p ( φ | d ) , the pos-terior probability distribution of φ given the measurements d . For instance, the measurement mean and variance of thephase are the first and second moments of this distribution.This can be calculated from the posterior of all the param-eters marginalized over the mode amplitudes: p ( φ | d ) = ∫ d N A p ( φ, A | d ) , (4)where N is the length of the vector A .For flat priors on the initial parameters, Bayes’ Theoremstates that the posterior distribution for the parameters isproportional to the likelihood function, p ( φ, A | d ) ∝ p ( d | φ, A ) .Since the data is Gaussian-distributed, the latter is givenby: p ( d | φ, A ) ∝ exp (cid:20) − χ ( d , φ, A ) (cid:21) (5) χ ( d , φ, A ) = σ (cid:213) i (cid:34) d i − (cid:213) m A m V mi ( φ ) (cid:35) . (6)Combining the above, we have: p ( φ | d ) ∝ ∫ d N A exp − σ (cid:213) i (cid:34) d i − (cid:213) m A m V mi ( φ ) (cid:35) (7)This is a multidimensional integral over the vector space of A , which would be prohibitively expensive to do numeri-cally. 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 doneanalytically. This same insight was used in Pennucci (2013,Chapter 2) where a similar integral appears over the fre-quency dependant template amplitudes in wide-band pulsartiming. However, while in wide-band timing this is a compu-tational convenience, here being able to perform this integralis essential due to the multi-modal nature of the likelihoodspace. Similarly, the integral appears in Lentati et al. (2015)when marginalizing over epoch dependant profile amplitudesin profile domain timing analysis.When performing the integral, there is a factor thatdepends on the expression (cid:213) i V mi ( φ ) V ni ( φ ) . This expression is independent of φ because the dependenceof V mi on φ is just a shift of the template modes, whichdisappears after summing over phase bins. Thus, there is noneed to carry out this part of the calculation. Note that thisis only true if the noise is uniform. The final expression is p ( φ | d ) ∝ exp − σ (cid:213) m (cid:213) i V mi ( φ ) d i (cid:213) j V mj ( φ ) d j , (8) MNRAS , 1–9 (2017) mproved Pulsar Timing Figure 4. Pulse profile fit residuals in the frequency do-main.
Real and imaginary parts of the Fourier components areshown. The red curve and blue dots are the model and the mea-sured flux, respectively. Green squares in lower panel are residu-als. Fitting for phase in the frequency domain allows for simpleinterpolation with time precision finer than the 1.024 ms cadenceof the dynamic spectra. or equivalently p ( φ | d ) ∝ exp (cid:20) − σ ( Vd ) T Vd (cid:21) . (9)which is proportional to p ( φ | d ) ∝ exp (cid:26) − χ (cid:2) d , φ, ˆ A ( d , φ ) (cid:3)(cid:27) , (10)where ˆ A ( d , φ ) is the linear-best-fit value for the templateamplitudes at fixed φ . Linear regression gives ˆ A ( d , φ ) = (cid:104) V ( φ ) V T ( φ ) (cid:105) − V ( φ ) d . (11)So succinctly, the procedure for estimating the time ofarrival is: (i) At fixed φ perform a linear fit for the mode amplitudes A , which we found in Section 3.1.(ii) Evaluate χ for these parameters and take e − χ / asthe likelihood.(iii) Repeat for a range of φ values, covering the region oflow χ that dominates the likelihood. (We repeat this stepfor a range of -20 to +20 phase bins from the fixed φ , withsteps of 0.02 phase bins.)(iv) Take the zeroth, first, and second moments to get thenormalization, phase estimate, and variance, respectively. InSection 3.3, we use the phase estimate and variance to findthe timing solution.While our procedure has been described in the profilespace, our actual analysis is performed in the Fourier do-main, as in Taylor (1992). We use only the first through 90 th harmonics to limit contamination from noise in the templatemodes and observed extraneous noise near 300 Hz. We use N = template modes in our fits. An example of a profile fitis shown in Figure 4.Finally, properly estimating the error on φ requires anestimate of the noise power σ . We use the value of σ thatresults in a reduced chi-squared of unity ( χ / ν , where ν isthe number of degrees of freedom) for the best fit (maximumlikelihood) parameters. We proceed to adjust the parameters of the pulsar timingmodel: φ E = φ + ∆ PP ( τ E − τ ) + (cid:219) P P ( τ E − τ ) + P (cid:20) − ∂τ E ∂α ∆ α − ∂τ E ∂δ ∆ δ (cid:21) (12)where, φ E is the barycentric phase of the pulse profile in-dexed by E , τ E is its barycentric time, τ is the refer-ence epoch, P is the period of the pulsar used for folding, ∂τ E / ∂α is the derivative of the barycentric time with respectto source right ascension, and ∂τ E / ∂δ is derivative of thebarycentric time with respect to source declination. Thereare also five free parameters in the timing model, includingthe period derivative (cid:219) P , period correction ∆ P , initial phaseoffset φ , right ascension correction ∆ α , and declination cor-rection ∆ δ . We proceed to adjust these parameters such thatthe above equation fits our phase measurements using stan-dard weighted least squares and extract best fit parameterswith uncertainties.One challenge to obtaining a timing solution with thisdata set is the large gap in our data. We have phase mea-surements for epochs over several months in 2011 and 2015but none in between. Compared to more uniformly sampleddata, this distribution of measurements weakly constrainsthe total number of pulsar rotations in the gap period. In-deed we find multiple χ minima corresponding to chang-ing the parameter ∆ P by integer multiples of ∼ × − s,which changes the number of rotations in the gap. However,the lowest of these minima is δ χ = . smaller than thenext lowest. Thus this minimum is strongly preferred overthe others, and our timing parameters are well constrained.Inspecting the timing residuals, we find one phase mea-surement that is an extreme, σ , outlier. Inspecting the full MNRAS000
Real and imaginary parts of the Fourier components areshown. The red curve and blue dots are the model and the mea-sured flux, respectively. Green squares in lower panel are residu-als. Fitting for phase in the frequency domain allows for simpleinterpolation with time precision finer than the 1.024 ms cadenceof the dynamic spectra. or equivalently p ( φ | d ) ∝ exp (cid:20) − σ ( Vd ) T Vd (cid:21) . (9)which is proportional to p ( φ | d ) ∝ exp (cid:26) − χ (cid:2) d , φ, ˆ A ( d , φ ) (cid:3)(cid:27) , (10)where ˆ A ( d , φ ) is the linear-best-fit value for the templateamplitudes at fixed φ . Linear regression gives ˆ A ( d , φ ) = (cid:104) V ( φ ) V T ( φ ) (cid:105) − V ( φ ) d . (11)So succinctly, the procedure for estimating the time ofarrival is: (i) At fixed φ perform a linear fit for the mode amplitudes A , which we found in Section 3.1.(ii) Evaluate χ for these parameters and take e − χ / asthe likelihood.(iii) Repeat for a range of φ values, covering the region oflow χ that dominates the likelihood. (We repeat this stepfor a range of -20 to +20 phase bins from the fixed φ , withsteps of 0.02 phase bins.)(iv) Take the zeroth, first, and second moments to get thenormalization, phase estimate, and variance, respectively. InSection 3.3, we use the phase estimate and variance to findthe timing solution.While our procedure has been described in the profilespace, our actual analysis is performed in the Fourier do-main, as in Taylor (1992). We use only the first through 90 th harmonics to limit contamination from noise in the templatemodes and observed extraneous noise near 300 Hz. We use N = template modes in our fits. An example of a profile fitis shown in Figure 4.Finally, properly estimating the error on φ requires anestimate of the noise power σ . We use the value of σ thatresults in a reduced chi-squared of unity ( χ / ν , where ν isthe number of degrees of freedom) for the best fit (maximumlikelihood) parameters. We proceed to adjust the parameters of the pulsar timingmodel: φ E = φ + ∆ PP ( τ E − τ ) + (cid:219) P P ( τ E − τ ) + P (cid:20) − ∂τ E ∂α ∆ α − ∂τ E ∂δ ∆ δ (cid:21) (12)where, φ E is the barycentric phase of the pulse profile in-dexed by E , τ E is its barycentric time, τ is the refer-ence epoch, P is the period of the pulsar used for folding, ∂τ E / ∂α is the derivative of the barycentric time with respectto source right ascension, and ∂τ E / ∂δ is derivative of thebarycentric time with respect to source declination. Thereare also five free parameters in the timing model, includingthe period derivative (cid:219) P , period correction ∆ P , initial phaseoffset φ , right ascension correction ∆ α , and declination cor-rection ∆ δ . We proceed to adjust these parameters such thatthe above equation fits our phase measurements using stan-dard weighted least squares and extract best fit parameterswith uncertainties.One challenge to obtaining a timing solution with thisdata set is the large gap in our data. We have phase mea-surements for epochs over several months in 2011 and 2015but none in between. Compared to more uniformly sampleddata, this distribution of measurements weakly constrainsthe total number of pulsar rotations in the gap period. In-deed we find multiple χ minima corresponding to chang-ing the parameter ∆ P by integer multiples of ∼ × − s,which changes the number of rotations in the gap. However,the lowest of these minima is δ χ = . smaller than thenext lowest. Thus this minimum is strongly preferred overthe others, and our timing parameters are well constrained.Inspecting the timing residuals, we find one phase mea-surement that is an extreme, σ , outlier. Inspecting the full MNRAS000 , 1–9 (2017)
H.-H. Lin et al. likelihood curve for the profile phase fit (as described in Sec-tion 3.2), we find the likelihood to be very non-Guassian,making such an outlier roughly . likely. This is not ter-ribly improbable given that we have 272 such data points,and is far more likely than a σ outlier for Gaussian data.Since our least-squares fit to the timing solution implicitlyassumes Gaussian-distributed phase errors, we exclude thisoutlying data point in the timing solution fit.The timing residuals are shown with the fitting resultsin Figure 5. The fit has a reduced chi-squared of χ / ν = . for 266 degrees of freedom, which is reasonable given thenon-Gaussian nature our phase measurements which wehave represented with 1- σ error bars. We inflate the un-certainties on the timing parameters by this factor.The parameters for the final timing solution forPSR J2139+0040 are given in Table 1. We find new rightascension and declination, period, period derivative, refer-ence time of arrival (derived from φ ), and also provide theepoch used for the timing model. We use the period and pe-riod derivative to estimate the magnetic field strength ( B ),the characteristic age ( t c ), and the spin-down luminosity( (cid:219) E )(Lorimer & Kramer 2004).To check for evidence of proper motion of the pulsar’ssky position, we fix all parameters other than the position,then separately fit the 2011 and 2015 data. We find the twooutput positions are consistent with each other and the solu-tion using the full data set, and as such we find no evidencefor proper motion.We used the TEMPO pulsar timing software to verify ourtiming solution. For this we converted our pairs of epoch-phase measurements to barycentric times of arrival then in-verted
PRESTO ’s bary command to convert to topocentrictimes of arrival. Feeding these and our timing solution pa-rameters into TEMPO we find agreement in the value of χ .Letting TEMPO refit our timing-model parameters we find nosignificant shifts.
We used the pointed data to determine DM. We align the 40folded profiles in phase based on our timing solution, thenstacked them into a single profile. We fit this profile for thedispersion using the first V mode as a template accordingto: d i f = A (cid:18) f
800 MHz (cid:19) β V i ( φ f ) (13) φ f ≡ φ + DM P .
808 s MHz pc / cm (cid:20) f − ( ) (cid:21) (14)Here, d i f is the frequency-dependant profile data, which arefit 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 φ and the dispersion phase delay.We subsequently calculate the dispersion-measure-based distance using Cordes-Lazio NE2001 Galactic FreeElectron Density Model (Cordes & Lazio 2002) website us-ing the new right ascension, declination, and DM. The re-sults of this analysis are given in Table 1 To assess the impact of PCA mode tracking we compared thetiming residuals using long and short integrations. We car-ried out these parallel analyses on sets of 16 individual 5.4 sprofiles from the pointed data. First, we stacked these pro-files into a single pulse profile with a total integration timeof 86 s and fit this profile with the technique described inSection 3.2. Second, we fit the 16 profiles individually, thencombined their resulting posterior distributions. We did thisfor the 40 sets of 16 profiles covering the full hour of pointeddata. The first technique is similar to the conventional tech-nique of taking long averages to smooth out pulse-to-pulsevariability. The second technique uses fine grained data al-lowing the PCA algorithm to track and compensate for thevariability.We find that using PCA mode tracking and fitting fine-grained 5.4 second averages results in a ∼ smaller un-certainty on average in the pulse time of arrival (calculatedfrom the second moment of the posterior as described inSection 3.2). While the improvement is substantial for thispulsar, it may not be generic. Pulsars are highly individual:they vary substantially in the degree and time scales modeswitching, so we anticipate the improvement available fromPCA mode tracking will also vary strongly from case to case. Pulsar timing observations often make use of profiles aver-aged for minutes. These long averages smooth out the effectof rapid mode switching. Since it is in general sub-optimalto average over a time-variable signal, we suggest that PCAMode Tracking could be applied broadly, by intentionallyusing short averages that display rapid mode switching—which the PCA then follows and compensates for. In par-ticular, it may be useful to test this technique using datafrom millisecond pulsars, since improvement of timing pre-cision of these pulsars allows more precise tests of generalrelativity, along with tighter constraints on long-wavelengthgravitational wave backgrounds.We have obtained a substantially improved timing so-lution for PSR J2139+0040, using untargeted observationsin conjunction with a ∼ (cid:219) P = ( . ± . ) × − (s s − ) is well below average, with ∼ of slow pulsars having a higher rate (Johnston &Karastergiou (2017)) as shown in Figure 6. This indicatesthat PSR J2139+0040 has an abnormally weak magneticfield.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 am-plitudes, we are able to fully sample the marginalized pos-terior for phase and calculate the distribution’s mean andstandard deviation. The subsequent fit to a timing modelachieves a reduced chi-squared of 1.26, which is reasonablegiven the non-Gaussian nature of our phase errors. A morecomplete treatment would be to use the phase measurement MNRAS , 1–9 (2017) mproved Pulsar Timing Figure 5. Timing residuals with error bars of WZ data from MJD 55694 to 57178.
The upper panel shows the residuals overfive years. The left lower panel shows detail of the residuals of 2011 WZ data. The denoted with a circle around in the left lower panel isan outlier resulting from the non-Gaussian nature of the error bars, which we exclude from the timing solution. The lower middle panelshows the 2015 WZ data. The lower right panel shows the 2015 pointed data. All data are shown offset from the best fit model indicatedwith the dashed red line. Right ascension, α +
21 : 39 : 25 . ( ) Declination, δ + ◦ (cid:48) . ( ) ” Period, P (s) 0.3124695464326(2)Period derivative, (cid:219) P (s s − ) . ( ) × − Dispersion measure, DM (pc cm − ) 31.585(8)Period epoch, τ (MJD) 56436.5Reference arrival time (Barycentric MJD, 900 MHz) 56436.506782091(13)Magnetic field strength, log [ B / G ] log [ t c / yr ] log [ (cid:219) E /( erg / s )] d (kpc) 1.7 Table 1. Timing solution and measured parameters of PSR J2139+0040. posteriors directly when fitting a timing solution, ratherthan reducing them to their first and second moments. Thiswould optimally extract the timing information from thephase measurements but would substantially complicate fit-ting a timing solution, since least-squares would no longerbe applicable.Our PCA mode tracking method addresses issues sim-ilar to those addressed via ”profile-domain pulsar timing”(Lentati et al. 2015; Lentati & Shannon 2015; Lentati et al.2017), particularly, what these authors refer to as “low-frequency stochasticity” and “phase-correlated stochastic-ity”. The profile-domain techniques uses pre-defined uncor- related shapelet components to contribute to the individualprofiles. In contrast, the PCA mode tracking method is non-parametric. There is no need to guess the shaplets, the PCAfinds them automatically. Compared to the profile-domainstrategy we anticipate that fewer components will be needed,yielding less degeneracy with the pulse phase and likely re-ducing timing error. The PCA mode tracking technique de-rives its mode templates from the data, and so for weakpulsars these waveforms may have noise that is absent whenusing the pre-defined templates of the profile domain tech-nique.A PCA technique was employed to characterize pul-
MNRAS000
MNRAS000 , 1–9 (2017)
H.-H. Lin et al.
Figure 6.
Periods and Period derivatives of known pul-sars(Manchester et al. 2005). Pulsar J2139+0040 lies at low (cid:219) P compared to other pulsars of similar period, but is not an ex-treme outlier. sar variability in Os(cid:32)lowski et al. (2011), where, rather thanfitting for the time of arrival and component amplitudes inhigh cadence data, the authors used the empirical correlationbetween the timing residuals and component amplitudes toretroactively correct for systematic timing error.Our results demonstrate the feasibility of timing pulsarscommensally with mapping surveys. In particular, upcominghydrogen intensity mapping surveys using MeerKAT (Fon-seca et al. 2017) and Phase 1 of the Square Kilometre Array(Santos et al. 2015) will map large fractions of the southernsky with thousands of hours of telescope time. The tele-scopes will operate as a collection of single dishes and willthus need to scan rapidly to overcome / f noise. There isthe potential to obtain pulsar timing solutions for free (interms of telescope time) using the data from these surveys inmuch the same way as we have done here. The challenge isin storing the data at the rapid cadence required for pulsarstudies, but, in single dish mode, data volumes are modestcompared to interferometric mode. Mapping surveys spenda 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 volumeof timing data.One of our goals in this work is to enable pulsar searches commensal with intensity mapping experiments. Having ob-tained a precise timing solution for a known pulsar in theGBTIM data, we have demonstrated that accumulation ofpulsar data over five years in few-pulse snippets can be ac-complished, although this may be more difficult for weakeror more erratic pulsars.This analysis was possible because the GBTIM dataused short 1 ms integrations. The intensity mapping datafrom CHIME (Bandura et al. 2014) and HIRAX (Newburghet al. 2016) will use longer integrations for their intensitymapping data, but these instrument will have additionaltransient-search backend hardware allowing high cadenceanalysis of the data streams. New search algorithms haverecently been proposed by Smith (2016) to reduce the cost of such long-term pulsar waveform assembly by several or-ders of magnitude. If these techniques could be employed atupcoming transit survey instruments such as CHIME andHIRAX the pay off could be a substantially increased rateof pulsar discovery. ACKNOWLEDGEMENTS
We thank Kendrick Smith, Ingrid Stairs, Maura McLaughlinand Alexander Roman for valuable discussions. K. W. M. issupported by the Canadian Institute for Theoretical Astro-physics National Fellows program. U.-L. P. acknowledgessupport from the Natural Sciences and Engineering Re-search Council of Canada. Research at Perimeter Instituteis supported by the Government of Canada through the De-partment of Innovation, Science and Economic DevelopmentCanada and by the Province of Ontario through the Ministryof Research, Innovation and Science. J. B. P. acknowledgessupport from NSF Award 1211777. Computations were per-formed on the GPC supercomputer at the SciNet HPC Con-sortium.
REFERENCES
Anderson S. B., 1993, PhD thesis, Caltech, http://resolver.caltech.edu/CaltechETD:etd-08282008-094151
Bandura K., et al., 2014, in Ground-based and Air-borne Telescopes V. p. 914522 ( arXiv:1406.2288 ),doi:10.1117/12.2054950Camilo F., Nice D. J., Shrauner J. A., Taylor J. H., 1996, ApJ,469, 819Chang T.-C., Pen U.-L., Bandura K., Peterson J. B., 2010, Na-ture, 466, 463Cordes J. M., Lazio T. J. W., 2002, preprint, ( arXiv:astro-ph/0207156 )DeBoer D. R., et al., 2017, PASP, 129, 045001Fonseca J., Maartens R., Santos M. G., 2017, MNRAS, 466, 2780Glazebrook K., et al., 2007, in Metcalfe N., Shanks T., eds, As-tronomical Society of the Pacific Conference Series Vol. 379,Cosmic Frontiers. p. 72 ( arXiv:astro-ph/0701876 )Hotan A. W., van Straten W., Manchester R. N., 2004, Publ.Astron. Soc. Australia, 21, 302Johnston S., Karastergiou A., 2017, MNRAS, 467, 3493Lentati L., Shannon R. M., 2015, MNRAS, 454, 1058Lentati L., Alexander P., Hobson M. P., 2015, MNRAS, 447, 2159Lentati L., et al., 2017, MNRAS, 466, 3706Lorimer D. R., Kramer M., 2004, Handbook of Pulsar AstronomyManchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ,129, 1993Masui K. W., et al., 2013, ApJ, 763, L20Masui K., et al., 2015, Nature, 528, 523Newburgh L. B., et al., 2016, in Ground-based and Air-borne Telescopes VI. p. 99065X ( arXiv:1607.02059 ),doi:10.1117/12.2234286Os(cid:32)lowski S., van Straten W., Hobbs G. B., Bailes M., DemorestP., 2011, MNRAS, 418, 1258Pennucci T. T., 2013, PhD thesis, University of Virginia (USA)Santos M., et al., 2015, Advancing Astrophysics with the SquareKilometre Array (AASKA14), p. 19Smith K. M., 2016, preprint, ( arXiv:1610.06831 )Switzer E. R., et al., 2013, MNRAS, 434, L46Taylor J. H., 1992, Philosophical Transactions of the Royal Soci-ety of London Series A, 341, 117 MNRAS , 1–9 (2017) mproved Pulsar Timing Yusifov I., K¨u¸c¨uk I., 2004, A&A, 422, 545This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000