© Y. J. Guo et al. 2021

1. Introduction

PSR J2222−0137 is a pulsar with a spin period (P) of 32.8 ms, which was discovered in the Green Bank Telescope (GBT) 350 MHz drift-scan pulsar survey (Boyles et al. 2013). It is in a binary system with an orbital period (P_b) of 2.44576 days and a projected semi-major axis of the pulsar’s orbit (x)¹ of 10.848 light-seconds (lt-s).

The small spin period derivative (5.8  ×  10⁻²⁰) implies that the pulsar was recycled by accretion of matter from its companion, during which event tidal torques would have circularised the orbit (Verbunt & Phinney 1995; Sepinsky et al. 2010). The fact that the orbit has a low eccentricity at present (e = 0.00038) implies that the companion has since become a white dwarf star (WD): Had it become instead a neutron star (NS), the associated supernova event would have caused a significant instantaneous mass loss and possibly a large kick that would in turn have, with very high probability, increased the eccentricity of the system by about three orders of magnitude (Tauris et al. 2017). The mass function of 0.229 M_⊙ implies that this WD is relatively massive. Optical observations have not detected the companion (Kaplan et al. 2014), implying that it is the coolest WD currently known.

As discussed by Cognard et al. (2017), this system has several characteristics that make it a unique gravitational laboratory. First, with a dispersion measure of 3.28 pc cm⁻³, it is one of the closest pulsars to the Solar System. This motivated a very long baseline interferometry (VLBI) astrometric campaign, from which Deller et al. (2013) obtained the most precise VLBI distance for any pulsar as well as precise values for the position and proper motion. Second, it is one of the very few systems where the orbital motion of the pulsar can be detected from VLBI astrometry, yielding a measurement of the longitude of ascending node, Ω.

Third, the edge-on orbit, the good timing precision, and the large mass of the companion allow a highly significant detection of the Shapiro delay, which was originally detected by Kaplan et al. (2014). From this effect, Cognard et al. (2017) obtained m_p  =  1.76(6) M_⊙ and m_c  =  1.293(25) M_⊙, where M_⊙ represents the solar mass parameter² and the numbers in parentheses represent, as in the remainder of the work, the 68% uncertainties in the last digit. This is the most massive double degenerate system known in our Galaxy³. Furthermore, it is also the largest NS birth mass known.

Fourth, the timing precision and the small but highly significant x ⋅ e product allow a measurement of the rate of advance of periastron ($\dot{\omega }$) for this system (Cognard et al. 2017). This allows precise and redundant measurements of the masses of the components of the system, and therefore a test of general relativity (GR).

Fifth, the timing precision and the large x/P_b ratio allow an unusually precise measurement of the variation in the orbital period (Ṗ_b). This can be compared to unusually precise theoretical predictions: The kinematic contributions to Ṗ_b owe their precision to the distance measurement; the predictions for the orbital decay according to GR and alternative gravity theories owe their precision to the well-measured masses.

Sixth, in this regard, the large difference in the compactness of the components of the system – a pulsar and a WD – is very important. Several alternative theories of gravity predict, for such systems, the emission of dipolar gravitational waves (DGWs) in addition to the quadrupolar gravitational waves predicted by GR (Eardley 1975; Damour & Esposito-Farèse 1992; Gérard & Wiaux 2002). This could potentially be detected in the measurement of Ṗ_b.

Finally, the system is exceptional even among the pulsar–WD systems that have been used to derive stringent limits on DGW emission, such as PSR J1738+0333 (Freire et al. 2012) and PSR J0348+0432 (Antoniadis et al. 2013). There are two reasons for this. Firstly, its mass estimates are more precise than for these two systems (Antoniadis et al. 2012, 2013); this is important for the interpretation of the measurements. Secondly, previous authors (Shibata et al. 2014) have made it clear that these tests should be carried out for a variety of NS masses in order to exclude strongly non-linear phenomena, such as spontaneous scalarisation (Damour & Esposito-Farèse 1993). Interestingly, the precise mass of PSR J2222−0137 places it in an intermediate, previous unexplored mass range. For this reason, even the relatively low-precision measurement of Ṗ_b by Cognard et al. (2017) has already provided useful constraints on alternative theories of gravity (Shao et al. 2017).

The dataset described by Cognard et al. (2017) ends in January 2017. Since then, regular pulsar timing observations with the 100 m Effelsberg radio telescope, the Lovell 76 m radio telescope, and the Nançay radio telescope have continued, using the same observing setups described by Cognard et al. (2017). The Effelsberg observations in particular were obtained in a set of orbital campaigns: To the two campaigns mentioned by Cognard et al. (2017), three more were added, which occurred in January 2018, August 2019, and October 2020. The last observation used in this work was taken in May 2021.

In addition, most of the pulse times of arrival (ToAs) derived from early discovery and follow-up data with the GBT, which were used by Boyles et al. (2013) and Kaplan et al. (2014), have now been included in our analysis. The addition of these ToAs significantly extends our timing baseline into the past, which now starts on June 23, 2009, and has a length of almost 12 years. This improves the precision of the measurements of proper motion, the rate of advance of periastron ($\dot{\omega }$), and, especially, the derivative of the orbital period, Ṗ_b.

Finally, we have observed the pulsar with the central beam of the 19-beam receiver of the Five hundred meter Aperture Spherical Telescope (FAST, Yao et al. 2021). These FAST data, taken at a frequency range between 1.0 and 1.5 GHz (40 MHz at each edge of the band is excised in the data reduction), provide the best polarimetric profile of PSR J2222−0137 to date, which is discussed below. The FAST observations will, in the near future, contribute greatly to improving the timing of this system.

In this work we address the proximate objective of this long-term timing project, which is to improve the precision of the physical parameters of the system, especially the post-Keplerian (PK) parameters (which quantify, in a theory-independent way, the observed relativistic effects, such as the aforementioned Shapiro delay, $\dot{\omega }$, and Ṗ_b). The ultimate objectives of the project – improved constraints on the nature of gravitational radiation and the behaviour of gravity for strongly self-gravitating systems – will be addressed in subsequent work.

The remainder of the paper is structured as follows: In Sect. 2 we revisit the VLBI astrometry of this pulsar using a coordinate convention for the pulsar’s orbit that is in full agreement with that used in pulsar timing. We also derive an absolute position for the pulsar with more realistic uncertainties. In Sect. 3 we present the new results on the polarimetry of the pulsar from a high S/N detection with FAST. In Sect. 4 we describe the processing of the radio timing data. In Sect. 5 we present the main timing results, with a detailed review of the different timing parameters and how they compare with previous estimates. In Sect. 6 we make a self-consistent estimate of the component masses and orbital orientation of the system and compare these with the VLBI results and the orientation of the pulsar obtained from the polarimetry. In Sect. 7 we list in detail the different contributions to Ṗ_b and estimate the observed excess variation relative to the GR prediction, which appears to be consistent with zero. Finally, in Sect. 8 we summarise our results and briefly point to further work on the implications of these timing results.

2. Re-analysis of the VLBI astrometry

2.1. Reference frame

Before describing the re-analysis of the VLBI data, we must first define the geometric parameters used in this paper. We use the ‘observer’s convention’, which is assumed for calculating all kinematic effects in the Damour-Deruelle-Kopeikin (DDK) orbital model in TEMPO and the T2 orbital model used in TEMPO2. This is different from the convention described by Damour & Taylor (1992) and used by Kramer et al. (2021) (see their Fig. 7).

To illustrate this convention, we refer the reader to Fig. 1. In this figure, the fundamental reference frame is centred at the centre of mass of the binary and has three axes, depicted in yellow: one towards the north (as seen by the observer), the second points to the east – both of these define the plane of the sky, in yellow – and the perpendicular direction towards the observer (n).

Fig. 1. Geometric parameters to be used throughout the paper according to the ‘observer’s convention’, which is used throughout the text. For the description of the fundamental reference frame (axes in yellow) and the orbital reference frame (axes in blue), see Sect. 2. For the orientation of the pulsar angular momentum, S, and its magnetic axis, μ, relative to this frame, see the detailed explanation in Sect. 3. In this figure, S is represented as being parallel to the direction of the orbital angular momentum, k. This is true for binaries such as PSR J2222−0137, where the pulsar was fully recycled by mass accretion from the companion and the companion then evolved into a WD. If the second-formed object in the system were a NS, then the kick from the supernova event that formed it might substantially change the orientation of the orbit (and thus the orientation of k) relative to S and greatly increase its eccentricity. The figure is by Vivek Venkatraman Krishnan.

The orbital plane is indicated in blue. It crosses the plane of the sky in the line of nodes. The ascending node is one of the two points where the pulsar, in its orbital motion, crosses the plane of the sky, it is the one where its distance from the observer is increasing. Its opposite point is the descending node. The orbital reference frame is defined by the three blue axes: the line pointing to the ascending node i, a perpendicular direction within the orbital plane, j, which defines superior conjunction, and finally a direction perpendicular to the orbital plane, k; which is parallel to the orbital angular momentum (not represented in the figure).

