A revisit of PSR J1909$-$3744 with 15-year high-precision timing

We combined the BON data presented in Desvignes et al. (2016) with the NUPPI data collected in the past eight years as described above, which delivered a new dataset with an overall timing baseline of nearly fifteen years. Whenever both BON and NUPPI data were available from the same observation, we kept only the NUPPI data. In total, the dataset contains 846 TOAs, 615 at L-band and 231 at S-band. This gave us an averaged observing cadence of more than one per week. We noted that in an earlier analysis using the NUPPI data of PSR J1909$-$3744 (Liu et al., 2014), the post-processing of the data was affected by an issue in the psrchive software package⁶⁶6This was induced by a lack of precision in polyco epoch stored in the header of PSRFITS files. The issue was fixed by an update of psrchive in late 2014.. Reprocessing the same dataset (which contains 30 epochs from MJD 56545 to 56592) now delivers timing residuals with a weighted rms of approximately 60 ns, more than a factor of 4 better than the previous value.

We used the temponest software package to perform timing analysis to our dataset. temponest is built based on tempo2 (Hobbs et al., 2006) and multinest (Feroz et al., 2009) software packages, to explore the parameter space of a non-linear pulsar timing model with Bayesian inference (Lentati et al., 2014). The fitted parameters included in the timing model are described in Table 1 (noted as “measured parameter”). To describe the orbital motion of the pulsar in the binary system, we used the “ELL1” timing model (Lange et al., 2001) as the orbit is known to have a very small eccentricity. Additionally, we included a set of parameters to characterize the noise in the dataset. The white noise parameters are the error factor ‘EFAC’ ($E_{\rm f}$) and an additional error added in quadrature ‘EQUAD’ ($E_{\rm q}$) which relate to the measurement uncertainty ($\sigma_{\rm r}$) of a TOA as:

For each observing system a pair of such parameters was included to capture its white-noise properties to yield a proper weight in the data combination. To describe the long-term red process in the dataset, we included two models to account for the chromatic component caused by dispersion measure (DM) variation and the monochromatic part of the signal, respectively. In both models, the signal is assumed to be a stationary, stochastic process with a power-law spectrum in the form of

where $S(f)$, $A$, $\gamma$ $f_{\rm r}$ are the power spectral density as a function of frequency $f$, the spectral amplitude, the spectral index and the reference frequency (set to 1 yr^-1), respectively. The constant $C_{\rm 0}$ is equal to $1$ for the chromatic term and is $12\pi^{2}$ for the monochromatic term. The spectrum had a low-frequency cutoff equal to the inverse of the data span (approximately 14.6 yr), and was sampled with integer multiples of the lowest frequency up to one over fourteen days. To align the data from different systems, we included an arbitrary time offset (known as JUMP) between the reference system and each of the rest systems. Both the timing and noise parameters were simultaneously fitted for during the temponest analysis, while the JUMPs were analytically marginalised. All model parameters were sampled with uniform priors, except for EQUAD and amplitudes of the red processes for which a log-uniform prior was used.

In Fig. 1 we present the timing residuals of the dataset obtained from the analysis with temponest, both before and after the subtraction of the red processes in the data. A long-term evolution of the residuals is prominent across the entire time span of the dataset, leading to an weighted rms of over 500 ns. We note that a similar feature in the PSR J1909$-$3744 timing residuals was also seen in the recent analysis presented by Arzoumanian et al. (2018a). Subtraction of the measured red noise models successfully whitened the residuals, leading to a weighted residual rms of 103 ns over a timescale of nearly fifteen years. The L-band data from NUPPI itself has a weighted residual rms of 86 ns. The posterior distribution of the red-noise parameters are presented in Fig. 2, where the logarithmic amplitude and spectral index of the achromatic red noise were measured to be $\log(A_{\rm red})=-13.92^{+0.07}_{-0.07}$ and $\gamma_{\rm red}=2.51^{+0.30}_{-0.27}$, highly consistent with published values by Caballero et al. (2016) and Arzoumanian et al. (2018a).

