A 34-Year Timing Solution of the Redback Millisecond Pulsar Terzan 5A

Before timing Ter5A, we removed the effects of the eclipse from observations because the time-variable ionized gas can cause unpredictable time delays in the measured pulses. To remove most of these effects, we cut all TOAs corresponding to an orbital phase between 0.0 and 0.5, where the pulsar was most likely to be eclipsed. We further manually inspected the data to remove any remaining large-scale effects of the eclipse. A summary of the phase-cleaned observations is presented in Table 1.

Figure 1 shows an example of what eclipses look like in the timing data and what the timing residuals look like after cleaning TOAs. As the pulsar moves behind a clump of ionized gas, this gas increases the effective dispersion measure (DM), causing pulse delays. The stochastic variability of the wind from the companion makes timing delays unpredictable. While every orbit has some eclipse around phase 0.25, the width and effect of this eclipse varies significantly between observations. Several observations show additional micro-eclipses or hours-long total disappearances. This is one reason that long-term timing of this system was previously prohibitively difficult.

We searched each observation of the GBT and Parkes data for pulsations using the spider_twister⁴⁴4https://github.com/alex88ridolfi/SPIDER_TWISTER package, which searches over orbital phase by performing trial folds to determine a most probable T0. With an adequate local orbit, we then re-folded the observations with one minute integrations to determine TOAs, which could then be used to determine a much more precise T0 via timing software (e.g. Tempo or PINT). For the GBT data, the folding and TOA determination were performed with the prepfold and get_TOAs.py commands, respectively, from PRESTO⁵⁵5https://github.com/scottransom/presto (Ransom, 2001, 2011).

While this produced good starting parameters for individual days, we needed to use a model that accommodated time-varying T0s to determine a long-term timing solution. We independently used two different techniques to determine the same long-term timing solution.

Using PINT, we fit long-term timing parameters through a bootstrapping process, starting with initial parameters from published solutions (spin frequency, frequency derivative, position, DM, orbital period, and semimajor axis; e.g.  Nice et al., 1990, 2000; Urquhart et al., 2020), and phase-connecting the pulsar in 2$-$3 year segments.

We followed a four-step phase-connection process. First, we determined orbital phases (i.e. the best T0) for each observation (see section 3.1). Second, we phase-connected overlapping two-year chunks of data, fitting for frequency, first and second frequency derivatives, and position. Third, we verified that the number of pulse rotations between each TOA in the overlapping sections were identical for each phase-connected section. Fourth, we joined all the sections together using the rotational counts between all the TOAs.

This process produced a phase-connected timing solution. We next held the T0 values fixed while fitting for long-term parameters: orbital period, spin frequency, spin frequency derivatives one through five, DM, and DM derivative. PINT struggled with fitting the $10000+$ TOAs with the binary_piecewise model, given the large number of free parameters, which is one of the reasons we pursued an additional “isolation” fitting technique (see section 3.3). To accommodate non-GBT data, we manually fit for the JUMPs between datasets. We did this by fixing all long-term parameters and then isolating portions of the data that had two different instruments and fitting for JUMPs in these isolated portions in a pair-wise round-robin fashion until instrumental offsets were smaller than offsets due to T0 measurement errors and much smaller than overall timing noise.

In order to nail down the short-term orbital variations within the long-term solution, we assumed that the orbital period was constant over roughly two to four-week periods, and used multiple observations within such periods to constrain the piecewise constant T0 for those data. This was more important for the data before the year 2000, which had fewer and less precise TOAs.

The orbital phase residuals between the measured and predicted T0s, assuming a constant orbital period (see sections 4.3 and 4.4), showed a stochastic-like pattern with a general quadratic trend. We used a least-squares fit to determine a best-fit orbital period derivative from the T0s, and also performed a Gaussian Process Regression which allowed us to predict and refine T0 measurements that were poorly determined in the first iteration.

