Algorithmic Pulsar Timer for Binaries

Even though the initial timing ephemeris gives residuals that appear random, there are in fact patterns within each cluster (i.e., single observation or closely spaced sets of observations). PINT has the infrastructure to extract valuable information from these patterns by fitting each cluster for a different arbitrary time offset. More specifically, PINT will find the parameters, including time offsets, such that the reduced χ² (${\chi }_{\nu }^{2}$, where ν represents the degrees of freedom) of the model is minimized. Allowing PINT to fit a cluster for an arbitrary time offset is called applying a JUMP to that cluster. The arbitrary time offset itself may be called a JUMP as well. The reason this method obtains a better solution is that the JUMPs account for unknown phase counts between clusters, while the main timing parameters change to account for the phase behavior within each and every cluster. This is particularly useful for finding accurate binary parameters, as these can have large effects on short timescales. When JUMPs are applied to multiple clusters, information on the pulsar phase between the clusters is lost. Thus, the JUMPs have to be sequentially removed to obtain a full solution.

To initialize a model, each cluster except one should receive a separate JUMP. At this point, fitting the TOAs for several parameters, like spin frequency (f or F0) and a few binary parameters should give a very low ${\chi }_{\nu }^{2}$. This ${\chi }_{\nu }^{2}$ should be very close to unity, and a ${\chi }_{\nu }^{2}\gtrsim 2$ should be cause for alarm—the initial ephemeris is likely too erroneous (e.g., an inaccurate estimate of the initial binary parameters) or the TOA errors have been underestimated. The initial ${\chi }_{\nu }^{2}$ will be referred to as the base ${\chi }_{\nu }^{2}$, or ${\chi }_{\nu ,\ \mathrm{base}}^{2}$.

In this section, we will describe why it is important to rank clusters in terms of importance, and the way that we rank them. When initializing a model, not every cluster can be JUMPed—there needs to be an anchor cluster. This first non-JUMPed cluster will be referred to as the starting cluster. As described in Section 5, the starting cluster can have a large impact on whether APTB can discover a solution. It is therefore prudent to find the best starting cluster before actually attempting any solution, as running APTB for each and every possible starting cluster would be computationally expensive and likely unnecessary.

We thus adopt a scoring system discussed in Section 4 and Equation (2) of PR22, with only a minor change. This scheme scores each TOA based on how close it is to all other TOAs. A dense group of TOAs will generally all have high scores, while isolated TOAs will ideally be given low scores. Our change to the method described in PR22 is that we apply an index, α, to the weighting scheme such that the score of the ith TOA, n_i , reads

$$\begin{eqnarray}&&{n}_{i}=\displaystyle \sum _{j\ne i}\displaystyle \frac{1}{| {t}_{i}-{t}_{j}{| }^{\alpha }},\end{eqnarray} \tag{ 1 }$$

where t_j refers to the arrival time of the jth TOA and t_i refers to the arrival time of the TOA being scored. A cluster's score is the highest score of its constituent TOAs.

APTB currently uses α = 0.3, although α's best value has not been extensively tested. Moreover, a "best value" is likely highly dependent on the TOA data. APTB uses α = 0.3, and not α = 1 as in APT, because the scoring system of APT often overweights particularly dense clusters that are not close enough to other clusters to allow phase connection to proceed.

When observing a pulsar, it is necessary to already have good estimates of the relevant parameters. Thus, a single observation should correctly identify the relative pulse numbers of its TOAs. In other words, a cluster should be phase-connected internally. In the case that a cluster is not phase-connected internally, APTB can often correct for the TOAs that break phase connection by tracking the phase during the cluster. The arbitrary definition of zero phase can cause the model to confuse residuals that cross the ±0.5 rotation boundary. This makes these TOAs appear to lack phase connection when they otherwise would be phase connected. APTB corrects these instances by applying appropriate phase wraps on individual TOAs until phase connection within every cluster holds. Applying a phase wrap changes the assigned phase of a single or group of TOAs by an integer step. The correction is done quickly via the NumPy unwrap routine.

