Measuring pulse times of arrival from broadband pulsar observations
K. Liu, G. Desvignes, I. Cognard, B. W. Stappers, J. P. W. Verbiest, K. J. Lee, D. J. Champion, M. Kramer, P. C. C. Freire, R. Karuppusamy
aa r X i v : . [ a s t r o - ph . GA ] J u l Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 5 November 2018 (MN L A TEX style file v2.2)
Measuring pulse times of arrival from broadband pulsarobservations
K. Liu, , G. Desvignes, I. Cognard, , B. W. Stappers, J. P. W. Verbiest, , K. J. Lee, , D. J. Champion, M. Kramer, , P. C. C. Freire, and R. Karuppusamy, Station de radioastronomie de Nan¸cay, Observatoire de Paris, CNRS/INSU, F-18330 Nan¸cay, France Laboratoire de Physique et Chimie de l’Environnement et de l’Espace LPC2E CNRS-Universit´e d’Orl´eans, F-45071 Orl´eans Cedex 02, France University of Manchester, Jodrell Bank Centre of Astrophysics, Alan-Turing Building, Manchester M13 9PL, UK Max-Planck-Institut f¨ur Radioastronomie, Auf dem H¨ugel 69, D-53121 Bonn, Germany Fakult¨at f¨ur Physik, Universit¨at Bielefeld, Postfach 100131, 33501 Bielefeld, Germany KIAA, Peking University, Beijing 100871, P.R. China
ABSTRACT
In recent years, instrumentation enabling pulsar observations with unprecedentedlyhigh fractional bandwidth has been under development which can be used to substan-tially improve the precision of pulsar timing experiments. The traditional template-matching method used to calculate pulse times-of-arrival (ToAs), may not functioneffectively on these broadband data due to a variety of effects such as diffractivescintillation in the interstellar medium, profile variation as a function of frequency,dispersion measure (DM) evolution and so forth. In this paper, we describe the chan-nelised Discrete Fourier Transform method that can greatly mitigate the influence ofthe aforementioned effects when measuring ToAs from broadband timing data. Themethod is tested on simulated data, and its potential in improving timing precision isshown. We further apply the method to PSR J1909 − Key words: methods: data analysis — pulsars: general — pulsars: individual(PSR J1909 − High-precision pulsar timing is the essential tool in thecurrently ongoing gravity tests with binary pulsars (e.g.Freire et al. 2012; Antoniadis et al. 2013), and gravita-tional wave detection with pulsar timing arrays (e.g.Hellings & Downs 1983; Foster & Backer 1990). Develop-ments in observing backends based on polyphase filterbank techniques now allow pulsar observations with largefractional bandwidth (e.g. & /
3, see Ransom et al. 2009;Stappers et al. 2011; Cognard et al. 2013). Broadband sys-tems such as the Ultra-broadband (UBB) receiver developedat the Effelsberg radio telescope which achieves an instanta-neous frequency coverage from 600 MHz to 3 GHz are beingimplemented (Karuppusamy et al. 2014). Extending the ob-serving bandwidth greatly increases the signal-to-noise ratio(S/N) of the observed pulsar profiles which can be translatedinto a decrease in measurement uncertainties of pulse times- of-arrival (ToAs, Downs & Reichley 1983). This is likely toimprove the resulting pulsar timing precision, under the cir-cumstance that the actual timing residuals are dominatedby white noise (Cordes & Shannon 2010; Liu et al. 2011).Typically, determination of ToAs is achieved by firstcorrecting for any frequency-dependent delays caused by dis-persion in the interstellar medium, then averaging the dataacross all frequency channels and finally cross-correlatingthe resulting pulse profile with a template (Taylor 1992).The template can either be a high-S/N profile obtained di-rectly from observations or an analytic model of the pulse(e.g. Kramer et al. 1999). However, applying this approachto broadband pulsar timing data may be problematic for anumber of reasons as follows: • The pulse profile shape is not constant across the ob-serving bandwidth. This can be due to intrinsic profilevariation across frequency, frequency-dependent scattering c (cid:13) K. Liu et al. timescales, instrumental channelisation and so forth. Whenthe observed flux density at different frequencies variesin time due to interstellar scintillation, the shape of thefrequency-averaged profile will be changed depending onwhich part of the band scintillates up and which part scintil-lates down. In turn, this variation in profile shape will affectthe cross-correlation and therefore the derived ToAs (e.g.Liu et al. 2011). • In many cases the pulsar dispersion measure (DM) hasbeen witnessed to evolve in time (e.g. Keith et al. 2013).Hence, applying a constant but not perfectly correct DMvalue for dedispersion would smear the profile shape differ-ently as the DM values vary, and thus shift the ToAs. Usingan incorrect DM can also induce misalignment of profilesfrom different frequencies, which would mimic profile varia-tion across the observing band and enhance the scintillationeffect mentioned above. • Variations in the radio frequency interference (RFI) en-vironment can change the usable frequency channels in dif-ferent observations. This effectively changes the frequencyrange from which the frequency-averaged profile is com-posed of due to the frequency-dependence of the pulse pro-file shape. Similar variations in the frequencies used (andtherefore in the resulting profile shape) could be caused byinstrumental failure or slight changes to the observing sys-tem.Consequently, an alternative approach that utilises thedata’s frequency information to avoid the corrupting effectslisted above, is required. Without frequency-averaging, anyfrequency-dependent changes in brightness will no longer re-sult in changes in the overall pulse shape. DM variations canalso be taken into account effectively if the data cover a largeenough frequency range.Besides the current standard template-matching schemein Taylor (1992), there are also alternative methods tar-geting different issues. Hotan et al. (2005) applied Gaus-sian interpolation when cross-correlating profiles with lowS/Ns. van Straten (2006) used full Stokes informationrather than only the total intensity for ToA determination.Os lowski et al. (2011) decomposed profiles in principal com-ponents to improve ToA measurements when the timingresiduals are dominated by pulse phase jitter. In this pa-per we develop a method to measure ToAs based on timingdata with frequency information.The structure of this paper is as follows. In Section 2we describe our method to generate ToAs with frequency-resolved data. In Section 3 we present results from testsbased on simulated data. Application of the method to realdata is shown in Section 4. We conclude with a brief discus-sion in Section 5.
In the traditional approach, the observed pulse profile P ( t ),is described by a one dimensional array modelled by P ( t ) = a + bS ( t − ∆ τ ) + n ( t ) , (1) Dispersion measure is defined as the integrated electron densitybetween the Earth and the pulsar. where S ( t ) is the standard template obtained from previousobservations, a is the baseline difference, b is a scaling fac-tor, ∆ τ is the phase difference between the profile and thetemplate, and n ( t ) is a noise component. The actual fittingroutine is often carried out in the Fourier domain, whereafter a discrete Fourier transform (DFT), the model can bewritten as (Taylor 1992) P k e iθ k = aNδ k + bS k e i ( φ k + k ∆ τ ) + G k , k = 0 , ..., N − , (2)where P k and S k are the amplitudes of the complex Fouriercoefficients, θ k and φ k are the phases, δ k = 1 ( k = 0) or 0( k = 0), N is the number of frequency bins, and G k rep-resents random noise equal to the Fourier transform of thesampled noise in the time-domain profile, n ( t ). The best es-timated value of ∆ τ , can then be found by minimising thegoodness-of-fit statistic χ ( b, τ ) = N − X k =1 (cid:12)(cid:12)(cid:12)(cid:12) P k − bS k e i ( φ k − θ k + k ∆ τ ) σ k (cid:12)(cid:12)(cid:12)(cid:12) , (3)where σ k is the root-mean-square intensity of the noise atfrequency k .If the frequency resolution is kept, the data arrays haveto be expanded into a second dimension and in this casethe dispersive delay of different frequency bands also needsto be taken into account. Mostly, the delay is found to ful-fill the scaling of t d ∝ DM /f , where f is the observingfrequency (Lorimer & Kramer 2005) . Accordingly, one canexpress the two-dimensional model as: P ( f, t ) = a ( f ) + b ( f ) S ( f, t − ∆ τ − D × ∆DM /f ) + n ( f, t ) , (4)where D is the dispersion constant and ∆DM is the dif-ference in dispersion measure between the profile and thetemplate. Note that
D × ∆DM /f vanishes when f → ∞ .Hence, ∆ τ can be directly related to the ToA at infinitefrequency, which is obtained after correcting the dispersiondelay.Similar to the treatment in the one-dimensional case,one can carry out a DFT individually on each frequencyband, and perform a two-dimensional fit over the entirebandwidth for ∆ τ and ∆DM. The approach will be laterreferred to as the “channelised DFT method”. After thetransformations, the model in the frequency domain can bewritten as P j,k e iθ j,k = a j Nδ k + b j S j,k e φ j,k + kτ j + G j,k , (5)where k = 0 , ..., N − j = 1 , ..., N b ; N is the number ofbins, N b is the number of frequency bands, and τ j = D × ∆ DMf j πP + ∆ τ. (6)Here τ j is in units of radians and P is the pulsar’s rotationalperiod. The data of different frequencies are assumed to befrom a single observation, which means they are aligned in Further propagation effects in the interstellar medium (e.g.scattering, as shown in Bhat et al. 2004) could introduce otherscalings, too, albeit at a less significant level. A quantitative sum-mary of such is not within the scope of this paper, but can befound in e.g. Lorimer & Kramer (2005). In this paper we use the definition of
D ≡ /K where K ≡ . × − MHz − cm pc s − .c (cid:13)000
D ≡ /K where K ≡ . × − MHz − cm pc s − .c (cid:13)000 , 000–000 -D template-matching time by default and thus ∆ τ is constant over the entire band.Following Taylor (1992), a j is immediately obtained from a j = P j, − b j S j, N . (7)The best estimates of ∆ τ and ∆DM can then be obtainedby minimising the function χ ( b i , ∆ τ, ∆DM) = N b X j =1 N − X k =1 (cid:12)(cid:12)(cid:12)(cid:12) P j,k − b j S j,k e i ( φ j,k − θ j,k + kτ j ) σ j,k (cid:12)(cid:12)(cid:12)(cid:12) . (8)Here the same as in Taylor (1992), the likelihood estimatoris calculated using all informative frequency harmonics . Inorder to perform the fitting, we modified the Levenberg-Marquardt (L-M) routine in Press et al. (1992) to adopta model of complex numbers, and recalculated the curva-ture matrix in this case. More details can be found in Ap-pendix A. The 1- σ errors of the fitted parameters are de-termined from the covariance matrix which is obtained bytaking the inverse of the curvature matrix. A similar ap-proach based on the same likelihood estimator as in Eq. (8),using however a different fitting routine, has been developedby Pennucci et al. (2014).It can be seen from the expression of χ in Eq. (8), thatfor data with a bandwidth that does not result in a signif-icant difference between f and f N b , ∆ τ and ∆DM wouldbe highly correlated. Following the calculation in Lee et al.(2014), in Fig. 1 we plot the correlation coefficient ( ρ ) be-tween the two parameters for a given observing bandwidthand lower bound of observing frequency. Clearly, ρ > . < ρ less than 0.7. Nevertheless, when the observing band-width is not enough to be sensitive to ∆DM fit, one canalways choose to fit solely for ∆ τ . Under this circumstance,the likelihood estimator is simplified into χ ( b i , ∆ τ ) = N b X j =1 N − X k =1 (cid:12)(cid:12)(cid:12)(cid:12) P j,k − b j S j,k e i ( φ j,k − θ j,k + k ∆ τ ) σ j,k (cid:12)(cid:12)(cid:12)(cid:12) . (9) In this section, the channelised DFT method is tested basedon simulated data. Here we carried out two types of tests:one to assess the functionality of the algorithm in general,and the other to evaluate the potential improvement to ToAmeasurements when applied in a realistic scenario. In Taylor (1992), the calculation of χ includes frequency binsonly to N/ . . . . . . . . . . . . . . . f low (GHz) B ( G H z ) Figure 1.
Correlation coefficient ( ρ ) between ∆ τ and ∆DM fordifferent bandwidths and frequency ranges. Here f low is the fre-quency at the lower edge of the bandwidth, ν is the total band-width. The Effelsberg UBB setting is represented by the star. Ideally, a functioning template-matching approach is ableto measure the offsets in phase and DM between the profileand the template, as well as their 1- σ uncertainties. In thefollowing tests, we simulate data with a bandwidth of 500MHz, between 1.2 and 1.7 GHz. A Gaussian shape with nofrequency dependence is assumed for the two-dimensionaltemplate. Individual observations are simulated by addingwhite noise to the standard profile, until a given S/N isreached. We hereby define the S/N as the ratio of the peakamplitude of the pulse and the rms in the off-pulse region,and use this definition throughout this paper. In order to show that the method measures the parametersconsistently with the calculated errors, we simulated 2 × profiles with randomly distributed offsets in phase and DMwith respect to the two-dimensional template. Fig. 2 showshistograms of measured ∆ τ (top) and ∆DM (bottom) valuesfrom the channelised DFT method, after subtraction of theinput offsets (∆ τ i and ∆DM i , separately). In both cases, thedistribution is well described by a Gaussian function andthe rms is consistent with expected from the measurementuncertainty. The intrinsic uncertainty of a ToA is expected to be onlyrelated to the S/N and the profile shape, and to scale in-versely with the S/N (Downs & Reichley 1983). In orderto show that the method measures the phase offset withthe expected uncertainty, we simulated profiles with differ-ent S/Ns and with different numbers of frequency channels. c (cid:13) , 000–000 K. Liu et al. N u m be r ∆τ - ∆τ i (bin) 0 500 1000 1500 2000-20 -15 -10 -5 0 5 10 15 20 N u m be r ∆ DM- ∆ DM i (10 -4 cm -3 pc) Figure 2.
Histograms of measured ∆ τ (top) and ∆DM (bottom)from 2 × profiles after subtraction of the input values (∆ τ i and∆DM i , separately). The profiles were simulated with S/N=100,16 frequency channels, and randomly distributed ∆ τ and ∆DMwith respect to the template. The distributions are well describedby Gaussian functions, with reduced χ of fit close to unity (Top:0.96; Bottom: 0.92). The rms in the two cases are 0.176 bins and4 . × − cm − pc, respectively, consistent with the measure-ment uncertainties which are 0.177 bins, 4 . × − cm − pc, re-spectively. Fig. 3 shows the measured phase uncertainties σ ∆ τ obtainedfrom these profiles. The values were compared with expectedmeasurement uncertainties calculated from the radiometerequation as in Downs & Reichley (1983). Clearly, the uncer-tainties for data of identical S/N all fully agree with the ex-pectation, and the number of channels has no effect (as is tobe expected for a pulse profile with no frequency evolution).Note that for this test the fit of ∆DM has been switchedoff, as the strong correlation (corresponding to ρ = 0 . τ would greatly worsen themeasurement uncertainty. In this case it is not expected thatthe resulting ToA uncertainties would be the same as calcu-lated from the radiometer equation. Fig. 4 shows the scalingof σ ∆ τ with S/N, as well as the influence of simultaneous∆DM determination. It can be seen that the scaling followsan inverse trend in either case, with or without simultaneous∆DM fitting. When both ∆DM and ∆ τ are fitted, due tothe high degree of correlation between these two parameters σ ∆ τ ( b i n ) N chan S/N=100S/N=300S/N=1000
Figure 3. σ ∆ τ obtained from the channelised DFT method forprofiles of different S/N values and numbers of channels. The linesrepresent the expected values from the radiometer equation fordifferent S/Ns. The fractional differences between the precisionsand the expectations are all within 2.5%. σ ∆ τ ( b i n ) S/N No ∆ DM fitWith ∆ DM fit
Figure 4. σ ∆ τ obtained from the channelised DFT method forprofiles with different S/Ns. Both cases are presented where ∆DMis simultaneously fitted or not. the measured uncertainties increase substantially (here by afactor of 6). It is known that in the low S/N regime (e.g. S/N close tounity) the traditional template-matching technique may failto recover the true phase offset within the measured un-certainty (e.g. Liu 2012). Fig. 5 demonstrates that this isalso true for the channelised DFT method. For this test, wegenerated pulse profiles with the number of frequency chan-nels, N chan equal to 8, 16 and 32 respectively, and variedthe S/N per frequency channel, S/N i , between 1 and 4. Foreach combination of N chan –S/N i , we simulated 10 profilesand determined how many of these resulted in measurementsconsistent with the input value for the phase offset (definedby falling into the range decided by the measurement uncer- c (cid:13) , 000–000 -D template-matching σ ∆ τ c on f i den c e l e v e l S/N i N chan =8N chan =16N chan =32N chan =1, 1D1- σ Figure 5.
Confidence level of σ ∆ τ by the channelised DFTmethod in the low S/N regime, based on profiles of different num-bers of frequency channels ( N chan ) and S/N values in individualchannels (S/N i ). The result achieved with the traditional method(1D) is also shown for comparison. tainty), while fitting for ∆ τ and ∆DM simultaneously. Thisfraction is shown on the y-axis of Fig. 5. Since the uncertain-ties are 1- σ values, one would expect all points in this plotto fall around a value of 0.68. This is clearly not the case forS/N i values less than 2.0, and the deviations from 0.68 occurat the same S/N i value for different channel numbers. Thesame situation was seen when investigating the consistencyof the measured ∆DM values. Statistics of ∆ τ measurementsbased on the traditional technique show deviation from 1- σ at the same S/N i value. Therefore, for more reliable mea-surements in this case, either frequency channels need to becombined to increase the S/N i , or an alternative approachhas to be considered (as in e.g. Hotan et al. 2005). As shown in Liu et al. (2011), using a template which is notentirely noise free may also limit the functionality of thetemplate-matching technique. This issue is demonstrated inFig. 6, where we evaluated the accuracy of the measureduncertainties using the same approach as the one followedfor Fig. 5, but this time with standard templates of S/N=10 and different N chan values. It can be seen that the σ ∆ τ valuesprove to be reliable, but only until the simulated profilesreach a S/N that is within an order of magnitude of theS/N of the template profile. Therefore, the algorithm worksas expected when the template profile has a sufficiently highS/N and the observations have a significantly lower S/N.Also in this case, the impact of the number of frequencychannels is limited, though fewer frequency channels (i.e.with higher S/N i values) do make for slightly more reliablemeasurement uncertainties. Again, the ∆DM measurementswere seen to be influenced by the noise of the template inthe same manner. As discussed in Section 1, applying the traditional template-matching approach to broadband pulsar timing data would σ ∆ τ c on f i den c e l e v e l S/NS/N=10 , N chan =8S/N=10 , N chan =16S/N=10 , N chan =321- σ Figure 6.
Confidence level of σ ∆ τ from the channelised DFTmethod, based on templates with a S/N= 10 and different num-bers of channels. A m p li t ude Phase 500 MHz1400 MHz3000 MHz
Figure 7.
Examples of frequency-dependent template shape atthree frequencies. Note that when later used to create profiles, thetemplate amplitudes from different frequencies were modulatedbased on a spectral index of -1.8. not properly produce ToAs in some circumstances. Below wedemonstrate improvements that are achieved using the chan-nelised DFT method. Here we simulate profiles and conductthe fitting based on a frequency-dependent template modeland examples can be found shown in Fig. 7.
When the DM of a pulsar varies between observations, us-ing a constant DM value for dedispersion will introduce abias in the measured ToAs. As discussed in Section 2, whenthe observing bandwidth is sufficiently large, use of a two-dimensional template can take the influence of variable DMsinto consideration when calculating the ToAs. For demon-stration purposes, we simulated 100 profiles with Gaussianlydistributed DM variations with a standard deviation of Note that the observed DMs are more likely to show gradualvariations (e.g. Keith et al. 2013). Nevertheless, our simulation isc (cid:13) , 000–000
K. Liu et al. ∆ τ ( b i n )
20 40 60 80 100 N TOA -1.5-1-0.5 0 0.5 1 1.5 0 10 20 30 40 50 60 70 80 90 100 ∆ D M - ∆ D M i ( - c m - p c ) N TOA
Figure 8.
Top:Fitted ∆ τ values of 100 profiles whose DM offsetsfrom the template follow a Gaussian distribution with standarddeviation 2 × − cm − pc. The assumed observing bandwidthis 500 MHz, from 1.2 GHz to 1.7 GHz. The input phase offset is10 bins. Here we used a S/N of 10 and kept 16 frequency bandswhen fitting for both ∆ τ and ∆DM. The fit is carried out bothwithout (left panel) and with (right panel) a fit for ∆DM. Themeasurements including a ∆DM fit have a reduced χ of 1.01,while those without result in a reduced χ of order ∼ i . The corresponding reduced χ is 1.02. × − cm − pc. The simulated data have a bandwidth of500 MHz between 1.2 and 1.7 GHz. The offsets in phase weremeasured using the channelised DFT method, both withand without a simultaneous fit for ∆DM. From the resultsin Fig. 8 (top plot, left panel), it is clear that the obtained∆ τ values are scattered beyond the tolerance of the mea-surement uncertainties if the variable DM values are notaccounted for. On the contrary, fitting for both parametersresults in ∆ τ and ∆DM values consistent with the input,and with reduced χ values close to unity for both measure-ment sets. enough to demonstrate the improvement by performing a simul-taneous fit for ∆ τ and ∆DM. t ( s ) Figure 9.
Example simulated dynamic spectrum with length 1 hrand frequency coverage from 500 MHz up to 3 GHz. The time andfrequency resolution are 10 s and 5 MHz, respectively. The scaleis given on the right hand side.
As mentioned in Section 1, the profile shape variability re-lated to diffractive scintillation would be most significant forbroadband pulsar timing data. In order to demonstrate theimprovement in timing precision with the channelised DFTmethod in this case, we simulated data covering a frequencyrange from 500 MHz up to 3 GHz. In total we created 1-hrobservations on 100 epochs. For each epoch, to better sim-ulate real observations, we generated a dynamic spectrumassuming a scintillation timescale of 20 min and frequencyscale of 50 MHz , respectively. An example of such a spec-trum can be found in Fig. 9. To create data at a given epoch,we firstly generated profiles for every 10 s and 5 MHz, basedon the template shape in Fig. 7. Here we assumed a pro-file S/N of 20 at 1.4 GHz based on 5 MHz bandwidth anda 1-hr integration. Then the amplitudes of the profiles wereweighted differently with respect to the dynamic spectrum.Next the profiles were integrated over the whole 1-hr ses-sion to be used for ToA measurements. When performingthe channelised DFT method we kept a frequency resolu-tion of 50 MHz.The results are summarised in Fig. 10. It can be seenthat using the channelised DFT method in both modes (withand without ∆DM fit) succeeds in obtaining measurementresiduals with reduced χ close to unity. The factor of nearly2 difference in the measurement precision is due to the cor-relation between ∆ τ and ∆DM which still corresponds to acorrelation coefficient of approximately 0.65. Simply apply-ing the traditional template-matching approach to the dataafter averaging over all frequency channels, leads to preci-sion more than two orders of magnitude worse and a reduced χ value of approximately 600. Note that when generating the dynamic spectra we did not con-sider the dependency of the scintillation timescale and frequencyupon observing frequency. However, the simulation is sufficientfor demonstrating the potential improvement from the new algo-rithm. c (cid:13) , 000–000 -D template-matching -1-0.8-0.6-0.4-0.200.20.40.60.81 0 20 40 60 80 R e s i dua l ( µ s )
20 40 60 80 20 40 60 80 100 N TOA
Figure 10.
Measurements of ∆ τ from 100 simulated profiles,which include a shape dependency with observing frequency andflux variation due to diffractive scintillation. Applying the chan-nelised DFT method without a ∆DM fit (left panel) leads tomeasurement residuals of RMS 47 ns and reduced χ of 1.07. In-cluding a simultaneous ∆DM fit (middle panel) reduces the pre-cision to 90 ns while the corresponding measurement reduced χ remains close to unity (0.97). The application of the traditionaltemplate-matching approach to the data after averaging over allfrequency channels (right panel) results in a RMS of 4.3 µ s and areduced χ of approximately 600. Not all of the data points aretherefore shown in the plot. Timing observations of millisecond pulsars have been con-ducted regularly at the Nan¸cay Radio Telescope (NRT).The legacy Berkeley-Orleans-Nan¸cay (BON) backend hasbeen in operation for nearly ten years, capable of produc-ing 128 MHz data which are coherently dedispersed online(Cognard & Theureau 2006). The Nan¸cay Ultimate PulsarProcessing Instrument (NUPPI), which started in late 2011,is a baseband recording system using a Reconfigurable OpenArchitecture Computing Hardware (ROACH) FPGA boarddeveloped by the CASPER group and GPUs. Here theanalog-to-digital converters firstly sample and digitise thesignal over a 512 MHz band at the Nyquist rate with dualpolarisations in 8-bit. Then a polyphase filter bank is per-formed to channelise the band into 128 channels each of4 MHz width. Next the channels are packetised into foursub-bands and sent to four GPU clusters individually via10 GbE links for online coherent dedispersion, at a DM of10.3940 cm − pc for PSR J1909 − − ∼
50 MHz (e.g. Cordes & Lazio 2002), and thus lessthan the observing bandwidth. Accordingly, we chose datacollected from ∼
30 epochs between MJD 56545 and 56592,with central frequency at 1488 MHz. The data were cali- http://casper.berkeley.edu/ -3 -2 -1
0 5 10 15 20 25 30 F R Epoch number142015481676
Figure 11.
Relative observed flux densities of the top three128 MHz sub-bands (centred at 1420, 1548 and 1676 MHz, respec-tively) with respect to the bottom one (centred at 1292 MHz) forall epochs. brated for polarisation with the common single-axis model(e.g. Ord et al. 2004), with reference to a noise diode posi-tioned at 45 degrees to the linear feed probes. The psrzap and pazi programs ( psrchive ’s RFI zapper) were used toclean any RFI. We formed a two-dimensional analytic tem-plate by fitting Gaussian components to the integration(e.g. Kramer 1994) on MJD 56592 when the source was thebrightest among our selected epochs. With the remainingdata we generated integrations of 18 −
25 min length, with anephemeris determined from the data produced by the legacyBON backend with a baseline of 8 yr. Most of the prepro-cessing was conducted with the psrchive software package(Hotan et al. 2004).
The scintillation has been significant within our observingfrequency band. In Fig. 11 we divide the whole bandwidthinto four 128 MHz sub-bands, and show the observed fluxdensities of the top three divided by that of the bottomone for all epochs. The ratios are seen to vary greatly fordifferent epochs, mostly within the scale of 10 − − − tempo2 software package(Hobbs et al. 2006) to calculate the timing residuals basedon the aforementioned ephemeris, without fitting for anyparameters. When generating ToAs with the traditionalmethod, we dedispersed the data with DM values both re-tained from online coherent dedispersion and determined byfitting for a constant DM with multi-frequency ToAs fromthe same dataset.A brief summary of the timing results can be found inTable 1. It can been seen that application of the channelisedDFT method (without fitting for ∆DM) achieves the low-est weighted timing precision and improves the ToA mea-surement uncertainties by ∼
20% in median. Using the DMmeasured from the dataset itself (10 . − pc) leads to asimilar timing precision, while retaining the DM used for co- c (cid:13) , 000–000 K. Liu et al. -3-2-1 0 1 2 1.3 1.4 1.5 1.6 1.7 O ff s e t ( µ s ) Frequency (GHz) DM=10.3940DM=10.3916DM=10.3905
Figure 12.
Phase offsets of the observed profile on MJD 56592at different frequencies based on the DM values as in Table 1.The solid line represents the best-fit quadratic curve when usingDM = 10 . − pc. The value of 10 . − pc is the best-estimated DM with the offsets. herent dedispersion (10 . − pc) results in significantlyworse rms and more systematics. This can be explained bythe fact that a deviation of DM from the true value may in-duce additional profile variation across the observing bandand enhance the scintillation effect. For demonstration, inFig. 12 we plot the phase offsets of the observed profile onMJD 56592 at different frequencies, given the DM valuesas in Table 1. Clearly, the DM value corresponding to lowertiming rms gives significantly less differential phase and thusprofile variation between frequencies. It is also interesting tonotice, that after subtracting all quadratic components theoffsets are mostly around zero. This indicates that the in-trinsic profile variation may not be significant within our ob-serving bandwidth, unless it actually appears as a quadraticdrift against frequency.Applying the channelised DFT method with ∆DM fitresults in worse rms residuals, which is expected due to thehigh correlation between ∆DM and ∆ τ in the data as in-dicated in Fig. 1. Nevertheless, in this case including ∆DMinto the template-matching fit may not be necessary as DMvariations are not significant within the observing time. InFig. 13 we plot the ∆DM measurements at each epoch,achieved by both the channelised DFT method and ToAsfrom multiple frequencies as in e.g. Keith et al. (2013). Herewe divided the whole bandwidth into eight 64-MHz sub-bands to create multi-frequency ToAs. It can be seen thatthe two methods lead to consistent measurements, both invalues and uncertainties. Within the observing time there isalso no evidence of gradual DM variations above the detec-tion threshold. In this paper, we have described the channelised DFTmethod that can conduct ToA measurements based on tim-ing data with frequency information. The functionality ofthe channelised DFT method has been tested, and thepotential for improvements to timing precision has beendemonstrated using simulated data. Furthermore, we have -25-20-15-10-5 0 5 10 15 20-20 -15 -10 -5 0 5 10 15 20 25 ∆ D M ( - c m - p c ) MJD-565652D template-matchingToA fitting
Figure 13.
Measured ∆DM on each epoch with the two-dimensional template-matching technique and the normal methodusing multi-frequency ToAs. The measurements from template-matching have rms and uncertainty median of 4 . × − cm − pcand 3 . × − cm − pc, respectively, while those obtained basedon multi-frequency ToAs achieve 6 . × − cm − pc and 2 . × − cm − pc. applied the method to timing data of PSR J1909 − ∼
20% in median. Our approach has also beenshown to achieve measurements of ∆DM on epochs consis-tent with methods based on multi-frequency ToAs (as in e.g.Keith et al. 2013).It is known that the observed dispersion delay can de-viate from the expected ∝ f − law coming from the coldplasma assumption (e.g. Armstrong et al. 1995; Keith et al.2013). This phenomenon can potentially be taken into ac-count with the channelised DFT method, by simply includ-ing a theoretical model in the merit function, i.e., Eq. (8),and enabling more parameter fittings. For this purpose,broadband observations would be highly requested to breakthe degeneracy between fitted parameters and to achieveenough sensitivity to measure the deviation.Note that application of the channelised DFT methodmight not achieve optimal modelling of DM variations as themeasurements use information only from a single observingsession without any interpolations over epochs like in otherwork (Kaspi et al. 1994; Keith et al. 2013; Lee et al. 2014;Lentati et al. 2014). Nevertheless, if the variation timescaleis well above the intervals between observing sessions, themethod can potentially be extended to combine data fromneighbouring epochs and perform a global fit for a uniqueoffset in DM, so as to increase the accuracy of determina-tion. Besides, when processing broadband timing data, thecurrent interpolation approaches do not fully use the in-formation that data from different frequencies are from asimultaneous observation. Therefore, it may also be worthattempting to combine these two types of method, to si-multaneously obtain timing residuals and DM modelling di-rectly from timing data. c (cid:13) , 000–000 -D template-matching Table 1.
Statistical results of ToAs obtained by the traditional template-matching approach (1D) with different DM values for dedis-persion, and the channelised DFT method without (2D ) and with ∆DM fit (2D ). The DM of 10.3916 cm − pc was achieved by fittingfor a constant DM with multi-frequency ToAs from the selected NUPPI dataset.MJD 1D DM (cm − pc) 10.3940 10.3916 10.3940 10.3940Max. σ ToA (ns) 934 912 839 4094Min. σ ToA (ns) 39 39 30 220Median σ ToA (ns) 149 149 114 1014Weighted rms (ns) 610 275 247 786Reduced χ ACKNOWLEDGEMENTS
We thank D. Stinebring for help with simulations of the dy-namic spectra. We are also grateful to the anonymous ref-eree who provided constructive suggestions to significantlyimprove the paper. K. Liu is supported by the ERC Ad-vanced Grant “LEAP”, Grant Agreement Number 227947(PI M. Kramer). The Nan¸cay Radio Observatory is oper-ated by the Paris Observatory, associated with the FrenchCentre National de la Recherche Scientifique.
REFERENCES
Antoniadis J. et al., 2013, Science, 340, 448Armstrong J. W., Rickett B. J., Spangler S. R., 1995, ApJ,443, 209Bhat N. D. R., Cordes J. M., Camilo F., Nice D. J., LorimerD. R., 2004, ApJ, 605, 759Cognard I., Theureau G., 2006, in IAU Joint Discussion,Vol. 2, IAU Joint DiscussionCognard I., Theureau G., Guillemot L., Liu K., Lassus A.,Desvignes G., 2013, in Cambresy L., Martins F., Nuss E.,Palacios A., ed, SF2A-2013: Proceedings of the Annualmeeting of the French Society of Astronomy and Astro-physics, p. 327Cordes J. M., Lazio T. J. W., 2002, astro-ph/0207156Cordes J. M., Shannon R. M., 2010, astro-ph/1107.3086Downs G. S., Reichley P. E., 1983, ApJS, 53, 169Foster R. S., Backer D. C., 1990, ApJ, 361, 300Freire P. C. C. et al., 2012, MNRAS, 423, 3328Hellings R. W., Downs G. S., 1983, ApJ, 265, L39Hobbs G. B., Edwards R. T., Manchester R. N., 2006, MN-RAS, 369, 655Hotan A. W., Bailes M., Ord S. M., 2005, MNRAS, 362,1267Hotan A. W., van Straten W., Manchester R. N., 2004,Proc. Astr. Soc. Aust., 21, 302Karuppusamy R. et al., 2014, in prep.Kaspi V. M., Taylor J. H., Ryba M., 1994, ApJ, 428, 713Keith M. J. et al., 2013, MNRAS, 429, 2161Kramer M., 1994, A&AS, 107, 527Kramer M., Xilouris K. M., Camilo F., Nice D., Lange C.,Backer D. C., Doroshenko O., 1999, ApJ, 520, 324Lee K. J. et al., 2014, MNRAS in press, astro-ph/1404.2084Lentati L., Alexander P., Hobson M. P., Feroz F., vanHaasteren R., Lee K. J., Shannon R. M., 2014, MNRAS,437, 3004Liu K., 2012, Ph.D. thesis, University of Manchester Liu K., Verbiest J. P. W., Kramer M., Stappers B. W., vanStraten W., Cordes J. M., 2011, MNRAS, 1442Lorimer D. R., Kramer M., 2005, Handbook of Pulsar As-tronomy. Cambridge University PressOrd S. M., van Straten W., Hotan A. W., Bailes M., 2004,MNRAS, 352, 804Os lowski S., van Straten W., Hobbs G. B., Bailes M., De-morest P., 2011, MNRAS, 418, 1258Pennucci T. T., Demorest P. B., Ransom S. M., 2014, astro-ph/1402.1672Press W. H., Teukolsky S. A., Vetterling W. T., FlanneryB. P., 1992, Numerical Recipes: The Art of Scientific Com-puting, 2 nd edition. Cambridge University Press, Cam-bridgeRansom S. M., Demorest P., Ford J., McCullough R., RayJ., DuPlain R., Brandt P., 2009, in American Astronomi-cal Society Meeting Abstracts, Vol. 214, American Astro-nomical Society Meeting Abstracts 214, p. 605.08Stappers B. W. et al., 2011, A&AS, 530, A80Taylor J. H., 1992, Phil. Trans. Roy. Soc. A, 341, 117van Straten W., 2006, ApJ, 642, 1004 APPENDIX A: DERIVATIVE CALCULATIONSFOR L-M ROUTINE
The L-M data modelling routine requires the input of agradient of χ and the curvature matrix. The calculationwas done in Press et al. (1992) for a model of real num-bers in a one-dimensional array, and here we derive the ex-pressions concerning complex numbers and two dimensions.With the measurement y j,k = m j,k + n j,k i , and the model y ( x j,k ,~a ) = m ( x j,k ,~a ) + n ( x j,k ,~a ) i , we have χ ( ~a ) = X j,k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m ( x j,k ,~a ) + n ( x j,k ,~a ) i − m j,k − n j,k iσ j,k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = X j,k σ j,k [ m ( x j,k ,~a ) + m j,k − m ( x j,k ,~a ) m j,k + n ( x j,k ,~a ) + n j,k − n ( x j,k ,~a ) n j,k ] , (A1)where ~a is an array containing the fitted parameters and i = √−
1. Hence, the gradient of χ is expressed by ∂χ ∂α p = X j,k σ j,k (cid:20) ( m ( x j,k ,~a ) − m j,k ) ∂m ( x j,k ,~a ) ∂α p + ( n ( x j,k ,~a ) − n j,k ) ∂n ( x j,k ,~a ) ∂α p (cid:21) , (A2) c (cid:13) , 000–000 K. Liu et al. and the second derivative matrix (Hessian matrix) is writtenas ∂χ ∂α p ∂α q = X j,k σ j,k [ m ( x j,k ,~a ) ∂m ( x j,k ,~a ) ∂α p ∂α q + n ( x j,k ,~a ) ∂n ( x j,k ,~a ) ∂α p ∂α q + ∂m ( x j,k ,~a ) ∂α p ∂m ( x j,k ,~a ) ∂α q + ∂n ( x j,k ,~a ) ∂α p ∂n ( x j,k ,~a ) ∂α q − m j,k ∂m ( x j,k ,~a ) ∂α p ∂α q − n j,k ∂n ( x j,k ,~a ) ∂α p ∂α q ] . (A3)The actual fitting routine uses the curvature matrix definedas α pq ≡ ∂χ ∂α p ∂α q . After eliminating the second derivativeterms so as to stabilise the iterations (explained in Chap-ter 15.5, Press et al. 1992), we then have α pq = X j,k σ j,k (cid:20) ∂m ( x j,k ,~a ) ∂α p ∂m ( x j,k ,~a ) ∂α q + ∂n ( x j,k ,~a ) ∂α p ∂n ( x j,k ,~a ) ∂α q (cid:21) . (A4)Note that in our model ~a = ( b , ..., b N b , ∆ τ, ∆DM) , (A5) m j,k = P j,k cos θ j,k , (A6) n j,k = P j,k sin θ j,k , (A7) m ( x j,k ,~a ) = b j S j,k cos( φ j,k + τ j ) , (A8) n ( x j,k ,~a ) = b j S j,k sin( φ j,k + τ j ) , (A9)then χ becomes χ = X j,k P j,k + b j S j,k − b j P j,k S j,k cos ( φ j,k − θ j,k + kτ j ) σ j,k , (A10)and the gradient of χ is written as ∂χ ∂α p = X i b p S p,k − P p,k S p,k cos( φ p,k − θ p,k + kτ p ) σ p,k , p N b X j,k kb j P j,k S j,k sin( φ j,k − θ j,k + kτ j ) σ j,k , p = N b + 1 X j,k kD b j P j,k S j,k sin( φ j,k − θ j,k + kτ j ) σ j,k f j , p = N b + 2(A11)and the curvature matrix is α pq = X k S p,k σ l , p, q N b , p = q , p, q N b , p = q , p N b , N b q N b + 2 X j,k k b j S j,k σ j , p = q = N b + 1 X j,k k D b j S j,k σ j f j , p = N b + 1 , q = N b + 2 X j,k (cid:18) kD b j S j,k σ j f j (cid:19) , p = q = N b + 2 . (A12) c (cid:13)000