The full orientation of the orbital plane relative to the reference is defined by two angles. The first is Ω. This is a ‘position angle’ (PA); these are measured in the plane of the sky, starting from the north-pointing axis and then increasing anti-clockwise, as seen from the observer. Ω is the PA of the ascending node (the PA of the descending node is given by Ω ± 180 deg). The second angle, the orbital inclination (i), is the angle between k and n. We can therefore see that, for i < 90 deg, the line of sight (LOS) component of the orbital angular momentum would point towards the Earth.

2.2. Re-analysis of the differential astrometry

In the analysis of Deller et al. (2013), it was assumed that the vector from the centre of mass to the pulsar (the pulsar’s ‘position vector’) follows the same convention as the angular momenta: its LOS component is positive when it points to us, the same sense as n in Fig. 1. This would imply, as described there, that orbital longitudes are measured from the descending node. However, and unlike the Damour & Taylor (1992) convention, the observer’s convention is not coherent in this respect: in pulsar timing, the LOS component of the position vector – the geometric component of the quantity measured most directly in pulsar timing, the time delays of the radio pulses – is always positive if it points away from us. This means that the orbital longitudes are always measured from the ascending node.

This prompted us to perform a re-analysis of the VLBI data, this time using a convention fully in agreement with the convention used in pulsar timing. We make use of the code binary_pulsar_MCMC.py⁴, which infers the pulsar reference position, proper motion, parallax, and unknown orbital parameters in a Bayesian fashion using the measured VLBI positions and uncertainties. This contrasts with Deller et al. (2013), in which the astrometric and orbital parameters and their uncertainties were estimated using bootstrap sampling and linear least squares fitting. The refined timing ephemeris also allowed us to place updated prior ranges on the orbital inclination during the fitting process: The inclination was restricted to 85.25 ± 0.25 deg or 94.75 ± 0.25 deg, whereas in Deller et al. 2013 only two values were trialled (86.9 and 93.1 deg). We do not apply any constraints based on the value of ẋ measured by pulsar timing (described in subsequent sections).

The numerical values of the resulting parameters and their symbols, to be used in the remainder of the paper, are presented in Table 1 together with the earlier estimates, their difference, and the significance of the change; the position coordinates, α′ and δ′, are calculated for the same epoch as Deller et al. (2013), and assuming the same position for the in-beam calibrator, FIRST J222201−013236 (hereafter J2222−0132); the results refer to the position of the pulsar as seen from the Solar System barycentre. Figure 2 shows a correlation plot for these parameters.

Fig. 2. Correlation plots for the astrometric parameters measured by VLBI listed in Table 1. All parameters have well-defined uncertainties, except for i: Systems with inclinations of 85.25 and 94.75 deg will have nearly identical motions in the plane of the sky. For this reason, the orbital inclinations were sampled for two intervals that correspond to the best estimate of sin i from timing: 85.0 − 85.5 deg and 94.5 − 95.0 deg. We present the offsets (Δ) relative to our best values of α′ and δ′ in the first column of Table 1.

By far the most significant change is that of Ω. Within measurement uncertainties, it is about 180 deg offset from the estimate by Deller et al. (2013). This is consistent with the idea that, basically, there was an exchange between descending and ascending node. When doing such an exchange, we find that the orbital offset of the pulsar relative to the centre of mass in δ′ is nearly symmetric under that change; for that reason, δ′ and μ_δ are 0.5 and 0.1-σ consistent with the values published by Deller et al. (2013).

However, this is not exactly true regarding α′: the orbital offset in α′ changes somewhat under the inversion (and to a smaller extent because of the slightly lower orbital inclination we derive compared to the earlier value by Kaplan et al. 2014). For this reason, it changes by about −1.8σ. It is also for this reason that μ_α has the second most significant change, about −0.9σ. The parallax changes by about −0.7σ; from its new value we obtain a new distance estimate of 268.0(1.2) pc.

2.3. Absolute VLBI position

The pulsar position provided by Deller et al. (2013) can be considered a relative position (α′ and δ′) with respect to the assumed position for the primary in-beam calibrator, J2222−0132. In order to compare with the timing position, an absolute VLBI position is required.

The absolute barycentric VLBI position of PSR J2222−0137 shown in Table 1 was estimated using the approach described by Ding et al. (2020). Our values of α and δ are different from the α′ and δ′ provided by Deller et al. (2013) for two reasons: (1) The absolute position of J2218−0335, the primary out-of-beam calibrator, is updated to the most recent estimate⁵, an update that directly shifts the absolute position of PSR J2222−0137 by the same amount, and (2) our estimate of the relative position of J2222−0132 with respect to J2218−0335 is refined using the archival VLBI data.

To obtain the uncertainty of α and δ, we have to take into account three additional contributions on top of the uncertainty in the relative separation between PSR J2222−0137 and J2222−0132: (1) the uncertainty of the absolute position of J2218−0335, (2) the uncertainty of the relative position of J2222−0132 with respect to J2218−0335 and (3) the unknown frequency-dependent source structure (‘core shift’) of J2218−0335 between the higher frequencies at which its reference position is defined and the lower frequencies at which it is observed here. In Deller et al. (2013), the first term (with a value of ∼0.1 mas) was assumed to dominate, but we show here that this is not the case. Using the scatter in the per-epoch positions of J2222−0132 obtained via phase referencing to J2218−0335 without self-calibration, we conservatively determined that the second term contributes errors of 0.6 mas in right ascension and 2.3 mas in declination (using a weighted mean of these per-epoch positions would yield a smaller uncertainty on this mean absolute position, but one that would depend sensitively on the assumed input uncertainties). For the core shift contribution, we adopted 0.8 mas in each direction, which is the median core shift between 1.5 GHz and 8 GHz reported by Sokolovsky et al. (2011), and added it in quadrature with other uncertainties.

The uncertainties of α and δ are ∼50 times larger than the uncertainties of α′ and δ′; this quantifies the difference in precision of in-beam astrometry and absolute astrometry for this system. Future multi-frequency observations targeting J2222−0132 could considerably reduce the uncertainty in α and δ. Regardless, we note that the estimation of α and δ is independent from (and has no impact on) α′, δ′, and the other parameters in Table 1, which are derived purely from in-beam astrometry.

In all our subsequent analysis, when referring to the VLBI astrometry, we use the re-derived values in Table 1.

3. Pulse polarisation results

We now refer again to Fig. 1 in order to define the angles used to describe the geometry of the pulsar. We now place the pulsar at the origin, in order to compare the directions of the orbital and pulsar vectors. The angle between the magnetic axis of the pulsar (μ) and the spin axis (S) is known as α_p. The closest approach of μ to the LOS (n) is an angle known as the impact parameter, β; this has to be smaller than the total radius of the pulsar emission cone (ρ), otherwise the pulsar beam does not intersect our LOS. Thus, the angle between S and n is known as ζ, and is given by α_p + β. The polarisation angle ψ is the PA of the projection of the spin axis of the pulsar in the plane of the sky, S_p. If the spin axis of the pulsar is aligned with the orbital angular momentum, then ζ = i and ψ = Ω + 90 deg.

In Fig. 3 we present a high-sensitivity (S/N = 23 000), high-resolution pulse profile of PSR J2222−0137 obtained with FAST at a central frequency of 1250 MHz, the integration time is 1735 s. The observing setup and data reduction is similar to that described by Yao et al. (2021). This profile is shown in the PSR/IEEE convention (van Straten et al. 2010).

Fig. 3. Polarisation profile as observed with FAST at a central frequency of 1250 MHz, averaged over a bandwidth of 420 MHz. Top panel: total intensity (black), the linear intensity (red), and circular polarisation (blue) as a function of longitude. The inset shows an enlarged version of the main pulse. Middle panel: low intensity level of the profile across the full rotational period, revealing the existence of an interpulse that is separated from the main pulse by about half a period. Bottom panel: values of the PA of the linear emission as a function of the longitude. An RVM has been fitted to the black PA values. The result is shown as a solid black line, while the dashed line indicates the RVM solution separated by 90 deg.

The quality of the profile not only allows the resolution of a number of features in the total, linear, and circular intensity but also reveals a number of faint profile components. Firstly, the main pulse shows polarised low-intensity components on both its leading and trailing side. Secondly, and most crucially, the sensitivity of FAST also uncovers a faint interpulse component separated from the main pulse by about half a period. The detection of these faint components is not a mere curiosity; it reveals the large-scale dipolar structure of the magnetic field of the pulsar, which is crucial for a determination of its geometry.

Inspecting the profile, the sense reversal of Stokes V (circular polarisation) underneath the main peak suggests that this longitude range can be identified with the location of the fiducial plane that is defined by the spin vector, the magnetic axis and the direction of the observer (see e.g., Lorimer & Kramer 2005). At the same time, the sudden drops in linearly polarised intensity, combined with the very rapid variation in V resolved in the FAST profile, reveal the existence of orthogonal polarisation modes (OPMs). These are clearly responsible for the large variation in the PA of the linear polarisation at these longitudes, which we mark in grey in Fig. 3.

