Exploring the effect of periastron advance in small-eccentricity binary pulsars

ABSTRACT

Short-orbital period small-eccentricity binary pulsars can, in principle, experience substantial advance of periastron. We explore the possibility of measuring this effect by implementing a timing model, ELL1k, in the popular tempo2 pulsar timing package. True secular variations in the Laplace–Lagrange parameters, present in our ELL1k model, can lead to measurable timing residuals while pursuing decade-long timing campaigns using the existing ELL1 timing model of Lange et al. (2001), especially for binaries exhibiting significant periastron advance. We also list the main differences between our approach and various implementations of the ELL1 model present in both tempo and tempo2 packages. Detailed tempo2 simulations suggest the possibility of constraining the apsidal motion constant of pulsar companions in certain observed binary pulsars with minuscule eccentricities such as PSR J1719−1438. Fortunately, the ELL1k timing model does not pose any challenges to the on-going Pulsar Timing Array campaigns that routinely employ the ELL1 timing model.

gravitation, binaries: general, pulsars: general

1 INTRODUCTION

Presently, binary pulsars provide the most accurate laboratories to test general relativity (GR) in quasi-stationary strong field regime (Wex 2014). This is mainly because of the technique of pulsar timing which requires an accurate prescription for determining the pulse phase as a function of time. In the case of binary pulsars like PSR B1913+16, this technique essentially provided an ideal clock for probing the nature of relativistic gravity (Taylor, Fowler & McCulloch 1979). Pulsar timing demonstrated the decay of orbital period in neutron star binaries, which provided the first observational evidence for the existence of gravitational waves (Taylor 1993). Additionally, the on-going Pulsar Timing Array (PTA) experiments aim to detect low-frequency gravitational waves in the 10⁻⁹–10⁻¹⁰ Hz frequency range (Verbiest et al. 2016). These efforts require accurate timing of millisecond pulsars (MSPs) due to their exquisite rotational stability.

Roughly, 10 per cent of the over 2600 currently known pulsars exist in binary systems with companions in various stages of stellar evolution (Manchester et al. 2005). Accurate timing of pulsars in such systems requires prescriptions to model delays in the times of arrival (TOAs) of pulses in an inertial barycentric frame due to the orbital motion of the pulsar and its companion (Blandford & Teukolsky 1976; Edwards, Hobbs & Manchester 2006). For double neutron star binaries in eccentric orbits, a number of relativistic orbital effects contribute to such delays (Damour & Deruelle 1986). This ensures that the binary orbit is specified by a number of post-Keplerian parameters in addition to the regular Keplerian parameters, namely, the orbital period P_a, argument of periastron ω, eccentricity e, time of periastron passage T₀, and the projected semimajor axis x. Timing of binary pulsars allow us to probe quasi-stationary strong gravitational fields due to our ability to measure their post-Keplerian parameters. However, timing of MSPs in binaries requires special care as a good fraction of them are in near-circular orbits.

For pulsar binaries with tiny orbital eccentricities, the TOAs do not prominently depend on Keplerian parameters ω and T₀. This results in large uncertainties in the usual χ² estimation of ω and T₀. This prompted Norbert Wex to describe such orbits in terms of the first and second Laplace–Lagrange (LL) parameters and certain time of ascending node passage |$T_{\rm ascnode }$| (Lange et al. 2001; Edwards et al. 2006). These parameters, namely |$\epsilon _1=e\, \sin \omega$|⁠, |$\epsilon _2=e\, \cos \omega$|⁠, and |$T_{\rm ascnode }=T_{0} - \frac{\omega }{n}$|⁠, where n is the mean motion, replace the regular Keplerian parameters e, ω, and T₀. Clearly, the use of rectangular components of the eccentricity vector, given by |$e\, \cos \omega$| and |$e\, \sin \omega$|⁠, to represent the periastron of the elliptical orbit is influenced by its use in celestial mechanics (Pannekoek 1948).

The ELL1 timing model, detailed in Lange et al. (2001), incorporates the effects of Rømer and Shapiro delays while neglecting the Einstein delay contributions. Additionally, linear-in-time variations of the sidereal orbital period, Laplace–Lagrange parameters and projected semimajor axis were introduced to model secular changes in these parameters. In the tempo2 implementation of this timing model, one employs ε₁₀, ε₂₀, |$\dot{\epsilon }_1$|⁠, |$\dot{\epsilon }_2$|⁠, and |$T_{\rm ascnode }$| as fitting parameters in the place of e₀, ω₀, |$\dot{e}$|⁠, |$\dot{\omega }$|⁠, and T₀ employed in the DD timing model for eccentric binaries (Damour & Deruelle 1986; Edwards et al. 2006). This is what one gathers from a detailed study of the appendix in Lange et al. (2001) and its tempo2 implementation (Edwards et al. 2006). It should be noted that the ELL1 model incorporates only the first order terms in orbital eccentricity. Very recently, it was pointed out by Zhu et al. (2018) that higher order eccentricity contributions to the Rømer delay should be relevant while timing nearly circular wide-orbit binary pulsars. These corrections were implemented in the ELL1+ model, an extension of ELL1. However, TOA measurement errors are expected to be much larger than the higher order e corrections to the Rømer delay for short orbital period binaries: systems of our current interest.

In this paper, we explore the implications of restricting the temporal variations of LL parameters to be linear-in-time, pursued in Lange et al. (2001) and implemented in tempo2 (Edwards et al. 2006). It turns out that the time evolutions of rectangular components of the orbital eccentricity vector are more general. This is because the rate of periastron advance per orbit need not be negligible for compact pulsar binaries with tiny orbital eccentricities as evident from equations (1–3) in Willems, Vecchio & Kalogera (2008). The more general time evolutions for the LL parameters can have observational implications due to the following reasons. These days it is indeed possible to measure TOAs with ∼100 ns uncertainties due to the advent of real-time coherent de-dispersion (Hankins & Rajkowski 1987) and availability of large collecting area telescopes like the Giant Metre-wave Radio Telescope equipped with 400 MHz wide-band receivers (Reddy et al. 2017). The TOA precision is likely to improve further with the advent of telescopes like the Five hundred metre Aperture Spherical Telescope (FAST; Nan et al. 2011) and the Square Kilometre Array (SKA; Combes 2015). Secondly, the secular systematics in timing residuals are expected to increase with observation span for timing campaigns lasting decades, such as those employed in various PTA experiments (Verbiest et al. 2016). This is relevant for us as many PTA millisecond pulsars are in small-eccentricity binary orbits and are timed by the ELL1 timing model. The final consideration is based on the possibility that the present investigation may allow us to extract additional astrophysical information from certain binary pulsar systems like PSR J1719−1438 (Bailes et al. 2011). This is because the classical contributions to the periastron advance may dominate over its general relativistic counterpart in a number of recently discovered MSP binaries. If this is indeed the case, the long-term timing of such tiny-eccentricity pulsar binaries may allow us to determine the apsidal motion constant of pulsar companions in such systems.

After completion of our numerical tempo2 experiments, we came to know that Norbert Wex had implemented a similar extension to his ELL1 timing model in tempo (Taylor et al. 2013). Unfortunately, no public documentation exists for this implementation. A close inspection reveals that our timing model and Wex’s tempo implementation¹ that extends his ELL1 timing model, detailed in Lange et al. (2001), differs by a term. Our tempo2 simulations reveal that this additional term is crucial for small-eccentricity binary pulsars experiencing substantial periastron advance and therefore should not be ignored.

In Section 2, we briefly summarize the ELL1 timing model. How we incorporate ‘exact’ temporal evolutions of the LL parameters and why long-term monitoring will be required to distinguish such evolutions from the present tempo2 evolutions for ε₁ and ε₂ are explained in Section 3. Why PSR J1719−1438 should be the most promising source to observe the effects of periastron advance is explained in Section 4. This section also details how we adapted the tempo2 software² (Edwards et al. 2006; Hobbs, Edwards & Manchester 2006) to explore observational consequences of the general time evolutions for ε₁ and ε₂. In Section 5, we summarize our results and discuss their implications. Expressions required to incorporate the ELL1k model in tempo2 are given in Appendix  A and a brief comparison of different timing models for near-circular binaries is given in Appendix  B.

2 THE ELL1 TIMING MODEL OF NORBERT WEX

In the ELL1 timing model, the barycentric arrival time of a pulse, t_obs, is related to its time of emission, t_em, by

\begin{eqnarray*} t_{\text{obs}}=t_{\text{em}} + \Delta _{\mathrm{ R}}(t_{\text{em}}) \, , \end{eqnarray*}

(1)

where Δ_R(t_em) is the Rømer delay due to the orbital motion of the pulsar around its companion. The standard expression for the Rømer delay associated with an eccentric Keplerian orbit is given by (Blandford & Teukolsky 1976)