Timing precision of the brightest MSPs could be limited by the intrinsic variability in the pulsar signal itself which is known as “jitter noise” (e.g., Cordes & Downs, 1985; Liu et al., 2011). To estimate such contribution in our data, we selected three epoch observations when the pulsar signal was the strongest due to interstellar scintillation, and measured the jitter noise following the approach described in Liu et al. (2012). For L-band, this gave us $\sigma_{\rm J,30min}=14.1\pm 1.7$ ns, where $\sigma_{\rm J,30min}$ is the amount of jitter noise with 30-min integration time. This result is well in agreement with previously published values (Shannon et al., 2014b; Lam et al., 2016, 2019). We did not measure any significant jitter noise in our S-band observation, but obtained a 95% upper limit of 78 ns for $\sigma_{\rm J,30min}$ at 2.5 GHz. The limit is in line with the prediction in Lam et al. (2019) that jitter noise of PSR J1909$-$3744 should be around 10 ns at this frequency. To explore the impact of jitter noise on our noise modelling, we incorporated our measurement into the L-band data, i.e., quadratically adding the corresponding jitter noise to the TOA errors before performing the timing analysis with temponest. By doing so, the white noise description here becomes $\sigma=\sqrt{E^{2}_{\rm q}+E^{2}_{\rm f}[\sigma^{2}_{\rm r}+\sigma^{2}_{\rm J}(\tau)]}$, where $\sigma_{\rm J}$ is calculated based on the integration time ($\tau$) of the observation as $\sigma_{\rm J}=\sigma_{\rm J,30min}/\sqrt{\tau/{\rm 30\,min}}$. In principle, the inclusion of jitter noise in the noise model would be fully covered by EQUAD if the integration times of the observations were identical (i.e., it is then a constant addition to all TOA errors just as EQUAD). This, however, is not the case here because the integration times of our observations range from 20 to 80 mins, and thus the jitter noise corresponding to each TOA can differ by up to a factor of two. Using this modified L-band dataset instead in the timing analysis gave a very minimal (well below 10%) change in EQUAD values of the L-band systems, while the EFAC values remained around unity. This is expected as the jitter noise with a typical 30-min integration time is well below the timing residual rms of the L-band data. In fact, jitter noise is dominant over the radiometer noise only in less than 1% of our observations.

To further verify our results, we carried out a separate investigation of the noise properties using the Enterprise software package which is also widely used in pulsar timing array data analysis (Ellis et al., 2019). We built the noise model consisting of the timing model, white noise and both the chromatic and achromatic red noise components. Those are the same components used in temponest, while the implementation and application are somewhat different. In detail, we performed a fully Bayesian analysis using a parallel tempering Markov Chain Monte Carlo sampler PTMCMC (Ellis & van Haasteren, 2017), adopting the same priors as in the temponest analysis but marginalising over the timing parameter errors (van Haasteren & Levin, 2012). The white noise is characterised by ‘EFAC’ and ‘EQUAD’, with the same description as in temponest (cf. Eq. 1). Both red-noise components were modelled as a stationary Gaussian process with a power-law power spectral density parameterised by an amplitude (at frequency $1/\mathrm{year}$) and a spectral index, the same as in Eq. 2 except for that here the constant $C_{\rm 0}$ is equal to $12\pi^{2}$ for both chromatic and monochromatic terms. The inferred posteriors, as shown in Fig. 2, are largely consistent with the results from temponest, with somewhat more pronounced tail in the distribution for the achromatic red noise distribution. This difference might be attributed to the use of different samplers (MULTINEST vs PTMCMC) in the multidimensional parameter space.

In addition, we noted that the irregularly distributed and unevenly numbered observations between L-band and S-band in the dataset may have an impact on the measurement of DM variation, in particular on separating its contribution to the overall red process from the achromatic component. To examine this potential issue and its potential impact on our timing results, we carried out the same timing analysis to a different version of the dataset where the sub-band TOAs were used for L-band observations with NUPPI. This offers more regular multi-frequency coverage and as shown in Fig. 2, thus provides better constraints on DM model parameters. The achromatic red noise parameters becomes less constrained, probably due to the drop of the best-estimated value of the amplitude. The spectral index is shifted to a higher value as expected since the effect from an “imperfect” DM modelling exhibits as a relatively shallow noise source. On the other hand, applying the noise model measurements obtained with the sub-band TOAs to subtract the red process in the original dataset resulted in very little difference in the residuals and the overall weighted rms. This suggests that our noise modelling with the broad-band TOAs should still be able to effectively extract the red processes in the data.

It is worth noting that the experimental analysis on incorporating jitter noise measurement in the TOAs and using sub-band TOAs in the dataset have both yielded highly consistent measurement in all timing parameters, compared with those from using the broad-band TOAs. This is demonstrated in Fig. 3, where the differences in measured timing parameters are shown to be all consistent with zero within 1-$\sigma$ confidence interval. We note that in Alam et al. (2020a) and Alam et al. (2020b) the analysis using broad-band and sub-band TOAs has also led to highly consistent timing results. Thus, for the analysis in the rest of the paper we will simply quote the timing results obtained from the broad-band dataset. A more detailed noise analysis using different forms of datasets will be presented in an upcoming paper (Chalumeau et al. in prep.).