Overall, however, the PA swing under the main pulse shows a positive slope, which in terms of the ‘rotating vector model’ (RVM; Radhakrishnan & Cooke 1969) implies a negative value of β. This is confirmed by a blind fit of the RVM to the black PA values, which results in the solid black line shown in the bottom panel of Fig. 3. The grey PA points corresponding to the rapid changes in PA and V are ignored here, even though including those PA values does not change the overall result. Instead, the presence of the interpulse with well-defined PA values provides important leverage and helps to break the usually existing co-variances between the RVM parameters (see e.g., Johnston & Kramer 2019).

We determine the geometry of the pulsar by fitting the RVM model using a Bayesian optimisation method as described by Johnston & Kramer (2019) and Kramer et al. (2021). Using uniform priors, we derive ζ, β, as well as Φ₀ (the location of the aforementioned fiducial plane relative to an arbitrary reference longitude on the neutron star; see Fig. 3) and the absolute PA defined by the linear polarisation at Φ₀, ψ₀. The code we used allows for the possible existence of OPMs, and indeed a number of PA values, especially those of the interpulse, follow an RVM swing (dashed line) that is separated from the main pulse swing (solid line) by 90 deg. This observation, that main pulse and interpulse emit in different modes, is quite common (Johnston & Kramer 2019).

The mean and 99% percentiles of the posterior distribution for these quantities are well defined and symmetric. They are ψ₀ = 40(1) deg, Φ₀ = 188.8(6) deg, ζ = 84.0(8) deg, and β = −7.2(6) deg; thus, α = 77(1) deg. We refer to Kramer et al. (2021) for a critical discussion of the assumptions and reliability of RVM fits in this context.

Apart from ζ, the second parameter that defines the 3D orientation of the pulsar spin is ψ. To determine it, we first remark that ψ₀ is measured at 1250 MHz (i.e., after it has been affected by Faraday rotation in the interstellar medium). Taking the rotation measure and its uncertainties into account, we obtain the de-rotated PA ${\psi }_0^{\prime }=30(10)\mathrm{deg}$. Because of the aforementioned OPMs, ${\psi }_0^{\prime }$ can differ from ψ by 0, ±90 deg. Thus ψ can have the following values: −60(10),30(10) and 120(10) deg.

Regarding the fidelity of polarisation measurements with FAST, we conducted tests on the centre beam of the FAST 19-beam system by observing some bright millisecond pulsars and solving for the full Muller matrix based on their known polarisation profiles. We found the cross-coupling between the two polarisations to be smaller than 1%. When tracking an object such as PSR J2222−0137, the FAST telescope measures and maintains the orientation of its 19-beam receiver at a fixed angle on the sky with a precision better than 0.1 degrees. Therefore, the systematic errors in the polarimetry of the 19-beam observing system should be negligible in our observations.

We now compare the polarisation profile from FAST with that presented by Cognard et al. (2017), which was taken with the Nançay Ultimate Pulsar Processing Instrument (NUPPI). There, we see the same modest degree of linear and circular polarisation; however, both the PA swing and the sign of the circular polarisation differ. We identify this change with a swap in the previous NUPPI data in both Stokes Q and V. Indeed, this was discovered before the FAST polarimetric data became available by a comparison of the NUPPI polarimetric profiles of other pulsars with previously published pulse profiles (see Guillemot et al., in prep.). New NUPPI data agree with the data presented here, albeit with lower sensitivity. Even though the PA swing presented by Cognard et al. (2017) showed the opposite sense to that shown here, an estimate of the geometry led to the same results (with much larger uncertainty), the reason is that instead of the PSR/IEEE convention, that work used the RVM convention, which reversed the sign of β a second time.

Finally, we can also compare the FAST profile to the recently published MeerKAT data (Kramer et al. 2021). The time resolution of the profile shown here is a factor of ∼20 better, so that the MeerKAT PA swing is smeared out in comparison and therefore appears to be different. We confirmed the consistency and correctness of both results by smearing the FAST data to the MeerKAT resolution.

4. Processing of the timing data

The timing observations used for this project are summarised in Table 2. As in Cognard et al. (2017), the reduction of our timing data (mainly radio frequency inference mitigation and polarisation calibration) was performed using the PSRCHIVE package⁶ (Hotan et al. 2004). For the observations and data analysis of the GBT data, we refer the reader to Kaplan et al. (2014). Here we used the 820 MHz and the L-band ToAs only, as these are the most extensive and useful datasets; the addition of the smaller P- and S-band ToA sets does not change any parameters noticeably, only increasing the complexity of the analysis. In what follows, we describe the improvements of the data analysis.

4.1. Derivation of the pulse times of arrival

As mentioned in Sect. 3, we now have a better understanding of the polarisation characteristics of the Nançay telescope. In order to take full advantage of these new polarimetric profiles, we used the matrix template matching (MTM) method implemented in the PAT routine of PSRCHIVE (van Straten 2006) to derive ToAs from the Nançay, but also from the Effelsberg data. In addition to the total intensity, the MTM method exploits the timing information available in the polarisation of the pulsar signal, by modelling the transformation between two polarised light curves in the Fourier domain.

This method seems to overestimate the ToA uncertainty to some extent (we had to multiply the ToA uncertainties from these two datasets by numbers smaller than 1 in order to obtain a reduced χ² of 1.0, see Table 2). However, it helped achieve a reduction of the root-mean-square (rms) of the timing residuals (which are the ToA minus the model prediction for its rotation number): The current weighted residual rms of 2.8 μs is significantly better than the global weighted rms of 3.4 μs reported by Cognard et al. (2017).

4.2. Dispersion measure model

Another change is the use of the piecewise-constant function (‘DMX’ model, Demorest et al. 2013) to describe the dispersion measure (DM) variations. In Cognard et al. (2017), a simple model using the DM and its first derivative was used. As we describe later in the paper, there are DM variations on relatively short timescales (tens of days), which are not captured by any simple model with a few DM derivatives and are large enough to influence our measurements of parameters that have long-term time signatures, like the spin period, its derivative, the position and especially the proper motion. For this reason, we used the DMX model. Our DMX model does not fit for a DM offset for the earlier GBT data since for those a single ToA was produced for the whole band.

4.3. Timing analysis and orbital models

The timing analysis is performed with TEMPO⁷, using the latest available version, 13.103. The telescope ToAs are first converted to the Bureau International des Poids et Mésures 2019 (BIPM2019) timescale, and then converted to the Solar System barycentre using a) the latest information on the Earth rotation (the Universal Time) provided by the International Earth Rotation Service and b) the Jet Propulsion Laboratory’s DE440 solar system ephemeris (Park et al. 2021). The resulting timing parameters are presented in Barycentric Dynamical Time (TDB). We use three orbital models to describe the orbital motion of the pulsar and the propagation of the radio signals to the Earth, all based on the Damour & Deruelle (DD) timing model (Damour & Deruelle 1986).

The first, the general relativity version of the DD model (DDGR), is a theory-dependent model that assumes the validity of GR and fits directly for the total mass of the system (M) and the companion mass (M_c). This model does not do a fully coherent analysis of the kinematic information.

The second, DDK, is a theory-independent model that takes into account the kinematic effects described by Kopeikin (1995, 1996) and implemented in TEMPO by Ingrid H. Stairs. We used this particular model as the basis of a self-consistent Bayesian analysis of the system, which also assumes the validity of GR.

The third, ELL1H+, is a theory-independent model based on the low-eccentricity approximation of the DD model known as ELL1 (Lange et al. 2001). Instead of the Keplerian orbital parameters of time of passage through periastron (T₀), longitude of periastron (ω), and orbital eccentricity (e) used in the DDGR and DDK models, the ELL1 and ELL1+ models use the times of ascending node (T_asc), and the Laplace-Lagrange parameters ϵ₁ = esinω and ϵ₂ = e cos ω. For the ELL1H+ model, we implemented the orthometric parameterisation of the Shapiro delay using its exact expression, Eq. (31) of Freire & Wex (2010).

The latest versions of the third model (referred to with the + sign) include extra terms for the expansion of the Rømer delay in orbital harmonics of order xe² (Zhu et al. 2019); these were implemented in TEMPO distributions 13.102 and later. Thus, and unlike the original ELL1 model, the new ELL1+ and ELL1H+ models can describe the orbit of PSR J2222−0137 well: The neglected xe³ terms of the Rømer delay are of the order of 0.6 ns, a quantity that is small in comparison with our timing precision. The advantage of these models is the avoidance of the strong correlation between T₀ and ω observed in the DD-like models. Also, by re-defining P_b as the time between passages through ascending node, we avoid its strong correlation with $\dot{\omega }$ seen in the DD-like models. Finally, by using the orthometric parameterisation, we avoid the large correlation between the Shapiro delay parameters r and s seen in the DD, DDK, and ELL1+ models (−0.86 in our timing); the correlation between the orthometric parameters is −0.54.

This lack of correlations has practical advantages: as we see below, for PSR J2222−0137, the T_asc and P_b measured in the ELL1H+ model are, respectively, ∼5400 and ∼4040 times more precise than the T₀ and P_b measured by the DDGR and DDK models. This has a consequence: we can state these parameters to their actual precision and still retain enough accuracy in the description of the orbital motion to do the timing. The DDGR and DDK values of T₀ and P_b in Table 4, also stated to their own uncertainties, are not precise enough for this purpose⁸.