\begin{eqnarray*} \Delta _{\mathrm{ R}}=\frac{x}{c}\left[\left(\cos u - e\right)\sin \omega + \sqrt{1-e^{2}}\sin u \, \cos \omega \right]\, , \end{eqnarray*}

(2)

where x = a_psin ι is the semimajor axis of the pulsar orbit projected on to the line of sight and ι is the orbital inclination. The Keplerian parameters are the eccentric anomaly u, orbital eccentricity e, and the argument of periastron ω. The semimajor axis of the pulsar orbit, a_p, is related to the semimajor axis of the relative orbit a by the standard relation |$a_\mathrm{ p}=\frac{m_\mathrm{ c}}{M} a$| where the total mass of the binary M = m_p + m_c is given in terms of the pulsar mass m_p and the companion mass m_c. The relativistic version of equation (2) may be found in Damour & Deruelle (1986).

It is possible to obtain an explicit analytic expression for the temporal evolution of Δ_R in the small eccentricity limit. This is due to the following standard expressions for sin u and cos u in terms of the mean anomaly l (Ch. 2 Brouwer & Clemence 1961)

\begin{eqnarray*} \cos u &=&\frac{-e}{2}+2\sum _{s=1}^{\infty }\frac{1}{s}\text{J}_{s}^{^{\prime }}(se)\cos sl \nonumber \\ &=&\cos l + \frac{e}{2}(\cos 2l-1)+\mathcal {O}(e^2)\, , \\ \sin u &=&\frac{2}{e}\sum _{s=1}^{\infty }\frac{1}{s}\text{J}_{s}(se)\sin sl \nonumber \\ &=&\sin l + \frac{e}{2} \sin 2l + \mathcal {O}(e^2) \, , \end{eqnarray*}

(3),(4)

where l is defined in terms of the mean motion n and the epoch of periastron passage T₀ as l = n(t − T₀) while J_s(se) and |$\text{J}_{s}^{^{\prime }}(se)$| stand for the Bessel function of the first kind and its derivative. General relativistic corrections to the above expressions, accurate to third post-Newtonian order, are available in Boetzel et al. (2017). The mean motion n is related to the orbital (anomalistic) period by n = 2π/P_a. Wex introduced a certain phase Φ = l + ω such that equation (2) with the help of equations (3) and (4) becomes (Lange et al. 2001)

\begin{eqnarray*} \Delta _{\mathrm{ R}}=\frac{x}{c} \left[\sin \Phi +\frac{e}{2}\left(\sin 2\Phi \:\cos \omega -\left(\cos 2\Phi +3\right)\sin \omega \right)\right] \end{eqnarray*}

(5)

The above expression, which is accurate to linear order in e, takes a simpler form in terms of the Laplace–Lagrange parameters ε₁ = e sin ω and ε₂ = e cos ω, and it reads

\begin{eqnarray*} \Delta _{\mathrm{ R}}=\frac{x}{c}\left[\sin \Phi +\frac{1}{2}\biggl (\epsilon _{2}\sin 2\Phi -\epsilon _{1}\left(\cos 2\Phi +3\right)\biggr)\right]. \end{eqnarray*}

(6)

A comparison with equation (A6) in Lange et al. (2001) reveals the presence of an additional term |$\frac{-3}{2} \frac{x}{c} \epsilon _1$| which was neglected in the ELL1 timing model. We observe that this term is absent in both tempo2 and tempo implementations of Wex’s timing model for small-eccentricity binary pulsars. To obtain explicit temporal evolution for Φ in terms of well-defined parameters, Wex introduced certain time of ascending node|$T_{\rm ascnode }$| to be

\begin{eqnarray*} T_{\rm ascnode }=T_{0}-\frac{\omega }{n}\, , \end{eqnarray*}

(7)

such that the time evolution for the phase is given by

\begin{eqnarray*} \Phi = l+\omega =n\left(t-T_{0}+\frac{\omega }{n} \right)=n \left(t- T_{\rm ascnode } \right) \, . \end{eqnarray*}

(8)

The introduction of |$T_{\rm ascnode }$| ensures that the usual Keplerian parameters, e, ω, and T₀ are replaced by parameters ε₁, ε₂, and |$T_{\rm ascnode }$|⁠, respectively, which are more appropriate for binaries having tiny orbital eccentricities.

When advance of periastron is present, it is convenient to introduce certain sidereal angular frequency |$n_\mathrm{ b}=n + \dot{\omega }$| such that the secular evolution of the phase becomes

\begin{eqnarray*} \Phi = n_{\mathrm{ b}}\left(t-T_{\rm ascnode }\right)+\frac{1}{2}\dot{n}_{\mathrm{ b}}\left(t-T_{\rm ascnode }\right)^{2} \, , \end{eqnarray*}

(9)

where we have also incorporated the orbital frequency derivative |$\dot{n}_\mathrm{ b}=\dot{n}$| due to the dissipative evolution of the orbit. The secular variations to x, ε₁, and ε₂ are provided in the tempo2 implementation of ELL1 by the following relations (Lange et al. 2001)

\begin{eqnarray*} x &=& x_0 + \dot{x} \, \left(t-T_{\rm ascnode }\right) \, , \nonumber \\ \epsilon _1 &=& \epsilon _{10} + \dot{ \epsilon _1} \, \left(t-T_{\rm ascnode }\right) \, , \nonumber \\ \epsilon _2 &=& \epsilon _{20} + \dot{ \epsilon _2} \, \left(t-T_{\rm ascnode }\right)\, , \end{eqnarray*}

(10)

where |$T_{\rm ascnode }$| is now defined as |$T_{\rm ascnode }=T_0 - \omega _0/n_b$|⁠, and ω₀, x₀, ε₁₀, and ε₂₀ are the values of the corresponding parameters at |$t=T_{\rm ascnode }$|⁠.

This ensures that ELL1 model in tempo2 employs ε₁₀, ε₂₀, |$\dot{\epsilon }_1$|⁠, |$\dot{\epsilon }_2$|⁠, and |$T_{\rm ascnode }$| as fittable parameters in addition to x₀, and |$\dot{x}$| if required. We note, in passing, that Einstein delay is not relevant for these systems while the Shapiro delay expression of Lange et al. (2001) is not altered by our considerations.

We would like to state explicitly that a prescription for the general temporal evolutions of the LL parameters in terms of |$\dot{e}$| and |$\dot{\omega }$| is available in the tempo software. Unfortunately, this prescription is missing in the currently popular tempo2 software, and as noted earlier no public documentation exists for the tempo2 implementation that generalizes the above-detailed ELL1 timing model. In what follows, we detail how we independently generalized linear-in-time evolutions for ε₁ and ε₂ while including the above mentioned additional term in the expression for the Rømer delay.

3 GENERALIZING THE LINEAR-IN-TIME EVOLUTIONS FOR ε₁ AND ε₂

In this section, we describe a simple prescription to model the ‘exact’ time evolutions for ε₁ and ε₂ due to linear-in-time evolutions for e and ω with the aim of extracting the effect of |$\dot{\omega }$| from binary pulsar timing. It turns out that the effect of periastron advance is implicitly present in the ELL1 timing model and this forces the conservative phase evolution to be |$\Phi=n_\mathrm{ b} \left(t - T_{\rm ascnode } \right)$| where |$n_\mathrm{ b}=n + \dot{\omega }$|⁠. However, it is fairly difficult to constrain |$\dot{\omega }$| from timing observations. This is because the times of ascending node passages are provided in terms of |$P_\mathrm{ b}=2\, \pi /n_\mathrm{ b}$| such that the epoch of Nth ascending node passage reads

\begin{eqnarray*} T_{\rm ascnode }^{(N)}=T_{\rm ascnode }^{(0)} + N\, P_{\mathrm{ b}} + \frac{1}{2}N^{2}P_{\mathrm{ b}}\dot{P}_{\mathrm{ b}}\, , \end{eqnarray*}

(11)

where we imposed linear-in-time evolution for n_b. Therefore, one may extract the time derivative of P_b by measuring |$T_{\rm ascnode }$| values at widely separated epochs as done, for example, in the case of the accreting millisecond X-ray pulsar SAX J1808.4−3658 (Patruno et al. 2012). In terms of P_b, time evolution of Φ is given by

\begin{eqnarray*} \Phi = 2\pi \left[\frac{\tau }{P_{\mathrm{ b}}}-\frac{1}{2}\dot{P}_{\mathrm{ b}}\left(\frac{\tau }{P_{\mathrm{ b}}}\right)^{2}\right]\, , \end{eqnarray*}

(12)

where we have defined |$\tau =t-T_{\rm ascnode }$| and neglected higher derivatives of P_b.

It is straightforward to show that the general temporal evolutions for the Laplace–Lagrange parameters are given by

\begin{eqnarray*} \epsilon _{1} (t) &=& (1+\xi \tau)\left(\epsilon _{10}\cos \dot{\omega }\tau +\epsilon _{20}\sin \dot{\omega }\tau \right), \nonumber \\ \epsilon _{2} (t) &=& (1+\xi \tau)\left(\epsilon _{20}\cos \dot{\omega }\tau -\epsilon _{10}\sin \dot{\omega }\tau \right), \end{eqnarray*}