The results of the timing parameters from the analysis are summarised in Table 1, where the uncertainties of the measurements are 1-$\sigma$ Bayesian credible interval obtained from one-dimensional marginalised posterior distribution. The posterior distributions of a subset of the timing and noise parameters from the temponest analysis can be found in Fig. 4. In Table 2 we compared a subset of timing parameters from our analysis with a few recent publications. It can be seen that all of our measurements have a broad consistency with previous analysis, and that our work has led to a clear improvement in measurement precision of $\dot{P}_{\rm b}$ and $\dot{x}$. We note that Alam et al. (2020a) reported a covariance between $\dot{x}$ and the red-noise terms. This, however, is not obvious from the posteriors shown in Fig. 4, possibly because the longer timing baseline here helps to mitigate the degeneracy between $\dot{x}$ and the red-noise parameters. Using the measurements of the Shapiro delay parameters, $s$ and $r$, the mass of the pulsar and the WD are derived to be $m_{\rm p}=1.492\pm 0.014$ $M_{\odot}$ and $m_{\rm c}=0.209\pm 0.001$ $M_{\odot}$, respectively, as demonstrated in Fig. 5. These new values resolves the slight inconsistency reported in Desvignes et al. (2016) from other works. Fig. 5 also shows the independent constraint on the mass ratio ($q\equiv m_{\rm p}/m_{\rm c}=7.0\pm 0.5$) reported in Antoniadis (2013) based on optical observations of the WD, which is consistent with the values determined from Shapiro delay, albeit with a larger uncertainty.

The orbital eccentricity can be calculated as $e=\sqrt{\kappa^{2}+\eta^{2}}=(1.15\pm 0.07)\times 10^{-7}$, which still makes PSR J1909$-$3744  the most circular system in all binary pulsars known by far (Manchester et al., 2005). Using the mass measurements, the relativistic periastron advance of the system is calculated to be 0.14 deg/yr (Lorimer & Kramer, 2005). Accordingly, the eccentricity vector has rotated by only 2 deg within the fifteen years of our timing baseline, which, if unmodelled, would cause an apparent variation in the $\hat{x}$, $\hat{y}$ components of the eccentricity vector by less than $10^{-10}$. This is well below their measurement precision shown in Table 1. Thus, the periastron advance is still not measurable with our dataset.

The observed secular variation of orbital period in binary pulsars can include contributions from multiple effects (Damour & Taylor, 1991; Lorimer & Kramer, 2005). As for the case of PSR J1909$-$3744 where the pulsar has negligible tidal interaction with the WD companion, the most relevant effects can be summarised below:

The first term, $\dot{P}_{\rm b}^{\rm Shk}$, denotes the contribution from the transverse relative motion of the system to the solar system barycentre (SSB) which is known as the Shklovskii effect. Using the measured proper motion of the pulsar ($\mu_{\rm\alpha}$, $\mu_{\rm\delta}$) and its distance ($d$) derived from timing parallax, this term is estimated to be

where $c$ is the speed of light. The second term, $\dot{P}_{\rm b}^{\rm Gal}$, is caused by the differential acceleration from the Galactic potential. As shown in Damour & Taylor (1991), this effect can be calculated based on astrometric measurements and a Galactic potential model. Using our timing results and the model provided in McMillan (2017), we have⁷⁷7The error in this calculation only reflects the error in the parallax. Accounting for a realistic error in the Galactic potential, one expects an error roughly an order of magnitude larger, which however, is still about an order of magnitude smaller than the error in $\dot{P}_{\rm b}^{\rm Shk}$.

The third term, $\dot{P}_{\rm b}^{\rm GR}$, represents the contribution from gravitational wave damping in GR, and following Peters (1964), is found to be⁸⁸8Whenever being in combination with $T_{\odot}$, masses are understood in solar units.

where $T_{\odot}=(\mathcal{GM})_{\odot}^{\mathrm{N}}/c^{3}$, and $(\mathcal{GM})_{\odot}^{\mathrm{N}}$ is the solar mass parameter as defined by the IAU (Mamajek et al., 2015). The last term, $\dot{P}^{\dot{m}}_{\rm b}$, is arising from the mass loss of the system. In the case of PSR J1909$-$3744, we consider the loss of rotational energy of the pulsar to be the dominant mass-loss effect, given that there is no expectation of significant stellar wind from a cooled WD with strong surface gravity like the companion of PSR J1909$-$3744 (e.g., Unglaub, 2008). Following Damour & Taylor (1991), this can be worked out as⁹⁹9The uncertainty here only includes the measurement error of the parameters, except $I_{\rm p}$ whose value is assumed. It should be noted that Eq. (7) is an approximation for slowly rotating neutron stars, whose precision nevertheless is better than 10% for PSR J1909$-$3744.

Here $I_{\rm p}$ and $\dot{P}$ denote the moment of inertia and the intrinsic period derivative of the pulsar, respectively. To obtain the estimate above, we used $I_{\rm p}=1.4\times 10^{45}$ g cm², a typical value derived using our mass measurement and different equations of state that are in agreement with all pulsar and LIGO/Virgo constraints (Lattimer & Prakash, 2001; Capano et al., 2020). We also recovered the intrinsic spin period derivative, by subtracting the contribution from the Shkloskii effect and the Galactic potential from the observed value as performed in the case of the orbital period derivative. These calculations are summarized in Table 1 which shows that the intrinsic $\dot{P}$ of the pulsar is in fact more than a factor of five smaller than the observed value.

