Issue 
A&A
Volume 654, October 2021



Article Number  A16  
Number of page(s)  17  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202141450  
Published online  01 October 2021 
PSR J2222−0137
I. Improved physical parameters for the system
^{1}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
email: yjguo@mpifrbonn.mpg.de
^{2}
Laboratoire de Physique et Chimie de l’Environnement et de l’Espace, Université d’Orléans/CNRS, 45071 Orléans Cedex 02, France
^{3}
Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, 18330 Nançay, France
^{4}
CAS Key Laboratory of FAST, NAOC, Chinese Academy of Sciences, Beijing 100101, PR China
^{5}
Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
^{6}
Centre for Astrophysics and Supercomputing, Swinburne University of Technology John St, Hawthorn, VIC 3122, Australia
^{7}
ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Hawthorn, Australia
^{8}
Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of WisconsinMilwaukee, PO Box 413 Milwaukee, WI 53201, USA
^{9}
Jodrell Bank Center for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK
^{10}
School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, PR China
^{11}
National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA
^{12}
LUTH, Observatoire de Paris, PSL Research University, CNRS, Université Paris Diderot, Sorbonne Paris Cité, 92195 Meudon, France
Received:
1
June
2021
Accepted:
17
July
2021
Context. The PSR J2222−0137 binary system has a set of features that make it a unique laboratory for tests of gravity theories.
Aims. To fully exploit the system’s potential for these tests, we aim to improve the measurements of its physical parameters, spin and orbital orientation, and postKeplerian parameters, which quantify the observed relativistic effects.
Methods. We describe an improved analysis of archival very long baseline interferometry (VLBI) data, which uses a coordinate convention in full agreement with that used in timing. We have also obtained much improved polarimetry of the pulsar with the Five hundred meter Aperture Spherical Telescope (FAST). We provide an improved analysis of significantly extended timing datasets taken with the Effelsberg, Nançay, and Lovell radio telescopes; this also includes previous timing data from the Green Bank Telescope.
Results. From the VLBI analysis, we have obtained a new estimate of the position angle of the ascending node, Ω = 189_{−18}^{+19} deg (all uncertainties are 68% confidence limits), and a new reference position for the pulsar with an improved and more conservative uncertainty estimate. The FAST polarimetric results, and in particular the detection of an interpulse, yield much improved estimates for the spin geometry of the pulsar, in particular an inclination of the spin axis of the pulsar of ∼84 deg. From the timing, we obtain a new ∼1% test of general relativity (GR) from the agreement of the Shapiro delay parameters and the rate of advance of periastron. Assuming GR in a selfconsistent analysis of all effects, we obtain much improved masses: 1.831(10) M_{⊙} for the pulsar and 1.319(4) M_{⊙} for the white dwarf companion; the total mass, 3.150(14) M_{⊙}, confirms this as the most massive double degenerate binary known in the Galaxy. This analysis also yields the orbital orientation; in particular, the orbital inclination is 85.27(4) deg – indicating a close alignment between the spin of the pulsar and the orbital angular momentum – and Ω = 187.7(5.7) deg, which matches our new VLBI estimate. Finally, the timing also yields a precise measurement of the variation in the orbital period, Ṗ_{b} = 0.251(8) × 10^{−12} ss^{−1}; this is consistent with the expected variation in the Doppler factor plus the orbital decay caused by the emission of gravitational waves predicted by GR. This agreement introduces stringent constraints on the emission of dipolar gravitational waves.
Key words: binaries: close / gravitational waves / pulsars: general / pulsars: individual: J2222−0137 / stars: neutron / white dwarfs
© Y. J. Guo et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
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 driftscan 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 semimajor axis of the pulsar’s orbit (x)^{1} of 10.848 lightseconds (lts).
The small spin period derivative (5.8 × 10^{−20}) 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^{−3}, 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 edgeon 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^{2} 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^{3}. 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 () 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 wellmeasured 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 & EspositoFarè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 nonlinear phenomena, such as spontaneous scalarisation (Damour & EspositoFarè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 lowprecision 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 followup 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 (), and, especially, the derivative of the orbital period, Ṗ_{b}.
Finally, we have observed the pulsar with the central beam of the 19beam 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 longterm timing project, which is to improve the precision of the physical parameters of the system, especially the postKeplerian (PK) parameters (which quantify, in a theoryindependent way, the observed relativistic effects, such as the aforementioned Shapiro delay, , and Ṗ_{b}). The ultimate objectives of the project – improved constraints on the nature of gravitational radiation and the behaviour of gravity for strongly selfgravitating 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 selfconsistent 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. Reanalysis of the VLBI astrometry
2.1. Reference frame
Before describing the reanalysis 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 DamourDeruelleKopeikin (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 secondformed 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 northpointing axis and then increasing anticlockwise, 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. Reanalysis 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 reanalysis 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^{4}, 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 inbeam 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 welldefined 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. 
Barycentric astrometric parameters derived from our reanalysis of the VLBI data on PSR J2222−0137, compared to the orbitalmotioncorrected parameters estimated by Deller et al. (2013) and the difference between the former and the latter.
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 inbeam 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 outofbeam calibrator, is updated to the most recent estimate^{5}, 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 frequencydependent 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 perepoch positions of J2222−0132 obtained via phase referencing to J2218−0335 without selfcalibration, 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 perepoch 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 inbeam astrometry and absolute astrometry for this system. Future multifrequency 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 inbeam astrometry.
In all our subsequent analysis, when referring to the VLBI astrometry, we use the rederived 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 highsensitivity (S/N = 23 000), highresolution 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 lowintensity 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 largescale 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 welldefined PA values provides important leverage and helps to break the usually existing covariances 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 Φ_{0} (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 Φ_{0}, ψ_{0}. 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 ψ_{0} = 40(1) deg, Φ_{0} = 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 ψ_{0} 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 derotated PA . Because of the aforementioned OPMs, 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 19beam system by observing some bright millisecond pulsars and solving for the full Muller matrix based on their known polarisation profiles. We found the crosscoupling 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 19beam 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 19beam 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^{6} (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 Lband ToAs only, as these are the most extensive and useful datasets; the addition of the smaller P and Sband 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.
Observations of J2222−0137 and data reduction parameters.
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 χ^{2} of 1.0, see Table 2). However, it helped achieve a reduction of the rootmeansquare (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 piecewiseconstant 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 longterm 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^{7}, 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 theorydependent 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 theoryindependent 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 selfconsistent Bayesian analysis of the system, which also assumes the validity of GR.
The third, ELL1H+, is a theoryindependent model based on the loweccentricity 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_{0}), 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 LaplaceLagrange parameters ϵ_{1} = esinω and ϵ_{2} = 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^{2} (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^{3} 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_{0} and ω observed in the DDlike models. Also, by redefining P_{b} as the time between passages through ascending node, we avoid its strong correlation with seen in the DDlike 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_{0} 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_{0} and P_{b} in Table 4, also stated to their own uncertainties, are not precise enough for this purpose^{8}.
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}.
Parameters for the PSR J2222−0137 binary system.
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 χ^{2} from the first to the second TEMPO iterations  all stayed at consistent values close to that of our best solution, χ^{2} ∼ 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 (frequencydependent) parameters, which describe nondispersive 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 selfconsistent Bayesian grid obtained with the DDK model, all derived assuming the VLBI proper motion.
Orbital parameters for PSR J2222−0137.
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 χ^{2} 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 χ^{2} of 1 for each dataset.
Fig. 4. Twelveyear 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); NCYS – Nançay at the S band; NCYL – Nançay at the L band; GBT820 – GBT at 800 MHz; and GBT1500 – 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, 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 χ^{2} of 10799.74. The proper motion is μ_{α} = 44.61(3) mas yr^{−1} and μ_{δ} = −5.26(5) mas yr^{−1}. Compared to the VLBI proper motion in the third column, the differences are Δμ_{α} = −0.09(5) mas yr^{−1} and Δμ_{δ} = 0.43(9) mas yr^{−1}, 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 shortterm 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 χ^{2}, 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^{−1} and Δμ_{δ} = 0.23(14) mas yr^{−1} (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 χ^{2} 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_{3} 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 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 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 χ^{2} 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 χ^{2} 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 . This is ∼6 times more precise than the measurement obtained by Cognard et al. (2017), ; 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 . Thus, our new measurement agrees with within 1σ. Therefore, the same applies to the total mass of the system derived from both methods: assuming GR, we obtain from 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. Massmass diagram of J2222−0137. In the main panel on the left we display the cos i − M_{c} plane, and on the right we display the M_{p} − M_{c} 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 h_{3}, ς, and 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), M_{p} (top right), and M_{c} (right), with the median value and the 1σ and 2σ confidence intervals indicated as vertical lines. 
This statement relies on the fact that the is relativistic. As discussed by Cognard et al. (2017), the largest additional contribution to , which is caused by the proper motion of the system, is of the order of a few times 10^{−16} deg yr^{−1} (i.e., about 100 times smaller than the current measurement uncertainty). Therefore, the assumption that is relativistic is fully warranted. This also means that, in the near future, improved measurements of 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 with the Shapiro delay and obtain selfconsistent 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 χ^{2} for that fit is nearly identical to that of the best ELL1H+ model, indicating again the selfconsistency of the relativistic effects.
5.4. Variation in the projected semimajor 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 semimajor axis, ẋ: When they assume the VLBI proper motion, they obtain ẋ = 3.5(30) × 10^{−15} lts s^{−1}, our new value is −7.76(48)×10^{−15} lts s^{−1}; 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):
Assuming the values of Ω and Θ_{μ} in Table 1 and of x, μ and i in Table 3, we obtain ẋ_{k} = ±6.12 × 10^{−15} lts s^{−1}, 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^{−15} lts s^{−1}, 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, , 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 spinorbit coupling, and ẋ_{SO}, the change of the orbital inclination also due to the spinorbit coupling. The spinorbit 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 spindown 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.
Contributions to ẋ of PSR J2222−0137.
The change in the orbital inclination due to spinorbit 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 LenseThirring 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^{9}), ẋ_{SOc} is also dominated by the LenseThirring effect. Both contributions are negligible^{10}. We want to add that, due to tidal torques during the Rochelobe overflow episode, one would expect an alignment of the spin axis, resulting in θ_{c} ≃ 0, further suppressing any spinorbit 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 nonlinear 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 , see detailed discussion by Ridolfi et al. (2019). From their Eq. (25), we obtain:
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 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 χ^{2} 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.
Details for the grid regions.
Furthermore, as we can see in Table 4, a DDK model based on the best i, Ω combination has a lower χ^{2} 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 selfconsistent 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 selfconsistent χ^{2} 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 (, γ, 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 χ^{2}, 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 χ^{2} 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 iM_{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 (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 selfconsistent 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 smallscale deviations from a perfect largescale 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^{−12} s s^{−1}. 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 reanalysis of VLBI data, we obtain:
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),
where β ≡ (d/R_{0}) 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)
In this calculation, we adopt the Galactic parameters in Gravity Collaboration (2021), where the distance from the Sun to the Galactic centre is R_{0} = 8.275(34) kpc and the Galactic circular velocity at the location of the Sun is Θ_{0} = 240.5(41) km s^{−1}^{11}. Assuming a 10% uncertainty in the vertical acceleration, we get
Subtracting these two terms from Ṗ_{b,obs} we obtain the ‘intrinsic’ variation in the orbital period:
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^{−12} s s^{−1}. 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}:
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 scalartensor theories of Sect. 5 of Cognard et al. (2017), we find α_{p} − α_{0}< 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} < 0.003 (95% C.L.) from Solar system experiments (Bertotti et al. 2003; EspositoFarèse 2006).
This tight limit on dipolar radiation, in combination with the large and welldetermined mass of PSR J2222−0137, makes this system an ideal laboratory for certain nonlinear aspects of strongfield gravity, like spontaneous scalarisation (Damour & EspositoFarè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^{−12} s s^{−1} 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 modelindependent ELL1H+ and DDK models is 0.004 × 10^{−12} s s^{−1}, 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.
Different contributions to Ṗ_{b}, in units of 10^{−12} s s^{−1}.
Finally, we note that the solar height z_{0} is ignored in the estimates made in Table 7. In Fig. 8, we show the variation in Ṗ_{b,Gal} as a function of z_{0} for the different Galactic models in Table 7. The differences of ∼10^{−15} s s^{−1} are also significantly smaller than the uncertainty of Ṗ_{b,obs}, but comparable to the differences between models. As an example, if we use z_{0} = 20.8 pc (Bennett & Bovy 2019), the Ṗ_{b,Gal} predicted by the analytical model becomes −0.0151(14)×10^{−12} s s^{−1}, 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 (z_{0}), 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_{0} will eventually limit the precision of Ṗ_{b,int} and Ṗ_{b,xs}.
8. Summary and perspectives
In this paper we present the results of a 12year timing of PSR J2222−0137, combining data from the Effelsberg, Nançay, and Lovell radio telescopes with early GBT data. Furthermore, we have reanalysed the astrometric VLBI data. Finally, we have also obtained polarimetric data from FAST.
The reanalysis 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 signaltonoise 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 shortterm 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_{3} and ς) and one for the rate of advance of periastron, . 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 semimajor 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 χ^{2}) to the data assuming only the changes in x expected from the geometry of the system.
From a selfconsistent 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 Rochelobe 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^{−12} s s^{−1}. 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^{−12} s s^{−1}, 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^{−20} (Zhu et al. 2019), and for PSR J2222−0137 the number is ∼3.6 × 10^{−20}; 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 (). 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.
This is an exact quantity defined in SI units as , which is similar to the precisely known product of Newton’s gravitational constant G and the mass of the Sun (Prša et al. 2016).
If GW190425 were a double NS system before its merger, then its mass was likely larger than that of PSR J2222−0137, about 3.4 M_{⊙}, (see Abbott et al. 2020).
Unlike in the J1141−6545 system (Venkatraman Krishnan et al. 2020), there was no substantial mass transfer to the WD that could have led to a significant spinup (Cognard et al. 2017). On the contrary, the tidal torques during the Rochelobe overflow phase most likely led to a (near) synchronisation with the orbital period (Tauris, priv. comm.).
Relevant numbers for the slowly rotating (ONeMg) WD companion were taken from Boshkayev et al. (2017).
Θ_{0} has been calculated based on the new R_{0} from Gravity Collaboration (2021) and the new proper motion measurements for Sgr A^{*} in Reid & Brunthaler (2020). For V_{⊙} we adopted the value used by Gravity Collaboration (2019). An analysis of the Solar motion with respect to nearby stars based on the Gaia Early Data Release 3 suggests a somewhat lower value for V_{⊙} (Gaia Collaboration 2021), which however we did not account for, since that 4 km s^{−1} shift is irrelevant for our analysis.
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 100m telescope of the MPIfR (MaxPlanckInstitut 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 49088816]. W.Z. is supported by the CASMPG 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
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, ApJ, 892, L3 [Google Scholar]
 Antoniadis, J., van Kerkwijk, M. H., Koester, D., et al. 2012, MNRAS, 423, 3316 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448 [Google Scholar]
 Arzoumanian, Z., Joshi, K., Rasio, F. A., & Thorsett, S. E. 1996, in IAU Colloq. 160: Pulsars: Problems and Progress, eds. S. Johnston, M. A. Walker, & M. Bailes, ASP Conf. Ser., 105, 525 [Google Scholar]
 Barker, B. M., & O’Connell, R. F. 1975, Phys. Rev. D, 12, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
 Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton University Press) [Google Scholar]
 Boshkayev, K., Quevedo, H., & Zhami, B. 2017, MNRAS, 464, 4349 [NASA ADS] [CrossRef] [Google Scholar]
 Boyles, J., Lynch, R. S., Ransom, S. M., et al. 2013, ApJ, 763, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Cognard, I., Freire, P. C. C., Guillemot, L., et al. 2017, ApJ, 844, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T., & Deruelle, N. 1986, Ann. Inst. Henri Poincaré Phys. Théor, 44, 263 [Google Scholar]
 Damour, T., & EspositoFarèse, G. 1992, Class. Quantum Gravity, 9, 2093 [CrossRef] [Google Scholar]
 Damour, T., & EspositoFarèse, G. 1993, Phys. Rev. Lett., 70, 2220 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T., & Taylor, J. H. 1991, ApJ, 366, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840 [NASA ADS] [CrossRef] [Google Scholar]
 Deller, A. T., Boyles, J., Lorimer, D. R., et al. 2013, ApJ, 770, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., et al. 2013, ApJ, 762, 94 [Google Scholar]
 Ding, H., Deller, A. T., Freire, P., et al. 2020, ApJ, 896, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Eardley, D. M. 1975, ApJ, 196, L59 [NASA ADS] [CrossRef] [Google Scholar]
 EspositoFarèse, G. 2006, The Tenth Marcel Grossmann Meeting. On recent developments in theoretical and experimental general relativity, gravitation and relativistic field theories, 647 [Google Scholar]
 Freire, P. C. C., & Wex, N. 2010, MNRAS, 409, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Freire, P. C. C., Wex, N., EspositoFarèse, G., et al. 2012, MNRAS, 423, 3328 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Smart, R. L., et al.) 2021, A&A, 649, A6 [EDP Sciences] [Google Scholar]
 Gérard, J. M., & Wiaux, Y. 2002, Phys. Rev. D, 66, 024040 [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2021, A&A, 647, A59 [CrossRef] [EDP Sciences] [Google Scholar]
 Holmberg, J., & Flynn, C. 2004, MNRAS, 352, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302 [Google Scholar]
 Johnston, S., & Kramer, M. 2019, MNRAS, 490, 4565 [NASA ADS] [CrossRef] [Google Scholar]
 Kaplan, D. L., Boyles, J., Dunlap, B. H., et al. 2014, ApJ, 789, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Kopeikin, S. M. 1995, ApJ, 439, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Kopeikin, S. M. 1996, ApJ, 467, L93 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, M., Stairs, I. H., Venkatraman Krishnan, V., et al. 2021, MNRAS, 504, 2094 [CrossRef] [Google Scholar]
 Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Lazaridis, K., Wex, N., Jessner, A., et al. 2009, MNRAS, 400, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Lorimer, D. R., & Kramer, M. 2004, Handbook of Pulsar Astronomy (Cambridge University Press), 4 [Google Scholar]
 Lorimer, D. R., & Kramer, M. 2005, Handbook of Pulsar Astronomy (Cambridge, England: Cambridge University Press) [Google Scholar]
 McKee, J. W., Freire, P. C. C., Berezina, M., et al. 2020, MNRAS, 499, 4082 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, P. J. 2017, MNRAS, 465, 76 [Google Scholar]
 Nice, D. J., & Taylor, J. H. 1995, ApJ, 441, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, AJ, 161, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Peters, P. C. 1964, Phys. Rev., 136, 1224 [Google Scholar]
 Piffl, T., Binney, J., McMillan, P. J., et al. 2014, MNRAS, 445, 3133 [NASA ADS] [CrossRef] [Google Scholar]
 Prša, A., Harmanec, P., Torres, G., et al. 2016, AJ, 152, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Radhakrishnan, V., & Cooke, D. J. 1969, ApJ, 3, 225 [Google Scholar]
 Reid, M. J., & Brunthaler, A. 2020, ApJ, 892, 39 [Google Scholar]
 Ridolfi, A., Freire, P. C. C., Gupta, Y., & Ransom, S. M. 2019, MNRAS, 490, 3860 [NASA ADS] [CrossRef] [Google Scholar]
 Sepinsky, J. F., Willems, B., Kalogera, V., & Rasio, F. A. 2010, ApJ, 724, 546 [NASA ADS] [CrossRef] [Google Scholar]
 Shao, L., & Wex, N. 2016, Sci. China Phys. Mech. Astron., 59 [CrossRef] [Google Scholar]
 Shao, L., Sennett, N., Buonanno, A., Kramer, M., & Wex, N. 2017, Phys. Rev. X, 7, 041025 [Google Scholar]
 Shibata, M., Taniguchi, K., Okawa, H., & Buonanno, A. 2014, Phys. Rev. D, 89, 084005 [CrossRef] [Google Scholar]
 Sokolovsky, K. V., Kovalev, Y. Y., Pushkarev, A. B., & Lobanov, A. P. 2011, A&A, 532, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Splaver, E. M., Nice, D. J., Arzoumanian, Z., et al. 2002, ApJ, 581, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Stovall, K., Freire, P. C. C., Antoniadis, J., et al. 2019, ApJ, 870, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170 [NASA ADS] [CrossRef] [Google Scholar]
 van Straten, W. 2006, ApJ, 642, 1004 [NASA ADS] [CrossRef] [Google Scholar]
 van Straten, W., Manchester, R. N., Johnston, S., & Reynolds, J. E. 2010, PASA, 27, 104 [NASA ADS] [Google Scholar]
 Venkatraman Krishnan, V., Bailes, M., van Straten, W., et al. 2020, Science, 367, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Verbunt, F., & Phinney, E. S. 1995, A&A, 296, 709 [NASA ADS] [Google Scholar]
 Yao, J., Zhu, W., Manchester, R. N., et al. 2021, Nat. Astron., 5, 788 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, W. W., Desvignes, G., Wex, N., et al. 2019, MNRAS, 482, 3249 [CrossRef] [Google Scholar]
All Tables
Barycentric astrometric parameters derived from our reanalysis of the VLBI data on PSR J2222−0137, compared to the orbitalmotioncorrected parameters estimated by Deller et al. (2013) and the difference between the former and the latter.
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 secondformed 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 welldefined 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. Twelveyear 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); NCYS – Nançay at the S band; NCYL – Nançay at the L band; GBT820 – GBT at 800 MHz; and GBT1500 – 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. Massmass diagram of J2222−0137. In the main panel on the left we display the cos i − M_{c} plane, and on the right we display the M_{p} − M_{c} 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 h_{3}, ς, and 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), M_{p} (top right), and M_{c} (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 (z_{0}), 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.