The initial timing ephemeris will usually not provide phase connection within every individual cluster nor between most of the clusters. Once every cluster has been corrected to be phase connected within itself, it is necessary to manually change the assigned pulse numbers for each cluster of TOAs so that at least the relative rotation count between two clusters can be compared. The goal of phase connection is to correct the relative rotation count between each and every cluster. The process of manually changing the pulse numbers of an entire cluster from those predicted by the (yet-to-be-correct) model is known as checking for phase wraps. See Table 1 for an example. Once phase wraps are used on individual clusters to phase connect the clusters within themselves, phase wraps are then only evaluated between JUMPed, or individually phase-connected, sections of the TOAs in order to map a gap (see Section 2.4).

Checking for phase wraps between clusters can be a computationally expensive process. It can be unclear how many possible phase wraps should be checked. Furthermore, if too many phase wraps are checked, the total number checked for every neighboring pair of TOAs grows quickly. See Section 3.3 for a more detailed explanation of what we call the phase wrap search space, but in short, if an algorithm is not careful, it may be in the process of checking more than 10²⁴ phase wraps, and even at a fraction of a second per phase wrap, this would take on the order of 10¹⁵ yr.

Table 1 shows an example where a cluster needs a single phase wrap removed to be correct. Unfortunately, we do not know this a priori. To attempt to determine if phase wraps are needed, APTB uses a technique referred to as mapping a gap (FR18). Below is a description of the gap-mapping technique used by APTB.

Let ${\chi }_{\nu }^{2}(n)$ be the ${\chi }_{\nu }^{2}$ after applying a phase wrap of n. If a phase wrap of m minimizes ${\chi }_{\nu }^{2}(n)$, then a Taylor series of ${\chi }_{\nu }^{2}(n)$ centered around m has no linear n term, and so the n² term is the leading nonconstant term. Therefore, the ${\chi }_{\nu }^{2}(n)$ has a quadratic dependence on the phase wrap when n is close to m (see also Figure 5 of FR18). Finding the vertex in phase wrap–${\chi }_{\nu }^{2}(n)$ space gives m, which represents the best phase wrap according to the best, but not necessarily fully correct, model. The parabola can be defined by sampling the ${\chi }_{\nu }^{2}(n)$ for three different phase wraps (i.e., by sampling three phase wrap-${\chi }_{\nu }^{2}(n)$ points) and the phase wrap of the vertex is given by the closest integer (i.e., round[m]) to

$$\begin{eqnarray}&&m=\displaystyle \frac{b}{2}\displaystyle \frac{{\chi }_{\nu }^{2}(-b)-{\chi }_{\nu }^{2}(b)}{{\chi }_{\nu }^{2}(b)+{\chi }_{\nu }^{2}(-b)-2{\chi }_{\nu }^{2}(0)}.\end{eqnarray} \tag{ 2 }$$

This is the formula for calculating the vertex x-coordinate of the parabola f(x) = a₂ x² + a₁ x + a₀ by sampling the values of f(b), f(−b), and f(0). APTB uses b = 5, as do FR18, though this can lead to an error if b = 5 gives a nonphysical parameter value, like a negative value of the projected semimajor axis ($a\sin i$). If this occurs, b = 4 is attempted, repeating lower b values until either the vertex is successfully calculated or b = 0. If the b = 0 iteration is reached, we cannot solve for the vertex of the parabola and the model is likely too poor to successfully phase connect, and so APTB will exit. A maximum b = 5 is chosen because b ≲ 3 would let small variations in the ${\chi }_{\nu }^{2}$, such as from TOA noise, greatly affect m in Equation (2). Additionally, b ≳ 7 would likely involve sampling phase wraps far from the phase wrap minimum—as the minimum is commonly at a phase wrap of order unity—weakening the quadratic-dependence assumption.

As for DRACULA, the gap-mapping method can only be successful if the TOAs are accurate to within their predicted error. If even a few TOAs are entirely erroneous, APTB will be fitting the parameters to compensate for fictitious differences in phase. Even worse, incorrect phase wraps will be treated as correct, and any model stemming from a model with even a single incorrect phase wrap can never be entirely phase connected when additional correct TOAs are included. Before phase connection is complete, it is nearly impossible to determine outliers, because the phase-connection process is sensitive enough to require the assumption of correct data. Therefore, it is prudent that APTB is given correct TOAs that are as accurate as possible, with special attention given to removing TOAs with high predicted error. APTB cannot currently address the cases where some TOAs are in error. A potential improvement to the algorithm would be to have branches (Section 3) investigate the removal of TOAs or entire clusters.