The analysis above clearly shows that the observed $\dot{P}_{\rm b}$ is dominated by the Shklovskii effect, next to which is the contribution from Galactic potential. Though the expected contribution from gravitational wave damping is well above the measurement uncertainty of the observed $\dot{P}_{\rm b}$, it is still below the estimation error of the Shklovskii contribution which mostly comes from the measurement uncertainty in parallax distance. Thus, with the current measurement precision $\dot{P}_{\rm b}^{\rm GW}$ is not separable from $\dot{P}^{\rm obs}_{\rm b}$ and we expect to be able to do so only when the precision in distance measurement is improved by at least a factor of 4 to 5. This could possibly be obtained by adding more existing data, extending the timing baseline, or using timing data from more sensitivity telescope (e.g., Bailes et al., 2018). Subtracting all aforementioned contributions from $\dot{P}^{\rm obs}_{\rm b}$ gives $\dot{P}_{\rm b}^{\rm Exs}=(-1.7\pm 7.8)\times 10^{-15}$ which is consistent with zero. This allows us to calculate the kinematic distance ($d_{\rm k}$) of the system with $\dot{P}^{\rm obs}_{\rm b}$ by directly assuming a zero excessive $\dot{P}_{\rm b}$ (Bell & Bailes, 1996), which gives $d_{\rm k}=1.157(3)$ kpc, a highly consistent and more precise distance estimate.

The observed secular variation of projected semi-major axis in binary pulsars can also be induced by a series of effects as listed below (Lorimer & Kramer, 2005):

The first term, $\dot{x}^{\rm Shk}$, stands for the contribution from the Shklovskii effect, and using our timing astrometric measurements, is estimated to be:

The second term, $\dot{x}^{\rm Gal}$, is the contribution from the differential Galactic potential, and following the same approach as for $\dot{P}_{\rm b}^{\rm Gal}$, is found to be

The third term, $\dot{x}^{\rm GW}$, is caused by gravitational wave damping. Following (Peters, 1964) and using our mass measurements, it is estimated to be

The fourth and fifth term, $d\epsilon_{\rm A}/dt$¹⁰¹⁰10Here $\epsilon$ is defined as in Eq. (2.6) of Damour & Taylor (1992). and $\dot{x}^{\rm SO}$, denote contribution from geodetic precession of the pulsar spin axis and the precession of the orbital plane due to spin-orbit coupling effect (classical and/or Lense-Thirring), respectively. Following Damour & Taylor (1991), $d\epsilon_{\rm A}/dt$ is of order $\Omega^{\rm geo}\cdot P/P_{\rm b}\sim 10^{-19}$. The spin-orbit precession induced by the companion’s rotation has so far been measured in a small number of different binary pulsar systems (Kaspi et al., 1996; Shannon et al., 2014a; Venkatraman Krishnan et al., 2020). However, spin-orbit contributions to $\dot{x}$ vanish when the spin axes of the pulsar and companion are aligned with the orbital angular momentum (e.g., Damour & Taylor, 1992; Wex, 1998), which is expected to be the case for fully recycled binary pulsars such as PSR J1909$-$3744 from the binary evolution aspect (e.g., Tauris, 2012). Additionally, the pulsar is in a relatively wide orbit with the WD companion, which will further suppress the spin-orbit coupling effects by a few orders of magnitude compared with the case of PSR J1141$-$6545 studied in (Venkatraman Krishnan et al., 2020).

Therefore, it can be concluded that the contribution from the relative motion between the pulsar and the SSB, $\dot{x}^{\rm PM}$, should be dominating in $\dot{x}^{\rm obs}$. Following the derivations presented in Arzoumanian et al. (1996) and Kopeikin (1996), this effect can be expressed as¹¹¹¹11Note that there is a conversion between the $\Omega$ and $i$ here and those in Kopeikin (1996): $\Omega=\pi/2-\Omega_{\rm K96}$, $i=\pi-i_{\rm K96}$.