4.4. Template alignment

Another important improvement in the data analysis was the use of pulse profile templates that have consistent phase definitions for all our datasets. In this way, the ToAs refer to a consistent longitude of the NS.

Normally, when combining data from different telescopes, the different datasets are not entirely consistent, because of different delays in the signal paths of the different observing systems. In order to take that into account, an arbitrary time offset Δt between datasets is fitted; this is done in TEMPO by bracketing the ToAs from a particular observing system with two JUMP statements. In a first iteration, TEMPO assumes these are phase offsets, unless they are provided by the ephemeris. Thus, no estimates of Δt are taken into account at this stage. At the end of that first iteration, these phase offsets are converted into the values of Δts, which are then added to the ToAs of the different datasets in subsequent iterations. These values should accurately characterise the different delays between the different observing systems.

However, if the templates have no consistent phase definitions, these Δts will be biased. We imagine two datasets taken at the same observing times, where the second has an extra signal delay relative to the first given by Δt. The ToAs for both datasets are derived with two different templates, where the second has a phase difference of −0.5 < Δϕ < 0.5 relative to the first. At the end of the first TEMPO iteration, the total observed phase difference between the two datasets is converted into a time offset for the second dataset (relative to the first) given by Δt′=Δt + ΔϕP (where P is the spin period of the pulsar). In subsequent iterations this biased estimate is added to the ToAs of the second dataset.

This has implications for our measurement of the orbital motion. If we measure T_asc for both datasets separately, without providing any time offsets (as in the first TEMPO iteration), we find that the second T_asc will differ from the first by ΔT_asc = Δt. If we do a second iteration where we use the previously determined Δt′s, we find that ΔT_asc = ΔϕP.

This is not a problem if the uncertainty of this measurement, δT_asc, is larger than ΔT_asc. If δT_asc < Δt, then the orbital phases of the two datasets are not consistent at the first iteration, a normal situation that gets fixed in following iterations with the Δt estimates. However, if δT_asc < ΔϕP, then this also happens for the second and following iterations. In particular, if ΔϕP > Δt, we even have a degradation of the quality of the TEMPO fit between the first iteration (which assumes all differences are phase offsets, including Δt) and later iterations (which assume all differences are time offsets, including ΔϕP). This does provide an easy way of diagnosing the problem, and is indeed is how we identified it in the first place.

For PSR J2222−0137, a template misalignment of 1/3 in spin phase – typical of the misalignments still present in the analysis of Cognard et al. (2017) – yields ΔT_asc = ΔϕP ∼ 10 ms. For δT_asc we have a value of 0.17 ms (see Table 3). Therefore, for this pulsar it is important to align the profiles to a phase precision of at least δT_asc/P  =  0.005. One can alternatively make TEMPO run with a single iteration, but this is only safe to do if we know in advance that, for all datasets, Δt < δT_asc.

Fixing this problem resulted in a major improvement in our timing of this pulsar. First, there was no longer a degradation by ∼1000 of the χ² from the first to the second TEMPO iterations - all stayed at consistent values close to that of our best solution, χ² ∼ 10628. Furthermore, once this alignment was done, we no longer needed to introduce added errors in quadrature to the different ToA datasets (this is done in TEMPO with EQUAD statements; their numerical values are in Table 1 of Cognard et al. 2017); this has contributed to the decrease in the residual rms mentioned above. Finally, the FD (frequency-dependent) parameters, which describe non-dispersive variations in the ToAs with radio frequency (see Table 2 of Cognard et al. 2017) also became unnecessary. The excellent quality of the first iteration implies that a single iteration without predetermined time offsets is fine and that, therefore, all our datasets have Δt < δT_asc.

A final improvement in the timing analysis will be described below, when we analyse the Shapiro delay and the reasons for the lower mass values derived by Kaplan et al. (2014).

5. Timing results

The timing parameters resulting from the different timing models listed above are given in Tables 3 and 4. In Table 3, we present the parameters of the ELL1H+ model, either fitting for position and proper motion or assuming the VLBI proper motion derived in Sect. 2. Table 4 shows the orbital parameters derived with the DDGR and DDK models, and the results of the self-consistent Bayesian grid obtained with the DDK model, all derived assuming the VLBI proper motion.

Figure 4 shows the residuals obtained with the ELL1H+ model where we fit for position and proper motion. For the 10772 ToAs used in our analysis, we obtain a weighted rms residual of 2.781 μs and a reduced χ² of 1.0051. The rms residuals for the individual observing systems are presented in Table 2, where we also listed the multiplication factor (this is done in TEMPO with EFAC statements) for the ToA uncertainties in order to achieve a reduced χ² of 1 for each dataset.

Fig. 4. Twelve-year timing data of PSR J2222−0137, for the ELL1H+ timing model where we used the DMX model and fit for position and proper motion (middle column in Table 3). Top panel: DM offsets relative to the reference DM (3.2805 cm−3 pc) as a function of date. Middle and bottom panels: timing residuals as a function of date (middle) and orbital phase (bottom, measured from ascending node). The residuals are displayed with different offsets for each instrument: EFF – Effelsberg; JB – Lovell telescope (Jodrell Bank); NCY-S – Nançay at the S band; NCY-L – Nançay at the L band; GBT-820 – GBT at 800 MHz; and GBT-1500 – GBT at 1500 MHz.

In what follows, we call the reader’s attention to the more important timing parameters. We also make some comparisons between our results and those presented by Cognard et al. (2017). Since the systematic issues discussed in Sect. 4.4 were still present in that earlier analysis, some of the differences in the values reported by both works are significant, particularly on the ẋ parameter.

5.1. Astrometric parameters

All positions reported in Table 3 are barycentric positions measured at MJD = 55743, making them directly comparable with the VLBI positions in Table 1. Below α and δ, we list their offsets relative to the VLBI absolute values in Table 1. Unlike the analysis presented by Cognard et al. (2017), these are consistent with zero for all cases. The primary reason is the more realistic uncertainty estimates for the absolute position presented in Sect. 2. Because the timing position is more precise, we fit for it from now on.

First, as Cognard et al. (2017), we use a simple description of the variation in the DM, which employs a small number of DM time derivatives. If we fit for parallax, proper motion, and position simultaneously, we obtain a timing parallax of 3.710(79) mas, which is in near perfect agreement with the VLBI parallax, $\varpi =3.{730}_{-0.016}^{+0.015}$ mas. Our timing parallax is ∼4 times more precise than that of Cognard et al. (2017), but its uncertainty is still ∼6 times larger than the VLBI parallax. For this reason, we use the VLBI parallax from now on.

In the first column of Table 3 we fit for the proper motion and DM derivatives. In this case, we get a χ² of 10799.74. The proper motion is μ_α = 44.61(3) mas yr⁻¹ and μ_δ = −5.26(5) mas yr⁻¹. Compared to the VLBI proper motion in the third column, the differences are Δμ_α = −0.09(5) mas yr⁻¹ and Δμ_δ = 0.43(9) mas yr⁻¹, which are ∼2σ and ∼5σ significant, respectively.

We have found that there are DM variations on short timescales which could affect the astrometric parameters. In order to take these short-term DM variations into account, we use the DMX model, fitting for a DM offset for ToAs within gaps of 60 days, the resulting DM offsets are depicted in the top panel of Fig. 4. This interval was chosen and adhered to before a detailed consistency analysis of all the PK parameters. The results are presented in the second column of Table 3, where we still fit for proper motion. This causes a very significant decrease in the χ², to 10609.25, but also causes (predictably) a degradation in the precision of all other timing parameters, especially the position and proper motion. The difference of this proper motion to the VLBI values is Δμ_α = −0.04(6) mas yr⁻¹ and Δμ_δ = 0.23(14) mas yr⁻¹ (i.e., they are 2σ consistent).

Finally, in the third column, we also use the DMX model but assume the more precise VLBI proper motion, as in all subsequent discussions. In this case, the χ² increases to 10629.20. This causes changes in the remaining parameters within their 1σ uncertainties, which is expected from the consistency of the proper motion. Nevertheless, a more precise VLBI proper motion will be important to help with future timing of this system.

5.2. Shapiro delay

The Shapiro delay was first detected by Kaplan et al. (2014), who used it to obtain M_p = 1.20(14) M_⊙ and M_c = 1.05(6) M_⊙. However, Cognard et al. (2017) found improved and significantly larger masses: M_p = 1.76(6) M_⊙ and M_c = 1.293(25) M_⊙.

In this work, we obtain an unusually precise measurement of the Shapiro delay: h₃ is 187-σ significant (i.e., an uncertainty of about 27 ns, which is 0.53% of the measured value); ς is measured with a relative uncertainty of 0.152%. Without the Shapiro delay, the residuals would have large trends (see Fig. 5). Assuming GR, we obtain: M_p = 1.81(3) M_⊙, M_c = 1.312(9) M_⊙ and i = 85.30(9) deg or i′ = 94.70(9) deg, respectively. The total system mass is M = 3.12(3) M_⊙. These results are 1-σ consistent with the mass measurement in Cognard et al. (2017) but more than twice as precise.