(13)

where |$\xi =\dot{e}/e_0$| and the subscript 0 stands for parameter values at τ = 0 such that

\begin{eqnarray*} \epsilon _{10} &=& e_0 \sin \omega _{0} \, ,\nonumber \\ \epsilon _{20} &=& e_0 \cos \omega _{0} \, . \end{eqnarray*}

(14)

Note that we let ω and e to vary linear-in-time, namely |$\omega=\omega _0 + \dot{\omega }\, \tau$| and |$e=e_0 + \dot{e} \tau$|⁠, to obtain the above time evolutions for ε₁ and ε₂. Further, close inspection of de/dt expressions due to dominant order gravitational radiation reaction effects and tidal dissipation allows us to impose ξ = 0 for these binaries (Peters 1964; Zahn 1978). For example, |$\dot{e}$| contributions, due to gravitational wave emission, appear at (⁠|$v$|/c)⁵ order and are small compared to the periastron advance contribution which occurs at (⁠|$v$|/c)² order where |$v$| is the orbital velocity. This ensures that the time-scale ξ⁻¹ associated with |$\dot{e}$| will be substantially higher than the time-scale relevant for the periastron advance. However, |$\dot{e}$| contributions should not be neglected while testing strong equivalence principle using small-eccentricity long-orbital period binary pulsars (Damour & Schäfer 1991). Neglecting ξ, equation (13) becomes

\begin{eqnarray*} \epsilon _{1} (t) &=& \epsilon _{10}\cos \dot{\omega }\tau +\epsilon _{20}\sin \dot{\omega }\tau , \nonumber \\ \epsilon _{2} (t) &=& \epsilon _{20}\cos \dot{\omega }\tau -\epsilon _{10}\sin \dot{\omega }\tau . \end{eqnarray*}

(15)

Indeed, for sufficiently small values of |$\dot{\omega }\tau$|⁠, the equation (15) can be written as

\begin{eqnarray*} \epsilon _{1} (t) &=& \epsilon _{10} + \epsilon _{20} \,\, \dot{\omega }\tau \, ,\nonumber \\ \epsilon _{2} (t) &=& \epsilon _{20} - \epsilon _{10} \,\, \dot{\omega }\tau \, . \end{eqnarray*}

(16)

The above equations are identical to the expressions for ε₁(t) and ε₂(t) present in the ELL1 model. This is because it is straightforward to show that |$\dot{\omega }$| is given either by |$\dot{\epsilon }_1/ \epsilon _{20}$| or |$-\dot{\epsilon }_2 / \epsilon _{10}$| with the help of equations (A14) and (A15) of Lange et al. (2001) while equating |$\dot{e}=0$|⁠.

The above arguments demonstrate that the ELL1 timing model, detailed in Lange et al. (2001), accounts only for linear-in-time variations of the LL parameters. This approximation may not be appropriate while pursuing timing campaigns spanning decades such that |$\tau \gtrsim \tau _{\dot{\omega }} \equiv {2\pi }/{\dot{\omega }}$|⁠. Additionally, we have included the |$\frac{-3}{2}\, \frac{x}{c}\, \epsilon _1$| term in our expression for the Rømer delay. In contrast, this term was neglected in Lange et al. (2001) as well as its tempo2 (Edwards et al. 2006) and tempo2 implementations. This term turned out to be crucial while invoking equation (15) for evolving the LL parameters of systems with significant periastron advance.

However, we would like to emphasize that the ELL1 model should be quite appropriate to fit the observed TOAs when the span of the timing campaign is much smaller than |$\tau _{\dot{\omega }}$|⁠. To demonstrate this, we write the barycentric arrival time of a pulse as

\begin{eqnarray*} t_{\text{obs}}=t_{\text{em}}+\Delta _{\mathrm{ R}}\, . \end{eqnarray*}

(17)

For nature’s perfect clocks, we have

\begin{eqnarray*} t_{\text{em}}=t_{\text{em}}^{(0)}+\frac{N}{f}\, , \end{eqnarray*}

(18)

where |$t_{\text{em}}^{(0)}$| is the time of emission of a ‘reference’ pulse, f is the pulsar rotation frequency, and N is some integer while ignoring any variations in pulsar frequency f. The expression for the Rømer delay, namely equation (6), may be written with the help of equation (16) as

\begin{eqnarray*} \Delta _{\mathrm{ R}} &=& \frac{x}{c}\left[\sin \Phi +\frac{1}{2}\left(\left(\epsilon _{20}+\dot{\epsilon }_{2}\tau \right)\sin 2\Phi -\left(\epsilon _{10}+\dot{\epsilon }_{1}\tau \right)\cos 2\Phi \right)\right] \nonumber \\ &&- \frac{3}{2}\frac{x}{c}\epsilon _{1} \nonumber \\ &=& \Delta _{\mathrm{ R}}^{({\tt family\, ELL1})}-\frac{3}{2}\frac{x}{c}\epsilon _{1}\, , \end{eqnarray*}

(19)

where we used the relations |$\dot{\epsilon }_1=\epsilon _{20}\dot{\omega }$| and |$\dot{\epsilon }_2=-\epsilon _{10}\dot{\omega }$|⁠. In the above equation, |$\Delta _{\mathrm{ R}}^{({\tt {ELL1}})}$| stands for the Rømer delay expressions present in the ELL1 model (Lange et al. 2001) that incorporate only linear-in-time variations for the Laplace–Lagrange parameters. The expression for the Rømer delay now becomes

\begin{eqnarray*} \Delta _{\mathrm{ R}} = \Delta _{\mathrm{ R}}^{({\tt {ELL1}})}-\frac{3}{2}\frac{x}{c}\epsilon _{10}-\frac{3}{2}\frac{x}{c}\epsilon _{20}\dot{\omega }\left(t_{\text{em}}^{(0)}+\frac{N}{f}-T_{\rm ascnode }\right)\, , \end{eqnarray*}

(20)

as |$\tau =t-T_{\rm ascnode }$|⁠. This allows us to write the barycentric arrival times of pulses as

\begin{eqnarray*} t_{\text{obs}} &=& t_{\text{em}}^{(0)}+\frac{N}{f}+\Delta _{\mathrm{ R}}^{({\tt {ELL1}})}-\frac{3}{2}\frac{x}{c}\epsilon _{10}-\frac{3}{2}\frac{x}{c}\epsilon _{20}\dot{\omega }\left(t_{\text{em}}^{(0)}+\frac{N}{f}-T_{\rm ascnode }\right)\nonumber \\ &=& t_{\text{em}}^{(0)}+\frac{N}{f}\left(1-\frac{3}{2}\frac{x}{c}\epsilon _{20}\dot{\omega }\right)+\Delta _{\mathrm{ R}}^{({\tt {ELL1}})} \nonumber \\ && -\,\frac{3}{2}\frac{x}{c}\left(\epsilon _{10}+\epsilon _{20}\dot{\omega }\left(t_{\text{em}}^{(0)}-T_{\rm ascnode }\right)\right) . \nonumber \\ \end{eqnarray*}

(21)

A close inspection reveals that the last term in equation (21) is a constant over time and can be neglected. Additionally, the effect of the ‘extra term’ can be absorbed into the pulsar frequency by redefining it as

\begin{eqnarray*} f^{\prime }=f\left(1-\frac{3}{2} \frac{x}{c} \epsilon _{20} \dot{\omega }\right)^{-1}. \end{eqnarray*}

(22)

The above conclusion holds even if we include |$\skew6\dot{f}$| in equation (18). A straightforward calculation allows us to re-write equation (21) as

\begin{eqnarray*} t_{\text{obs}}=t_{\text{em}}^{(0)}+\frac{N}{f^{\prime }}-\frac{\skew6\dot{f}^{\prime }N^{2}}{2f^{\prime 3}}+\Delta _{\mathrm{ R}}^{(\tt {ELL1})}\, , \end{eqnarray*}

(23)

where we neglected certain constant terms and

\begin{eqnarray*} \skew6\dot{f}^{\prime }=\skew6\dot{f}\left(1-\frac{3}{2}\frac{x}{c}\epsilon _{20}\dot{\omega }\right)^{-2}. \end{eqnarray*}

(24)

The above arguments demonstrate that the ELL1 timing model should be sufficient to analyse timing data, provided |$\tau \ll \tau _{\dot{\omega }}$| for ensuring the validity of equation (16). However, the model may not be appropriate while pursuing long-term timing campaigns of short orbital period pulsar binaries that can accommodate, in principle, non-negligible periastron advance. It should be noted that the higher order corrections to equation (16) may be partially absorbed into higher derivatives of f, namely |$\skew6\dot{f}$|⁠, |$\ddot{f}$| etc., provided |$\tau \ll \tau _{\dot{\omega }}$|⁠. This may be relevant during the coherent timing of pulsar binaries like SAX J1808.4−3658 (Patruno et al. 2012). However, it may not be desirable to treat temporal variations in the LL parameters to be linear-in-time while pursuing long-term timing campaigns of short-period binary pulsars. This is what we pursue in the next section.