The result is a phase-connected timing solution with relatively well-behaved TOAs in individual observations and relatively smooth long-term behavior, as well as a time series of orbital phase wander as determined by the measured values of T0 (Figures 3 and 5).

Besides this piecewise technique, and as part of a related project on long-term timing of other globular cluster redbacks (Corcoran et al, in prep), we also phase-connected the system without breaking it into chunks by removing the effects of the binary and effectively “isolating” the pulsar. We used PINT and the accurately known T0 values measured over short timescales (see Table LABEL:tab:T0s_all), and the long-term $P_{\mathrm{orb}}$, and $\dot{P}_{\mathrm{orb}}$ values to subtract off the Roemer delays from the binary orbit from the TOAs for each observation, resulting in a much simpler dataset: an “isolated” pulsar with far fewer parameters to fit.

This isolation technique effectively de-couples the orbital variations from the long-term spin behavior of the pulsar. Since the arrival times are topocentric (i.e. in the reference frame of the observatories), but the binary Roemer delays need to be computed in the barycentric frame, we had to perform a first-order correction to the predicted Roemer delays to compensate for the difference in those frames. We tested this technique with simulated binary TOAs and were able to to remove all orbital effects to better than 100 ns. It is likely to be useful in other long-term redback timing efforts, and can be implemented in PINT with just a few lines of code.

Ultimately we used the short-term timing position measured with the GBT Spigot data (see Figure 8 and Table 3) rather than a long-term timing-fitted position in our final model, and Gaia-measured proper motion values (Vasiliev & Baumgardt, 2021) for the cluster rather than the timing-fitted proper motion of Ter5A (see section 4.5). The long-term timing parameters are presented in Table 5.

If we use the obtained pulse numbers but do not accommodate for changing T0, and fit for spin frequency, and the first and second spin frequency derivatives, we do obtain a phase-connected timing solution with a smooth overall trend, but within any individual observation the pulse phase varies by 5$-$10%. This solution is shown in Figure 2 and in Table 5.

We fit the data through five spin frequency derivatives to show the extent of the un-modeled spin-noise present in the data. Figure 3 shows the evolution of fit with successive derivatives. All plots shown have been fit for DM, DM derivative, frequency, and the number of frequency derivatives indicated.

Ter5A is a complex system and over a 34-year timespan the DM changes significantly. As a result, the DM is not well fit for much of the dataset. However, compared to the timing noise systematics, the variations due to the DM are small (i.e. few percent of a rotation for the F1-based solution and up to $\sim$10% of a rotation for the F5-based solution) and not prohibitive to finding a timing solution. In principle, we could fit for the DM on each observation, but simultaneously fitting for $\sim$ 300 additional free parameters becomes computationally difficult, and in addition, on many days the DM changes significantly on minute to hour timescales. We acknowledge this issue and leave it for future work.

There is an apparent abrupt timing feature around the year 2000 (i.e. MJD 51700). This is noted in Nice et al. (2000) and corroborated by our observations. While it is unfortunate that this feature occurs during the least-well sampled period of our data, the unambiguous phase connection of the Parkes and Jodrell Bank observations during that time period suggests this feature is real and corresponds to some physical change in the system. We note that whatever this physical change is, it appears strongly only once the second frequency derivative has been subtracted. We also note that several other timing features of smaller amplitude appear later in the data (e.g. around MJDs 55000, 56000, and 58000).

We measure a spin frequency derivative (i.e. F1) of $\sim$1.35$\times$10^-16 Hz s^-1, depending on the exact spin model used. This positive value is consistent with other measured values (e.g.  Lyne et al., 1990; Nice & Thorsett, 1992; Prager et al., 2017). The intrinsic spin-down rate of the pulsar must be negative. The positive measured first spin derivative value is likely due to the globular cluster accelerating the pulsar toward the observer from the far side of the cluster. There are additional mechanisms that contribute to apparent positive F1 values, and they usually affect the orbit and spin equally. We discuss this in the next two sections.