Fig. 5. Residuals as a function of the orbital phase for the same ephemeris as Fig. 4. Top panel: we display the residuals predicted by the model but without taking into account the Shapiro delay. Lower panel: we see the same while fitting for other Keplerian parameters. As we see, some of the Shapiro delay is absorbed by this fit, but some still remains; this remnant represents the ‘measurable’ part of the Shapiro delay; this is given, in the limit of perfect orbital sampling, by Eq. (31) of Freire & Wex (2010). In this plot, the Nançay residuals are displayed in black, the Effelsberg residuals in red, and all others in grey, all without error bars for the sake of clarity.

The addition of the early GBT data allows an investigation of the reasons for the low masses derived by Kaplan et al. (2014). As it turns out, this is not caused by the correlation between Shapiro delay and $\dot{\omega }$ in the GBT data, as suggested by Cognard et al. (2017), although that correlation, already identified by Kaplan et al. (2014), is real. Our much longer timing baseline, with its far better constrained timing parameters, has helped identify a set of six ToAs in the GBT 820 MHz data (taken on June 28, 2009) that have extra delays of 30–40 μs (i.e., ∼1/1000^th of a full rotation). The causes for these extra delays have not been found, but they are highly significant, since they are systematic and much larger than the uncertainties of those ToAs. If we exclude those ToAs, then the Shapiro delay and $\dot{\omega }$ obtained with the GBT data alone are in 1-σ agreement with the values obtained by Cognard et al. (2017). If we do not exclude them, then the masses we obtain with the GBT dataset are in near agreement with the values obtained by Kaplan et al. (2014).

The exclusion of these six ToAs has caused a significant decrease in the reduced χ² associated with the GBT data: Kaplan et al. (2014) needed to increase the uncertainty estimates of their ToAs by a factor of 2.7 in order to achieve a reduced χ² of 1. With the exclusion of those six ToAs, we can achieve the same using factors of 1.50 and 1.72 only (see the EFAC factors in Table 2); and did this with a timing solution that is strongly constrained by nine years of subsequent data. These smaller EFAC factors are much more commonly found in the timing of recycled pulsars. We are therefore confident that we have found the reason for the much lower mass estimates reported by Kaplan et al. (2014).

5.3. Advance of periastron

From our timing we derive ${\dot{\omega }}_{\mathrm{obs}}=0.09607(48)deg{\mathrm{yr}}^{-1}$. This is ∼6 times more precise than the measurement obtained by Cognard et al. (2017), ${\dot{\omega }}_{\mathrm{obs}}=0.1033(29)deg{\mathrm{yr}}^{-1}$; however, the difference is ∼2.5-σ significant. This is caused by the problem mentioned in Sect. 4, the use of inconsistent definitions of spin phase for the different datasets.

For the masses obtained from the Shapiro delay, GR predicts ${\dot{\omega }}_{\mathrm{GR}}=0.09576(67)deg{\mathrm{yr}}^{-1}$. Thus, our new measurement agrees with ${\dot{\omega }}_{\mathrm{GR}}$ within 1σ. Therefore, the same applies to the total mass of the system derived from both methods: assuming GR, we obtain from ${\dot{\omega }}_{\mathrm{obs}}$M  = 3.139(16) M_⊙. This represents a successful and precise (∼1%) test of GR. This is illustrated in Fig. 6 by the fact that all mass constraints from the different PK parameters intersect in the same regions of the diagrams.

Fig. 6. Mass-mass diagram of J2222−0137. In the main panel on the left we display the cos i − Mc plane, and on the right we display the Mp − Mc plane; the grey region is excluded by the constraint sin i ≤ 1. The black contours indicate the 68.3% and 95.4% confidence region derived from a 3D χ2 map of the cos i − Ω − M plane using the DDK orbital model combined with GR equations. The constraints (according to GR) from h3, ς, and $\dot{\omega }$ in the ELL1H+ model are shown in blue, green, and red lines, respectively, where the solid and dotted lines indicate the nominal and ±1σ measurements. The fact that they all meet in the same regions of these planes represents a ∼1% test of GR, which the theory passes. The side panels display the 1D PDFs for cos i (top left), Mp (top right), and Mc (right), with the median value and the 1σ and 2σ confidence intervals indicated as vertical lines.

This statement relies on the fact that the ${\dot{\omega }}_{\mathrm{obs}}$ is relativistic. As discussed by Cognard et al. (2017), the largest additional contribution to ${\dot{\omega }}_{\mathrm{obs}}$, which is caused by the proper motion of the system, is of the order of a few times 10⁻¹⁶ deg yr⁻¹ (i.e., about 100 times smaller than the current measurement uncertainty). Therefore, the assumption that ${\dot{\omega }}_{\mathrm{obs}}$ is relativistic is fully warranted. This also means that, in the near future, improved measurements of $\dot{\omega }$ will translate directly in a better constrained M which, as we can see in the right panel of Fig. 6, will also yield improved component masses.

We can use the DDGR model to combine $\dot{\omega }$ with the Shapiro delay and obtain self-consistent and more precise mass measurements; doing this we obtain M = 3.135(19) M_⊙, M_p = 1.820(14) M_⊙ and M_c = 1.315(6) M_⊙. These values are 1.5 times more precise but slightly larger than the value derived from Shapiro delay only. The χ² for that fit is nearly identical to that of the best ELL1H+ model, indicating again the self-consistency of the relativistic effects.

5.4. Variation in the projected semi-major axis of the pulsar’s orbit

Relative to Cognard et al. (2017), the orbital parameter that has had the most significant change was the rate of change of the projected semi-major axis, ẋ: When they assume the VLBI proper motion, they obtain ẋ = 3.5(30) × 10⁻¹⁵ lt-s s⁻¹, our new value is −7.76(48)×10⁻¹⁵ lt-s s⁻¹; the change is 3.8-σ significant. The uncertainty has decreased by a factor of six.

The observed value of ẋ is dominated by the secular change of the orbital inclination caused by the proper motion (Arzoumanian et al. 1996; Kopeikin 1996):

$$\begin{array}{c}\hfill {\dot{x}}_{\mathrm{PM}}=x\mu cotisin({\mathrm{\Theta }}_{\mu }-\mathrm{\Omega }).\end{array}$$(1)

Assuming the values of Ω and Θ_μ in Table 1 and of x, μ and i in Table 3, we obtain ẋ_k = ±6.12 × 10⁻¹⁵ lt-s s⁻¹, where the sign depends on whether i′ = 94.70(9) deg or i  =  85.30(9) deg. The magnitude of this effect is very close to maximal, since there is an angle of very nearly 90 deg between Θ_μ and Ω.

In Fig. 7 we present the constraints on the orbital inclination and Ω derived from this measurement. The ẋ curves do not intersect the constraint on cos i, but they come close in two regions, Ω  ≃  7 deg, i = 94.7 deg, and Ω  ≃  187 deg, i = 85.3 deg (Regions 1 and 2 in Table 6). This lack of an intersection reflects the fact that the difference between the most negative ẋ possible and the observed value is −1.64(48) × 10⁻¹⁵ lt-s s⁻¹, a difference that is 3.4σ significant. This is a robust result in all our fits: the ẋ is not strongly correlated with any other timing parameter.

Fig. 7. Orbital orientation constraints for PSR J2222−0137. In the main panel we display the cos i − Ω plane, where randomly aligned orbits will have a priori constant probability density. The black contours show the 68.3% and 95.4% confidence region derived from a 3D χ2 map of cos i − Ω − M space using the DDK model with the additional assumption that GR is the correct theory of gravity. The solid red lines indicate the regions that are consistent with the nominal and ±1σ measurements of ẋ, and the dashed vertical lines indicate the position of ς obtained in the ELL1H+ model. The dashed horizontal lines show the value of Ω from VLBI observations. The side panels display the 1D PDFs for cos i (top) and Ω (right). The solution with positive cos i has 99.66% of the total probability.

We now look into possible causes for this discrepancy. In Lorimer & Kramer (2004), there is an extensive list of effects that can contribute to the observed ẋ, these are: ẋ_GW, the variation in x caused by the orbital decay due to the emission of gravitational waves, ${\dot{x}}_{\dot{\mathrm{D}}}$, which is caused by the variation in the Doppler shift of the system D (analogous to the kinematic effect on Ṗ_b discussed in Sect. 7), ẋ_A, which is caused by a change in the aberration due to spin precession, which is itself caused by spin-orbit coupling, and ẋ_SO, the change of the orbital inclination also due to the spin-orbit coupling. The spin-orbit coupling has several classical and relativistic terms caused by the rotations of the pulsar and of the companion. Finally, the ẋ_ṁ is caused by mass loss.

