PSR J2222--0137. I. Improved physical parameters for the system

The PSR J2222-0137 binary system is a unique laboratory for testing gravity theories. To fully exploit its potential for the tests, we aim to improve the measurements of its physical parameters: spin, orbital orientation, and post-Keplerian parameters which quantify the observed relativistic effects. We present improved analysis of archival VLBI data, using a coordinate convention in full agreement with that used in timing. We also obtain much improved polarimetry with FAST. We provide an analysis of significantly extended timing data taken with Effelsberg, Nancay, Lovell and Green Bank telescopes. From VLBI analysis we obtain a new estimate of the position angle of ascending node, Omega=189(19) deg, and a new position of the pulsar with more conservative uncertainty. The FAST polarimetry and in particular the detection of an interpulse, yield much improved estimate for the spin geometry of the pulsar, in particular an inclination of the spin axis of 84 deg. From the timing we obtain a new 1% test of general relativity (GR) from the agreement of the Shapiro delay and the advance rate of periastron. Assuming GR in a self-consistent analysis of all effects, we obtain much improved mass: 1.831(10) M_sun for the pulsar and 1.319(4) M_sun for the companion; the total mass, 3.150(14) M_sun confirms it as the most massive double degenerate binary known in the Galaxy. This analysis also yields the orbital orientation: the orbital inclination is 85.27(4) deg, indicating a close alignment between the spin of the pulsar and the orbital angular momentum; Omega = 188(6) deg, matching our VLBI result. We also obtain precise value of the orbital period derivative, 0.251(8)e-12 s s^-1, consistent with the expected variation of Doppler factor plus the orbital decay caused by emission of gravitational wave (GW) predicted by GR. This agreement introduces stringent constraint on the emission of dipolar GW.