From Figure 3, it is clear that spin variations are not correlated with orbital variations. X-ray activity also does not appear to correlate with either of these variations, nor do the days in which the pulsar was undetected in radio observations.

The total fractional sensitivity of the pulsar’s spin (median spin phase error divided by the total pulse phase covered by the dataset, T$\times$F0) is $4.4\times 10^{-15}$ while the total fractional sensitivity of the pulsar’s orbit (median T0 phase error divided by the total orbital phase covered by the dataset, T$\times P_{\mathrm{orb}}^{-1}$) is $9.3\times 10^{-11}$, so the spin residuals are substantially more sensitive to torques and other effects on the system than the T0 residuals.

We examine the spectrum of the spin residuals by taking a Fourier decomposition of the F0$-$F1 solution (Figure 3, top panel) using a single “median residual” TOA from each observation. We calculated the powers of the components using the get_wavex_amps() and get_wavex_freqs() functions of pint.utils and determined a best-fit line using scipy.curve_fit(). This gives a steep power-law-like spectrum with a spectral index of $-5.54(7)$, as shown in Figure 4. The error represents the span between the maximum and minimum slope values estimated by multi-variate gaussian sampling using the fit covariance matrix. The first two harmonics are the most significant and dominate the measured spectral index.

In the absence of external effects, such as that of accelerations from the globular cluster, tidal effects in the bloated companion star, gravitational wave radiation from the compact orbit, and mass transfer or loss from the companion, we can expect T0 will evolve according to

where T0₀ is some reference measured T0, $n$ is an integer number of orbits that have elapsed between T0₀ and T0_n, and $P_{\rm orb}$ is the orbital period of the binary (1.82 hours). When we compare our measured T0 values with those predicted by this equation using a precise estimate for $P_{\mathrm{orb}}$ (see Table 5), we find that the measured values vary across $\sim$10 seconds, as shown in Figure 5. There is an apparent quadratic trend corresponding to an orbital period derivative ($\dot{P}_{\mathrm{orb}}$). Adjusting for the effects of $\dot{P}_{\mathrm{orb}}$, as described in the next section, results in a random process varying across $\sim$5 seconds, also shown in Figure 5.

We also performed a Gaussian Process Regression on the measured T0 values (see Figure 5) to estimate how smooth the variations might be in time and to better estimate T0 values on days where it could not be measured well. On several days, micro-eclipses or significant ionized gas effects during the observation cause strong systematic offsets in the T0 measurements due to unmodeled DM changes. Our final values are presented in Table LABEL:tab:T0s_all.

We determined the Power Spectral Density (PSD) (see Figure 6) of these $\dot{P}_{\mathrm{orb}}$ subtracted orbital variations by using scipy.curve_fit() to take a best fit line of the portion of the PSD between the first bin ($8.0\times 10^{-5}$ d^-1) and where the spectrum turns into white noise ($10^{-3}$ d^-1). The first bin is partially covariant with $\dot{P}_{\mathrm{orb}}$ and therefore has a reduced power. This produces a power-law with a spectral index $\gamma$ of $-0.9(2)$.

Three other works have assessed the spectral index of T0 variations in redbacks. From an analysis of Fermi-LAT observations of three redback pulsars, Thongmeearkom et al. (2024) found orbital phase wander spectral indices constrained to $\gamma<-2.4$, $\gamma=-3.81^{+0.32}_{-0.48}$, and $\gamma<-5.4$. Corcoran et al. (in prep) surveys five redbacks and measures spectral indices of $\gamma=-0.7(2)$, $-0.7(2)$, $-0.9(1)$, $-1.1(2)$, and $-1.7(2)$. Clark et al. (2021)’s analysis concluded there were too few significant frequency bins for a meaningful measurement. Our measurement of $\gamma=-0.9(2)$ is shallower than the indices observed in Thongmeearkom et al. (2024) and consistent with those observed in Corcoran et al. (in prep).