4 ELL1k MODEL AND ITS OBSERVATIONAL IMPLICATIONS

To probe the observational consequences of our equations (6) and (15), we implemented a modified version of the ELL1 timing model in tempo2. This model, available in the latest version of tempo2, is referred to as ELL1k due to the use of k to represent the dimensionless fractional periastron advance per orbit (Damour & Schäfer 1988). The relevant expressions, required to implement the ELL1k timing model, are listed in Appendix  A and they do contain contributions from the above mentioned additional term. To probe the relevance of our improvements, we simulate binary pulsar TOAs using our ELL1k model for different observation spans and white timing noise amplitudes using the fake plug-in of tempo2. The resulting TOAs are fitted using the ELL1 model. We explore the goodness of the fit by varying |$\dot{\omega }$| values, observation spans, and the amplitude of the white timing noise. If the fit exhibits systematic variations in the residuals, we conclude that our ELL1k model should be relevant for analysing TOAs from such binary pulsar configurations. Additionally, we explore the implications of red (i.e. non-Gaussian) timing noise in these simulations and our observations are summarized in Section 4.3. We begin by describing why the long-term timing of a recently discovered system should be interesting from our point of view.

4.1 Diamond planet-pulsar binary as a possible test bed for our ELL1k timing model

The recently discovered PSR J1719−1438 is a unique binary millisecond pulsar with an ultra-low mass/planetary companion having an orbital period of about 2.18 h and eccentricity of ∼8 × 10⁻⁴ (Bailes et al. 2011). We list in Table 1 a few relevant orbital parameters, extracted from Ng et al. (2014). Accurate timing enabled Bailes et al. (2011) to demonstrate that the mass of the companion should be greater than |$0.0011292\, \mathrm{M}_{\odot }$| (about 1 Jupiter mass) while its radius is constrained to be less than 28432.9 km (about 40 per cent of Jupiter radius).

Table 1.

Orbital parameters of PSR J1719−1438, obtained using the ELL1 timing model, listed in Ng et al. (2014).

Parameter	Unit	Value
Pb	d	0.09 070 629
x	lt-s	0.0 018 212
|$T_{\rm ascnode }$|	MJD	552 35.516 505
ε1		−7 × 10−4
ε2		4 × 10−4

It turns out that this binary system should exhibit, in principle, the advance of periastron. For PSR J1719−1438, the advance of periastron should arise due to the dominant order general relativistic, static tidal, and classical spin-orbit (SO) interactions. The relevant expressions for |$\dot{\omega }_{\text{GR}}$|⁠, |$\dot{\omega }_{\text{tidal}}$|⁠, and |$\dot{\omega }_{\text{SO}}$| are given by (see e.g. Willems et al. 2008)

\begin{eqnarray*} \dot{\omega }_{\text{GR}} &=& \frac{3}{1-e^{2}}\left(\frac{GM}{c^{3}}\right)^{2/3}\left(\frac{2\pi }{P_{\mathrm{ b}}}\right)^{5/3}\, , \\ \dot{\omega }_{\text{tidal}} &=& 15\left(\frac{R_{\mathrm{ c}}}{a}\right)^{5}\left(\frac{m_{\mathrm{ p}}}{m_{\mathrm{ c}}}\right)\frac{2\pi }{P_{\mathrm{ b}}}k_{2}\left(\frac{1+\frac{3}{2}e^{2}+\frac{1}{8}e^{4}}{\left(1-e^{2}\right)^{5}}\right) \, , \\ \dot{\omega }_{\text{SO}} &=& \frac{2\pi }{P_{\mathrm{ b}}}\left(\frac{R_{\mathrm{ c}}}{a}\right)^{5}\left(\frac{M}{m_{c}}\right)\left(\frac{P_{\mathrm{ b}}}{P_{\mathrm{ c}}}\right)^{2}\frac{k_{2}}{\left(1-e^{2}\right)^{2}}\, , \end{eqnarray*}

(25a),(25b),(25c)

where R_c, P_c and k₂ are respectively the radius, rotational period, and quadrupolar tidal Love number (also known as quadrupolar apsidal motion constant) of the companion.

It should be obvious that general relativistic contributions are well-constrained by the existing observations (Ng et al. 2014). However, additional assumptions will be required to estimate classical contributions to |$\dot{\omega }$|⁠. We may estimate a value for the semimajor axis of the relative orbit a to be 2.1789239 lt-s by assuming an edge-on orbit and a neutron star mass |$m_\mathrm{ p}=1.35 \, \mathrm{M}_{\odot }$| (Manchester et al. 2005). Unfortunately, it is rather difficult to provide a firm estimate for k₂ as we do not know the exact nature of the companion. It is reasonable to expect that the companion’s k₂ may lie between the following two possible ranges. In case the companion is a low-mass white dwarf, k₂ can be between ∼0.01 to ∼0.1 (Valsecchi et al. 2012). Possible k₂ values can be from ∼0.001 to ∼0.01 for main-sequence stars (Stothers 1974) while k₂ value can be as high as ∼0.5 for the Jovian planets (Wahl, Hubbard & Militzer 2016). We choose a rather conservative estimate for the companion’s apsidal motion constant and let k₂ = 0.01, and impose co-rotation of the companion with its orbital motion. This leads to

\begin{eqnarray*} \dot{\omega }_{\text{GR}} &\simeq& 13.3^{\circ }\text{ yr}^{ -1}\, ,\nonumber \\ \dot{\omega }_{\text{tidal}} &\simeq& 40.5^{\circ }\text{ yr}^{ -1}\, ,\nonumber \\ \dot{\omega }_{\text{SO}} &\simeq& 2.7^{\circ }\text{ yr}^{ -1}\, , \end{eqnarray*}

(26)

and the total |$\dot{\omega } \simeq 56.5^{\circ }\text{ yr}^{ -1}$| for PSR J1719−1438.

It is important to note that our |$\dot{\omega }_{\text{tidal}}$| and |$\dot{\omega }_{\text{SO}}$| estimates are extremely sensitive to the ratio R_c/a, as evident from equations (25b) and (25c). For example, if the actual R_c value is 50 per cent of the current upper bound, |$\dot{\omega }_{\text{tidal}}$| and |$\dot{\omega }_{\text{SO}}$| will reduce by a factor of 2⁵ = 32 from the above listed values. This implies that our possible estimates for |$\dot{\omega }_{\text{tidal}}$| and |$\dot{\omega }_{\text{SO}}$| should be treated with a high degree of caution. However, our |$\dot{\omega }_{\text{GR}}$| estimate should be fairly accurate due to its crucial dependence on the accurately measured P_b value. Additionally, the inferred total mass of the binary system is expected to be close to the typical pulsar mass, namely |$1.35\, \mathrm{M}_{\odot }$| as observations strongly suggest the presence of an ultra-low mass companion in PSR J1719−1438 (Bailes et al. 2011). Therefore, we employ our |$\dot{\omega }_{\text{GR}}$| estimate while exploring the observational implications of the ELL1k model for the PSR J1719−1438 binary system.

We display in three figures our simulations for the above binary that spans 2 yr (Fig. 1a–c), 20 yr (Fig. 2a–c), and 100 yr (Fig. 3a–c). For these plots, we simulated one observation each in every 45 d and chose a white timing noise amplitude of 100 ns while assuming the absence of any red timing noise. Timing residuals arise as we employ the standard ELL1 timing model to analyse the TOAs generated with the ELL1k model. Plots in Fig. 1(a)–(c) demonstrate that it will not be possible to identify the effects of |$\dot{\omega }$| if the total duration of the observation is too short compared to the periastron advance time-scale, namely |$\tau _{\dot{\omega }}=2\, \pi / \dot{\omega }$|⁠, which in this case is ∼27 yr while assuming only the GR contribution. Clearly, any systematics arising from our proposed modifications are absorbed in f as was shown in Section 3. However, the effect of |$\dot{\omega }$| manifests as systematic variations in the post-fit residuals when span of the observations is comparable to |$\tau _{\dot{\omega }}$| as evident from plots in Fig. 2(a)–(c). This point is further emphasized in Fig. 3(a)–(c) which show the case where the total observation duration is much greater than the periastron advance time-scale for the PSR J1719−1438 binary system.