Going down all acceptable branches is like trying to produce a tree-like structure that is actively forming with every decision made. The tree's structure is only revealed by exploring it. APTB creates the tree structure by mapping each gap and testing phase wraps. Phase wraps that fail the pruning condition result in that branch being cut. Phase wraps that pass are ranked based on lowest to highest ${\chi }_{\nu }^{2}$. The best phase wrap (lowest ${\chi }_{\nu }^{2}$) in the short term is explored first, with its own branches being explored next. Thus, this is a depth-first search (see Figure 2).

Zoom In Zoom Out Reset image size

While exploring several branches is often needed to phase connect pulsars with limited data, it is important to know when a dead end is reached. In pulsar timing, an infinite number of options are available, so we have to be clever in deciding which branches, or phase wraps, are in fact dead ends. Thus, a branch must be pruned when it becomes clear that it cannot lead to the correct solution.

When it is time for a cluster to be unJUMPed, the gap-mapping technique is applied. Again, gap mapping allows APTB to know a reasonably good phase wrap to start with. The ${\chi }_{\nu }^{2}$ is calculated for the vertex point as well as for the phase wraps of ±1 from the vertex. The minimum ${\chi }_{\nu }^{2}$ of these three phase wraps, now referred to as PW_calc, must be lower than the pruning condition. The pruning condition, PC, used by APTB is ${PC}={\chi }_{\nu ,\ \mathrm{base}}^{2}+1$. If PW_calc ≥ PC, then the F-test for model comparison is done on the current model versus the current model plus a new parameter. The F-test provides a way to compare two models that have different degrees of freedom. We refer to Section 3 of PR22 for a brief discussion on the F-test for model comparison. This model is given one more chance and the gap-mapping technique is applied one more time, with PW_calc calculated again. If PW_calc ≥ PC still, this branch is pruned. This model is given a second chance because in some instances evaluating the F-test with different parameters should be done before mapping the gap.

When a branch is created, certain parameters are checked to see if they are physically relevant. For instance, if $a\sin i$ (the projected semimajor axis of the orbit) is negative, then the branch is pruned. Optionally, a branch may be pruned if $\dot{f}$ (F1) becomes positive, though this is not recommended, as F1 could become negative in a later model stemming from this branch. However, a pulsar can be observed to have a non-intrinsic positive F1 if it has a negative radial acceleration. This is common for pulsars in globular clusters, due to the gravitational potential of the clusters. In fact, two of the pulsars discussed in Section 6, PSRs J1748−2446aq and J1748−2446at, have a positive F1. Also optionally, if the R.A. or decl. leaves a set boundary, the branch is pruned. No other pruning mechanisms are implemented currently.

We pause the explanation of the algorithm to elaborate on why and when APTB may succeed or fail. The models that APTB explores based on its different phase wrap decisions all comprise the phase wrap search space (PWSS). The size of this space can be the dominant factor in determining whether APTB will succeed or not. In the case of Figure 2, the size of the PWSS is seven models. Generally speaking, the size of the PWSS, S, is

$$\begin{eqnarray}&&S=\displaystyle \sum _{d=0}^{{N}_{\mathrm{clusters}}}W(d),\end{eqnarray} \tag{ 3 }$$

where N_clusters is the number of clusters and W(d) is the number of models, or width, of the solution tree at a depth d. The depth refers to the number of clusters already phase connected for a given model, minus one. In Figure 2, W(0) = 1 and W(1) =W(2) = 3. If APTB takes τ seconds on average per model, then the total runtime will be τ S.

The nature of W(d) marks the boundary between a solvable and unsolvable pulsar using our algorithm. We will start with just one example where APTB would fail. If every gap has two plausible phase wraps (two branches), then W(d) = 2^d , implying S ≥ 2^d . For the case of PSR J1748−2446aq (Section 6), τ ∼ 10 s and N_clusters ∼ 80. Therefore, the total runtime would be on the order of 10 · 2⁸⁰ s ∼ 10¹⁷ yr. Clearly, APTB would be unsuccessful for PSR J1748−2446aq if it needed to check two branches for every gap.