We compare our results to the Applegate model (Applegate, 1992), which predicts that spider pulsar systems should show periodic modulations in orbital period due to cycles of magnetic activity and changing quadrupole moment of the companion star. We see no evidence for such periodicity on the scale of up to $\sim$20 years.

The orbital period change due to gravitational wave emission (Peters & Mathews, 1963) can be expressed as a function of the component masses $m_{1}$ and $m_{2}$ and the well-known orbital parameters (e.g., Taylor & Weisberg 1982; Freire & Wex 2024):

where $e$ is the eccentricity, $P_{\mathrm{orb}}$ is the orbital period, $M$ is the total mass of the system, and $T_{\odot}\equiv GM_{\odot}/c^{3}=4.925490947\,\mu$s. For Ter5A, the eccentricity is zero, simplifying Eqn 2 to

The observed pulsar spin-down to spin period ratio $\left(\dot{P}/P\right)_{\mathrm{obs}}$ can be converted to an acceleration:

$A_{\mathrm{obs}}$ is the sum of several accelerations:

where $A_{\mathrm{int}}=c\left(\dot{P}/P\right)_{\mathrm{int}}$ is the intrinsic spin-down rate of the pulsar, $A_{\mathrm{GC}}$ is the contribution due to acceleration of the pulsar by the globular cluster, and $A_{\mathrm{gal}}$ is the contribution due to the acceleration of the pulsar by the Galactic gravitational potential. $A_{\mathrm{Shk}}$ is the apparent acceleration due to the Shklovskii effect which is given by

where $\mu$ is the proper motion of the pulsar and $d$ is the distance to the pulsar. Of these, $A_{\mathrm{GC}}$ is the dominant term since the pulsar is observed to have a positive frequency derivative.

Orbital period accelerations follow a similar form to Eqn. 5:

$A_{\mathrm{orb,meas}}=c\times\left(\dot{P}_{\mathrm{orb}}/P_{\mathrm{orb}}\right)_{\mathrm{meas}}$ is the observed orbital acceleration. $A_{\mathrm{GR}}$ is an acceleration-like term representing the quantity $c\times\left(\dot{P}_{\mathrm{orb}}/P_{\mathrm{orb}}\right)_{\mathrm{GR}}$, the GR-predicted contraction due to gravitational wave emission. $A_{\mathrm{RB}}$ is the apparent acceleration due to variations in the companion in the redback system. $A_{\mathrm{RB}}$ is not directly measurable, but can be constrained by the span of reasonable $A_{\mathrm{GR}}$ values, in combination with the estimates of the other acceleration parameters.

By combining Eqns 5 and 7, we can eliminate the dependency on $A_{\mathrm{gal}}$, $A_{\mathrm{GC}}$, and $A_{\mathrm{Shk}}$ (for a breakdown of the predicted contributions of these terms, see Table 2):

The quantity $A_{\mathrm{obs}}=c\left(\dot{P}/P\right)_{\mathrm{obs}}$ has a 7% fractional error ($-4.62(34)\times 10^{-10}$ m s^-2) purely due to measurement uncertainty. From fitting linear (i.e. $P_{\mathrm{orb}}$) and quadratic (i.e. $\dot{P}_{\mathrm{orb}}$) terms to the raw T0 residuals (Figure 5, top) we measure $\dot{P}_{\mathrm{orb}}=-2.5(3)\times 10^{-13}$.

Assuming a reasonable age for the pulsar of $5.0\times 10^{9}$ yr, within a factor of ten, and a surface magnetic field strength $B$ of $8\times 10^{8}$ G based on the mean of the 16 closest-period pulsars in the ATNF catalogue⁶⁶6https://www.atnf.csiro.au/research/pulsar/psrcat/, excluding globular cluster pulsars, we can predict $\left(\dot{P}/P\right)_{\mathrm{int}}$. Varying $B$ across a reasonable range of $6\times 10^{8}$ to $10^{9}$ G, we find $\left(\dot{P}/P\right)_{\mathrm{int}}=4.7^{+2.6}_{-2.0}\times 10^{-18}$, resulting in an intrinsic acceleration of $A_{\mathrm{int}}=14.0^{+7.8}_{-6.0}\times 10^{-10}$ m s^-2.