It can be seen that one can derive the ascending node of the binary system, $\Omega$, once $\dot{x}^{\rm PM}$ and $i$ are measured. This is shown in Fig. 6, where we plot the contribution to $\dot{x}/x$ from the proper motion for PSR J1909$-$3744 as a function of $\Omega$, in conjunction with the measured $\dot{x}/x$. Accordingly, we have identified four possible solutions considering the ambiguity in $i$ ($\sin{i}$ equals to the shape of Shapiro delay, $s$): $\Omega=217\pm 5$, $352\pm 5$ deg when $i<90$ deg, and $\Omega=37\pm 5$, $172\pm 5$ deg when $i>90$ deg. We noted that Reardon et al. (2016) reported only one possible solution of such and Alam et al. (2020a) reported two. In principle, this ambiguity can be resolved if the annual-orbital parallax of the binary system is measurable, via directly fitting for the “Kopeikin” (KOM and KIN) parameters in the tempo2 software package as e.g., has been carried out in the PSR J1713+0747 system (see e.g. Desvignes et al., 2016). However, with the formula in Kopeikin (1995), we estimate the scale of the feature in the timing data induced by annual-orbital parallax to be $\sim 6$ ns, well below any of the published timing precision of PSR J1909$-$3744. This is mostly due to the fact that PSR J1909$-$3744 has a small and edge-on orbit with a distance on order of kpc. In practice, the Kopeikin parameters can also be highly correlated with some of the Keplerian parameters (e.g., the orbital period), which would further suppress the subtractable feature from the fit. Thus, our current timing precision does not allow to distinguish among the four possible solutions by modelling the annual-orbital parallax of the system. This has been verified by our attempt to switch on the fit for the Kopeikin parameters in our analysis, where all of the four solutions can be found, resulting in the same post-fit residual rms and returning the same goodness-of-fit. We however note that studying the orbital variations in the interstellar scintillation of the pulsar could potentially provide additional constraint over orbital inclination and help to resolve the ambiguity in our solutions (e.g., Lyne, 1984; Reardon et al., 2019).

The orbital eccentricity of PSR J1909$-$3744  is a record low among all known binary pulsars with a value of just $e=1.15\times 10^{-7}$. Given a potential age of from 0.5 to 7 Gyr of this binary pulsar (inferred from the WD cooling age, as discussed in Section 3.2), it is of interest to investigate whether we can use it as a probe of the stellar density in the Galactic disk. It is anticipated that orbital eccentricities might be induced by dynamical interactions with field stars, similar to the examples of resulting high eccentricities seen among some binary pulsars in dense stellar environments like globular clusters. The question is if for an assumed old system like PSR J1909$-$3744  (depending on the exact WD cooling models discussed in Section 3.2), the column density of field stars within its cross-section radius, accumulated from its motion in the disk (Section 3.3) over its lifetime of several Gyr, may reach a critical level.

Rasio & Heggie (1995) investigated resulting eccentricities of binary MSPs induced via such dynamical interactions, and they derived the following expression (their Eq. 6):

where $\eta\equiv t_{9}\,n_{4}/v_{10}$, and $t_{9}$ is the age in units of Gyr, $n_{4}$ is the stellar density (all assumed to be of mass $1\;M_{\odot}$) in units of $10^{4}\;{\rm pc}^{-3}$, and $v_{10}$ being the one-dimensional velocity dispersion in units of $10\;{\rm km\,s}^{-1}$. Given $P_{\rm b}=1.53\;{\rm days}$ and $e=1.15\times 10^{-7}$, it yields $\eta\simeq 0.50$. Assuming $v_{10}\simeq 5$ (based on the relative velocity between PSR J1909$-$3744  and local field stars), and $t_{9}\simeq 10$ as an extreme upper limit, we find that $n_{4}\la 0.25$. I.e. the critical number density needed to explain the small observed eccentricity of $1.15\times 10^{-7}$ for PSR J1909$-$3744  based on dynamical interactions is about $2\,500\;{\rm stars\;pc}^{-3}$. This value is four orders of magnitude larger that the number density of stars in the solar neighborhood, $n_{\odot}=0.12\;{\rm pc}^{-3}$ (e.g., Holmberg & Flynn, 2000).

Therefore, we conclude that it is not surprising that PSR J1909$-$3744  was able to retain its small eccentricity although potentially plowing through the Galactic disk for possibly up to several Gyr. We also note that this eccentricity is expected to be a residual value from its formation process via mass transfer from the progenitor star of the WD (Phinney, 1992).

Short orbital period ($P_{\rm b}\lesssim 1$ d) pulsar-WD systems have turned out to provide some of the most stringent limits on dipolar gravitational radiation, a prediction by many alternative gravity theories (Wex, 2014). With its orbital period of approximately 1.5 days, PSR J1909$-$3744 is rather slow compared to the systems like PSR J1738+0333 (Freire et al., 2012) or PSR J0348+0432 (Antoniadis et al., 2013). However, its outstanding timing precision and the corresponding precise measurement of $\dot{P}_{\rm b}$ makes it still interesting for a test of dipolar gravitational wave damping. Moreover, while for PSR J1738+0333 and PSR J0348+0432 optical measurements and the modelling of WD spectra were required to obtain the masses of pulsar and companion, for PSR J1909$-$3744 we can estimate these masses directly from the high precision timing observations due to the presence of a prominent Shapiro delay.

Here we discuss our dipolar radiation test within the framework of Bergmann-Wagoner theories. Bergmann-Wagoner theories represent the most general mono-scalar–tensor theories that are at most quadratic in the derivatives of the fields in their action (Will, 2018). Quite a number of well known mono-scalar-tensor theories belong to this class, like Jordan-Fierz-Brans-Dicke (JFBD) gravity (Jordan, 1955; Fierz, 1956; Brans & Dicke, 1961), DEF gravity (Damour & Esposito-Farese, 1993, 1996), MO gravity (Mendes & Ortiz, 2016), $f(R)$ gravity (De Felice & Tsujikawa, 2010), and massive Brans-Dicke gravity (Alsing et al., 2012).