APTB can be successful due to most pulsars following a different W(d). The initial model is precise enough to phase connect TOAs within an observation but not between observations (i.e., clusters). Therefore, W(1) can be quite high, and W(2) through roughly W(5) can become exponentially higher. Eventually, a model has phase connected enough clusters to be much more precise in its predictive power, and so it admits exponentially fewer possible branches, significantly lowering the width at higher depths. This all suggests we adopt

$$\begin{eqnarray}W(d)\sim \left\{\begin{array}{lr}{r}_{1}^{d}, & \mathrm{if}\ d\leqslant {d}_{0}\\ {r}_{1}^{{d}_{0}}{r}_{2}^{(d-{d}_{0})}, & \mathrm{if}\ d\gt {d}_{0}\end{array}\right\},\end{eqnarray} \tag{ 4 }$$

where r₁ and r₂ are the growth and decay rates, respectively, and dictate W(d) based on the cutoff depth, d₀. The values of r₁, r₂, and d₀ will determine if APTB can solve a pulsar. Non-integer values of W(d) can be interpreted as a characteristic value. Assuming W(d ≲ d₀) ≫ W(d ≫ d₀), we can simplify further to $S\sim {\sum }_{d=0}^{{d}_{0}}W(d)\sim {r}_{1}^{{d}_{0}}$. Therefore, if either r₁ or d₀ is too large, S will also be too large. Here, r₁ represents how quickly the model will branch until a depth of d₀ is reached. Pulsars with bad data (e.g., due to erroneous TOAs or low cluster density) will initially accept many phase wraps (increasing r₁) and could continue to accept them for longer (increasing d₀).

While both human and computer-based timing rely on the search space not becoming too large, a computer is able to probe a much larger search space. This is helpful in the cases where the cutoff depth allows for S no larger than on the order of hundreds of models to a few thousand models, but not when the search space balloons to tens of thousands of models or more. A human could take days to search S ∼ 100 models, but our algorithm can take minutes or hours and is fully automated. Furthermore, APTB would be more systematic than a human, searching for phase wraps in the same way every time. When the search space is in the tens of models, then both humans and APTB would find no issue, but APTB would still likely phase connect faster.

This is all to say that, while APTB can solve trickier pulsars than a human could solve, APTB can still fail even when functioning as intended. Thinking of our algorithm with the PWSS in mind can give a first-level analysis of the limits of our algorithm.

Most binary MSPs have nearly circular orbits, where is it difficult to measure the longitude of periastron (ω) and epoch of periastron passage (T₀). This is because the TOAs do not clearly indicate the location of periastron, and thus the epoch of periastron is not clearly indicated either (Lange et al. 2001). The ELL1 model was designed to address this difficulty by measuring the orbital phase with respect to the epoch of ascending node (T_asc) rather than T₀. Only intended for nearly circular orbits, the ELL1 model is a Newtonian binary model that neglects eccentricity (e) terms of order e² or higher, although most timing software have improved the accuracy of the ELL1 model by adding e² or even e³ terms. Therefore, this model describes residuals that follow a slightly perturbed sinusoid as a function of the orbital phase. Part of APTB's initialization is to not start fitting for ${\epsilon }_{1}=e\sin \omega$ and ${\epsilon }_{2}=e\cos \omega$ immediately. ₁ and ₂ can only be fit after the unJUMPed clusters span 5 times the binary orbital period (P_b ), though this can be changed by the user. In contrast, the model should immediately begin fitting for P_b , T_asc, and $a\sin i$, in order to correct for orbital effects within every cluster. Also, for this reason, these parameters should be known to a satisfactory precision beforehand. Notably, when the F-test processes are conducted (Section 3.2), both ₁ and ₂ are tested together. If they significantly help the model, they are added to the model together. In many cases, the eccentricity is negligible, so assuming ₁ = ₂ = 0 is often sufficient for phase connection.

It should be noted that, while using the ELL1 model, APTB does admit post-Keplerian parameters such as the time derivative of the binary period ($\dot{{P}_{b}}$). These are very rarely required for phase connection on short timescales, though, and they are usually included after phase connection has already been achieved. We briefly discuss the use of $\dot{{P}_{b}}$ in the timing of PSR J1748−2446ar in Section 6. Post-Keplerian parameters are only handled if the user explicitly includes them in the initial parameter file.