A 1.5 M_⊙ pulsar with a 0.09 M_⊙ companion star are reasonable estimates for the masses of the system components. MSPs are typically more massive than canonical pulsars, and the system is highly likely to have a large inclination due to the strong eclipses, meaning 0.09 M_⊙ is a likely companion mass. This yields the prediction shown in Table 2 for a predicted $\dot{P}_{\mathrm{orb}}$ of $-2.34^{+0.22}_{-0.17}\times 10^{-13}$, which is consistent with our measured value of $-2.5(3)\times 10^{-13}$. The error range on the predicted value comes from the range of likely surface magnetic field strengths and galactic accelerations.

Strader et al. (2019) measured the masses of 10 redbacks and redback candidates, finding an average mass for redback neutron stars of $1.78\pm 0.09$ M_⊙ with a dispersion of $\sigma=0.21\pm 0.09$ M_⊙. We therefore examine the viability of a 1.8 M_⊙ pulsar. A 1.8 M_⊙ pulsar with a 0.09 M_⊙ companion predicts a $\dot{P}_{\mathrm{orb}}$ of $-2.59^{+0.17}_{-0.24}\times 10^{-13}$, consistent with our measured value of $-2.5(3)\times 10^{-13}$. For the lowest reasonable values of magnetic field and galactic acceleration, a 1.9$-$2.0 M_⊙ pulsar is still consistent with our measured value of $\dot{P}_{\mathrm{orb}}$. Although this higher mass pulsar cannot be definitively ruled out, it seems unlikely that it could have reached this mass while only being spun up to an 11.56 ms pulse period.

We therefore cannot make strong constraints on Ter5A’s mass; it is most likely in the range of 1.5$-$1.8 M_⊙. As a less-recycled redback, it is likely also a less massive one. If we do indeed have a minimum-likely-mass system, the contribution of $\dot{P}_{\mathrm{RB}}$ is no greater in magnitude than $-6\times 10^{-14}$, but we do not need any $\dot{P}_{\rm RB}$ to explain our results. We conclude that $\dot{P}_{\mathrm{orb}}$ is dominated by gravitational wave radiation.

We note that the time for a binary to coalesce due to GR is given by (e.g., Shapiro & Teukolsky, 1983):

where $r$ is the total separation between the two objects. For Ter5A, if the companion was a compact star that would lose no mass, the coalescence time would be $\sim$250 Myr. In reality, though, the companion is not compact and more complex binary interactions would occur as the orbit shrinks, dramatically changing this GR-driven evolution.

We also note that we do not expect a measurable value of the contraction of the semimajor axis in timing residuals due to the contraction of the orbit from gravitational wave emission. Over the time period of our observations, $\dot{P}_{\mathrm{orb}}$ predicts a change in period that, via Kepler’s third law, corresponds to a fractional change of $2.3\times 10^{-8}$. The fractional error of our timing semi-major axis measurement is $6.4\times 10^{-6}$, so the change in semimajor axis is significantly smaller than our errors and is currently undetectable via timing.

In theory, the relativistic precession of the system (i.e., Susobhanan et al., 2018), predicted to be $\sim$20 ^∘/yr, should be easily measurable. However, in practice, the eccentricity is so close to zero that this is not possible — the precession is entirely covariant with the orbital period.