In Bergmann-Wagoner theories, the field equations for the (physical) “Jordan-frame metric” $g_{\mu\nu}$ and the scalar field $\phi$ can be derived from the following action

where $G_{0}$ denotes the fundamental (“bare”) gravitational constant, $g$ is the determinant of the metric $g_{\mu\nu}$, and $R$ the curvature scalar. The two functions $\omega(\phi)$ and $U(\phi)$ denote the coupling function and the scalar potential, respectively. Here we will assume that the influence of $U(\phi)$ is negligible on the relevant scales (cf. Alsing et al., 2012). $S_{\rm mat}$ is the action of the matter fields $\psi_{\rm mat}^{A}$, which couple universally to the spacetime metric $g_{\mu\nu}$.

The Newtonian gravitational constant, as measured in a Cavendish-type experiment, is given by

where $\phi_{0}$ is the cosmological background field and $\zeta\equiv 1/(2\omega(\phi_{0})+4)$. A further quantity, which we will need below, is the sensitivity

which accounts for the dependence of the mass of the $i$-th body (for a fixed number of baryons) on a change in its ambient scalar field (Will, 2018).

For our dipolar radiation test with PSR J1909$-$3744 we can utilize three post-Keplerian parameters. These are the range

and shape

of the Shapiro delay caused by the WD companion, and a potential (intrinsic) change of the orbital period due to scalar dipolar radiation

where $\tilde{m}_{i}\equiv G_{\rm N}m_{i}/c^{3}$. In Eq. (22) we have neglected the dependence on the eccentricity, since the PSR J1909$-$3744 orbit is practically circular. The sensitivity of the WD companion has been neglected as well, since generally $s_{\rm c}\ll s_{\rm p}$ (see e.g. Will, 2018, for details). In principle, $\dot{P}_{\rm b}$ also has additional multipole contributions, dominated by the tensorial quadrupole (Damour & Esposito-Farese, 1992a; Will, 2018). They are all well below the current measurement precision (cf. Eqs. 6) and (23)) and therefore can be ignored in the following calculations.

While $|\zeta|\lesssim 10^{-5}$ due to Solar system experiments (Bertotti et al., 2003), one cannot a priori assume that $\zeta s_{\rm p}$ and $\zeta s_{\rm p}^{2}$ are small as well. In fact, within Damour–Esposito-Farèse gravity it has been found that $\zeta s_{\rm p}^{2}$ can be of order unity even if $\zeta\rightarrow 0$. In this regime of the so-called “spontaneous scalarization” the strong-field gravity of a neutron star can lead to very large sensitivities, $s_{\rm p}$.

For a given theory, i.e. a given $\omega(\phi)$, and a given equation of state for neutron star matter, one can calculate the sensitivity $s_{\rm p}$, and use the three post-Keplerian parameters, i.e. Eqs. (20–22), with the two unknown masses as a consistency test (see e.g. Damour & Taylor, 1992). Here we will follow a more generic approach by leaving $\omega(\phi)$ unspecified and only making use of the Solar system constraint on $|\zeta|$. For a given $\zeta$, Eqs. (20–22) can then be solved for $\tilde{m}_{\rm p}$, $\tilde{m}_{\rm c}$, and $s_{\rm p}$.

The numerical values for the Shapiro range in Eq. (20) and the Shapiro shape (21) can directly be taken from Table 1: $r=1.029\pm 0.005\,\mu{\rm s}$ and $s=\sin i=0.998005\pm 0.000065$. In Eq. (22) we need the intrinsic change of the orbital period, i.e. we need to subtract Shklovskii ($\dot{P}_{\rm b}^{\rm Shk}$) and Galactic ($\dot{P}_{\rm b}^{\rm Gal}$) contributions from the observed orbital period derivative $\dot{P}_{\rm b}$ in Table 1 (see e.g. Lorimer & Kramer, 2005). Using the parallax distance $d_{\pi}$ from Table 1 and the Galactic potential of McMillan (2017) we find for the intrinsic change of the orbital period

It is worth mentioning that as discussed in Section 3.1.2, the observed secular change in orbital period is predominantly coming from the Shklovskii contribution. There is no measurable intrinsic orbital period change, as one can see from Eq. (23), and the Galactic correction is only $2.7\pm 0.1\,{\rm fs\,s^{-1}}$ and therefore smaller than the error in the Shklovskii correction. For comparison, with the masses from Table 1 one finds for the orbital period change as predicted by the GR’s quadrupole formula $\dot{P}_{\rm b}^{\rm GR}=-2.79\pm 0.03\,{\rm fs\,s^{-1}}$ (see Eq. (6)), significantly smaller than the error in $\dot{P}_{\rm b}^{\rm int}$.