Though less common, some binary pulsars are elliptical enough to require binary models that make no assumptions about the eccentricity. One such model is the BT model, a Newtonian model. If the eccentricity is important for phase-connection purposes, this binary model should be used. APTB will begin fitting for all major binary parameters immediately, and therefore no testing via the F-test is done on these parameters. The base parameters of the BT model are P_b , $a\sin i$, T₀, e, and ω. Like in the ELL1 model, APTB can also handle post-Keplerian parameters, provided the user includes them in the initial parameter file.

Another noncircular orbit model is the DD model, a post-Newtonian approximation relativistic model that is theory-independent. This model is intended for pulsars with large eccentricities and measurable secular and periodic post-Newtonian effects. Similarly to the BT model, all major parameters will be fit for immediately. These parameters are P_b , $a\sin i$, T₀, e, ω, and the time derivative of the longitude of periastron ($\dot{\omega }$).

The model is initialized by applying JUMPs to every cluster except one. APTB then ensures that each cluster is phase-connected internally. The highest-scoring cluster, as defined by Equation (1), is chosen as the starting cluster and is the only cluster not JUMPed. If the ${\chi }_{\nu }^{2}$ is above three at this point, there may be a problem that APTB cannot fix. APTB will attempt to fix this by phase connecting and fitting two more times, but this is likely to fail. Either the starting timing model is too poor, the TOAs are inaccurate, or the TOA uncertainties are underestimated. In any case, the user must do something on their end before proceeding. This can include fixing one of the above problems or overriding this warning (see Section 5.2).

Provided the model has been properly initialized, the main algorithm can begin. The starting cluster is the first member of the unJUMPed clusters. The cluster closest in time to the unJUMPed clusters is unJUMPed and therefore added to this group. The number of pulsar rotations between the closest cluster and the other members is not known, so the phase wraps between the closest cluster and the previously phase-connected data are checked according to the gap-mapping methodology (Section 2.4). Gap mapping (ideally) finds the minimum ${\chi }_{\nu }^{2}$ for any particular phase wrap, but another phase wrap may indeed prove to be the correct one. For this reason, APTB incrementally checks the phase wraps 1, 2, 3,...,N_below below the best phase wrap until the model fails the pruning condition at a phase wrap of N_below + 1 below the best phase wrap. The same is repeated for 1, 2, 3,...,N_above above the best phase wrap, with the model failing the pruning condition at the N_above + 1 phase wrap above the best phase wrap. The phase wraps between N_below and N_above, inclusive, are acceptable models that represent branches on the solution tree.

All acceptable solutions are stored at this point, and attached to the parent model that created them. The next used model is determined by the branch search hierarchy (see Section 3.1 and Figure 2). This next model is then fit for new parameters, if they significantly improve the model. As in PR22, an F-test value of 0.005 or lower is deemed a significant improvement (i.e., the probability that the model fit improvement could be due to noise is <0.005). Not every nonfit parameter is tested. The R.A. and decl. can be tested after the unJUMPed clusters span at least T_R.A. and 1.3 · T_R.A., respectively, where

$$\begin{eqnarray}&&{T}_{{\rm{R}}.{\rm{A}}.}=\left(-\displaystyle \frac{30}{700}\cdot \displaystyle \frac{f}{\mathrm{Hz}}+40\right)\mathrm{days}.\end{eqnarray} \tag{ 5 }$$

The functional form of T_R.A. is such that T_R.A. depends linearly on the pulsar's spin frequency, T_R.A.(f = 700 Hz) = 10 days, and the maximum cluster span until the R.A. is tested is 40 days (i.e., ${\mathtt{\max }}[{T}_{{\rm{R}}.{\rm{A}}.}]$ = 40 days). The linear dependence condition is motivated by the fact that the sky position becomes important quickly for pulsars with higher TOA precision, which is usually directly proportional to the spin frequency. The three requirements for Equation (5) have not been derived from first principles, but rather from experience phase connecting many pulsars manually. F1 can be tested after the unJUMPed clusters span is long enough that a typical F1 would cause the residuals to exceed ± 0.35 in phase. A typical F1 is estimated based on the pulsar's spin frequency and rough placement in the P–$\dot{P}$ diagram (see, e.g., Figure 6.3 of Condon & Ransom 2016).