Fitting for position over the entire dataset yields a value that is entirely inconsistent with VLA imaging positions and prior timing-based positions. The source of this discrepancy is two sided. Eclipse variability during long integrations can cause systematic positional offsets in VLA imaging (for further discussion of positional uncertainties, see Fruchter & Goss, 2000). On the other hand, long-term pulsar timing systematics also produce massive uncertainties and can strongly affect the position signal in the timing data. However, fitting for position over just 2-3 years of data at a time produces reasonable positions that are both self-consistent and consistent with the proper motion measured by Gaia (Vasiliev & Baumgardt, 2021). While we cannot get extremely precise positions with this method, short-term timing is good enough to allow multi-wavelength follow-ups, with the caveat that Ter5A’s very low ecliptic latitude ($\beta=-1.4^{\circ}$) means declination measurements via timing are highly imprecise.

Fitting for proper motion over the entire dataset likewise produces nonsensical results; again, long-term timing systematics skew the proper motion measurements. Timing over the entire dataset measures a proper motion in the RA direction of $-1.53(11)$ mas/yr and in the Dec direction of $-24.5(3.2)$ mas/yr. The dispersion in the proper motion of Terzan 5 stars is $\sim$0.5 mas/yr (Taylor et al., 2022; Gaia Collaboration et al., 2021), so while the RA-proper motion is plausible in comparison with Gaia-measured values ($-1.53(11)$ mas/yr compared to $-1.989(68)$ mas/yr) (Vasiliev & Baumgardt, 2021), the declination-proper motion radically differs ($-24.5(3.2)$ mas/yr compared to $-5.243(66$) mas/yr). The source of this discrepancy is likely a combination of the effects of low ecliptic latitude on declination precision and covariance with other measured parameters.

The large timing declination proper motion is also unphysical. If Ter5A is indeed a cluster member, its elliptical orbit about the cluster core would never take it on the measured trajectory. It cannot even be sensibly ejected in this direction. While in principle these results could suggest that Ter5A is not actually a cluster member, this is highly unlikely. The observed positive spin frequency derivative practically guarantees cluster membership.

However, if we instead fit for the position of the pulsar over just a few years of well-sampled, smoothly behaved data at a time, we get results in the RA-direction that are consistent with Gaia’s proper motion measurements, as shown in Figure 8. However, the declination results are too ridden with systematic errors to draw any conclusions.

We compare our timing positions to VLA positions in Table 3 and Figure 8. The details of the first two VLA positions are described in Fruchter & Goss (2000) and Urquhart et al. (2020). The last position, from Urquhart et al. (in prep), is obtained by propagating pulsar timing positions and Gaia proper motions from the 55000 MJD epoch to 59669 MJD. Continuum images were shifted by matching to 20 timed pulsar positions, though Ter5A was excluded from the frame shift calculation. The new positional uncertainties come from the quadrature sum of the rms associated with the PSR frame shift and the thermal fluctuations in the continuum position fitting. With Ter5A being such a bright source, we are dominated by the PSR frame shift uncertainties.

As a result, we use the SPIGOT timing position and the proper motion measured by Gaia (Vasiliev & Baumgardt, 2021) for our long-term timing solution. These are the parameters that produced the residuals in Figure 3.

Ter5A has an associated variable X-ray source, and it is one of the hardest-spectrum and X-ray brightest MSPs in Terzan 5 (Bogdanov et al., 2021; Bahramian et al., 2020). This source might indicate intermittent accretion. Notably, Ter5A has been observed to have high- and low- X-ray emission states, with approximately 40% of its X-ray counts measured during two comparatively short observations where it was $\sim$9 times brighter than during other observations (Bogdanov et al., 2021). In the “high”-flux state, Bogdanov et al. (2021) reports an unabsorbed 0.5-8 keV X-ray luminosity of $L_{X}=8.6(2.1)\times 10^{31}$ ergs/s, assuming a distance of 5.9 kpc to the cluster. In the “low”-flux, or quiescent state, this luminosity is $L_{X}<1.3\times 10^{31}$ ergs/s with a photon index of $1.5\pm 2.0$ (Bogdanov et al., 2021), which places it in the redback-black widow overlap region in $L_{X}$–photon index space (see Figure 10 of Urquhart et al. (2020)). The quiescent $L_{X}$ is lower than that of a typical redback, probably due to a lower-than-typical spindown luminosity.