The inferences, gleaned from these figures, are clearly consistent with Section 3 conclusions. The linear variation in the timing residuals, present in Fig. 1(b) is fully expected due to the additional term in our equation (19). In comparison, the periodic nature of the timing residuals present in Fig. 3(a)–(c) naturally arise from the ‘exact’ temporal evolutions of the Laplace–Lagrange parameters, namely equation (15). These simulations open up the tantalizing possibility of measuring the effects of |$\dot{\omega }$| in the PSR J1719−1438 binary system, provided high-cadence 10 yr timing data is available. It is possible that the effects of |$\dot{\omega }$| may barely be detectable now as the system was discovered in 2011 (Bailes et al. 2011), especially if there exists significant contributions from the tidal and classical spin-orbit interactions. Obviously, this requires the availability of high-precision TOAs and the absence of significant red timing noise. The present simulations should provide sufficient motivation to pursue long-term, high-cadence timing of PSR J1719−1438. The cadence of an observation campaign for measuring |$\dot{\omega }$| should be such that it samples a cycle of the periastron (with period |$\tau _{\dot{\omega }}$|⁠) at a sufficient rate. We have checked that the systematic timing residuals, visible in Figs 2(a)–(c) and 3(a)–(c), persist even when the cadence is as sparse as one observation per year. Finally, we would like to note that the 100 ns TOA uncertainties, invoked in our simulations, are usually achievable by the existing PTA telescopes.

Figure 1.

tempo2 simulations of PSR J1719−1438 spanning 2 yr duration. For these simulations, TOAs were generated using the fake plug-in of tempo2 with the ELL1k timing model for parameters listed in Table 1 and we let |$\dot{\omega }$| take its GR value, namely 13.3° yr⁻¹. The top panel (a) shows the pre-fit timing residuals to our data where the fit was done using the ELL1 model for Table 1 parameters. The middle panel (b) shows the post-fit residuals after fitting ε₁, ε₂, |${\dot{\epsilon }}_{1}$|⁠, and |${\dot{\epsilon }}_{2}$| of the ELL1 model. The bottom panel (c) shows the post-fit residuals after additionally fitting for the pulsar spin frequency f. The bottom panel shows that |$\dot{\omega }$| induced variations can be absorbed into the unknown pulsar frequency if |$\tau \ll \tau _{\dot{\omega }}$|⁠.

We observe that the extra term, namely −3xε₁/2c in equation (6), contributed significantly to the long-term oscillatory behaviour of the timing residuals as visible in Figs 2(a) and 3(a). In the absence of this extra term, the resulting timing residuals turned out to be smaller in magnitudes and their long time-scale variations were similar to what we display in Fig. 4. This suggests that the additional term should be relevant while extracting the effect of periastron advance from the timing of binary pulsars like PSR J1719−1438. Additionally, the extra term forces the observed spin frequency to be different from its intrinsic one as evident from equation (22). This may have implications while searching for continuous gravitational waves from pulsars in binaries (Watts et al. 2008). This is because continuous gravitational waves from such systems are expected at 2f rather than at 2f′ and further investigations will be required to quantify its relevance.

An additional point that requires further investigation is related to the inclusion of Shapiro delay. In general, the Shapiro delay is degenerate with the Rømer delay (Lange et al. 2001; Freire & Wex 2010). This ensures that the measured x and ε₁₀ values are different from their intrinsic values. However, our extra term in the Rømer delay can, in principle, break this degeneracy, provided the observation span is sufficiently large and the system exhibits significant periastron advance. This interesting possibility may be relevant for systems like PSR J0348+0432 (Antoniadis et al. 2013) and demands further studies. This may require incorporating |$\dot{\omega }$| in an exact manner in the ELL1H model of Freire & Wex (2010), which expresses Shapiro delay as a Fourier series in the orbital phase Φ for low eccentricity, low-to moderate-inclination binaries. In the next subsection, we explore possible binary pulsar systems where the |$\dot{\omega }$| effect may be measurable in the near future.

4.2 Are there additional ELL1k-relevant pulsar binaries?

A close inspection of equation (25) reveals that pulsar binaries with short orbital periods and non-compact companions should be interesting candidates for the ELL1k timing model. It turns out that there exist a number of tiny-e binaries which can, in principle, exhibit non-negligible |$\dot{\omega }$| as evident from the current ATNF pulsar catalogue³ (Manchester et al. 2005). In Table 2, we display the most promising ELL1k-relevant pulsar binaries present in the catalogue for which the GR-induced periastron advance time-scale is less than 500 yr. Unfortunately, it is rather difficult to provide plausible estimates for the |$\dot{\omega }_{\text{tidal}}$| and |$\dot{\omega }_{\text{SO}}$| contributions in many of these interesting candidates. This is expected due to our ignorance about the likely values for the companion mass m_c, its radius R_c, and the associated Love number k₂.

Table 2.

List of potential rotation-powered binary pulsars with tiny orbital eccentricities where the ELL1k timing model may become relevant. This list shows that there are six systems where |$\tau _{\dot{\omega }_{\text{GR}} }$| is less than 50 yr. There are few PTA pulsars in the list and clearly ELL1k, at present, may not be relevant for their timing (PTA pulsars are marked by the grey background). We use the following notations to denote various companions: UL=ultra-low mass, He=He white dwarf, CO=C–O white dwarf, MS=Main sequence, ??=Unknown. The companion radius arises from the Roche lobe radius arguments only for the first system and its R_c value is marked with the boldface font. The star (⋆) symbol marks the cases where the main-sequence mass-radius relation is used to estimate R_c value.

PSR J	Pb	x	e	mc	Comp.	a	Rc	|$\dot{\omega }_\text{GR}$|	|$\dot{\omega }_\text{tidal}$|	|$\dot{\omega }_\text{SO}$|	|$\dot{\omega }$|	|$\tau _{\dot{\omega }}$|	|$\tau _{\dot{\omega }_\text{GR}}$|
	(d)	(lt-s)		(10−3 M⊙)	Type	(103 km)	(103 km)	(° yr−1)	(° yr−1)	(° yr−1)	(° yr−1)	(yr)	(yr)
1719−1438	0.090706	0.00182	8.0E-4	1.1	UL	653	28	13.3	40.5	2.7	56.5	6.4	27.1
0636+5129	0.066551	0.00899	2.2E-5	6.9	UL	532	26	22.3	17.3	1.2	40.8	8.8	16.1
2339−0533	0.193098	0.61166	2.1E-4	257.2	MS	1146	204⋆	4.2	9.7	0.8	14.7	24.5	85.1
2215+5135	0.172502	0.46814	1.1E-5	207.9	MS	1052	167⋆	5.0	7.5	0.6	13.1	27.5	72.0
0348+0432	0.102424	0.14098	2.4E-6	83.9	He	723	18	11.3	0.0	0.0	11.3	31.8	31.9
0023+0923	0.138799	0.03484	2.4E-5	16.4	UL	871	24	6.6	0.2	0.0	6.8	53.2	54.7
0024−7204I	0.229792	0.03845	6.3E-5	12.9	UL	1218	25	2.8	0.0	0.0	2.9	125.6	127.0
1957+2516	0.238145	0.28335	2.8E-5	96.6	??	1272	17	2.8	0.0	0.0	2.8	129.5	129.5
0751+1807	0.263144	0.39662	3.3E-6	128.3	He	1370	16	2.4	0.0	0.0	2.4	150.8	150.8
1446−4701	0.277666	0.06401	2.1E-5	19.0	UL	1384	23	2.1	0.0	0.0	2.1	173.0	173.6
0610−2100	0.286016	0.07349	3.0E-5	21.4	UL	1412	23	2.0	0.0	0.0	2.0	181.7	182.1
1731−1847	0.311134	0.12016	2.9E-5	33.3	UL	1498	22	1.7	0.0	0.0	1.7	208.2	208.4
1952+2630	0.391879	2.79820	4.1E-5	925.7	CO	2062	3	1.6	0.0	0.0	1.6	219.6	219.6
1816+4510	0.360893	0.59541	7.8E-6	158.2	He	1702	15	1.4	0.0	0.0	1.4	251.9	251.9
1738+0333	0.354791	0.34343	3.4E-7	89.5	He	1657	18	1.4	0.0	0.0	1.4	252.5	252.6
1757−5322	0.453311	2.08653	4.0E-6	556.7	CO	2142	3	1.1	0.0	0.0	1.1	315.0	315.0
0024−7204U	0.429106	0.52695	1.5E-4	122.8	He	1895	16	1.1	0.0	0.0	1.1	341.5	341.5
2214+3000	0.416633	0.05908	8.2E-6	13.3	UL	1811	24	1.1	0.0	0.0	1.1	341.6	342.3
1431−4715	0.449739	0.55006	2.3E-5	124.3	He	1956	16	1.0	0.0	0.0	1.0	369.1	369.1

Timing of tiny eccentric binary pulsars usually yields an estimate for the binary mass function in the absence of any measurable Shapiro delay. Therefore, an estimate for the companion mass requires us to impose edge-on orbital orientation and let the pulsar mass be |$m_\mathrm{ p}=1.35\, \, \mathrm{M}_{\odot }$|⁠. This arises from the definition of the mass function