Currently, APTB has the capability to conduct an F-test and fit for $\ddot{f}$ (F2), but the user has to specify a minimum unJUMPed cluster time span, because the default is not to test for F2. We choose to omit F2 by default because F2 usually only becomes significant in time spans longer than are necessary for phase connection. Furthermore, typical F2 values vary widely, and estimating when F2 should be fit is difficult; a wrong guess may prevent the correct solution from being found, due to possible covariances with the other parameters. In most cases, F2 can be ignored until after a successful phase connection, and then manually fit, as we did for PSR J1748−2446aq (Section 6).

After the appropriate F-tests and new parameter additions, the model is then brought back to the beginning of the algorithm, as all model selection from the solution tree is done after checking for phase wraps. The first step is, again, to unJUMP the closest unJUMPed cluster, and so on.

Once every TOA has been unJUMPed, if ${\chi }_{\nu }^{2}\lt 10$, then the model is deemed a possible solution and saved. APTB then moves to the next starting cluster to see if another (or the same) solution arises.

In some cases, every branch is pruned before finding a solution. This means APTB has failed to find the solution with this particular starting cluster. The algorithm will attempt to run again with the second-best starting cluster and will keep selecting starting clusters until the top five (default) clusters by score have been attempted. While we hope to make APTB robust enough to find the solution regardless of a starting point, some assumptions break down when the model is too poor or the data are too sparse. When searching for correct phase wraps, APTB uses the quadratic dependence of the ${\chi }_{\nu }^{2}$ to efficiently find the minimum. If the minimum phase wrap is too large or the model is too poor, this quadratic dependence falls apart quickly. This could pin every early-on model into a local minimum rather than the global minimum. If this local minimum is too wide, as in the parabola covers a large number of phase wraps, or if there are too many local minima, then APTB may not find any solutions—or it may find a very large number of wrong solutions.

APTB will store all explored models on the user's local disk. These files are crucial in the case APTB is very close to finding a solution but failed due to unforeseen circumstances. For example, if the model requires F2 or the derivative of the binary orbital period, it may have pruned the correct branch. By finding this almost-complete model, the remaining few steps can be manually completed quickly. For more details on the files saved and how to use them, see the APTB User Guide.

APTB comes equipped with several ways to modify its behavior from the command line. We will describe the most useful modifications, with the rest being listed and briefly described in the APTB User Guide and the Python file on GitHub. By default, APTB will check for phase wraps and also use the solution tree structure. These can separately be turned off by including --no-check_phase_wraps and --no-branches, respectively, in the command line. However, it should be noted that branching only works if phase-wrap checking is on. If phase-wrap checking is turned off, then APTB will assume that the cluster that just had its JUMP removed has a phase wrap of zero with the unJUMPed clusters. In other words, at each step, APTB assumes the model correctly predicts the relative pulse numbers after removing a JUMP.

If phase-wrap checking is on but branching is not, then APTB will map the gap and find the phase wrap with the lowest ${\chi }_{\nu }^{2}$, but will not pursue nearby phase wraps that are also acceptable. For some pulsars, this is sufficient to time them successfully. While APTB is really meant to excel in the case where both phase-wrap checking and branching are required to time the pulsar, the algorithm can also easily and quickly time the less-challenging case.

When APTB finds a potential solution, there may be other unexplored branches in the solution tree that may yield other potential solutions. By default, APTB will explore these branches to rule out, or discover, any degenerate solutions. This can be turned off, but doing so can have similar effects to turning off branching, in that the correct solution may not be initially the best one. The advantage to stopping APTB from pursuing branches after finding a potential solution is that there may be a large number of potential solutions (tens or hundreds), which can consume much more computing time and storage than anticipated.

Some command line options affect the pruning condition. By default, if ${\chi }_{\nu ,\ \mathrm{base}}^{2}\gt 3$, then the algorithm will not proceed. This can be overridden if the user thinks APTB should try to time this pulsar anyway. The pruning condition by default is set to ${\chi }_{\nu ,\ \mathrm{base}}^{2}+1$, but users may specify their own values. When APTB can fit for R.A., decl., F1, and certain binary parameters like ₁ and ₂ can be adjusted manually. Finally, a rectangle in R.A.–decl. space can be specified. A model will be pruned if its astrometric position leaves this rectangle.