Bogdanov et al. (2021) suggest two possible interpretations of the data: 1) these luminosity spikes are short-lived burst or flare activity indicating two transitions from a radio pulsar state to an accretion disk-dominated state, or 2) flaring of the companion star. The spectral and variability properties suggest intrabinary shock radiation (Bogdanov et al., 2021). The observed flaring is inconsistent with observed transitional MSPs, but could be consistent with magnetic activity from the companion’s surface.

We do not have contemporaneous radio observations during the flared activity. However, a comparison of Ter5A’s X-ray quiet days (in which emission is at Chandra’s detection threshold and not much higher) and its flared days with spin frequency residuals and orbital variations is presented in the bottom-most panel of Figure 3. There is no apparent correlation between X-ray activity and changes in either the spin or orbital behavior of the pulsar.

We detect Ter5A in the radio almost every time we observe it. Though there are gaps in our observations, we see very little evidence that Ter5A is radio-silent for any extended period of time and therefore very little evidence that Ter5A spends very much time actively accreting. Even if it were radio-silent for extended periods of time, that is not in and of itself a smoking gun for transitional behavior. For example, J1653$-$0158 is a 1.97 ms binary pulsar in a 75-minute orbit that has only been detected in gamma-rays despite 26 radio searches on 8 different observatories over 11 years across all orbital phases and several wavelengths (Nieder et al., 2020). However, J1653$-$0158’s radio non-detections are attributed to it being eclipsed by ablated material the vast majority of the time. Radio disappearances could simply be due to excess gas in the system and not prolonged accretion-driven X-ray emission.

Ter5A’s behavior is also in stark contrast with the known redback and transitional MSP M28I. M28I also experiences X-ray flux enhancements of 1-2 orders of magnitude when it switches from pulsar- to disk- state (Vurgun et al., 2022). However, in Chandra data, M28I was detected as a pulsar three times from July to September of 2002, as a disk-accreting X-Ray Binary (XRB) on two observations four days apart in August of 2008, and as a pulsar again on three observations from May to November of 2015 (Vurgun et al., 2022). While data are sparse, odds are low that M28I was by chance observed on two consecutive observations as an XRB if its presence in that state is rare. Furthermore, Papitto et al. (2013) reports that M28I’s behavior is consistent with months-long X-ray outbursts. We see no evidence for months-long state changes in Ter5A. M28I is likely stable for much longer periods of time than Ter5A is. We conclude that Ter5A is likely not a transitional MSP. If it is switching between accreting and non-accreting states, these switches are rare and short in duration.

The Applegate model predicts orbital variability due to deformations in the companion star from magnetic activity. It could potentially explain the observed orbital variations if those variations have a clear periodicity.

We investigate if the Applegate model is energetically viable. We assume that the companion star is approximately the radius of its Roche lobe, using Nice et al. (1990)’s value of 0.2 $R_{\odot}$. Consulting a $P$-$\dot{P}$ diagram, we note the intrinsic spindown of this pulsar is likely of order $\dot{P}\sim 10^{-19}-10^{-20}$. We assume $\dot{P}=10^{-19}$ for the following calculations. The energy radiated by the pulsar, or its spindown luminosity, is given by $L_{\mathrm{spin}}=4\pi^{2}I\dot{P}P^{-3}$, where $I$ is taken to be $10^{45}$ g cm². Assuming isotropic emission, the companion captures an energy of $L_{\mathrm{capt}}=4.3\times 10^{31}$ ergs/s of this spindown energy, or 0.01 L_⊙. This energy blows mass off the companion at a maximum rate given by

