An algorithm for determining the rotation count of pulsars
MMNRAS , 1–13 (2015) Preprint 15 October 2018 Compiled using MNRAS L A TEX style file v3.0
An algorithm for determining the rotation count of pulsars
Paulo C. C. Freire (cid:63) and Alessandro Ridolfi Max-Planck-Institut f¨ur Radioastronomie, Auf dem H¨ugel 69, D-53121 Bonn, Germany
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present here a simple, systematic method for determining the correct global ro-tation count of a radio pulsar; an essential step for the derivation of an accuratephase-coherent ephemeris. We then build on this method by developing a new algo-rithm for determining the global rotational count for pulsars with sparse timing datasets. This makes it possible to obtain phase-coherent ephemerides for pulsars for whichthis has been impossible until now. As an example, we do this for PSR J0024 − Key words: methods: data analysis – (stars:) pulsars: general – (stars:) pulsars:individual (PSR J0024 − In the 50 years since the discovery of radio pulsars, theirstudy has been a major scientific success story. They are ex-tremely versatile and powerful tools for studying fundamen-tal physics, in particular the study of gravity, the fundamen-tal properties of space-time and of gravitational waves (seeWex 2014; Berti et al. 2015 and the many references therein).Measurements of large neutron star masses (e.g., Antoniadiset al. 2013; Fonseca et al. 2016) had profound implicationsfor the study of cold nuclear matter at supra-nuclear densi-ties (see ¨Ozel & Freire 2016 and references therein), a fun-damental question in nuclear physics. The orbital propertiesand abundance of double neutron star systems provided thefirst estimates of the rate of NS-NS mergers, the first “guar-anteed” source of events for ground-based gravitational wavedetectors. Furthermore pulsars can be used to detect verylow frequency gravitational waves directly (The NANOGravCollaboration et al. 2015; Desvignes et al. 2016; Reardonet al. 2016; Verbiest et al. 2016).Apart from these “fundamental physics” applications,radio pulsars are also superb astrophysical tools: they al-low a much deeper understanding of the late stages of theevolution of massive stars, including the supernova events(important for understanding the origin of the most usefulpulsars, the recycled pulsars, see e.g., Lorimer 2008; Tauriset al. 2017), the dynamics and history of globular clustersand the ionized interstellar medium. (cid:63)
E-mail: [email protected]
Most of these applications rely on a single, simple technique, pulsar timing . This technique achieves its full power for re-cycled pulsars. Put simply, this consists of the study of thetimes of arrival (TOAs) of the pulses (normally determinedby adding many individual pulses coherently in phase) atone or several telescopes, most often radio telescopes. TheseTOAs can be determined very accurately, in some cases withprecisions better than 100 ns, but more commonly a few µ s;they correspond to times where a particular longitude ofthe pulsar (normally close to the spin phase where the ob-served radio emission is at a maximum) is aligned towardsthe Earth. Therefore, in between any two radio pulses thepulsar has rotated an integer number of times.A good timing solution (henceforth a phase-coherentephemeris ) must be capable of accurately predicting theTOAs. This ephemeris consists of a specified mathematicaldescription with a few crucial free parameters that describethe spin of the pulsar plus the transformation between thereference frame of the pulsar and that of the receiving tele-scope. In the reference frame of most recycled pulsars, theirrotation can be described by a spin frequency ( ν ), plus asmall, but constant and negative spin frequency derivative( (cid:219) ν ). In the reference frame of the receiving radio telescope,these arrival times are affected by the motion of the pulsarrelative to the radio telescope. This has several components:the observatory’s motion caused by the Earth’s rotation, themotion of the Earth relative to the Solar System barycentre(SSB), and the motion of the pulsar relative to the SSB,which might be affected by the pulsar’s own orbital motion(if it happens to have any companion(s)). In order to correct © a r X i v : . [ a s t r o - ph . I M ] F e b Paulo C. C. Freire for the movement of the radio telescope relative to the SSB,we need to have an estimate of the position of the pulsarin the sky; i.e., the right ascension ( α ), declination ( δ ) and,for measurements spanning many years, the proper motionin these two coordinates ( µ α , µ δ ). Finally, if the pulsar is ina binary system, then its orbital motion can be normallyparameterized by 5 Keplerian parameters: in the Damour &Deruelle (1986) model (known as the DD model), for exam-ple, these parameters are the orbital period ( P b ), the semi-major axis of the pulsar’s orbit projected along the line ofsight ( x ), the orbital eccentricity ( e ), the longitude of peri-astron ( ω ) and the time of passage through periastron ( T ).In some cases, a few additional “post-Keplerian” parameterscan be measured, these are caused by geometric (Kopeikin1996) and relativistic (e.g., Damour & Taylor 1992) effects.In the case of “black widow” or “redback” pulsars, wherethere are Newtonian perturbations to the orbit, additionalparameters might be necessary (e.g., Shaifullah et al. 2016).It is clear from the number of parameters (and their mea-surement precision) that these phase-coherent ephemeridesprovide a wealth of scientific information.These ephemerides are derived from the TOAs using atiming program, like tempo , tempo2 (Hobbs et al. 2006;Edwards et al. 2006) or PINT ; in what follows we will beusing tempo because it allows a very simple implementationof the methods to be described. For each TOA, tempo willfirst correct it using tables where the local time standard atthe radio telescope is compared to a more stable time scale,e.g., the universal coordinated time (UTC). Then, the posi-tion of the radio telescope relative to the SSB is calculated,first by using Earth rotation tables to calculate where theradio telescope is relative to the Earth’s centre, then by us-ing a Solar System ephemeris (like DE 421, Folkner et al.2008 or DE 430, Folkner et al. 2014) to translate the latterto the SSB. This vector is then projected along the direc-tion to the pulsar (derived from α, δ, µ α and µ δ ) to calculatethe effect of the motion of the radio telescope relative tothe SSB on the TOAs. Then, if the pulsar is in a binary,a binary model (like the aformentioned DD model) is usedto subtract the time delays caused by the orbital motion. Itis only after this stage that tempo calculates the residuals :these are the TOAs minus the prediction of the model forthe corresponding rotation of the pulsar . The best-fit timingparameters are determined by varying the timing parame-ters in a way that minimizes the sum of the squares of theresiduals. In the last paragraph, we have placed the emphasis on onefundamental point, which is rarely discussed in the litera-ture: tempo can only determine a phase-coherent ephemerisif it is comparing a TOA with the model estimate for thatpulse, not the one before or after. This can only be achieved after the correct rotation count between any two TOAs (aset of integers) has been established. This is what we call a global rotation count .As we will see, when a new pulsar is discovered, the http://tempo.sourceforge.net/ https://github.com/nanograv/PINT rotation count within a single observation is generally wellknown. However, the rotation count between successive ob-servations is generally not known. As we will show, this issueis easy to solve for most pulsars, but it is hard to solve ifthere are no closely spaced detections of the pulsar, as in thecase of scintillating pulsars, or pulsars in eclipsing systems.In the case of very faint pulsars, the timing precision mightbe too poor for the determination of the rotation count evenfor closely spaced observations.The number of known faint pulsars has been growingsignificantly, but the time available at the major radio tele-scopes that can detect them remains unchanged, this impliesautomatically less time available to follow up each new pul-sar. This can be partly compensated by more sensitivity andlarger bandwidths, but also a more careful coordination ofobservations, and by dropping many pulsars that are scien-tifically less rewarding.However, as we will see below, even if a pulsar is seento have a high scientific value and its observation becomesa priority at a particular radio telescope, finding the correctrotation count might still be extremely difficult to achieve.Frequent eclipses or the effects of scintillation mean thateven a dense set of observations does not necessarily trans-late into a dense set of detections with a minimal usefulsignal-to-noise ratio (S/N). Thus, in some cases, a phase-coherent ephemeris has not been obtained to date. In Section 2 we present the observations and measurementsof TOAs for a recently discovered pulsar, PSR J0024 − tempo as a practical toolfor deriving results: we believe that the use of a this specificprogram on a specific example is helpful for the illustrationof the concepts presented in the paper. In Section 3, wepresent the concept of global rotation count in more detailand describe a simple, standard method for achieving it thatcan be applied to most pulsars; this assumes familiarity withthe tempo timing package. This standard method can de-termine global rotation counts for most pulsars but it cannotachieve this for 47 Tuc aa owing to the extreme sparsity ofits detections. In Section 4, we highlight a method for de-termining the global rotation count for such pulsars with aminimal number of tempo iterations, and present its imple-mentation as a UNIX shell script. In Section 5, we use thismethod to determine the rotation count of 47 Tuc aa, pre-senting a previously unavailable phase-coherent ephemerisfor this pulsar and discussing the significance of its param-eters. Finally, in Section 6, we summarize our results andhighlight some of the possibilities opened up by the methoddiscussed here.Throughout this work, we will assume basic familiar-ity with the concepts of pulsar signal analysis, in particulardedispersion, folding and the use of low-noise pulse profiletemplates to estimate TOAs; for a review see Lorimer &Kramer (2004). We will introduce some concepts as we goalong, whenever this happens the new concept is presentedin italic. MNRAS , 1–13 (2015) etermining the rotation count of pulsars The pulsar discussed in this work, 47 Tuc aa, is one of the25 radio millisecond pulsars (MSPs) known in the globu-lar cluster NGC 104, also known as 47 Tucanae (henceforth47 Tuc). All these pulsars are detected in the same data (thetelescope beam covers all their positions simultaneously),obtained during 519 observations of 47 Tuc carried out from1997 until 2013 (i.e., a time span of 16 years) with the 64-m Parkes radio telescope in New South Wales, Australia.These observations are part of a long-term project dedi-cated to these pulsars. The project’s scientific motivation,setup of the observations, data taking (with the AnalogueFilterbank, or AFB) and reduction and some of its resultsare described by Ridolfi et al. (2016), Freire et al. (2017) andtheir references.47 Tuc aa was discovered in this data set by Pan et al.(2016). With a spin period of only 1.845 ms, it is the fastest-spinning pulsar known in 47 Tuc. It also has the highest dis-persion measure (DM) for any known pulsar in the cluster, . ( ) pc cm − . The detections of this object are extremelyrare: in the 519 observations of the data set, the pulsar wasonly detected (in the sense of yielding at least one usableTOA) on 22 occasions; this represents a detection rate of4.2 per cent, and an average of 1.38 detections per year.In most of these detections the signal-to-noise ratio (S/N)is smaller than 6. Given this extreme sparsity, no phase-coherent ephemeris could be derived by Pan et al. (2016).The reason for the small number of detections is thesmall flux density of the pulsar combined with diffractivescintillation: like most pulsars in 47 Tuc, 47 Tuc aa has a fluxdensity that is well below the telescope’s sensitivity limit.Most of these pulsars are only detectable when scintillationamplifies their signal (see e.g., Camilo et al. 2000). 47 Tuc aais so faint that only very rarely is the scintillation amplifi-cation large enough to make the pulsar detectable, even inobservations lasting 8 hours. This means that although wecan choose when to observe the cluster and the pulsars in it,we cannot choose when to detect any particular pulsar. Thisis also one of the reasons why, despite the lack of improve-ment of the sensitivity, new pulsars are still being discoveredin this data set: for instance, 47 Tuc aa is not detected onthe vast majority ( ∼ ) of observations, and other evenfainter pulsars might be appearing even more infrequently.Although no phase-coherent ephemeris was presentedfor 47 Tuc aa in Pan et al. (2016), the few detections of thisobject show that, after correction for the Doppler shift dueto the Earth’s motion, its spin period does not change mea-surably with time. This implies that, unlike most MSPs in47 Tuc, this particular object is not part of a binary system.Apart from this, not much was known about it, in particularits location relative to the centre of the cluster, its propermotion or its spin period derivative. These parameters areexquisitely well determined for the other 22 MSPs in 47 Tucwith phase-coherent ephemerides (Ridolfi et al. 2016; Freireet al. 2017). Figure 1.
Best detection of 47 Tuc aa from an observation madeon 2007 August 2. In the main greyscale plot, we display the radiointensity at 20 cm (darker meaning larger intensity) as a functionof the spin phase (displayed twice for clarity) and time, which isdivided into 60 segments. As we can see, the pulsar signal (thevertical gray line) does not drift perceptibly in phase, indicatingthat the spin period used to fold these data is accurate. The topplot represents the integrated radio intensity (vertical axis) as afunction of spin phase, also displayed twice for clarity.
In Fig. 1, we display one of the best detections of 47 Tuc aa,from an observation on MJD = 54314 (2007 August 2). Inthis plot, we can see that the pulsations appear at a con-stant spin phase. This shows conclusively that the estimateof the spin frequency ( ν = . Hz), determined fromthe search procedure (Pan et al. 2016), is precise enough topredict the arrival times of the pulses within that observa-tion. If the spin period were in error, we should see a driftof the phase with time.This can be confirmed with tempo in the following way:(i) Group the 60 sub-integrations displayed in Fig. 1 into,for example, 6 longer sub-integrations and, using the tech-niques described in Ridolfi et al. (2016), make a single TOA
MNRAS000
MNRAS000 , 1–13 (2015)
Paulo C. C. Freire for each sub-integration. A TOA list suitable for use in tempo should look like this : MODE 17 1390.000 54314.7954878129250 16.660 0.000007 1390.000 54314.8211805391308 17.735 0.000007 1390.000 54314.8468749751720 15.921 0.000007 1390.000 54314.8725694346731 21.125 0.000007 1390.000 54314.8982638745009 19.083 0.000007 1390.000 54314.9222205617143 19.950 0.00000
The “MODE 1” flag indicates that tempo will make aweighted fit, with the weight of each TOA given by the in-verse of its uncertainty. The first column indicates the radiotelescope where the data were taken (in this case 7 identifiesParkes), the second column the central radio frequency inMHz, the third column is the TOA itself in MJD, the fourthcolumn is the TOA uncertainty in µ s and the last one indi-cates any previous DM corrections (none in this case).(ii) Make a simple ephemeris file (“47Tucaa.par”) contain-ing the best estimate of the spin frequency (we start fromthe estimate in Pan et al. 2016), no spin period derivative(this is generally very small), and the coordinates of the cen-tre of 47 Tuc (the pulsar is not likely to be more than . (cid:48) away from that centre, see Freire et al. 2017, except in onecase, 47 Tuc X, see Ridolfi et al. 2016) and the proper mo-tion of the cluster ( idem ). We also put in a parallax derivedfrom the known cluster distance (4.69 kpc, see discussion inFreire et al. 2017 and references therein); this has little ef-fect on the timing but it is a known quantity anyway. Thisephemeris looks like this: PSR J0024-7205AARAJ 00:24:05.67DECJ -72:04:52.62PMRA 5.16PMDEC -2.85PX 0.2132F0 541.893656 1F1 0PEPOCH 51600DM 24.971EPHEM DE421CLK UTC(NIST)UNITS TDBNITS 1
The flag “1” after F0 (the spin frequency, ν ) means thatwe are fitting for this quantity. The reference epoch (givenby the “PEPOCH” flag) should generally be the barycentrictime at which ν was measured. For a MSP with a smallvariation of the spin period, this is not so crucial, so we setPEPOCH to MJD = 51600 (2000 February 26) because it isthe reference epoch used for all other ephemerides in Freireet al. (2017).With these files, we then run tempo and look at the resid-uals. The pre-fit and post-fit residuals plots are displayed inFig. 2; they have reduced χ of 4.24 and 1.06 respectively.In that Figure, tempo confirms what was already ob-vious from Fig. 1: the ephemeris “47Tucaa.par” predicts thecorrect rotation count within the observation: all TOAs arearriving within a very small phase window (from 0.0, de-fined automatically as the phase of the first TOA, to 0.03, arange that is of the same order as the TOA uncertainties).This means that the rotation count (5933580) between the The program accepts a variety of formats. In this case we choseuse the “Princeton” format because of its simplicity.
Figure 2.
Pre-fit (top) and post-fit (bottom) residuals for the6 TOAs derived from the best observation of 47 Tuc aa, takenon 2007 August 2. The horizontal axis represents time and thevertical axis represents the residual, in spin phase units. 0.01 ofthe spin phase corresponds to 18.5 µ s, which is comparable to theTOA uncertainties. first and last TOAs from the barycentric spin frequency es-timated from the initial spin frequency is correct, i.e., thisgroup of TOAs is phase-connected (henceforth connected ).After fitting for ν , the situation improves further, with allresiduals even closer to 0.From this set of TOAs, we can already have a very pre-cise measurement of the spin frequency (or period, P ≡ / ν ).According to the tempo output (printed in file tempo.lis)the time elapsed between the first and last TOA is 3.04hours, more precisely δ T t = δ T b = P ) of0.0018453805301 s, exactly the value returned by tempo .The uncertainty of the spin period δ P should be the uncer-tainty of δ T b ( µ s ) divided by the number of rotations, i.e., δ P = . × − s . The uncertainty estimated by tempo isactually smaller ( . × − s ) since it does not only takeinto account the first and last TOAs, but the other 4 in-termediate ones as well. The parameters from this fit arepresented in row 1 of Table 1.For the remainder of this paper it is important to notice MNRAS , 1–13 (2015) etermining the rotation count of pulsars that the absolute phase of the pulses is immaterial, since tempo automatically assigns a spin phase of zero to the firstTOA. All that matters is the rotation count between TOAs.By minimizing the residual rms, all that tempo did was toadjust the spin frequency in order to remove the residualslope from this single observation. Such residual slopes areessential because they give us the initial constraints on theglobal rotation count and the ephemeris to be derived fromit. Although the rotation count within an observation is, as wehave seen, generally well known, this is not the case for therotation count between more widely spaced TOAs.We can see this by repeating the procedure in the pre-vious section for all available observations:(i) Make TOAs for all observations where the pulsar hasappeared. Whenever possible, make at least 3 TOAs - thisis important for obtaining residual slopes from each day’sdata.(ii) Use the same ephemeris as in the previous step,(iii) Run tempo and look at the residuals.For the ephemeris listed above, the pre-fit residuals aredisplayed in Fig. 3. These are scattered evenly through thewhole spin phase, from − . to . , with no discerniblepattern. This means that our first ephemeris is not pre-cise enough to predict the rotation count between all theseTOAs, otherwise all the residuals would be confined to anarrow range of phases, as in Fig. 2. Thus, the startingephemeris is a partial ephemeris: it is certainly not phase-coherent.This lack of precision stems from a variety of factors:a) Lack of precision of the spin frequency ν (this is preciseenough to predict TOAs for a few hours, but probably notenough to do so for many years) b) the possibility that thesky position of the pulsar is slightly different from the po-sition of the centre of 47 Tuc; this introduces small delaysdue to the Earth’s motion that are nonetheless larger thanone rotation of the pulsar and c) the unknown variation of ν with time. In what follows, we will assume at first that thiscan be described by a single spin frequency derivative, (cid:219) ν .In both plots of Fig. 3, we see that whenever multipleTOAs are derived from successive sub-integrations withinthe same observation, their residuals appear in (nearly) hor-izontal groups, except for two cases that appear near phase0.5 that will be discussed below. This means that thosegroups of TOAs are connected (see Section 3.1). The timeinterval between two such groups is a gap ; thus by definitionthe rotation count within a gap is not necessarily known atthis stage. The first step toward the determination of the rotation countis to make a more precise estimate of the spin frequencyand other parameters using all connected TOA groups, rightfrom the start. This has the great advantage that, for everystep of the work, the estimates of all timing parameters areconstrained by the residuals slopes from all connected TOAgroups. Since the initial ephemeris does not provide a reliablerotation count across all gaps, we cannot assume to knowit. We can remove this assumption in tempo by fitting anarbitrary time offset for each group of connected TOAs. Thiscan be achieved – in tempo – by bracketing the TOAs fromall observations but one (in this case the last) with JUMPstatements, as exemplified here:
MODE 1JUMP7 1390.000 51413.6357638815603 18.173 0.000007 1390.000 51413.6499999905185 23.038 0.00000JUMPJUMP7 1390.000 51413.6895833852018 14.824 0.000007 1390.000 51413.7381944448891 14.797 0.000007 1390.000 51413.7868055551721 18.587 0.000007 1390.000 51413.8354166713418 22.973 0.00000JUMPJUMP7 1390.000 51490.5017355131067 24.352 0.00000JUMP...(some TOAs not shown here)...C DON’T BRACKET WITH JUMPS7 1390.000 54816.4975696298174 31.740 0.00000
Then, as before, run tempo . The pre-fit residuals arethe same as displayed in the top of Fig 3.Running tempo at this stage, we obtain the residualsdisplayed at the bottom of Fig. 3. These have a very highreduced χ of 225.15. The reason for this (and this is some-thing we should always beware during this process) is thatthere are groups of residuals close to the spin phase of 0.5:individual residuals in such groups might appear at a spinphases near 0.5 and − . . This is happening for the two pairsof TOAs called “a” (at MJD = 51492) and “b” (at MJD= 52531) respectively. In this situation, the wrong rotationcount is being assumed within those observations, so tempo cannot produce a good estimate of the parameters. In orderto correct this, one should introduce or subtract (whicheveris suitable) an extra rotation in between the TOAs in eachgroup, so that they appear with the same phase in the pre-fit residual plot. This can be done with a PHASE +1/ − (...)JUMP7 1390.000 51492.5850727545785 31.229 0.00000PHASE -17 1390.000 51492.6440955080127 32.860 0.00000JUMPC OptionalPHASE +1(...) and (...)JUMP7 1390.000 52531.6708365991755 24.719 0.00000PHASE -17 1390.000 52531.7347222220090 20.900 0.000007 1390.000 52531.7968733372661 17.472 0.00000JUMPC OptionalPHASE +1(...) MNRAS000
Then, as before, run tempo . The pre-fit residuals arethe same as displayed in the top of Fig 3.Running tempo at this stage, we obtain the residualsdisplayed at the bottom of Fig. 3. These have a very highreduced χ of 225.15. The reason for this (and this is some-thing we should always beware during this process) is thatthere are groups of residuals close to the spin phase of 0.5:individual residuals in such groups might appear at a spinphases near 0.5 and − . . This is happening for the two pairsof TOAs called “a” (at MJD = 51492) and “b” (at MJD= 52531) respectively. In this situation, the wrong rotationcount is being assumed within those observations, so tempo cannot produce a good estimate of the parameters. In orderto correct this, one should introduce or subtract (whicheveris suitable) an extra rotation in between the TOAs in eachgroup, so that they appear with the same phase in the pre-fit residual plot. This can be done with a PHASE +1/ − (...)JUMP7 1390.000 51492.5850727545785 31.229 0.00000PHASE -17 1390.000 51492.6440955080127 32.860 0.00000JUMPC OptionalPHASE +1(...) and (...)JUMP7 1390.000 52531.6708365991755 24.719 0.00000PHASE -17 1390.000 52531.7347222220090 20.900 0.000007 1390.000 52531.7968733372661 17.472 0.00000JUMPC OptionalPHASE +1(...) MNRAS000 , 1–13 (2015)
Paulo C. C. Freire
Table 1.
Post-fit parameters for several iterations described in the text. *Parameter not being fit, derived from initial assumptions (seetext).Iteration ν (cid:219) ν α δ Parameters TOAsnumber (Hz) ( − Hz s − ) (hh:mm:ss.s) ( ◦ : (cid:48) : (cid:48)(cid:48) ) fitted fitted1 541.8936548(11) 0* 00:24 05.67* − ν − ν Initial groups + JUMPS3 541.89365651(5) 0* 00:24:07.1(8) − ν, α, δ Initial groups + JUMPS4 541.8936564(6) 12(6) 00:24:07.0(8) − ν, α, δ, (cid:219) ν Initial groups + JUMPS5 541.8936568(6) 4(4) 00:24:07.4(8) − ν, α, δ, (cid:219) ν − ν, α, δ, (cid:219) ν, µ α , µ δ All connected
The optional PHASE +1 statements were introduced inboth cases so that all residual phases continued being dis-played in the interval from − tempo converge on an ephemeriswith a reduced χ close to 1. The pre and post-fit residu-als as a function of epoch obtained after this correction arepresented in Fig. 4.In the lower plot of Fig. 4, all post-fit residuals appearwithin a narrow window of spin phases. This is mostly (butnot entirely) a consequence of the fact that we have fittedtime offsets for each group of connected TOAs. Solely re-moving these offsets does not eliminate the residual slopesfor each group of connected TOAs; this is done by fittingthe timing parameters.Let us now exemplify this. If we fit only for ν , we obtainthe spin frequency (presented in the second row in Table 1)that is more than twice as precise as the spin period estimateobtained from the best observation in Section 3.1. However,this ephemeris produced residuals with a reduced χ of 3.79– still with non-zero residual slopes within each day – so wemight need to fit for (an)other parameter(s).Because the total time span is larger than one year, it ispossible to fit for ν , α and δ , as the covariance between theseparameters can be broken. If we do so, we obtain residualswith a reduced χ of 1.54, which is a major improvement.The solution we obtain then is presented in the third row ofTable 1. We can see here that, although the spin frequencyis not quite as precise as in the previous column, we canalready pinpoint the location of the pulsar with a precisionof a few arcseconds, which is a major improvement over theprevious (and, as we can now see, correct) assumption thatthe pulsar is near the centre of 47 Tuc.If we fit for ν , α and δ and (cid:219) ν , then the reduced χ decreases to 1.40. The resulting best ephemeris is presentedin the fourth row of Table 1. Although the resultant (cid:219) P ≡− (cid:219) ν / ν is not very significant (− . ± . ) × − s s − , it iswell within the range of what one finds for other pulsars in47 Tuc. This implies that this parameter should be fit atthis stage. Our starting ephemeris file is then the same as inSection 3.1, only with a few extra fitting flags (1’s) addedafter α , δ and (cid:219) ν lines.To summarize: although we have not determined therotation count for a single gap, by using the JUMP state-ments we can use the information contained in the residualslopes within all connected TOA groups to derive a reason-ably precise position, period, and even a strong constraint on (cid:219) ν . These estimates cannot be derived from any individualgroup of TOAs. We now arrive at the most important step on the path toget the correct global rotation count for the pulsar. We lookat the TOA list and find closely spaced groups of connectedTOAs, i.e., short gaps. As an example, for the best detectionof 47 Tuc aa, there is an observation earlier that day wherethe pulsar was also detected. Here are the TOAs: (...)JUMP7 1390.000 54314.7031249750801 15.335 0.000007 1390.000 54314.7515008653730 25.434 0.00000JUMPJUMP7 1390.000 54314.7954878129250 16.660 0.000007 1390.000 54314.8211805391308 17.735 0.000007 1390.000 54314.8468749751720 15.921 0.000007 1390.000 54314.8725694346731 21.125 0.000007 1390.000 54314.8982638745009 19.083 0.000007 1390.000 54314.9222205617143 19.950 0.00000JUMP(...)
The question is now: is our original ephemeris (presented inSection 3.1) precise enough to predict the rotation count forthis gap or not?To find out, we comment out the inner pair of JUMPstatements and introduce a PHASE +N statement, where Nis an integer that determines the correction to the rotationcount predicted by the ephemeris. In what follows, such inte-gers (or combinations of k integers when discussing k gaps)are referred to as k -gap solutions , a 1-gap solution will bereferred to simply as a solution . (...)JUMP7 1390.000 54314.7031249750801 15.335 0.000007 1390.000 54314.7515008653730 25.434 0.00000C JUMPPHASE +0C JUMP7 1390.000 54314.7954878129250 16.660 0.000007 1390.000 54314.8211805391308 17.735 0.000007 1390.000 54314.8468749751720 15.921 0.000007 1390.000 54314.8725694346731 21.125 0.000007 1390.000 54314.8982638745009 19.083 0.000007 1390.000 54314.9222205617143 19.950 0.00000JUMP(...) Now we run tempo , still with the original ephemeris.
MNRAS , 1–13 (2015) etermining the rotation count of pulsars aa bbaa bb Figure 3.
Pre-fit (top) and post-fit (bottom) residuals for allTOAs for 47 Tuc aa as a function of epoch. The pre-fit residu-als are scattered evenly across the full spin phase space (between − . and . ), i.e., our starting ephemeris is not precise enoughto predict the phase even approximately. However, for TOAs de-rived from a single observation (i.e., those where we know thereis connection), the residuals appear in nearly horizontal groups,this means that the predicted spin period for each observation isaccurate. If a group of residuals appears near phase 0.5, then someresiduals in that group might wrap around and appear near phase − χ is 114.9. Assigning the integers − , − , − , 0, 1, 2 and 3 to N, weobtain the following values for the reduced χ : 2063.0, 914.1,226.9, 1.37, 237.5, 935.3 and 2094.8 (see Fig. 5).The quality of the fit is much higher for N = 0 (with areduced χ of 1.37) than for the second best fit ( − , wherethe reduced χ = . ). In such cases we can say with greatconfidence that there is a unique acceptable solution , “0”.This means that the rotation count of this gap is unambigu-ous , and in this case identical to the rotation count predictedby the initial ephemeris. Figure 4.
Pre-fit and post-fit residuals for all TOAs for 47 Tuc aaas a function of epoch, now with the correct rotation count forepochs “a” and “b”. Although the pre-fit residuals are scatteredevenly across the full spin phase space (between − . and . ),the post-fit residuals all appear near phase zero, with a reduced χ near 1.9. This is partially a consequence of the fact that weare fitting an arbitrary time offset to each observation. Flatteningthe residuals slopes within individual observations constrains thetiming parameters. In this situation we can eliminate the gap, i.e., declarethe groups on both sides to be mutually connected (in whatfollows, we simply describe this as connecting the gap ). Be-cause the gap has been connected, there is no need to rein-troduce its JUMP statements in TOA list, but whateverPHASE +N statement worked best should stay. This con-nection results in more tightly constrained timing parame-ters.In Fig. 5, we can see that the relation between the re-duced χ and N follows very closely a parabola, which isdepicted by the solid red line. In the sections that follow, wemake use of this fact. MNRAS000
Pre-fit and post-fit residuals for all TOAs for 47 Tuc aaas a function of epoch, now with the correct rotation count forepochs “a” and “b”. Although the pre-fit residuals are scatteredevenly across the full spin phase space (between − . and . ),the post-fit residuals all appear near phase zero, with a reduced χ near 1.9. This is partially a consequence of the fact that weare fitting an arbitrary time offset to each observation. Flatteningthe residuals slopes within individual observations constrains thetiming parameters. In this situation we can eliminate the gap, i.e., declarethe groups on both sides to be mutually connected (in whatfollows, we simply describe this as connecting the gap ). Be-cause the gap has been connected, there is no need to rein-troduce its JUMP statements in TOA list, but whateverPHASE +N statement worked best should stay. This con-nection results in more tightly constrained timing parame-ters.In Fig. 5, we can see that the relation between the re-duced χ and N follows very closely a parabola, which isdepicted by the solid red line. In the sections that follow, wemake use of this fact. MNRAS000 , 1–13 (2015)
Paulo C. C. Freire R edu c ed c h i Rotation number 0 500 1000 1500 2000 2500 3000 -3 -2 -1 0 1 2 3 R edu c ed c h i Rotation number
Figure 5.
Reduced χ as a function of the solution for the smallgap in MJD = 54314. The values derived by tempo (green crosses)follow closely a parabolic curve; the latter is represented by thesolid red line. After this stage, we look for other gaps shorter than 1 day.The two gaps we found (within MJD 51413 and 53164) alsohave unique solutions, in both cases 0, so both gaps canbe connected. Since we now have slightly larger connectedTOA groups, the post-fit timing parameters will be moreprecise. However, as we can see in the fifth row of Table 1,the improvement in this case is marginal.In order to proceed, we must try to find other gapswith unique solutions, since such gaps can be immediatelyconnected. The longer connected dataset will then result inan ephemeris with more tightly constrained parameters; thelatter will more likely make other gaps unambiguous. This isa runaway process, which means that, using this technique,we can, for most pulsars, determine unique solutions for allgaps, i.e., determine their global rotation count. This can bedone even if the detections appear to be very sparse at first– all that matters is that there is a good distribution of gaplengths .We will mention here a technical detail of this process,which will become apparent to anyone using this technique.Although the solutions for the shorter gaps will be relativelysmall integers, as the length of the gaps increases, the magni-tude of the solutions increases significantly. In order to miti-gate the problem, the initial ephemeris can be replaced withthe post-fit ephemeris calculated by tempo at any stage. Forinstance, we could have replaced the initial guess ephemerisin Section 3.1 with the better ephemeris derived in the pre-vious section. For this to work, all previous PHASE +Nstatements in the TOA list must be commented out, sincethe corresponding rotation count is already taken into ac-count by the new starting ephemeris. After such a change,however, we must always be careful about groups of TOAsat phases close to 0.5, as mentioned in Section 3.3.Another important point to keep in mind is that for thistechnique to succeed, three conditions must be met:a) The TOAs must be accurate. If one TOA is for somereason in error, it will derail the whole process. For this reason, we must be careful with the selection of the pulseprofiles while deriving TOAs.b) The TOAs should be consistent with each other (i.e.,preferably come from the same telescope and timing instru-ment), also the template used to derive them should be thesame andc) the pre-fit ephemeris must contain an appropriateset of timing parameters. For many isolated pulsars, α , δ , ν and (cid:219) ν are perfectly adequate; as we will see this is also thecase for 47 Tuc aa. For longer timescales, proper motion pa-rameters might be necessary. However, even among isolatedpulsars, particularly the younger ones, there are phenomenalike glitches and timing noise that can greatly complicatethe analysis.When sparsely detected pulsars are part of binary sys-tems (something that becomes immediately apparent fromthe variation of the spin period), we can make a prelimi-nary determination of the orbital parameters from the ob-served spin periods and their derivatives using the methodpresented by Freire et al. (2001a) ; this was used to findthe orbits of several pulsars in 47 Tuc. Those initial or-bital parameters can then be refined by fitting the observedDoppler-shifted spin period as a function of time to a Kep-lerian model (something that can be done using a programlike FITORBIT ).Phase-coherent ephemerides for such pulsars can be de-rived using the method we have just described above; how-ever, to achieve that we must exercise good judgment onthe binary parameters that must be fit. For a pulsar inwide, nearly circular orbit with a compact object, it is likelythat fitting for Keplerian parameters from the start will beenough to find phase connection. If the orbit is compact andeccentric, then fitting for the rate of advance of periastron( (cid:219) ω ) from the start might be necessary. In the case of eclips-ing binary pulsars, then unpredictable changes in the orbitalperiod must also be taken into account.Back to 47 Tuc aa, and before we continue, we note thatat this stage of the work the reduced χ of the residuals is1.40. Our previous experience timing other faint MSPs sug-gests this is likely to be caused by a slight under-estimate ofthe TOA uncertainties. We can fix this by re-scaling the lat-ter. In tempo , this is done by introducing a flag in the TOAlist called EFAC, followed by a numerical value. All TOAuncertainties are then multiplied by this factor. In what fol-lows, we use an EFAC of (cid:113) χ (cid:39) . . In this way thereduced χ becomes 1.0. This (slightly) increases the uncer-tainties of the timing parameters to more realistic values. Apart from the single-day gaps mentioned in Section 3.5,there are, in the case of 47 Tuc aa, no wider gaps with uniquesolutions, i.e., all remaining gaps are ambiguous . In suchcases it is not clear how to proceed.In what follows, we will address this problem. We startby presenting an automatic method for finding and listing The related
CIRCORBIT software can be found in
MNRAS , 1–13 (2015) etermining the rotation count of pulsars all acceptable solutions for each gap, this is called the gapmapper . By acceptable solutions we mean a solution havinga reduced χ < . , where that factor should not be too largeas to allow an exceedingly large number of solutions, but itshould not be so small as to exclude the correct solution,which might have a reduced χ larger than 1.0. As we willsee in Section 5.1, the eventual correct global solution has,using the four parameters being fit, a reduced χ = . (already using the EFAC factor of 1.185 mentioned above).Using a smaller χ threshold might have excluded the cor-rect solution. The limit of 2.0 is chosen because for mostMSPs the reduced χ in their published timing solutions issmaller than that. One tool we will need to do the connection is to learn howto find the best solutions for any gap automatically. An ef-ficient method to achieve this can be derived from Fig. 5 inSection 3.5: the reduced χ varies as a parabola. We can findthe minimum of a parabola if we sample it (i.e., use tempo to evaluate the reduced χ ) for three (not necessarily ac-ceptable) solutions. We generally choose these to be − b , 0and b , with b = being a good choice. In this case, then theminimum will be the integer M closest to m = b χ (− b ) − χ ( b ) χ ( b ) + χ (− b ) − χ ( ) . (1)Then, starting from M , we evaluate the integers immediatelyabove and below it to establish the quality of fits, keeping arecord of the acceptable solutions.Using this simple device, we have mapped all the gapsin our 47 Tuc aa data set. For each gap i we list in Table 2the best solution ( M i ) and the number of acceptable 1-gapsolutions ( n i ) around it (roughly one half of that number ofsolutions above and below M i ). One possible way of finding the timing solution for 47 Tuc aawould be to try all possible combinations of acceptable so-lutions for the 21 gaps. However, this would clearly be im-practical: the product of all n i for this data set is , or approximately . × . While mapping these gaps each tempo iteration took about 0.17 seconds; this means that,with the same computer, the whole process would take . × years. As we will see below, the vast majority of thesecalculations are unnecessary.We now describe the fundamental concept of thesolution-finding algorithm. Although in the previous stepwe are not able to find a single unambiguous gap, we canassume, for one of the gaps (say, number 2), that one of itstwo acceptable 1-gap solutions ( − Table 2.
Unconnected gaps in the 47 Tuc aa TOA data set andcorresponding best solution for each gap ( M i ) and number ofacceptable solutions ( n i ). i Gap M i n i −
14 582 51490 - 51492 0 23 51492 - 51582 −
53 864 51582 - 51589 −
75 51589 - 51690 +8 1156 51690 - 51809 +34 847 51809 - 52193 −
38 1928 52193 - 52262 −
729 52262 - 52307 −
27 4610 52307 - 52526 +33 16811 52526 - 52531 0 312 52531 - 52693 −
81 14613 52693 - 52925 +33 14014 52925 - 52930 −
515 52930 - 53152 −
70 17216 53152 - 53164 +5 1317 53164 - 53280 +22 4918 53280 - 53348 −
36 8019 53348 - 54314 −
128 69720 54314 - 54646 −
74 30421 54646 - 54816 −
37 168 k -gap solution stands out in terms of its quality.The reason why this works is because of what has beendescribed in Section 3: the assumption that there is phaseconnection for the first gap (in this case number 2) increasesthe precision of the ephemeris, which is derived from a largergroup of connected TOAs. This means that, when we at-tempt to connect a new gap (in this case 11), the numberof acceptable solutions decreases (in this case from three to two for each of the 2 assumptions).This means that there are acceptable 1-gap solutionsfor gaps 2 and 11 that do not provide an acceptable solu-tion when they are assumed for both gaps at the same time,i.e., they don’t work well together as a 2-gap solution. Sincethe reduced χ of this 2-gap solution is high, it is unlikelythat any k -gap (with k > ) solutions based on this 2-gapsolution will have a low reduced χ . We can therefore forgetabout this 2-gap solution when attempting to find the rota-tion count for the following gaps, i.e., the combination doesnot pass the χ sieve.As we will see in Section 5.1 for the specific case of47 Tuc aa, this sieve is extremely powerful because it dec-imates the number of possible solutions every time we testthe surviving k -gap solutions against gap k + . The effectmight appear to be small at every iteration, but it increasesexponentially with every stage. The number of acceptablesolutions increases fast at first, but as more gaps are tested,the number of acceptable solutions stabilizes and then startsdecreasing. Eventually one solution starts standing out interms of the quality of the fit. This solution is then gener-ally able to unambiguously connect all other gaps. MNRAS000
37 168 k -gap solution stands out in terms of its quality.The reason why this works is because of what has beendescribed in Section 3: the assumption that there is phaseconnection for the first gap (in this case number 2) increasesthe precision of the ephemeris, which is derived from a largergroup of connected TOAs. This means that, when we at-tempt to connect a new gap (in this case 11), the numberof acceptable solutions decreases (in this case from three to two for each of the 2 assumptions).This means that there are acceptable 1-gap solutionsfor gaps 2 and 11 that do not provide an acceptable solu-tion when they are assumed for both gaps at the same time,i.e., they don’t work well together as a 2-gap solution. Sincethe reduced χ of this 2-gap solution is high, it is unlikelythat any k -gap (with k > ) solutions based on this 2-gapsolution will have a low reduced χ . We can therefore forgetabout this 2-gap solution when attempting to find the rota-tion count for the following gaps, i.e., the combination doesnot pass the χ sieve.As we will see in Section 5.1 for the specific case of47 Tuc aa, this sieve is extremely powerful because it dec-imates the number of possible solutions every time we testthe surviving k -gap solutions against gap k + . The effectmight appear to be small at every iteration, but it increasesexponentially with every stage. The number of acceptablesolutions increases fast at first, but as more gaps are tested,the number of acceptable solutions stabilizes and then startsdecreasing. Eventually one solution starts standing out interms of the quality of the fit. This solution is then gener-ally able to unambiguously connect all other gaps. MNRAS000 , 1–13 (2015) Paulo C. C. Freire
The relative simplicity of the method outlined above al-lows a simple implementation as a UNIX shell script (calledsieve.sh); this is freely available in the Dracula Github repos-itory , complete with a short description of its usage. Theadvantage of such a script is that it can run on any UNIX orLINUX platform (which normally have sed and awk by de-fault), without the need for any special program other than tempo .However, there are still many improvements that couldbe made to speed up the process and make it more robust.One improvement that has already been partly imple-mented is that the script can be made to run much faster byrunning in shared memory, as in the case of our Github im-plementation. This avoids writing output files to hard disk,a process that causes disk wear and takes significantly moretime than the calculations themselves.Two further improvements have to do with limitationsof tempo . The first one is that the values of the reduced χ being reported have a limited range: if they are larger than999999.99 they are written by the program as a string ofasterisks. In such cases it is impossible to determine the po-sition of the minimum using eq. 1. Second, each time tempo is called it spends most of the time a) reading clock correc-tion files and correcting the TOAs, b) reading Earth rotationfiles and calculating the observatory position and c) readingSolar System ephemerides and calculating the Earth posi-tion. These operations only need to be done once. We notethat the PINT timing program does not have these limita-tions, and for that reason it will certainly allow much fasterdiscovery of rotation counts.A fourth improvement has to do with the script itself.In the current version the user decides which gap is to beconnected next. This is not necessarily the gap that yieldsthe smallest number of solutions. An obvious next step is tohave the script determine what is the best gap to connectnext. This minimization is likely to be very important, sinceany reduction in the number of solutions at each step hasexponential consequences when many steps are considered.A fifth improvement has to do with an analysis of thetiming solutions that are produced with every fit – the infor-mation they provide has not been used in the work above. Ifthe estimated positions are well outside the telescope beam,or if the pulsar has a negative (cid:219) P , then the solutions can besafely excluded. However, the latter condition must be usedwith care: in a globular cluster, many pulsars have negative (cid:219) P ; as we shall see this is the case for 47 Tuc aa.Another technique that should provide an extreme ac-celeration is to parallelize the computation of solutions forsuccessive gaps. This is conceptually very simple: while aversion of the script is still searching for solutions for gap A(using the gap mapper), it passes any acceptable 1-gap solu-tions it finds to a second version of the script working on gapB. This sorts (according to reduced χ ) the new 1-gap so-lution against all previous unprocessed 1-gap solutions thathave been passed to it from gap A and then works on thebest one (deleting it from the list of tasks to do). Assumingthat 1-gap solution, it will use the gap mapper to look for so-lutions for gap B. If it finds any acceptable 2-gap solutions, it https://github.com/pfreire163/Dracula Table 3.
Successive gaps (order given by k ) being tested for47 Tuc aa; see i numbers in Table 2. P - product of the n i forthe gaps being tested, which is the number of solutions we wouldhave to test if we were not sieving the acceptable solutions at eachstage. N k - number of k -gap solutions actually being tested ateach stage. Here we can see how this value becomes progressivelysmaller than P , demonstrating the power of the sieve technique.In the last column, we indicate the correct solution for these gaps(number 2 of the 30 solutions at step 10). These are small cor-rections to the rotation count predicted by our initial ephemeris. k i n i P N k Solution1 2 2 2 2 02 11 × × −
24 14 × −
15 16 ×
13 = 2730 26 66 17 ×
49 = 133770 36 317 18 ×
80 = 10701600 90 −
338 1 ×
58 = 620692800 135 09 3 ×
86 = 53379580800 84 − ×
115 = 6138651792000 32 25 will pass them to an analogous script using the gap mapperto find solutions for gap C. This script passes the acceptable3-gap solutions it finds to an analogous script working ongap D, etc. This saves a lot of time not only because of theparallelism - there is no need to wait for the conclusion ofthe discovery of all acceptable ( k − ) -gap solutions to startsearching for k -gap solutions - but also because the sortingby reduced χ at every stage focuses the processing on thebest possible solutions first at every stage. This should leadto a bare minimum of tempo iterations. Using the technique described in Sections 4.1 and 4.2,we were able to determine a phase-coherent ephemeris for47 Tuc aa with 1278 tempo iterations used in the gap map-per, plus an extra 476 tempo iterations that evaluated thequality of valid solutions around minima, plus a similar num-ber of iterations that evaluated non-valid combinations. Thefirst two sets of iterations took a total time of about 3 min-utes. In what follows, we give a detailed account of this pro-cess and describe the resulting phase-coherent ephemeris.
First, we ran the sieve.sh script looking for solutions for asingle gap (gap n.2 in Table 2), using the gap mapper. Thereare only two 1-gap solutions, ( − ) and (0). Then we ran thesame script again to find solutions for gap 11 based on thetwo 1-gaps solutions for gap 2. It found four acceptable 2-gap solutions: (0, 1), ( − , 0), ( − , − ) and (0, 0). Then, foreach of these the script searches for acceptable solutions forgap 4, finding 15 acceptable 3-gap solutions. The successivegaps being tested and the number N k of acceptable k -gapsolutions found are listed in Table 3.The number of k -gap solutions the script has to testat each stage becomes progressively insignificant compared MNRAS , 1–13 (2015) etermining the rotation count of pulsars to P , which is the number of all possible combinations ofintegers we would have to test if we were not sieving theacceptable solutions at each stage. This number, P , is theproduct of the number of 1-gap solutions ( n i ) for the gaps al-ready considered. This is a clear demonstration of the powerof sieving.After finding solutions for the tenth gap (number 5 inTable 3), the script starts running into a problem mentionedin Section 4.3: the values of reduced χ at one or more of thesample points used by the gap mapper ( − b , 0 and b ) startbecoming too large. This has a good implication: it meansthat the ephemerides resulting from the assumption of these10-gap solutions have a strong predictive power, otherwisethe values of reduced χ for “wrong” solutions would not beso large. This is the stage at which we should start looking inmore detail at the individual solutions. Fortunately at thisstage there are only 32 acceptable 10-gap solutions to lookat. We wrote a second script (test.sh, also found in theDracula Github repository) that, for any set of k -gap solu-tions, applies the PHASE statements contained in each ofthem and then uses tempo to derive an ephemeris that as-sumes that rotation count. Then the script switches to thelatter as a starting ephemeris, for this to work it removesall the PHASE statements from the previous TOA list sincethose are already taken into account by the new ephemeris. Ifall pre-fit residuals obtained with this ephemeris fall within arelatively narrow phase range, or if they have a clear, slowly-varying long-term pattern, then we know we have the correctsolution.Using this script to look at all 10-gap solutions for47 Tuc aa, we find that the ephemeris based on the sec-ond best (reduced χ = . , integers presented in the lastcolumn of Table 3) is able to predict all subsequent TOAsto less than 5% of a rotation (see Fig. 6). This behaviouris unique among the 32 acceptable 10-gap solutions. Thismeans that this is the correct solution. Adding these integersto the rotation count predicted by the original ephemeris weobtain the correct rotation count.We then use the ephemeris derived from the correct 10-gap solution as the new starting ephemeris. This implies wemust remove all previous PHASE statements. Furthermore,because we’re confident that this ephemeris predicts the cor-rect global rotation count, we can remove all 10 remaininggaps (i.e., we remove all JUMP statements from the remain-der of the TOA.tim file) - if we were not yet confident aboutthis, we could use the same sort of test as in Section 3.4 toassure ourselves that we indeed have unique “0” solutions forall remaining gaps. With all gaps removed, we use tempo to derive a global phase-coherent ephemeris, with much im-proved parameters compared to the initial ephemeris. Thisalmost eliminates the small residual trend we still see inFig. 6, with a post-fit reduced χ that is similar as in theprevious step (1.39), despite the fact that 10 degrees of free-dom (the JUMP statements) have been eliminated for thisfit. This is another strong indication that this is the correctsolution. Fitting for µ α , µ δ and the second spin frequencyderivative ( (cid:220) ν ), we obtain an even lower reduced χ , 1.33.The resulting ephemeris (with the EFAC updated to 1.367so as to yield a reduced χ of 1.0 and more conservative un-certainty estimates) is presented in Table 4. The parameteruncertainties were derived using a bootstrap method, similar Figure 6.
Pre-fit residuals using the ephemeris obtained withthe second-best 10-gap solution. The narrow range of phases forall the TOAs outside these 10 gaps indicates that this particularephemeris is based on the correct rotation count.
Table 4.
Timing solution for PSR J0024 − Observation and data reduction parametersFitting program . . . . . . . . . . . . . . . . . . . . . . TEMPOTime Units . . . . . . . . . . . . . . . . . . . . . . . . . . . TDBSolar system ephemeris . . . . . . . . . . . . . . . DE421Reference Epoch (MJD) . . . . . . . . . . . . . . 51600Span of Timing Data (MJD) . . . . . . . . . . 51143 - 54816Number of TOAs . . . . . . . . . . . . . . . . . . . . . 49EFAC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.367RMS Residual ( µ s ) . . . . . . . . . . . . . . . . . . . 26.4Reduced χ . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.00Timing parametersRight Ascension, α (J2000) . . . . . . . . . . . 00:24:07.2783(8)Declination, δ (J2000) . . . . . . . . . . . . . . . . − ν (Hz) . . . . . . . . . . . . . . . . 541.89365494899(14)First derivative of ν , (cid:219) ν ( − Hz s − ) . . . +1.34754(35)Second derivative of ν , (cid:220) ν ( − Hz s − ) . +7.0(2.9)Proper motion in α , µ α (mas yr − ) . . . 4.6(8)Proper motion in δ , µ δ (mas yr − ) . . . . − cm − ) . . . 24.941(7)Derived parametersSpin Period, P (ms) . . . . . . . . . . . . . . . . . . 1.8453805296800(6)Spin Period Derivative, (cid:219) P ( − s s − ) − α , θ α ( (cid:48) ) 0.123Angular offset from centre in δ , θ δ ( (cid:48) ) − θ ⊥ ( (cid:48) ) 0.465 θ ⊥ (cluster core radii) . . . . . . . . . . . . . . . . . 1.34Projected offset from centre, R ⊥ (pc) . 0.63 to what is used by Freire et al. (2017). The TOA residualscalculated with this are depicted in Fig. 7. The comparisonof some of the parameters with those of the early partialephemerides is presented in Table 1. The pulsar ephemeris includes a very precise position. Asshown in Fig. 4 of Freire et al. (2017), the pulsar is located, in
MNRAS000
MNRAS000 , 1–13 (2015) Paulo C. C. Freire
Figure 7.
TOA residuals for 47 Tuc aa. Notice the sparseness of the data set. All residuals appear within a narrow range of spin phases.This happens without any JUMP statements, indicating that the ephemeris is correct and can predict all times of arrival. No trends inthe residuals are visible, this is an indication that the ephemeris gives a good description of the TOAs measured for the pulsar. projection, almost due South from the centre of the cluster,at a distance of 0.465 arcminutes, or 1.34 core radii.Unlike virtually all pulsars in the Galactic disk, this pul-sar (and many others in globular clusters) has a negative (cid:219) P .For a rotationally powered pulsar, the intrinsic (cid:219) P is alwayspositive; this means that, for 47 Tuc aa, this quantity is dom-inated by the contribution from a negative acceleration ofthe pulsar in the gravitational field of 47 Tuc (see discussionin Freire et al. 2017). A negative acceleration means that theline-of-sight acceleration of the pulsar is in the direction ofthe Earth; this implies in turn that the pulsar is in the farhalf of the cluster. This is expected from the relatively highDM of 47 Tuc aa: previous observations have shown that thepulsars in the far side of the cluster (those with negative (cid:219) P )have larger DMs than the remaining pulsars. This happensbecause there is a cloud of ionized gas in the centre of thisglobular cluster (Freire et al. 2001b). The magnitude of thisnegative acceleration can be explained by the model of thegravitational field of the cluster presented by Freire et al.(2017), see Fig. 6 of that paper.We also have a marginal detection of (cid:220) ν ; this is causedby a change in the line-of-sight acceleration in the gravita-tional field of the cluster that results from the change of theposition of the pulsar in the cluster; these can also be causedby the motion relative to nearby stars. The derived jerk isdisplayed graphically with the jerks of other pulsars in Fig.7 of Freire et al. (2017).Since the TOA data set covers a total of 9 years, wecan also measure the proper motion; this is also listed inTable 4. This proper motion is roughly consistent with theproper motion of the remaining pulsars presented by Freireet al. (2017), but it has a lower precision given the smallnumber of TOAs.The precise location provided by the phase-coherentephemeris has allowed a detection of the pulsar in X-rays(Bhattacharya et al. 2017), as for the other 22 MSPs in47 Tuc with such ephemerides. Its X-ray flux is the smallestfor any MSP in this cluster. This somewhat surprising findis used to put an upper bound for amplitude of r-mode os- cillations in this pulsar as α < . × − and to constrainthe shape of the r-mode instability window. In the first part of this paper we have described a simpletechnique for determining the correct global rotation countof any pulsar. This technique allows for a simple determina-tion of phase-coherent ephemerides for the vast majority ofpulsars. However, for some pulsars the sparseness and/or lowprecision of the TOAs does not allow such a simple deriva-tion of the global rotation count.To solve such cases, we developed an algorithm thattests solutions for k gaps at the same time. Assuming eachof the acceptable solutions for k − gaps, the algorithmsearches (using the “gap mapper”) for acceptable solutionsfor gap number k ; if these exist then together with the as-sumed ( k − ) -gap solution they make new k -gap solution(s). k is increased until a single k -gap solution becomes clearlysuperior to the rest. The algorithm naturally eliminates so-lutions with a bad reduced χ at every step. This is whatwe call “sieving”.This should not be described as a “brute-force” tech-nique since estimating the suitability of all combinationsof acceptable 1-gap solutions would clearly be impractical.Sieving makes the problem not only tractable, but relativelycheap from a computational point of view. Implementingthe hierarchical parallel mode described in Section 4.3 andadding automatization should further reduce the amountof time needed to find the correct rotation count with thismethod; this will likely be necessary for binary pulsars withsparse data sets and/or poor timing precision. We have alsosuggested changes to tempo that might also produce furthergains in speed and reliability; these issues have already beenaddressed in the development of PINT , which will likely bea better vehicle for future development.We have implemented this technique with a simple shellscript, which we have made freely available. We have demon-
MNRAS , 1–13 (2015) etermining the rotation count of pulsars strated the technique (and its implementation) by find-ing the correct global rotation count and resultant phase-coherent ephemeris for 47 Tuc aa, an isolated MSP in theglobular cluster 47 Tuc with a very sparse set of detec-tions. The scientific implications of the parameters in thisephemeris have already been presented in previous studies.The automatic determination of the phase-coherentephemeris for a MSP with such a sparse set of detectionsimplies that this process might be achievable with signifi-cantly smaller amounts of telescope time than necessary un-til now. This will likely become a more pressing issue oncehigh-sensitivity instruments like the SKA come online: thisinstrument will be finding many thousands of faint pulsars(Keane 2017; Levin et al. 2017) that no other telescope willbe able to detect, which means that SKA time for follow-upwill be limited for each pulsar. ACKNOWLEDGEMENTS
We would like to thank our collaborators in the 47 Tuctiming project (Michael Kramer, Dick Manchester, AndrewLyne, Fernando Camilo, Dunc Lorimer, Christine Jordan,John Sarkissian and Nichi D’Amico) for collecting and help-ing to analyze the amazing AFB data set on 47 Tuc over theyears. We would also like to thank Norbert Wex and An-drew Cameron for suggestions and comments on the script,and Erik Madsen for his contributions to the Github im-plementation. Both authors gratefully acknowledge finan-cial support by the European Research Council for the ERCStarting grant BEACON under contract No. 279702, andcontinued support from the Max Planck Society. Finally,we acknowledge the referee, Dr. Scott Ransom, for valuablesuggestions that have significantly improved this work.
REFERENCES
Antoniadis J., et al., 2013, Science, 340, 448Berti E., et al., 2015, Classical and Quantum Gravity, 32, 243001Bhattacharya S., Heinke C. O., Chugunov A. I., Freire P. C. C.,Ridolfi A., Bogdanov S., 2017, MNRAS, 472, 3706Camilo F., Lorimer D. R., Freire P., Lyne A. G., ManchesterR. N., 2000, ApJ, 535, 975Damour T., Deruelle N., 1986, Ann. Inst. Henri Poincar´ePhys. Th´eor., Vol. 44, No. 3, p. 263 - 292, 44, 263Damour T., Taylor J. H., 1992, Phys. Rev. D, 45, 1840Desvignes G., et al., 2016, MNRAS, 458, 3341Edwards R. T., Hobbs G. B., Manchester R. N., 2006, MNRAS,372, 1549Folkner W. M., Williams J. G., Boggs D. H., 2008, Technical Re-port IOM 343R-08-003, JPL Planetary and Lunar EphemerisDE421. NASA Jet Propulsion LaboratoryFolkner W. M., Williams J. G., Boggs D. H., Park R. S., KuchynkaP., 2014, Interplanetary Network Progress Report, 196, 1Fonseca E., et al., 2016, ApJ, 832, 167Freire P. C., Kramer M., Lyne A. G., 2001a, MNRAS, 322, 885Freire P. C., Kramer M., Lyne A. G., Camilo F., ManchesterR. N., D’Amico N., 2001b, ApJ, 557, L105Freire P. C. C., et al., 2017, MNRAS, 471, 857Hobbs G. B., Edwards R. T., Manchester R. N., 2006, MNRAS,369, 655Keane E. F., 2017, preprint, ( arXiv:1711.01910 )Kopeikin S. M., 1996, ApJ, 467, L93Levin L., et al., 2017, preprint, ( arXiv:1712.01008 ) Lorimer D. R., 2008, Living Reviews in Relativity, 11Lorimer D. R., Kramer M., 2004, Handbook of Pulsar Astronomy¨Ozel F., Freire P., 2016, ARA&A, 54, 401Pan Z., Hobbs G., Li D., Ridolfi A., Wang P., Freire P., 2016,MNRAS, 459, L26Reardon D. J., et al., 2016, MNRAS, 455, 1751Ridolfi A., et al., 2016, MNRAS, 462, 2918Shaifullah G., et al., 2016, MNRAS, 462, 1029Tauris T. M., et al., 2017, preprint, ( arXiv:1706.09438 )The NANOGrav Collaboration et al., 2015, ApJ, 813, 65Verbiest J. P. W., et al., 2016, MNRAS, 458, 1267Wex N., 2014, preprint, ( arXiv:1402.5594 )This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000