Introduction
PSR J2222−0137 is a pulsar with spin period (P) of 32.8 ms, discovered in the Green Bank Telescope (GBT) 350 MHz driftscan pulsar survey . 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 light-seconds (lt-s).
1. 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 and also precise values for the position and proper motion. 2. 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, Ω. 3. The edge-on orbit, the good timing precision and the large mass of the companion allow a highly significant detection of the Shapiro delay, which was originally detected by Kaplan et al. (2014). From this effect, Cognard et al. (2017) obtained m p = 1.76(6) M and m c = 1.293(25) M , where M represents the solar mass parameter 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. 4. 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 measurement of the masses of the components of the system, and therefore a test of GR. 5. The timing precision and the large x/P b ratio allow an unusually precise measurement of the variation of 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 general relativity (GR) and alternative gravity theories owe their precision to the wellmeasured masses. 6. 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 (DGW), in addition to the quadrupolar gravitational waves predicted by GR (Eardley 1975;Damour & Esposito-Farèse 1992;Gérard & Wiaux 2002). This could be detected in the measurement ofṖ b . 7. 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). First, because its mass estimates are more precise than for the latter systems (Antoniadis et al. 2012(Antoniadis et al. , 2013; this is important for the interpretation of the measurements. Second, previous authors (Shibata et al. 2014) have made it clear that these tests should be carried out for a variety of NS masses in order to exclude strongly non-linear phenomena like spontaneous scalarization (Damour & Esposito-Farèse 1993). Interestingly, the precise mass of PSR J2222−0137 places it in an intermediate, previous unexplored mass range. For this reason, even the relatively low-precision measurement ofṖ b by Cognard et al. (2017) has already provided useful constraints on alternative theories of gravity (Shao et al. 2017).
The data set 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 have been obtained in a set of orbital campaigns: To the two campaigns mentioned by Cognard et al. (2017) three more were added, which happened in 2018 January, 2019 August and 2020 October. The last observation used in this work was taken in 2021 May.
In addition, most of the pulse times of arrival (ToAs) derived from early discovery and follow up data with the GBT, which were used by Boyles et al. (2013) and Kaplan et al. (2014), have now been included in our analysis. The addition of these ToAs significantly extends our timing baseline to the past, which now starts in 2009 June 23 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 19-beam receiver of the Five hundred meter Aperture Spherical Telescope (FAST, Yao et al. 2021). These FAST data, taken at a frequency range between 1.0 and 1.5 GHz (40MHz at each edge of the band is excised in data reduction), provide the best polarimetric profile of PSR J2222−0137 to date, which will be discussed below. The FAST observations will, in the near future, contribute greatly to improved timing of this system.
In this work, we will address the proximate objective of this long-term timing project, which is to improve the precision of the physical parameters of the system, especially the post-Keplerian parameters (which quantify, in a theory-independent way, the observed relativistic effects, like the aforementioned Shapiro delay,ω andṖ b ). The ultimate objectives of the projectimproved constraints on the nature of gravitational radiation and the behaviour of gravity for strongly self-gravitating systemswill be addressed in subsequent work.
The remainder of the paper is structured as follows: In section 2, we will re-visit the Very Long Baseline Interferometry (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 section 3, we present the new results on the polarimetry of the pulsar from a high S/N detection with FAST. In Section 4, we describe the processing of the radio timing data. In Section 5, we present the main timing results, with a detailed review of the different timing parameters, and how they compare with previous estimates. In section 6, we make a self-consistent estimate of the component masses and orbital orientation of the system, and compare these with the VLBI results and the orientation of the pulsar obtained from the polarimetry. In section 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 section 8, we summarize our results and briefly point to further work on the implications of these timing results. . 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 section 2. For the orientation of the pulsar angular momentum, S, and its magnetic axis, µ relative to this frame, see detailed explanation in section 3. In this Figure, S is represented as being parallel to the direction of the orbital angular momentum k. This is true for binaries like PSR J2222−0137 where the pulsar was fully recycled by mass accretion from the companion and the companion then evolved to a WD. If the second-formed object in the system is a NS, then the kick from the supernova event that formed it might substantially change the orientation of the orbit (thus of k) relative to S, and greatly increase its eccentricity. Figure by Vivek Venkatraman Krishnan.

Reference frame
Before describing the re-analysis of the VLBI data, we must first define the geometric parameters used in this paper. We will use the "observer's convention", which is assumed for calculating all kinematic effects in the 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 North (as seen by the observer), the second points to the East -both of these define the plane of the sky, in yellowand the perpendicular direction towards the observer (n).
The orbital plane is indicated in blue. It crosses the plane of the sky in the line of nodes. The ascending node is one of the two points where the pulsar, in its orbital motion, crosses the plane of the sky, it is the one where its distance from the observer is increasing. Its opposite point is the descending node. The orbital reference frame is defined by the three blue axes: the line pointing to the ascending node i, a perpendicular direction within the orbital plane, j, which defines superior conjunction, and finally a direction perpendicular to the orbital plane, k; which is parallel to the orbital angular momentum (not represented in the Figure).
The full orientation of the orbital plane relative to the reference is defined by two angles. The first is Ω. This is a "position angle" (PA); these are measured in the plane of the sky, starting from the north-pointing axis and then increasing anti-clockwise, as seen from the observer. Ω is the PA of the ascending node (the PA of the descending node is given by Ω ± 180 deg). The second angle, the orbital inclination (i), is the angle between k and n. We can therefore see that, for i < 90 deg, the line of sight (los) component of the orbital angular momentum would point towards the Earth.

Re-analysis of the differential astrometry
In the analysis of Deller et al. (2013), it was assumed that the vector from the centre of mass to the pulsar (the pulsar's "position vector") follows the same convention as the angular momenta: its los component is positive when it points to us, the same sense as n in Fig. 1. This would imply, as described there, that orbital longitudes are measured from the descending node. However, and unlike the Damour & Taylor (1992) convention, the observer's convention is not coherent in this respect: in pulsar timing, the los component of the position vector -the geometric component of the quantity measured most directly in pulsar timing, the time delays of the radio pulses -is always positive if it points away from us. This means that the orbital longitudes are always measured from the ascending node. This prompted us to perform a re-analysis of the VLBI data, this time using a convention fully in agreement with the convention used in pulsar timing. We make use of the code binary_pulsar_MCMC.py 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 trialed (86.9 and 93.1 deg). We do not apply any constraints based on the value ofẋ measured by pulsar timing (described in subsequent sections).
The numerical values of the resulting parameters and their symbols, to be used in the remainder of the paper, are presented in Table 1 together with the earlier estimates, their difference, and the significance of the change; the position coordinates, α and δ , are calculated for the same epoch as Deller et al. (2013), and assuming the same position for the in-beam calibrator, FIRST J222201−013236 (hereafter J2222−0132); the results refer to the position of the pulsar as seen from the Solar System barycentre. Fig. 2 shows a correlation plot for these parameters.
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). Table 1. Barycentric astrometric parameters derived from our re-analysis of the VLBI data on PSR J2222−0137, compared to the orbital-motion corrected parameters estimated by Deller et al. (2013), and the difference between the former to the latter. PA means "position angle". The errors on the difference are those of our uncertainty estimate. The estimates of α and δ do not take into account the newly derived position for the in-beam calibrator (J2222−0132) and its related uncertainties; these are taken into account in the estimate of the absolute VLBI position, α and δ (see section 2.3).

Measured Parameters
This work 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.

Absolute VLBI position
The pulsar position provided by Deller et al. (2013) can be considered a relative position (α and δ ) with respect to the assumed position for the primary in-beam calibrator, J2222−0132. In order to compare with the timing position, an absolute VLBI position is required.
The absolute barycentric VLBI position of PSR J2222−0137 shown in Table 1 was estimated using the approach described by Ding et al. (2020). Our values of α and δ are different from the α and δ provided by Deller et al. (2013) for two reasons: 1) The absolute position of J2218−0335, the primary out-of-beam calibrator, is updated to the most recent estimate 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 frequency-dependent source structure ("core shift") of J2218−0335 between the higher frequencies at which its reference position is defined and the lower frequencies at which it is observed here. In Deller et al. (2013), the first term (with a value of ∼0.1 mas) was assumed to dominate, but we show here that this is not the case. Using the scatter in the per-epoch positions of J2222−0132 obtained via phase referencing to J2218−0335 without self-calibration, we conservatively determined that the second term contributes errors of 0.6 mas in right ascension and 2.3 mas in declination (using a weighted mean of these per-epoch positions would yield a smaller uncertainty on this mean absolute position, but one that would depend sensitively on the assumed input uncertainties). For the core shift contribution, we adopted 0.8 mas in each direction, which is the median core shift between 1.5 GHz and 8 GHz reported by Sokolovsky et al. (2011), and added it in quadrature with other uncertainties.
The uncertainties of α and δ are ∼50 times larger than the uncertainties of α and δ ; this quantifies the difference in precision of in-beam astrometry and absolute astrometry for this system. Future multi-frequency observations targeting J2222−0132 could considerably reduce the uncertainty in α and δ. Regardless, we note that the estimation of α and δ is independent from (and has no impact on) α , δ and other parameters in Table 1; which are derived purely from in-beam astrometry.
In all our subsequent analysis, when referring to the VLBI astrometry, we will use the re-derived values in Table 1. 5 http://astrogeo.org/sol/rfc/rfc_2021b/rfc_2021b_cat. html

Pulse Polarisation Results
We will now refer again to Figure 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 line of sight. 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 Figure 3, we present a high-sensitivity (S/N = 23,000), high-resolution pulse profile of PSR J2222−0137 obtained with FAST at a central frequency of 1250 MHz, the integration time is 1735 s. The observing setup and data reduction is similar to that described by Yao et al. (2021). This profile is shown in the PSR/IEEE convention (van Straten et al. 2010).
The quality of the profile allows not only the resolution of a number of features in the total, linear and circular intensity, but also reveals a number of faint profile components. Firstly, the main pulse shows polarised low-intensity components on both its leading and trailing side. Secondly, and most crucially, the sensitivity of FAST also uncovers a faint interpulse component separated from the main pulse by about half a period. The detection of these faint components is not a mere curiosity; it reveals the large-scale dipolar structure of the magnetic field of the pulsar, which is crucial for a determination of its geometry.
Inspecting the profile, the sense reversal of Stokes V (circular polarisation) underneath the main peak suggests that this longitude range can be identified with the location of the fiducial plane that is defined by the spin vector, the magnetic axis and the direction of the observer (see e.g., Lorimer & Kramer 2005). At the same time, the sudden drops in linearly polarised intensity, combined with the very rapid variation in V resolved in the FAST profile, reveal the existence of orthogonal polarisation modes (OPMs). These are clearly responsible for the large variation of the PA of the linear polarisation at these longitudes, which we mark in gray in Figure 3.
Overall, however, the PA swing under the main pulse shows a positive slope, which in terms of the "Rotating Vector Model" (RVM, Radhakrishnan & Cooke 1969) implies a negative value of β. This is confirmed by a blind fit of the RVM to the black PA values, which results in the solid black line shown in the bottom panel of Fig. 3. The grey PA points corresponding to the rapid changes in PA and V are ignored here, even though including those PA values does not change the overall result. Instead, the presence of the interpulse with well defined PA values provides important leverage and helps to break the usually existing 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 Figure 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, are following a 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).
Apart from ζ, the second parameter that defines the 3-D 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 ISM. Taking the rotation measure and its uncertainties into account, we obtain the de-rotated PA ψ 0 = 30(10) deg. Because of the aforementioned OPMs, ψ 0 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 have conducted tests on the center beam of the FAST 19 beam system by observing some bright MSPs and solving for the full Muller matrix based on their known polarisation profiles. We found that the cross-coupling between the two polarisations to be smaller than 1%. When tracking an object such as PSR J2222-0137, the FAST telescope measures and maintains the orientation of its 19 beam receiver at a fixed angle on the sky with a precision better than 0.1 degrees. Therefore, the systematic errors in the polarimetry of the 19-beam observing system should be negligible in our observations.
We now compare the polarisation profile from FAST with that presented by Cognard et al. (2017), which was taken with the Nançay Ultimate Pulsar Processing Instrument (NUPPI). There, we see the same modest degree of linear and circular polarisation, however, the PA swing, as well as 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 have confirmed the consistency and correctness of both results by smearing the FAST data to the MeerKAT resolution.

Processing of the timing data
The timing observations used for this project are summarized in Table 2. As in Cognard et al. (2017), the reduction of our timing data (mainly RFI 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 use the 820-MHz and the L-band ToAs only, as these are the most extensive and useful data sets; as the addition of the smaller P and S-band ToA sets does not change any parameters noticeably, only increasing the complexity of the analysis. In what follows, we describe mostly the improvements of the data analysis.

Derivation of the pulse times of arrival
As mentioned in Section 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 modeling 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 data sets 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 root-mean-square (rms) of the timing residuals (which are the ToA minus the model prediction for its rotation number): the current weighted residual rms of 2.8 µs is significantly better than the global weighted rms of 3.4 µs reported by Cognard et al. (2017).

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

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 BIMP2019 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 "DD" timing model of Damour & Deruelle (1986): 1. DDGR -a theory-dependent model that assumes the validity of GR and fits directly for the total mass of the system (M) and the companion mass (M c ). This model does not do a fully coherent analysis of the kinematic information. 2. DDK -a theory-independent model that takes into account the kinematic effects described by Kopeikin (1995Kopeikin ( , 1996, implemented in tempo by Ingrid H. Stairs. We use this particular model as the basis of a self-consistent Bayesian analysis of the system, which also assumes the validity of GR. 7 https://sourceforge.net/projects/tempo/ 3. ELL1H+ -a theory-independent 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 Laplace-Langrange parameters 1 = e sin ω and 2 = e cos ω. For the ELL1H+ model, we implemented the orthometric parameterization of the Shapiro delay using its exact expression, eq. (31) of Freire & Wex (2010).
The latest versions of this 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 DD-like models. Also, by re-defining P b as the time between passages through ascending node, we avoid its strong correlation withω seen in the DD-like models. Finally, by using the orthometric parameterisation, we avoid the large correlation between the Shapiro delay parameters r and s seen in the DD, DDK and ELL1+ models (−0.86 in our timing); the correlation between the orthometric parameters is −0.54.
This lack of correlations has practical advantages: as we will 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 .

Template alignment
Another important improvement in the data analysis was the use of pulse profile templates that have consistent phase definitions for all our data sets. In this way, the ToAs refer to a consistent longitude of the neutron star.
A&A proofs: manuscript no. main  (2) 10.8480213 (2) 10.8480213 (2)  Notes. Timing parameters derived using tempo when fitting for proper motion and DM derivatives (left column), fitting for proper motion and using DMX model (middle column), fixing proper motion at VLBI value and using the DMX model (right column). In all columns we use the VLBI parallax (see text). The α and δ offsets are calculated relative to the VLBI positions (see Table 1). The binary parameters are derived using the ELL1H+ orbital model. Numbers in parentheses represent 1σ uncertainties in the last digits. The reference epoch is MJD = 58000, and the position epoch is MJD = 55743 for consistency with Deller et al. (2013) and the analysis in section 2.
Normally, when combining data from different telescopes, the different data sets 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 data sets 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 data sets in subsequent iterations. These values should accurately characterize the different delays between the different observing systems.
However, if the templates have no consistent phase definitions, these ∆ts will be biased. Let us imagine two data sets 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 data sets 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 data sets is converted into a time offset for the second data set (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 data set.
This has implications for our measurement of the orbital motion. If we measure T asc for both data sets separately, without  (10) 10.8480235 (2 providing any time offsets (as in the first tempo iteration), we will 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 will 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 data sets 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 will 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 data sets, ∆t < δT asc .
Fixing this problem resulted in a major improvement in our timing of this pulsar. First, there was no longer a degradation by ∼ 1000 of the χ 2 from the first to the second tempo iterationsall 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 data sets (these are the EQUAD parameters in Table  1 of Cognard et al. 2017); this has contributed to the decrease in the residual rms mentioned above. Finally, the FD parameters, which describe non-dispersive variations in the times of arrival 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 pre-determined time offsets is fine and that, therefore, all our data sets 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).

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 section 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. 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 (EFAC) for the ToA uncertainties in order to achieve a reduced χ 2 of 1 for each data set.
In what follows, we will call the reader's attention to the more important timing parameters. We will also make some comparisons between our results and those presented by Cognard et al. (2017). Since the systematic issues discussed in section 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.

Astrometric Parameters
All positions reported in Table 3 Table 3). Top panel: DM offsets relative to the reference DM (3.2805 cm −3 pc) as a function of date. Middle and bottom panel: timing residuals as a function of date (middle) and orbital phase (bottom, measured from ascending node). The residuals are displayed with different offsets for each instrument, EFF -Effelsberg, JB -Lovell telescope (Jodrell Bank), NCY-S -Nançay at S-band, and NCY-L -Nançay at L-band, GBT-820 -GBT at 800 MHz, GBT-1500 -GBT at 1500 MHz. 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 sec-tion 2. Because the timing position is more precise, we will fit for it from now on.
First, as Cognard et al. (2017), we use a simple description of the variation of 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  Fig. 4. In the top panel we display the residuals predicted by the model, but without taking into account the Shapiro delay. In the 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; the latter 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 clarity. of 3.710(79) mas, which is in near perfect agreement with the VLBI parallax, = 3.730 +0.015 −0.016 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 will use VLBI parallax from now on.
We have found that there are DM variations on short timescales which could affect the astrometric parameters. In order to take these short-term DM variations into account, we use the DMX model, fitting for a DM offset for TOAs within gaps of 60 days, the resulting DM offsets are depicted in the top panel of Fig. 4. This interval was chosen and adhered to before a detailed consistency analysis of all the PK parameters. The results are presented in the second column of Table 3, where we still fit for proper motion. This causes a very significant decrease in the χ 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.

Shapiro Delay
The Shapiro delay was first detected by Kaplan et al. (2014), which used it to obtain M p = 1.20 (14)  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 Figure 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.
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 2009 June 28) 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 data set 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 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).

Advance of Periastron
From our timing we deriveω obs = 0.09607(48) deg yr −1 . This is ∼6 times more precise than the measurement obtained by Cognard et al. (2017),ω obs = 0.1033(29) deg yr −1 ; however, the difference is ∼ 2.5-σ significant. This is caused by the problem mentioned in Section 4, the use of inconsistent definitions of spin phase for the different data sets.
For the masses obtained from the Shapiro delay, GR predictsω GR = 0.09576(67) deg yr −1 . Thus, our new measurement agrees withω GR within 1σ. Therefore, the same applies to the total mass of the system derived from both methods: assuming GR, we obtain fromω obs M = 3.139(16) M . This represents a successful and precise (∼ 1%) test of GR. This is illustrated in Figure 6 by the fact that all mass constraints from the different PK parameters intersect in the same regions of the diagrams.
This statement relies on the fact that theω obs is relativistic. As discussed by Cognard et al. (2017), the largest additional contribution toω obs , which is caused by the proper motion of the system, is of the order of a few times 10 −6 deg yr −1 , i.e., about 100 times smaller than the current measurement uncertainty. Therefore, the assumption thatω obs 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 Figure 6, will also yield improved component masses.
We can use the DDGR model to combineω with the Shapiro delay and obtain self-consistent and more precise mass measurements; doing this we obtain M = 3.135(19)M , M p = 1.820(14)M and M c = 1.315(6)M . These values are 1.5 times more precise but slightly larger than the value derived from Shapiro delay only. The χ 2 for that fit is nearly identical to that of the best ELL1H+ model, indicating again the self-consistency of the relativistic effects.

Variation of 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 semi-major axis,ẋ: When they assume the VLBI proper motion, they obtainẋ = 3.5(30) × 10 −15 lt-s s −1 , our new value is −7.76(48) × 10 −15 lt-s s −1 ; the change is 3.8-σ significant. The uncertainty has decreased by a factor of 6. 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 lt-s 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 Figure 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 lt-s s −1 , a difference that is 3.4-σ significant. This is a robust result in all our fits: thė x is not strongly correlated with any other timing parameter.
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 of x caused by the orbital decay due to the emission of gravitational waves,ẋḊ, which is caused by the variation of the Doppler shift of the system D (analogous to the kinematic effect onṖ b discussed in section 7),ẋ A , which is caused by a change in the aberration due to spin precession, which is itself caused by spin-orbit coupling, andẋ SO -the change of the orbital inclination also due to the spin-orbit coupling. The spin-orbit coupling has several classical and relativistic terms caused by the rotations of the pulsar and of the companion. Finally, theẋṁ is caused by mass loss.
We have calculated these terms systematically, the results are also in Table 5. For the mass loss contribution toẋ, we have assumed the spin-down energy of the pulsar. Theẋ A term depends on the misalignment angle between the angular momentum of the pulsar and the angular momentum of the orbit, θ p (not de- Notes. Theoretical expectations for the different contributions toẋ. All in units of 10 −15 lt-s s −1 . The last term will only contribute toẋ obs in orbital models that do not account for the Einstein delay, γ, which all our models do. picted 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 Sec. 3 to constrainẋ A (see eqs. (2.5a), (2.25b), (3.24), and (3.35) in Damour & Taylor (1992)). The conclusion is the same. The change in the orbital inclination due to spin-orbit coupling can be split into contributions from the pulsar and the companion (see e.g. Barker & O'Connell 1975;Damour & Taylor 1992), i.e.ẋ SO =ẋ SOp +ẋ SOc . Due to the compactness of the pulsar,ẋ SOp is clearly dominated by the Lense-Thirring effect, and depends on the orientation (Φ p , θ p ) and magnitude of the pulsar spin.ẋ SOc depends on the unknown spin period of the WD companion (P c ), as well as the unknown spin-orientation (Φ c , θ c ).
For reasonable values of P c (∼ hours 9 ),ẋ SOc is also dominated by the Lense-Thirring effect. Both contributions are negligible. 10 We want to add that, due to tidal torques during the Roche-lobe overflow episode, one would expect an alignment of the spin axis, resulting in θ c 0, further suppressing any spin-orbit contribution toẋ by the WD.
There is an extra term listed by Lorimer & Kramer (2005) that could be caused by a possible third component of the system. If such a component existed, the motion of the PSR J2222−0137 binary around the centre of mass of the triple would cause a non-linear variation of 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, 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 post-Keplerian 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: sin ω. (2) For PSR J2222−0137 GR predicts a small Einstein delay (γ = 4.545 µs), which is largely the result of the small eccentricity. This and the relatively smallω 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 .

Annual orbital parallax
An alternative explanation for the largeẋ has to do with the fact that, apart from the secular variation of 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 of 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.
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 will base all subsequent discussions on this particular orbital model.

A self-consistent estimate of the component masses and orbital orientation of the system
In order to better determine the uncertainties and correlations between the masses and orbital configuration, we have made a self-consistent χ 2 map with the DDK model, where we additionally assume the validity of GR. Since, as discussed in section 5.4, we expect no significant additional contributions toẋ, we assume that any variations of x are caused by variations of 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 post-Keplerian parameters (M c ,ω, γ, 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 3-D probability density function (pdf) for the cos i, Ω, M space (Splaver et al. 2002). This 3-D pdf is then projected onto several planes and axes: 2-D pdfs are calculated for the cos i-Ω and the derived cos i-M c and M c -M p planes, and 1-D pdfs are calculated for the Ω, cos i, M, and derived M c and M p axes. These 2-D pdfs are represented by the black contours in the main panels of Figures 6 and 7, and some of the 1-D 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 3-D function: Some features, like 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 2-D or 3-D 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 section 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 large-scale dipolar field (the grey points in Fig. 3).
Regarding our measurement of ψ b , it is near one of the three possible values for the PA of the pulsar spin, ψ = −60(10) = 300(10) deg (see section 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.

Variation of 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 will now discuss the implications of this measurement in more detail. First, we update the estimate of the contribution from Shklovskii effect. Using the distance and proper motion from our re-analysis of VLBI data, we obtain: 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 ≡ |d sin b| ≤ 1.5 kpc, the vertical component of Galactic acceleration K z can be approximated as (Holmberg & Flynn 2004;Lazaridis et al. 2009) K z (10 −9 cm s −2 ) 2.27 z kpc + 3.68(1 − e −4.3 z kpc ). Assuming a 10% uncertainty in the vertical acceleration, we geṫ P b,Gal = −0.0142(13) × 10 −12 s s −1 .
Subtracting these two terms fromṖ b,obs we obtain the "intrinsic" variation of 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 GWs 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 ob-servedṖ b : which agrees well with zero. This limits any additional effects beyond GR, like a variation of Newton's gravitational constant or the emission of dipolar gravitational waves predicted by some alternative theories of gravity. For instance, following the calculations for scalar-tensor theories of section 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). Note, α 0 < 0.003 (95% C.L.) from Solar system experiments (Bertotti et al. 2003;Esposito-Farèse 2006). This tight limit on dipolar radiation, in combination with the large and well determined mass of PSR J2222−0137, makes this system an ideal laboratory for certain non-linear aspects of strong-field gravity, like spontaneous scalarization (Damour & Esposito-Farèse 1993;Shao et al. 2017). This will be explored in detail in a forthcoming publication (Zhao et al., in prep.).
Let us 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 thė P b,obs obtained with the model-independent 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.
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 ofṖ b,Gal Table 7. Different contributions toṖ b , in units of 10 −12 s s −1 . Several different models for Galactic potential are used for comparison. The value z 0 = 0 was used for these calculations. 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. For now, none of these differences change the fact that thė P 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 .

Summary and perspectives
In this paper, we present the results of 12-year timing of PSR J2222−0137, combining the data from Effelsberg, Nançay and Lovell radio telescopes with early GBT data. Furthermore, we have re-analyzed the astrometric VLBI data. Finally, we have also obtained polarimetric data from FAST.
The re-analysis of the VLBI data confirms most of the results presented by Deller et al. (2013), except for Ω, which changed by ∼ 180 deg. This resulted from our use of conventions that are fully consistent with those used in pulsar timing. We have also calculated the absolute position, with more realistic uncertainty estimates. Because of these, there is no longer a significant disagreement with the timing position.
The very high signal-to-noise ratio of the FAST data yields polarimetry consistent with the (corrected) Nançay data, and it has allowed a detection of several faint emission regions, which include, importantly, an interpulse. This has allowed an unambiguous determination of the geometry of the pulsar, in particular a precise determination of its 3-D 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 data sets. 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 3 − 5σ level. If we use instead the DMX model, which can describe short-term DM changes, then the proper motion is consistent with the VLBI result, but with much larger uncertainties. Because of this, in our timing we use DMX model and fix 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 post-Keplerian 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 of the semi-major axis,ẋ, is larger (in magnitude) than expected, the difference to the expected value is 3.4-σ significant. It is likely that this is caused by the presence of effects like 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 self-consistent analysis that assumes the validity of GR and takes all kinematic effects into account, we obtain a large pulsar mass of 1.831(10) M and a companion WD mass of 1.319(4) M . This is the largest confirmed NS birth mass (Cognard et al. 2017). This is only one of two recycled pulsar / 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 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.