\begin{eqnarray*} f(m_{\mathrm{ c}})=\frac{4\pi ^{2}}{G}\frac{\left(a_{\mathrm{ p}}\sin \iota \right)^{3}}{P_{\mathrm{ b}}^{2}}=\frac{\left(m_{\mathrm{ c}}\sin \iota \right)^{3}}{M^{2}}\, . \end{eqnarray*}

(27)

Measurement of the Rømer delay provides an estimate for |$a_\mathrm{ p}\, \sin \iota$| and P_b which allows us to estimate m_c with the above assumptions. This is how we obtain the m_c values listed in Table 2.

To obtain an estimate for the companion’s radius, we need to specify its nature. If the companion is believed to be a white dwarf, we invoke the well-known white dwarf mass-radius relation to compute R_c, given m_c. Specifically, we employ the following numerical approximation to the white dwarf mass-radius relation (Carvalho, Marinho & Malheiro 2015)

\begin{eqnarray*} \frac{m}{\mathrm{M}_{\odot }}=\frac{\alpha R + \beta }{e^{\gamma R^2}+\delta }\, , \end{eqnarray*}

(28)

where α = 2.325 × 10⁻⁵ km⁻¹, β = 0.4617, γ = 7.277 × 10⁻⁹ km⁻², and δ = −0.644. A caveat is that the white dwarf mass-radius relation may not be appropriate for ultra-low mass companions, and may give unrealistically large radii. Influenced by Bailes et al. (2011), we place upper bounds on companion radii by arguing that the companion must be smaller than its Roche lobe. In such cases, we invoke the standard radius of the Roche lobe for ultra-low mass companion (Paczynski 1971)

\begin{eqnarray*} R_{\mathrm{ L}}=0.462 \, a\left(\frac{m_{2}}{M}\right)^{\frac{1}{3}}\, . \end{eqnarray*}

(29)

If the companion radius given by equation (28) is larger than the above upper bound, given by equation (29), we employ the lower R_c estimate in our calculations. We highlight such cases by using boldface font for R_c in Table 2. The third category involves the possibility that pulsar companions can be main-sequence stars. In this case, we employ the following empirical main-sequence mass-radius relation, relevant for low-mass stars, to estimate the companion radius (Demircan & Kahraman 1991).

\begin{eqnarray*} \frac{R}{\mathrm{R}_{\odot }} = 1.06 \left(\frac{m}{\mathrm{M}_{\odot }}\right)^{0.945} \, . \end{eqnarray*}

(30)

These cases are marked by the star (⋆) symbol in the R_c column of Table 2.

The third parameter, required to obtain |$\dot{\omega }_{\text{tidal}}$| and |$\dot{\omega }_{\text{SO}}$| estimates, is clearly the apsidal motion constant k₂. We let k₂ = 0.01 for compact degenerate companions and k₂ = 0.001 for main-sequence companions. Unfortunately, it is rather impossible to obtain an estimate for the rotational period of pulsar companions. Therefore, we choose P_c to be identical to the binary orbital period which requires tidal locking. Clearly, these assumptions ensure that the listed |$\dot{\omega }_{\text{tidal}}$| and |$\dot{\omega }_{\text{SO}}$| values in Table 2 should be taken as very rough estimates.

We gather from Table 2 that there are six pulsars for which |$\dot{\omega }$| may be observable within ∼50 yr of their discovery. The effect of periastron advance in these systems, namely PSRs J1719−1438 (Bailes et al. 2011), J0636+5129 (Stovall et al. 2014), J2339−0533 (Romani & Shaw 2011; Ray et al. 2014), J2215+5135 (Hessels et al. 2011), J0348+0432 (Antoniadis et al. 2013; Lynch et al. 2013), and J0023+0923 (Hessels et al. 2011), may be more predominant if they can have non-negligible |$\dot{\omega }_{\text{tidal}}$| and |$\dot{\omega }_{\text{SO}}$| contributions. In particular, J1719−1438 and J0636+5129 may well have shorter periastron advance time-scales (say <5 yr), provided they have significant tidal and spin-orbit contributions to |$\dot{\omega }$|⁠. Unfortunately, J2215+5135 is a redback pulsar for which it is difficult to obtain a coherent timing solution over long time spans (Abdo et al. 2013) and therefore the ELL1k model may not be relevant for this binary. High-cadence timing of these pulsars (except J2215+5135) should allow us to obtain observationally relevant bounds on |$\dot{\omega }$|⁠. For these pulsars, measurement of |$\dot{\omega }$| should lead to constraints on the companion’s radius and/or Love number.

A number of pulsars listed in Table 2 are employed in the present PTA experiments for detecting nanohertz gravitational waves (Verbiest et al. 2016). Fortunately, none of these pulsars have periastron advance time-scale less than 100 yr, and therefore we do not expect advance of periastron to be relevant for Pulsar Timing Arrays. In Table 2, PTA pulsars are marked with grey background. In what follows, we invoke a more realistic scenario where red noise is present in our simulated TOA measurements to probe the observational feasibility of our above conclusions.

4.3 On the effects of Red Timing Noise for our tempo2 simulations

In this section, we investigate how the presence of red timing noise (RN) affects the measurability of |$\dot{\omega }$|⁠. The RN refers to unexplained slow, stochastic wandering of TOAs observed in many pulsars (Arzoumanian et al. 1994). These modulations have a power-law spectrum with a low-frequency cut-off and their power spectral density takes the form (Lasky et al. 2015)

\begin{eqnarray*} \Phi _{\text{RN}}(f) = A\left(1+\frac{f^{2}}{f_{\mathrm{ c}}^{2}}\right)^{-{\rho }/{2}} \, , \end{eqnarray*}

(31)

where A is the spectral density amplitude, f_c is the low-frequency cut-off and ρ is the power-law index. The spectral power associated with the above spectral density is

\begin{eqnarray*} P_{\text{RN}}=\int _0^\infty \Phi _{\text{RN}}(f) \, \mathrm{ d}f=Af_{\mathrm{ c}}\frac{\sqrt{\pi }\Gamma \left(\frac{\rho -1}{2}\right)}{2\Gamma \left(\frac{\rho }{2}\right)} \end{eqnarray*}

(32)

when ρ > 1.

The long-time-scale temporal variations in timing residuals, visible in Fig. 2(c), can be mimicked by TOA measurements affected by RN. In Fig. 4, we demonstrate the consequence of RN while probing the implications of our ELL1k timing model. To obtain Fig. 4, we assume the presence of RN in the simulation data that provided Fig. 2(c). Thereafter, we apply harmonic whitening, detailed in Hobbs et al. (2004), to remove it. This ensures that the long-time-scale variation due to the ‘extra term’ in equation (6) will also be removed as evident from Fig. 4. Interestingly, there are traces of systematic variations present in the residuals of Fig. 4. These variations in timing residuals may be associated with ε₂ sin (2Φ) and ε₁ cos (2Φ) terms in equation (6) which vary in the orbital time-scale while the overall envelope reflects the periastron advance time-scale. These inferences force us to speculate a possible origin, in principle, for the observed red timing noise present in tiny-e binary pulsars. It will be interesting to explore what fraction of the red noise arises from the unaccounted effect of |$\dot{\omega }$| while determining their ephemeris.

Figure 2.

We repeat the tempo2 analysis of Fig. 1(a)–(c) for a time span of 20 yr. Fig. 2(c) shows that the general evolution of the LL parameters cannot be accounted for by the ELL1 timing model even after fitting for f and its time derivative |$\skew6\dot{f}$|⁠.

Figure 3.

We repeat the tempo2 simulations for PSR J1719−1438 for a time window of 100 yr. The three plots were obtained similarly to the Fig. 1(a)–(c). The effect of periastron advance is impossible to miss in the residuals while timing PSR J1719−1438 using the existing ELL1 timing model.

Figure 4.

Post-fit timing residuals after applying harmonic whitening to the simulated data responsible for Fig. 2(c) (Reduced χ² = 8.89).

We move on to quantify how accurately we can measure |$\dot{\omega }$| in the presence of RN. For this purpose, we injected several instances of red noise with different spectral parameters with the help of the simRedNoise plug-in of tempo2 into TOAs simulated using the fake plug-in. The simulations, as expected, were done for a PSR J1719−1438 like system while incorporating the effect of |$\dot{\omega }$| using the ELL1k timing model. Invoking the FITWAVES routine of tempo2 (Hobbs et al. 2004), we applied the harmonic whitening to the simulated TOAs in order to remove the red noise. Thereafter, we tried to recover the |$\dot{\omega }$| value by fitting the whitened data with the ELL1k timing model. We also tried to fit the whitened data with the ELL1 model to see if the systematic timing residuals due to our modifications are washed out by the red noise. The results of such tempo2 experiments are summarized in Table 3 where we also list the parameters of the injected red noise.

Table 3.