With the numbers for the three post-Keplerian parameters in hand, we can now use Eqs. (20–22) to calculate the mass parameters $\tilde{m}_{\rm p}$ and $\tilde{m}_{\rm c}$ and put generic (i.e. independent of the details of $\omega(\phi)$) constraints on combinations of $\zeta$ and $s_{\rm p}$. As it turns out, for the allowed range of $|\zeta|\lesssim 10^{-5}$ the constraints on $s_{\rm p}$ do come exclusively from Eq. (22). As a result, for all values of $\zeta$ compatible with Solar system experiments one has

In a comparison with the currently best limit on dipolar radiation obtained from PSR J1738+0333 (Freire et al., 2012; Zhu et al., 2019) one has to keep in mind that $\zeta=\kappa_{\rm D}/4$ (Will, 2018). Assuming Eq. (7) in Zhu et al. (2019) for the sensitivity, in order to make the limit comparable, one obtains from Eq. (24)

This limit is still about an order of magnitude weaker than the one from PSR J1738+0333. Nevertheless, the limit is qualitatively different as it is based on three post-Keplerian parameters, while the PSR J1738+0333 limit used optical observations and WD models to obtain the masses of pulsar and companion. As a result, the mass of PSR J1909$-$3744 is known with higher precision (better than 1%; see Tab. 1), which is important for strong field effects that depend critically on the neutron-star mass. In Fig. 12 we give a specific example to illustrate this.

High-precision timing with binary pulsars is a also powerful tool in constraining the gravitational Lorentz invariance violation via its orbital dynamics. This type of study is typically carried out using the parameterized post-Newtonian (PPN) formalism (Will, 2014; Wex, 2014; Shao & Wex, 2016; Will, 2018). Here we investigate the application of the PSR J1909$-$3744 system to such gravity experiments.

In the PPN framework, Lorentz invariance violation is described by three PPN parameters, $\alpha_{1}$, $\alpha_{2}$, and $\alpha_{3}$. While $\alpha_{1}$ and $\alpha_{2}$ only pertain to semi-conservative dynamics, $\alpha_{3}$ additionally introduces nonconservative effects, namely, it simultaneously introduces preferred-frame effects and violates the energy-momentum conservation laws (Will, 2014, 2018). The $\alpha_{3}$ parameter was severely bound, down to the level of $\sim 10^{-20}$, via the long-term high-precision timing of PSR J1713+0747 (Zhu et al., 2019). Therefore, here we only consider constraints of $\alpha_{1}$ and $\alpha_{2}$. Because neutron stars are strongly self-gravitating objects, in order to distinguish them from weak-field objects, in the following we use $\hat{\alpha}_{1}$ and $\hat{\alpha}_{2}$ respectively to denote the strong-field counterparts of $\alpha_{1}$ and $\alpha_{2}$.

As shown in Shao & Wex (2012), for a nearly circular binary orbit, $\hat{\alpha}_{2}$ introduces a precession of the orbital angular momentum around the direction of ${\bf w}$, where ${\bf w}$ is the absolute velocity of the binary with respect to a preferred frame. Usually, the frame wherein the cosmic microwave background (CMB) is isotropic is chosen. The $\hat{\alpha}_{2}$-induced precession causes the change of the orientation of the orbit with respect to the observer, notably that the orbital inclination angle $i$ is changing. This change introduces a nonzero time derivative of the projected semi-major axis, and can be bound via the timing parameter $\dot{x}$. The contribution is given by (Shao & Wex, 2012),

where $n_{b}\equiv 2\pi/P_{b}$, $\psi$ is the angle between ${\bf w}$ and the orbital norm $\hat{\bf k}$, and $\vartheta$ is the angle between $\mathbf{w}_{\perp}\equiv{\bf w}-\left({\bf w}\cdot\hat{\bf k}\right)\hat{\bf k}$ and the direction of ascending node. As discussed in Section 3.1, for PSR J1909$-$3744 the proper motion effect is able to fully account for the observed $\dot{x}$. Hence, there is no need to have a nonzero $\hat{\alpha}_{2}$ in order to reproduce the measured $\dot{x}$, which is in fact translated into a bound on $\hat{\alpha}_{2}$.

We follow the method developed in Shao & Wex (2012) to obtain the upper limit on $\hat{\alpha}_{2}$. We randomize orbital ascending node $\Omega$ uniformly in $\left[0,360^{\circ}\right)$, and give equal probabilities to the configuration $i\in\left[0,90^{\circ}\right)$ and the configuration $i\in\left[90,180^{\circ}\right]$. We subtract the contribution to $\dot{x}/x$ from the proper motion (12) and assign the remaining $\dot{x}/x$ to the $\hat{\alpha}_{2}$-induced precession (26). Such a treatment renders the test a probabilistic one (see Sec. 3 in Shao & Wex, 2012), thus the obtained results on $\hat{\alpha}_{2}$ should only be interpreted as upper bounds (see Fig. 13). In the test, the CMB frame is chosen to be the preferred frame, and the absolute 3-dimensional systematic velocity is obtained for the binary by combining radio timing and optical observations (Antoniadis, 2013). All uncertainties in relevant parameters are taken into account in the Monte Carlo simulations. The probability density for $\hat{\alpha}_{2}$ is given in Fig. 13, from where we have,