We calculated these terms systematically; the results are also given in Table 5. For the mass loss contribution to ẋ, we assumed the spin-down energy of the pulsar. The ẋ_A term depends on the misalignment angle between the angular momentum of the pulsar and the angular momentum of the orbit, θ_p (not depicted in Fig. 1) and on a related precession phase, Φ_p. Since the pulsar was recycled with matter that was necessarily orbiting in the orbital plane, the two angular momenta should be very nearly aligned (i.e., θ_p ≃ 0). Thus ẋ_A can be safely neglected. Independent of that, one can use the polarisation information from Sect. 3 to constrain ẋ_A (see Eqs. (2.5a), (2.25b), (3.24), and (3.35) in Damour & Taylor 1992). The conclusion is the same.

The change in the orbital inclination due to spin-orbit coupling can be split into contributions from the pulsar and the companion (see e.g., Barker & O’Connell 1975; Damour & Taylor 1992; i.e., ẋ_SO = ẋ_SOp + ẋ_SOc). Due to the compactness of the pulsar, ẋ_SOp is clearly dominated by the Lense-Thirring effect, and depends on the orientation (Φ_p, θ_p) and magnitude of the pulsar spin. ẋ_SOc depends on the unknown spin period of the WD companion (P_c), as well as the unknown spin orientation (Φ_c, θ_c). For reasonable values of P_c (∼ hours⁹), ẋ_SOc is also dominated by the Lense-Thirring effect. Both contributions are negligible¹⁰. We want to add that, due to tidal torques during the Roche-lobe overflow episode, one would expect an alignment of the spin axis, resulting in θ_c ≃ 0, further suppressing any spin-orbit contribution to ẋ by the WD.

There is an extra term listed by Lorimer & Kramer (2005) that could be caused by a possible third component of the system. If such a component existed, the motion of the PSR J2222−0137 binary around the centre of mass of the triple would cause a non-linear variation in the Doppler shift of the pulsar, which would be observable (at the very least) as higher derivatives of the spin frequency of the pulsar. We do not detect any such effects in the timing, and thus we have no evidence of any extra component of this system. For this reason, we do not consider any related contributions to ẋ.

In addition to the terms listed by Lorimer & Kramer (2005), there is a contribution to ẋ (ẋ_γ) that results from the correlation of this term with one of the PK parameters, the Einstein delay (γ), this is inevitable for systems with small values of $\dot{\omega }$, see detailed discussion by Ridolfi et al. (2019). From their Eq. (25), we obtain:

$$\begin{array}{c}\hfill {\dot{x}}_{\gamma }=-\frac{\gamma \dot{\omega }}{\sqrt{1-e^2}}sin\omega .\end{array}$$(2)

For PSR J2222−0137, GR predicts a small Einstein delay (γ = 4.545 μs), which is largely the result of the small eccentricity. This and the relatively small $\dot{\omega }$ in turn yield a ẋ_γ that is about 1/2 of the measurement uncertainty of ẋ_obs, see Table 5. This term will contribute to the ẋ_obs measured with the ELL1H+ model, but it will not contribute to the ẋ_obs measured using the DDGR model since the latter already takes the Einstein delay into consideration as a relativistic effect. The difference between the two ẋ measurements is mostly (but not entirely) due to ẋ_γ.

Thus, it is clear that all terms are much smaller than ẋ_PM, and none of them can explain the difference between ẋ_obs and ẋ_PM.

5.5. Annual orbital parallax

An alternative explanation for the large ẋ has to do with the fact that, apart from the secular variation in x, there is a yearly modulation of x (and ω) caused by the changing viewing angle of the pulsar’s orbit due to the Earth’s orbital motion (Kopeikin 1995); this effect is known as the annual orbital parallax. This is not taken into account by the DDGR and ELL1H+ models; if it is significant, it could potentially be absorbed into the secular ẋ estimated by those models.

In order to test this we use the DDK model, which apart from the secular effects on x and ω, also takes into account their annual variations. These are not fitted explicitly (via, for instance, the ẋ parameter); they are calculated internally from all the geometric parameters of the model, in particular i and Ω, which are fitted directly.

However, before proceeding, we must urge a note of caution related to the use of this model. In it, we always use the γ calculated by the DDGR model. This is important because, as discussed above, γ is correlated with ẋ. If we fail to include γ in the DDK model, it will find biased values of i and Ω that will account for secular variation in x that is different (by ẋ_γ) from the ẋ caused by the proper motion of the system.

Generally, when one has a good constraint on the orbital inclination and a measurement of ẋ, there are four possible degenerate combinations of Ω and i that can satisfy those constraints, these reduce to two if ẋ is at its maximum possible value. Detecting the annual orbital parallax can eliminate this degeneracy (see e.g., Stovall et al. 2019). In Table 6 we can see that the local χ² minimum at Ω ∼ 190 deg is lower than the minimum at Ω ∼ 0 deg. This difference is an indication that the annual orbital parallax is significant (for a precise quantification, see following section). This is not surprising given the relatively small distance to the Earth and large size of the pulsar’s orbit.

Furthermore, as we can see in Table 4, a DDK model based on the best i, Ω combination has a lower χ² than the ELL1H+ or DDGR models where we fit for a freely varying ẋ. Thus, by taking the annual orbital parallax into account, we can find a model that provides a better fit to the data that assumes no changes in x other than those expected from the geometry of the system.

It is therefore clear that the DDK model provides a superior description of the orbital geometry and motion of the system. For that reason, we base all subsequent discussions on this particular orbital model.

6. A self-consistent estimate of the component masses and orbital orientation of the system

In order to better determine the uncertainties and correlations between the masses and orbital configuration, we made a self-consistent χ² map with the DDK model, where we additionally assume the validity of GR. Since, as discussed in Sect. 5.4, we expect no significant additional contributions to ẋ, we assume that any variations in x are caused by variations in i (i.e., caused by the geometry of the system and its orientation relative to the Earth), all of which are automatically taken into account by the DDK model.

The process is described in detail by Stovall et al. (2019); briefly, for each point in a cos i, Ω and M grid, we hold Ω and i fixed in its corresponding DDK model; from these two values, the astrometric parameters and the orbital Keplerian parameters all kinematic effects on x and ω are estimated internally by the model. Other relevant PK parameters ($M_{\mathrm{c}},\dot{\omega }$, γ, but not ẋ_b, which is kept as a free parameter because of other kinematic effects) are derived by our script from the known mass function, i and M using the GR equations, and then used as fixed inputs to the DDK model for that point of the grid. We then run TEMPO to fit this DDK model to the data, allowing all other timing parameters to vary, and record the value of χ², which is assigned to the respective point in the grid. The two regions of the cos i, Ω, M space that we sampled are listed in Table 6; for areas outside these two regions, the quality of the fit is just too poor to contribute any significant probability.

The resulting 3D grids of χ² values are then used to calculate a Bayesian 3D probability density function (PDF) for the cos i, Ω, M space (Splaver et al. 2002). This 3D PDF is then projected onto several planes and axes: 2D PDFs are calculated for the cos i-Ω and the derived cos i-M_c and M_c-M_p planes, and 1D PDFs are calculated for the Ω, cos i, M, and derived M_c and M_p axes. These 2D PDFs are represented by the black contours in the main panels of Figs. 6 and 7, and some of the 1D PDFs are represented in the side panels of those figures; their medians and ±1σ uncertainties are presented in Table 4. These numerical values are valid, but do not capture the full complexity of the underlying 3D function: Some features, such as the positive correlation between cos i and M_c (left main panel of Fig. 6) and the very high correlation between M_p and M_c (right main panel of that same figure), are captured only by the 2D or 3D PDFs. The latter correlation implies that continued timing, which will keep improving the precision of $\dot{\omega }$ (and thus of M), will result in much improved measurements of the individual masses.

The overall values for M_c and M_p derived from this self-consistent approach are slightly larger and more precise than those derived by the DDGR and DDK models, but about 1σ consistent with them. They are also 1-σ consistent with the masses derived from Shapiro delay alone. With regards to the orbital orientation, we see that a solution in Region 2 is preferred, with a total probability of 99.66%. The solution in Region 1 has a total probability of 0.34%, the difference between the two regions reaches a statistical significance close to 3σ.

Two conclusions can be derived from this. First, the Ω obtained from our Bayesian analysis of the timing data agrees well within 1σ with the VLBI estimate in Table 1. Second, the small amount of probability for the solution in Region 1 means that our timing yields a ∼3σ detection of the annual orbital parallax.

The PA of the orbital angular momentum, ψ_b ≡ Ω + 90 deg =  280.4(57) deg, means that the orbital angular momentum points nearly westwards. This is nearly opposite to the (mostly eastwards) PA of the proper motion, Θ_μ (see Table 1).

We now discuss the alignment of the pulsar spin axis with the orbital angular momentum. If they are aligned, then i = ζ = ∼84 deg (see Sect. 3). Although the timing value we measured for i is close to ζ, it is not consistent: the difference between them is 1.24 deg, which is outside the 99% uncertainty range for ζ (0.8 deg). Taken at face value, this small difference suggests a minor misalignment between the spin axis of the pulsar and the orbital angular momentum. However, before jumping into that conclusion, we reiterate the fact that RVM fit has important systematic issues, one of them being that there are obviously small-scale deviations from a perfect large-scale dipolar field (the grey points in Fig. 3).