Results which probe the effect of red noise (RN) for a PSR J1719−1438 like system. The simulations were obtained by injecting several instances of red noise into TOAs generated using a white noise amplitude of 100 ns and |$\dot{\omega }$| value of 13.3° yr⁻¹. Total simulated duration is 28 yr with one observation in every 45 d. The listed results were obtained by fitting the simulated TOAs with different configurations – ELL1 without whitening, ELL1 after harmonic whitening, ELL1k without whitening and ELL1k after harmonic whitening using 15 harmonics. FW stands for harmonic whitening using FITWAVES. For fits using ELL1k, the measured |$\dot{\omega }$| values (⁠|$\dot{\omega }_{\text{fit}}$|⁠) and the 1σ uncertainty reported by tempo2 (⁠|$\sigma _{\dot{\omega }_{\text{fit}}}$|⁠) are also listed. |$\varepsilon _{\dot{\omega }}$| is the percentage error in |$\dot{\omega }_{\text{fit}}$| compared to the true value of |$\dot{\omega }$| computed as |$\left|\frac{\dot{\omega }_{\text{fit}}-\dot{\omega }_{\text{true}}}{\dot{\omega }_{\text{true}}}\right|\times 100 \%$|⁠. Values quoted for ELL1 and ELL1+FW were obtained after fitting for f, |$\skew6\dot{f}$|⁠, ε₁₀, ε₂₀, |$\dot{\epsilon }_{1}$|⁠, and |$\dot{\epsilon }_{2}$|⁠. Values quoted for ELL1k and ELL1k+FW were obtained after fitting for f, |$\skew6\dot{f}$|⁠, ε₁₀, ε₂₀, and |$\dot{\omega }$| with an initial guess value of 12° yr⁻¹ for |$\dot{\omega }$|⁠. A careful inspection of various columns of reduced χ² (χ²/dof) reveals the usefulness of the ELL1k model even in the presence of substantial red noise.

Injected RN parameters			Fitting Results
PRN	fc	ρ	ELL1	ELL1+FW	ELL1k				ELL1k+FW
			χ2/dof	χ2/dof	χ2/dof	|$\dot{\omega }_{\text{fit}}$|	|$\sigma _{\dot{\omega }_{\text{fit}}}$|	|$\varepsilon _{\dot{\omega }}$|	χ2/dof	|$\dot{\omega }_{\text{fit}}$|	|$\sigma _{\dot{\omega }_{\text{fit}}}$|	|$\varepsilon _{\dot{\omega }}$|
(s2 yr−2)	(yr−1)					(° yr−1)	(° yr−1)	(%)		(° yr−1)	(° yr−1)	(%)
0	–	–	67.4	15.5	0.9	13.41	0.06	0.82	0.9	13.39	0.07	0.67
1E-29	0.100	3.5	68.4	16.8	1.1	13.43	0.07	0.97	0.9	13.48	0.07	1.32
1E-29	0.037	3.5	60.9	16.1	1.1	13.17	0.07	0.95	1.0	13.27	0.08	0.25
1E-29	0.010	3.5	67.1	16.0	1.1	13.18	0.06	0.88	1.1	13.15	0.08	1.12
1E-28	0.100	3.5	59.3	16.9	4.8	13.12	0.14	1.35	1.1	13.33	0.08	0.23
1E-28	0.037	3.5	86.1	16.7	1.7	13.37	0.08	0.54	1.0	13.33	0.07	0.21
1E-28	0.010	3.5	73.5	16.4	1.2	13.34	0.07	0.31	1.2	13.24	0.08	0.47
1E-27	0.100	3.5	72.9	16.3	31.6	13.31	0.45	0.09	1.8	13.35	0.11	0.39
1E-27	0.037	3.5	121.2	16.2	9.0	14.30	0.17	7.55	1.2	13.39	0.08	0.69
1E-27	0.010	3.5	81.8	16.0	2.1	13.14	0.09	1.23	1.0	13.19	0.08	0.79
1E-26	0.100	3.5	748.6	22.6	224.8	20.74	0.44	55.90	8.0	13.38	0.22	0.62
1E-26	0.037	3.5	111.7	18.2	85.9	13.94	0.81	4.84	2.0	13.30	0.10	0.03
1E-26	0.010	3.5	87.3	16.8	7.0	13.27	0.19	0.22	1.0	13.39	0.08	0.69
1E-25	0.100	3.5	4043	67.9	3311	18.05	1.89	35.69	52.8	13.16	0.55	1.02
1E-25	0.037	3.5	1312	25.2	564.4	17.52	0.75	31.71	8.2	13.21	0.21	0.66
1E-25	0.010	3.5	122.1	17.2	48.7	18.28	0.68	37.43	1.2	13.28	0.08	0.16
1E-24	0.100	3.5	45449	863	23828	19.66	1.00	47.83	845	13.50	2.37	1.50
1E-24	0.037	3.5	12942	57.1	5979	19.19	0.87	44.28	41.9	13.21	0.49	0.66
1E-24	0.010	3.5	1201	18.2	568	17.05	0.93	28.17	2.9	13.26	0.13	0.30
1E-23	0.100	3.5	635501	9269	362386	24.40	0.81	83.49	9221	10.30	12.77	22.52
1E-23	0.037	3.5	193171	553	66253	24.48	0.47	84.08	529	13.82	1.47	3.94
1E-23	0.010	3.5	38419	50.1	8612	19.42	0.45	46.00	34.1	13.48	0.44	1.37

From Table 3 it is clear that ELL1k gives better reduced χ² values than ELL1 in all cases, indicating that red noise does not completely wash out the effects introduced by our modifications. However, when the red noise power is large (P_RN ≳ 10⁻²⁵ in Table 3), the systematic timing residuals such as those seen in Fig. 4 are no longer visually identifiable post-harmonic whitening (see Fig. 5). Therefore, the decision to employ the ELL1k model over ELL1 model should be made after estimating |$\tau _{\dot{\omega }}$| values rather than by visual inspection of timing residuals.

Figure 5.

Post-fit timing residuals after applying harmonic whitening to simulated TOAs injected with red noise having spectral parameters P_RN = 10⁻²⁵ s² yr⁻², f_c = 0.1 yr⁻¹, and ρ = 3.5 (Reduced χ² = 58.7). Total duration of our tempo2 simulations is 28 yr with one observation in every 45 d. Note the absence of systematic timing residuals such as those seen in Fig. 4.

It can be seen from Table 3 that the value of |$\dot{\omega }$| can be obtained with reasonable (⁠|${\sim } 1{{\ \rm per\ cent}}$|⁠) accuracy by applying harmonic whitening when the red noise amplitude is sufficiently small. It is still possible to obtain an estimate for |$\dot{\omega }$| even for higher red noise amplitudes albeit with poorer accuracies and higher reduced χ² values. Answer to the question how small an amplitude is ‘sufficiently small’ crucially depends on f_c as the accuracy of |$\dot{\omega }$| measurement reduces significantly for larger values of f_c. This is to be expected as harmonic whitening removes low-frequency noise more effectively than high-frequency noise.

A perhaps surprising observation is that applying harmonic whitening does not affect the accuracy of the |$\dot{\omega }$| measurement significantly. This is due to the fact that the timing residuals introduced by the exact secular evolution of LL parameters vary in two distinct time-scales (i.e. orbital and periastron advance) as we have noticed earlier in this section. Although the long-time-scale variation introduced by the ‘extra term’ is removed by FITWAVES, the contributions which arise from the general temporal evolutions of the LL parameters should allow us to obtain a |$\dot{\omega }$| estimate even in the presence of red noise. This point is further emphasized by the observation that the accuracy of |$\dot{\omega }$| measurement does not decrease even when the red noise frequency cut-off corresponds to the periastron advance time-scale itself (f_c = 0.037 yr|$^{-1} \sim \tau _{\dot{\omega }}^{-1}$|⁠).

5 SUMMARY AND DISCUSSIONS

We explored the possibility of measuring the effect of periastron advance in small-eccentricity binary pulsars. This was pursued by implementing a timing model, ELL1k, in tempo2 that essentially incorporated general temporal evolutions for the crucial LL parameters present in the ELL1 timing model. This present prescription should be relevant for binary pulsars that can, in principle, experience significant periastron advance. With the help of tempo2 simulations and using PSR J1719−1438 as an example, we probed the observational implications of the ELL1k timing model in comparison with the tempo2 implementation of Wex’s ELL1 model. This comparison provides significant timing residuals while simulating a decade-long timing of PSR J1719−1438 even while using the conservative |$\dot{\omega }_{\text{GR}}$| estimate for the periastron advance. Our tempo2 simulations suggest the possibility of measuring the effect of |$\dot{\omega }$| in PSR J1719−1438 by pursuing a high-cadence decade-long timing campaign. Additionally, we investigated the measurability of periastron advance in the presence of red timing noise in our tempo2 simulations for PSR J1719−1438. This is achieved by injecting several instances of red noise with varying spectral properties into the simulated TOAs. We concluded that |$\dot{\omega }$| can be measured to about 1 per cent accuracy provided the amplitude of the red noise is sufficiently small. Furthermore, the tolerable amplitude of red-noise increases with the red-noise time-scale (1/f_c). Interestingly, if the measured |$\dot{\omega }$| turns out to be larger than the expected |$\dot{\omega }_\text{GR} \sim$| 13° yr⁻¹, it should lead to constraints on the apsidal motion constant k₂ of its unique companion. We also compiled a list of binary pulsars with tiny orbital eccentricities where |$\dot{\omega }$| effects may become noticeable within a reasonable time-scale especially if the classical contributions to periastron advance dominates over their GR counterpart. More importantly, we showed that the linear-in-time evolutions of LL parameters do not introduce any significant timing residuals in the currently employed PTA pulsars with tiny orbital eccentricities.