These limits are weaker than the best existing limit from the spin precession of the Sun (Nordtvedt, 1987) and solitary pulsars (Shao et al., 2013), yet it still represents a new modest limit from a completely different system.

As for the PPN parameter $\hat{\alpha}_{1}$, Damour & Esposito-Farèse (1992b) elegantly picturised the orbital polarization phenomenon for near-circular binary orbits. In this picture, a nonzero $\hat{\alpha}_{1}$ induces a polarization of the eccentricity vector, ${\bf e}\equiv e\hat{\bf a}$ where $\hat{\bf a}$ is the direction to the periastron. The polarization is towards a direction in the orbital plane that is perpendicular to the absolute velocity. Consequently, the observed eccentricity vector is a vectorial superposition of a constant-length, relativistically rotating eccentricity, ${\bf e}_{R}(t)$, and a forced eccentricity, ${\bf e}_{F}$,

The forced eccentricity is given by (Damour & Esposito-Farèse, 1992b; Shao & Wex, 2012),

where $q\equiv m_{p}/m_{c}$, ${\cal V}_{O}\equiv\left({\cal G}Mn_{b}\right)^{1/3}$ with ${\cal G}$ the effective gravitational constant and $M\equiv m_{p}+m_{c}$. The ${\dot{\omega}_{\mathrm{PN}}}$ in Eq. (30) is the constantly rotating rate of ${\bf e}_{R}(t)$, which, in GR, is the well-known periastron advance rate.

The above picture was applied in Shao & Wex (2012) to PSRs J1012+5307 and J1738+0333, and the yet tightest limit on $\hat{\alpha}_{1}$ was achieved. The test used the observed eccentricity, decomposed into the Laplace parameters $\eta\equiv e\sin\omega$ and $\kappa\equiv e\cos\omega$ (Lange et al., 2001). It requires that the $\hat{\alpha}_{1}$-introduced time variations in $\eta$ and $\kappa$ to be consistent with their observed uncertainties (Shao & Wex, 2012). The best constraint, $\hat{\alpha}_{1}=-0.4_{-3.1}^{+3.7}\times 10^{-5}$ at the 95% confidence level, comes from PSR J1738+0333 (Freire et al., 2012).

To obtain constraint on $\hat{\alpha}_{1}$ with PSR J1909$-$3744, we repeated the same timing analysis as described in Section 3.1.1, with time derivatives of $\eta$ and $\kappa$ ($\dot{\eta}$ and $\dot{\kappa}$) directly included as fitted timing parameters. This allowed us to have high-precision measurements of both $\eta$ and $\kappa$ and their time derivatives, and enabled us to perform a more straightforward test of $\hat{\alpha}_{1}$ by directly using the time derivative of ${\bf e}(t)$. With a nonzero $\hat{\alpha}_{1}$, after averaging over an orbital timescale, we have (Shao & Wex, 2012)

where $\hat{\bf b}\equiv\hat{\bf k}\times\hat{\bf a}$, and $\mathbf{w}_{\perp}$ is the projection of ${\bf w}$ onto the orbital plane.

We set up Markov-chain Monte Carlo simulations which account for all the uncertainties. For each value of $\Omega$, we use the emcee package (Foreman-Mackey et al., 2013) to explore the posterior distribution of $\hat{\alpha}_{1}$ that is consistent with both $\dot{\eta}$ and $\dot{\kappa}$. From the posterior, we obtain the median, as well as the 1-$\sigma$, 2-$\sigma$, and 3-$\sigma$ enclosed regions for $\hat{\alpha}_{1}$. The result for the configuration $i<90^{\circ}$ is given in Fig. 14, while the other configuration with $i\geq 90^{\circ}$ is merely a mirrored case of Fig. 14 (see e.g. Figs. 6 & 7 in Shao & Wex, 2012). From the figure, we see that the loosest limit is from $\Omega\sim 135^{\circ}$,

The limit is slightly worse than that from PSR J1738+0333. If we assume that the $\hat{\alpha}_{2}$ limit from the spin precession of solitary pulsars (Nordtvedt, 1987; Shao et al., 2013) is applicable to PSR J1909$-$3744, we can use Eq. (12) to determine $\Omega$ (cf. Fig. 6). The corresponding allowed values of $\Omega$ are indicated by the gray bands in Fig. 14, from which we can obtain,

The limit is fractionally better than that from PSR J1738+0333. More importantly, our limit here is obtained via a more straightforward method, directly using Eq. (31), and sophisticated statistical analysis was carried out using the Bayesian inference.