Regarding our measurement of ψ_b, it is near one of the three possible values for the PA of the pulsar spin, ψ = −60(10) = 300(10) deg (see Sect. 3). The difference between them is −20(12) deg, where we have added their uncertainties in quadrature. This difference is not statistically significant, and consistent with a pulsar spin aligned with the orbital angular momentum.

7. Variation in the orbital period

The observed orbital period derivative obtained with the DDK model in Table 4 is Ṗ_b,obs = 0.2509(76) × 10⁻¹² s s⁻¹. This is ∼12 times more precise than the estimate made by Cognard et al. (2017).

We now discuss the implications of this measurement in more detail. First, we update the estimate of the contribution from Shklovskii effect. Using the distance and proper motion from our re-analysis of VLBI data, we obtain:

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b},\mathrm{Shk}}={\mu }^2\frac{d}{c}{P}_{\mathrm{b}}=0.2794(12)\times {10}^{-12}\mathrm{s}{\mathrm{s}}^{-1}.\end{array}$$(3)

The contribution of Galactic acceleration can be calculated with the analytical formulae provided by Damour & Taylor (1991), Nice & Taylor (1995), and Lazaridis et al. (2009),

$$\begin{array}{c}\hfill \frac{{\dot{P}}_{\mathrm{b},\mathrm{Gal}}}{P_{\mathrm{b}}}=-\frac{K_{\mathrm{z}}|sinb|}{c}-\frac{{\mathrm{\Theta }}_0^2}{R_{0}c}(cosl+\frac{\beta }{{\beta }^2+{sin}^{2}l})cosb,\end{array}$$(4)

where β ≡ (d/R₀) cos b − cos l. For Galactic height z ≡ |dsinb| ≤ 1.5 kpc, the vertical component of Galactic acceleration K_z can be approximated as (Holmberg & Flynn 2004; Lazaridis et al. 2009)

$$\begin{array}{c}\hfill K_z({10}^{-9}\mathrm{cm}{\mathrm{s}}^{-2})\simeq 2.27z_{\mathrm{kpc}}+3.68(1-e^{-4.3z_{\mathrm{kpc}}}).\end{array}$$(5)

In this calculation, we adopt the Galactic parameters in Gravity Collaboration (2021), where the distance from the Sun to the Galactic centre is R₀ = 8.275(34) kpc and the Galactic circular velocity at the location of the Sun is Θ₀ = 240.5(41) km s⁻¹¹¹. Assuming a 10% uncertainty in the vertical acceleration, we get

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b},\mathrm{Gal}}=-0.0142(13)\times {10}^{-12}\mathrm{s}{\mathrm{s}}^{-1}.\end{array}$$(6)

Subtracting these two terms from Ṗ_b,obs we obtain the ‘intrinsic’ variation in the orbital period:

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b},\mathrm{int}}={\dot{P}}_{\mathrm{b},\mathrm{obs}}-{\dot{P}}_{\mathrm{b},\mathrm{Shk}}-{\dot{P}}_{\mathrm{b},\mathrm{Gal}}=-0.0143(76)\times {10}^{-12}\mathrm{s}{\mathrm{s}}^{-1}.\end{array}$$(7)

An intrinsic Ṗ_b is expected originate from orbital decay of the system caused by the emission of gravitational waves, Ṗ_b,GR. Using the masses and orbital parameters derived from the DDGR model and the relation of Peters (1964), which provides the leading order estimate for the orbital decay caused by the emission of quadrupolar gravitational waves in GR, we obtain a slightly higher and more precise value than Cognard et al. (2017), Ṗ_b,GR = −0.00809(5) × 10⁻¹² s s⁻¹. This is 1-σ consistent with Ṗ_b,int and similar to its measurement precision.

Subtracting Ṗ_b,GR from Ṗ_b,int, we obtain the excess in the observed Ṗ_b:

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b},\mathrm{xs}}={\dot{P}}_{\mathrm{b},\mathrm{int}}-{\dot{P}}_{\mathrm{b},\mathrm{GR}}=-0.0063(76)\times {10}^{-12}\mathrm{s}{\mathrm{s}}^{-1},\end{array}$$(8)

which agrees well with zero. This limits any additional effects beyond GR, such as a variation in Newton’s gravitational constant or the emission of DGWs predicted by some alternative theories of gravity. For instance, following the calculations for scalar-tensor theories of Sect. 5 of Cognard et al. (2017), we find |α_p − α₀|< 0.005 (95% C.L.), which is a significant improvement compared to their Eq. (6), and comparable to the limits of Freire et al. (2012) and Antoniadis et al. (2013). It should be noted that we know α₀ < 0.003 (95% C.L.) from Solar system experiments (Bertotti et al. 2003; Esposito-Farèse 2006).

This tight limit on dipolar radiation, in combination with the large and well-determined mass of PSR J2222−0137, makes this system an ideal laboratory for certain non-linear aspects of strong-field gravity, like spontaneous scalarisation (Damour & Esposito-Farèse 1993; Shao et al. 2017). This will be explored in detail in a forthcoming publication (Zhao et al., in prep.).

We now discuss how robust this estimate is. First, as we can see from Table 3, this value depends on the DM model and the assumptions we make relative to the proper motion and position. The difference of ∼0.01 × 10⁻¹² s s⁻¹ is comparable to the uncertainty of Ṗ_b,obs. As shown in Table 4, the choice of orbital model has a smaller influence: the difference between the Ṗ_b,obs obtained with the model-independent ELL1H+ and DDK models is 0.004 × 10⁻¹² s s⁻¹, which is smaller than the 1σ uncertainty for Ṗ_b,obs.

Another possible source of uncertainty is the model used to estimate Ṗ_b,Gal. We now compare the predictions of different models of the gravitational field of the Galaxy, following the outline of the analysis of Zhu et al. (2019) for PSR J1713+0747; these are summarised in Table 7. We find that the differences between the predictions of Ṗ_b,Gal for different models are smaller than the current uncertainty of Ṗ_b,obs, but are larger than the estimated uncertainties of Ṗ_b,Gal according to each model.

Finally, we note that the solar height z₀ is ignored in the estimates made in Table 7. In Fig. 8, we show the variation in Ṗ_b,Gal as a function of z₀ for the different Galactic models in Table 7. The differences of ∼10⁻¹⁵ s s⁻¹ are also significantly smaller than the uncertainty of Ṗ_b,obs, but comparable to the differences between models. As an example, if we use z₀ = 20.8 pc (Bennett & Bovy 2019), the Ṗ_b,Gal predicted by the analytical model becomes −0.0151(14)×10⁻¹² s s⁻¹, a difference similar to the 10% uncertainty in the vertical acceleration of that model.

Fig. 8. Variation in Ṗb,Gal with the Galactic height of the Sun (z0), for the Galactic potential models listed in Table 7. The dashed vertical line on the right corresponds to the estimate of Bennett & Bovy (2019).

For now, none of these differences change the fact that the Ṗ_b,xs is 1σ consistent with zero. However, as the measurement of Ṗ_b,obs improves, these uncertainties in the Galactic model and z₀ will eventually limit the precision of Ṗ_b,int and Ṗ_b,xs.

8. Summary and perspectives

In this paper we present the results of a 12-year timing of PSR J2222−0137, combining data from the Effelsberg, Nançay, and Lovell radio telescopes with early GBT data. Furthermore, we have re-analysed the astrometric VLBI data. Finally, we have also obtained polarimetric data from FAST.

The re-analysis of the VLBI data confirms most of the results presented by Deller et al. (2013), except for Ω, which changed by ∼180 deg. This resulted from our use of conventions that are fully consistent with those used in pulsar timing. We have also calculated the absolute position, with more realistic uncertainty estimates. Because of these, there is no longer a significant disagreement with the timing position.

The very high signal-to-noise ratio of the FAST data yields polarimetry consistent with the (corrected) Nançay data, and it has allowed a detection of several faint emission regions, which include, importantly, an interpulse. This has allowed an unambiguous determination of the geometry of the pulsar, in particular a precise determination of its 3D orientation.

Regarding the timing, one of the most important things we have learned from this system is the great importance of consistent spin phase definitions for all the templates used to derive ToAs from the different datasets. Without this, we have no consistent measurements of the orbital motion of the pulsar. Fixing this issue has resulted in a very significant improvement in the quality of the timing relative to previous analyses.

If we use a few DM derivatives to model the DM variations, the proper motion shows discrepancies from the VLBI result at the 3 − 5σ level. If we use instead the DMX model, which can describe short-term DM changes, then the proper motion is consistent with the VLBI result, but with much larger uncertainties. Because of this, in our timing we used the DMX model and fixed the parallax and proper motion to the VLBI values.

Relative to previous work, our improved timing analysis results in a much more precise measurement of three PK parameters, two for the Shapiro delay (h₃ and ς) and one for the rate of advance of periastron, $\dot{\omega }$. The mutual agreement between the mass estimates obtained with these parameters within the framework of GR provides a successful ∼1% test of that theory.