The present effort should be interesting for the following scenarios. The first one involves the detection of sub-|$\mu$|Hz GWs using an ensemble of MSPs in a PTA experiment that demands post-fit timing residuals of the order of 10 nano-seconds. Unfortunately, this is not attainable for most of the current PTA MSPs and many of them are part of nearly circular binary systems. A possible source of higher post-fit timing residuals can be unmodelled systematic effects as evident from our plots. Therefore, the use of a more refined timing model can, in principle, lower post-fit residuals, especially for MSPs in compact orbits with small eccentricities. However, it is unlikely that the ELL1k model will lead to any improvements to the current PTA sensitivity and this is mainly due to the large |$\tau _{{\dot{\omega }}}$| values for the current list of PTA pulsars. However, it will be helpful to check the relevance of this timing model while including new MSPs into the existing list of PTA pulsars. The second scenario involves accreting millisecond X-ray pulsars like SAX J1808.4−3658. It should be interesting to explore the implications of this model while analysing timing data associated with its many observed outbursts. For such accretion-powered systems, the classical contributions to |$\dot{\omega }$| can be quite large and therefore it may be worthwhile to perform coherent timing of its outburst data by employing our timing model. This may, in principle, lead to an estimate for the k₂ value of its brown-dwarf companion and detailed analysis and simulations will be required to quantify our suggestion. Note that it will require LISA observations of eccentric Galactic binaries for estimating k₂ values of degenerate objects (Willems et al. 2008). Finally, this timing model should be relevant during the FAST-SKA era as the MSP population is expected to quadruple during this era (Levin et al. 2018). It is reasonable to expect that SKA will monitor dozens of short orbital period MSPs in nearly circular orbits and the ELL1k timing model will be required for the high-cadence timing of such systems.

ACKNOWLEDGEMENTS

We are grateful to Norbert Wex for his valuable comments and suggestions and thank the anonymous referee for her/his constructive comments which helped us to improve the manuscript. Additionally, we would like to express our gratitude to George Hobbs and Michael Keith for promptly incorporating our ELL1k timing model into the tempo2 pulsar timing package.

Footnotes

1

Available at https://sourceforge.net/p/tempo/tempo/ci/master/tree/src/bnryell1.f.

2

Available at https://bitbucket.org/psrsoft/tempo2.

3

Available at http://www.atnf.csiro.au/research/pulsar/psrcat/.

APPENDIX A: EXPRESSIONS REQUIRED TO IMPLEMENT ELL1k MODEL IN tempo2

We list below the partial derivatives of the Rømer delay, namely equation (6), with respect to the relevant binary parameters while neglecting terms of |$\mathcal {O}(e)$| and |$\mathcal {O}({\dot{P}}_b)$|⁠. These expressions are required to implement the ELL1k timing model in tempo2 and in Table A1 we provide a list of the binary parameters and their tempo2 notations.

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial P_\mathrm{ b}} = -\frac{x}{c}n_{\mathrm{ b}}\frac{\tau }{P_{\mathrm{ b}}}\cos \Phi \, , \end{eqnarray*}

(1a)

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial \dot{P}_\mathrm{ b}} = -\frac{1}{2}\frac{x}{c}n_{\mathrm{ b}}\frac{\tau ^{2}}{P_{\mathrm{ b}}}\cos \Phi \, , \end{eqnarray*}

(1b)

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial (x/c)} = \sin \Phi \, , \end{eqnarray*}

(1c)

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial (\dot{x}/c)} = \tau \sin \Phi \, , \end{eqnarray*}

(1d)

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial \epsilon _{10}} = -\frac{x}{2c}\left(\left(\cos 2\Phi +3\right)\cos \dot{\omega }\tau +\sin 2\Phi \,\,\sin \dot{\omega }\tau \right) \, , \end{eqnarray*}

(1e)

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial \epsilon _{20}} = -\frac{x}{2c}\left(\left(\cos 2\Phi +3\right)\sin \dot{\omega }\tau -\sin 2\Phi \,\,\cos \dot{\omega }\tau \right) \, , \end{eqnarray*}

(1f)

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial {\dot{\omega }}} = -\frac{x}{2c}\left(\epsilon _{2}\left(\cos 2\Phi +3\right)+\epsilon _{1}\sin 2\Phi \right)\tau \, , \end{eqnarray*}

(1g)

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial T_{\rm ascnode }} = -\frac{x}{c} n_{b} \cos \Phi \, . \end{eqnarray*}

(1h)

In addition, if onewishes to fit for the logarithmic derivative of eccentricity ξ, the corresponding partial derivative is given by

\begin{eqnarray*} \frac{\partial \Delta _\mathrm{ R}}{\partial \xi } = \frac{x}{2c}\left(\epsilon _{2}\sin 2\Phi -\epsilon _{1}\left(\cos 2\Phi +3\right)\right)\tau \, . \end{eqnarray*}

(1i)

Table A1.

Binary model parameters for ELL1k.

Parameter	tempo2 parameter name	Unit
Pb	PB	d
|$\dot{P}_{\mathrm{ b}}$|	PBDOT
x	A1	lt-s
|$\dot{x}$|	A1DOT	lt-s/s
ε10	EPS1
ε20	EPS2
|$\dot{\omega }$|	OMDOT	° yr−1
ξ		yr−1
|$T_{\rm ascnode }$|	TASC	MJD

Note that ξ is not present in the standard list of tempo2 parameters.

APPENDIX B: COMPARISON OF TIMING MODELS FOR LOW-ECCENTRICITY PULSAR BINARIES

With our proposed ELL1k model described in this work, there are now four different timing models for low-eccentricity binaries. We provide in Table B1 a brief comparison of these models.

Table B1.

Comparison of different timing models for low-eccentricity binary pulsars. Listed references provide further details on the various models and their parameters.

Binary Model	Model Parameters	Description
ELL1	Pb, |$\dot{P}_{\mathrm{ b}}$|⁠, x, |$\dot{x}$|⁠, ε10, ε20, |$\dot{\epsilon }_{1}$|⁠, |$\dot{\epsilon }_{2}$|⁠, |$T_{\rm ascnode }$|⁠, mc, sin ι	Model for low-eccentricity binaries. Uses Laplace–Lagrange parameters and |$T_{\rm ascnode }$| instead of usual Keplerian parameters e, ω, and T0. Incorporates |$\mathcal {O}(e)$| corrections of the Rømer and Shapiro delays. Reference : Lange et al. (2001)
ELL1H	Pb, |$\dot{P}_{\mathrm{ b}}$|⁠, x, |$\dot{x}$|⁠, ε10, ε20, |$\dot{\epsilon }_{1}$|⁠, |$\dot{\epsilon }_{2}$|⁠, |$T_{\rm ascnode }$|⁠, h3, h4, ς, Nharm	Model for low eccentricity, low-to-moderate inclination binaries with measurable Shapiro delay. Same as ELL1 except that this model uses a truncated Fourier expansion of Shapiro delay instead of the exact |$\mathcal {O}(e)$| expression, where hi is the amplitude of the ith harmonic, |$\varsigma =\frac{\sin \iota }{1+|\cos \iota |}$| and Nharm is the number of harmonics in the Fourier expansion. Reference : Freire & Wex (2010)
ELL1+	Pb, |$\dot{P}_{\mathrm{ b}}$|⁠, x, |$\dot{x}$|⁠, ε10, ε20, |$\dot{\epsilon }_{1}$|⁠, |$\dot{\epsilon }_{2}$|⁠, |$T_{\rm ascnode }$|⁠, mc, sin ι	Model for low-eccentricity, wide-orbit binaries where xe2/c is comparable to or larger than the timing precision. Incorporates |$\mathcal {O}(e^{2})$| terms in the Rømer delay. Reference : Zhu et al. (2018)
ELL1k	Pb, |$\dot{P}_{\mathrm{ b}}$|⁠, x, |$\dot{x}$|⁠, ε10, ε20, |$\dot{\omega }$|⁠, |$T_{\rm ascnode }$|⁠, mc, sin ι	Model for low-eccentricity, compact binaries showing significant advance of periastron. Incorporates an extra, slowly varying term in the Rømer delay that was neglected in ELL1. Reference : This work