which predicts a value of $1.9\times 10^{-12}$ M_⊙/yr. This agrees with predictions made by Shaham (1995) and results in an evaporation time much longer than a Hubble time.

The energy available due to tidal heating, denoted $L_{T}$ for tidal luminosity, is given by

according to Equation 28 of Applegate & Shaham (1994), where $t_{M}$ is the evaporation time in years of the companion and $a$ is the separation. This yields $L_{T}=1.7\times 10^{29}$ ergs/s, or $4.2\times 10^{-5}$ L_⊙.

The Applegate model predicts a change in transported angular momentum $\Delta{J}$ given by Equation 27 in Applegate (1992):

where $R_{c}$ is the radius of the companion, here taken to be the Roche lobe radius. The associated energy is given by

where $\Omega_{\mathrm{diff}}$ is the differential rotation between layers of the star. This is not measurable, but since existing constraints are of order $\sim$0.01 of the total orbital angular velocity of the system, we will take the generous assumption that it is 0.1$\Omega$ to test the model’s energetic viability. Noting that $\Delta{P}_{\mathrm{orb}}/\Delta{t}\sim\dot{P}_{\mathrm{orb}}$ for $\Delta t$ the length of the dataset, Eqn 13 can be rewritten to express the luminosity predicted by the Applegate model $\Delta L_{A}$:

This yields $\Delta L_{A}=3.5\times 10^{29}$ ergs/s, or $8.9\times 10^{-5}\,\mathrm{L}_{\odot}$. We conclude that the Applegate model is energetically viable - there is more than sufficient spindown energy ($0.01\,\mathrm{L}_{\odot})$, and an available tidal heating energy reservoir of the same order of magnitude ($4.2\times 10^{-5}\,\mathrm{L}_{\odot})$, to drive these changes.

While energetically viable, we observe no periodicity in the T0 measurements in Figure 5 to suggest that this process is taking place on a timescale even of a few decades. However, it is certainly possible that some other form of energy exchange is occurring between the orbit and the companion.

If any of the ablated material falls back on the pulsar, it should torque the pulsar and affect its spin as the second derivative of phase $\phi$:

This quantity is related to (but not equivalent to) the second derivative of the timing residuals. From Figure 9, it’s apparent that the slope of the first derivative is high around the odd feature at about MJD 51500, so that feature does indeed correspond to some type of torque. It is unclear whether this torque is intrinsic (such as internal seismic activity) or extrinsic (such as accretion or the effects of a nearby star), but similar sharp features also appear later in the data.

These sharp turnarounds are inconsistent with glitches, as the pulses after the turnarounds are showing up later, not earlier, as would be expected for the increase in spin frequency following a glitch. These events all show the pulsar slowing down in its rotation. One possible explanation is that gas interacting with the magnetosphere is hindering Ter5A’s rotation, but the mechanism driving these spin features is very unclear.

There remain contaminants in the data. In order to limit the number of TOAs and focus on phase connection, much of the radio frequency information available in the original data has been discarded. For an example and analysis of how DM can vary over the course of a single Ter5A observation, see Figure 10 of Bilous et al. (2019). This figure is an analysis of the first observation in the GUPPI panel in Figure 7, and shows a DM change of $\sim$0.6 pc cm^-3 in a systematic manner over the course of the observation. This is an example of a day in which DM variations – which we do not account for in our study – prohibit precise measurement of T0. Similar DM-related effects are likely present during most observations, but with smaller amplitudes.

Of the nine redback pulsars with long-term timing solutions published or in preparation (Thongmeearkom et al. (2024), Corcoran et al, in prep), Ter5A has by far the most timing noise. Six redbacks show effectively flat timing residuals when fit out to the second frequency derivative, and two show flat residuals when fit to the fourth frequency derivative (Thongmeearkom et al. (2024), Corcoran et al, in prep), while Ter5A shows systematics at the $\sim$1 ms level even when fit to the fifth frequency derivative.