The secular variation in the semi-major axis, ẋ, is larger (in magnitude) than expected; the difference with the expected value is 3.4σ significant. It is likely that this is caused by the presence of effects, such as the annual orbital parallax, which are not taken into account in the models that fit explicitly for ẋ. Indeed, a DDK model, which takes the annual orbital parallax into account, provides the best fit (with the lowest χ²) to the data assuming only the changes in x expected from the geometry of the system.

From a self-consistent analysis that assumes the validity of GR and takes all kinematic effects into account, we obtain a large pulsar mass of 1.831(10) M_⊙ and a companion WD mass of 1.319(4) M_⊙, which is the largest confirmed NS birth mass (Cognard et al. 2017). This is only one of two recycled pulsar and massive (> 0.6 M_⊙) WD binary systems with precisely measured masses, the other being PSR J2045+3633 (McKee et al. 2020). The total mass of the system is 3.150(14) M_⊙, confirming this as the most massive double degenerate binary known in the Galaxy. The resulting orbital orientation, which favours an inclination angle of 85.27 deg and Ω = 188 deg, is fully consistent with the VLBI astrometry. It is also consistent with the orientation of the pulsar spin derived from polarimetry, showing that, within experimental precision, the spin axis of the pulsar is aligned with the orbital angular momentum.

The relatively long spin period of PSR J2222−0137 means that not too much angular momentum was transferred in this case; thus, in principle, there could be a measurable misalignment. However, taking into account the characteristics of its current companion, we come to the conclusion that the pulsar was already ∼50 Myr old when the Roche-lobe overflow started; at that time it was likely much slower than it is now. This would imply that most of the spin seen today did originate from the recycling process. Thus, the observed alignment between the pulsar spin and the orbital angular momentum is to be expected.

We have also obtained a very precise measurement of the variation in the orbital period, Ṗ_b = 0.2509(76) × 10⁻¹² s s⁻¹. After subtracting the precise estimates for the kinematic effects, we find an intrinsic variation in the orbital period (Ṗ_b,int) that is consistent with the orbital decay caused by the emission of quadrupolar gravitational waves predicted by GR (Ṗ_b,GR). Subtracting Ṗ_b,GR from Ṗ_b,int we obtain an excess orbital decay, Ṗ_b,xs = −0.0063(76) × 10⁻¹² s s⁻¹, that is consistent with zero. This represents an important constraint on alternative theories of gravity, particularly since the mass of PSR J2222−0137 falls into a range that so far is poorly constrained in terms of phenomena such as spontaneous scalarisation (Shao & Wex 2016).

This system also has the potential to improve the constraints on the variation in Newton’s gravitational constant, G. Those constraints are proportional to the precision of Ṗ_b,xs/P_b. For PSR J1713+0747, this number is ∼2.5 × 10⁻²⁰ (Zhu et al. 2019), and for PSR J2222−0137 the number is ∼3.6 × 10⁻²⁰; these two values are comparable with each other despite the much shorter timing baseline for the latter pulsar. Indeed, the precision of Ṗ_b,xs is currently limited by the measurement of Ṗ_b,obs, and this decreases rapidly with time ($T^{-\frac{5}{2}}$). As we have seen, this might soon be limited by uncertainties in the Galactic gravitational potential; however, these are also expected to improve with the dynamical data provided by the Gaia mission.

The limits on alternative theories of gravity and on the variation in G that result from these measurements and the future prospects for improved measurements of the limits to be derived from this system will be discussed in greater detail in future publications.

Acknowledgments

We thank Thomas Tauris for discussions on the evolution of the WD companion of PSR J2222−0137, and Vivek Venkatraman Krishnan for Fig. 1. This study is partly based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from “Programme National de Cosmologie et Galaxies” (PNCG) of CNRS/INSU, France. D.L.K was supported by NSF Physics Frontiers Center award number 1430284. J.W. McKee is a CITA Postdoctoral Fellow: This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference #CITA 490888-16]. W.Z. is supported by the CAS-MPG LEGACY project, the National Key R&D Program of China No. 2017YFA0402600, the National SKA Program of China No. 2020SKA0120200 and the National Natural Science Foundation of China No. 11873067, No. 12041303.

References

All Tables

All Figures

Fig. 1. Geometric parameters to be used throughout the paper according to the ‘observer’s convention’, which is used throughout the text. For the description of the fundamental reference frame (axes in yellow) and the orbital reference frame (axes in blue), see Sect. 2. For the orientation of the pulsar angular momentum, S, and its magnetic axis, μ, relative to this frame, see the detailed explanation in Sect. 3. In this figure, S is represented as being parallel to the direction of the orbital angular momentum, k. This is true for binaries such as PSR J2222−0137, where the pulsar was fully recycled by mass accretion from the companion and the companion then evolved into a WD. If the second-formed object in the system were a NS, then the kick from the supernova event that formed it might substantially change the orientation of the orbit (and thus the orientation of k) relative to S and greatly increase its eccentricity. The figure is by Vivek Venkatraman Krishnan.

In the text

Fig. 2. Correlation plots for the astrometric parameters measured by VLBI listed in Table 1. All parameters have well-defined uncertainties, except for i: Systems with inclinations of 85.25 and 94.75 deg will have nearly identical motions in the plane of the sky. For this reason, the orbital inclinations were sampled for two intervals that correspond to the best estimate of sin i from timing: 85.0 − 85.5 deg and 94.5 − 95.0 deg. We present the offsets (Δ) relative to our best values of α′ and δ′ in the first column of Table 1.

In the text

Fig. 3. Polarisation profile as observed with FAST at a central frequency of 1250 MHz, averaged over a bandwidth of 420 MHz. Top panel: total intensity (black), the linear intensity (red), and circular polarisation (blue) as a function of longitude. The inset shows an enlarged version of the main pulse. Middle panel: low intensity level of the profile across the full rotational period, revealing the existence of an interpulse that is separated from the main pulse by about half a period. Bottom panel: values of the PA of the linear emission as a function of the longitude. An RVM has been fitted to the black PA values. The result is shown as a solid black line, while the dashed line indicates the RVM solution separated by 90 deg.

In the text

Fig. 4. Twelve-year timing data of PSR J2222−0137, for the ELL1H+ timing model where we used the DMX model and fit for position and proper motion (middle column in Table 3). Top panel: DM offsets relative to the reference DM (3.2805 cm−3 pc) as a function of date. Middle and bottom panels: timing residuals as a function of date (middle) and orbital phase (bottom, measured from ascending node). The residuals are displayed with different offsets for each instrument: EFF – Effelsberg; JB – Lovell telescope (Jodrell Bank); NCY-S – Nançay at the S band; NCY-L – Nançay at the L band; GBT-820 – GBT at 800 MHz; and GBT-1500 – GBT at 1500 MHz.

In the text

Fig. 5. Residuals as a function of the orbital phase for the same ephemeris as Fig. 4. Top panel: we display the residuals predicted by the model but without taking into account the Shapiro delay. Lower panel: we see the same while fitting for other Keplerian parameters. As we see, some of the Shapiro delay is absorbed by this fit, but some still remains; this remnant represents the ‘measurable’ part of the Shapiro delay; this is given, in the limit of perfect orbital sampling, by Eq. (31) of Freire & Wex (2010). In this plot, the Nançay residuals are displayed in black, the Effelsberg residuals in red, and all others in grey, all without error bars for the sake of clarity.

In the text

Fig. 6. Mass-mass diagram of J2222−0137. In the main panel on the left we display the cos i − Mc plane, and on the right we display the Mp − Mc plane; the grey region is excluded by the constraint sin i ≤ 1. The black contours indicate the 68.3% and 95.4% confidence region derived from a 3D χ2 map of the cos i − Ω − M plane using the DDK orbital model combined with GR equations. The constraints (according to GR) from h3, ς, and $\dot{\omega }$ in the ELL1H+ model are shown in blue, green, and red lines, respectively, where the solid and dotted lines indicate the nominal and ±1σ measurements. The fact that they all meet in the same regions of these planes represents a ∼1% test of GR, which the theory passes. The side panels display the 1D PDFs for cos i (top left), Mp (top right), and Mc (right), with the median value and the 1σ and 2σ confidence intervals indicated as vertical lines.

In the text

Fig. 7. Orbital orientation constraints for PSR J2222−0137. In the main panel we display the cos i − Ω plane, where randomly aligned orbits will have a priori constant probability density. The black contours show the 68.3% and 95.4% confidence region derived from a 3D χ2 map of cos i − Ω − M space using the DDK model with the additional assumption that GR is the correct theory of gravity. The solid red lines indicate the regions that are consistent with the nominal and ±1σ measurements of ẋ, and the dashed vertical lines indicate the position of ς obtained in the ELL1H+ model. The dashed horizontal lines show the value of Ω from VLBI observations. The side panels display the 1D PDFs for cos i (top) and Ω (right). The solution with positive cos i has 99.66% of the total probability.

In the text

	Fig. 8. Variation in Ṗb,Gal with the Galactic height of the Sun (z0), for the Galactic potential models listed in Table 7. The dashed vertical line on the right corresponds to the estimate of Bennett & Bovy (2019).
In the text