Strong-Field Gravity Tests with the Double Pulsar

M. Kramer, I.H. Stairs, R.N. Manchester, N. Wex, A.T. Deller, W.A. Coles,
M. Ali, M. Burgay, F. Camilo, I. Cognard, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

M. Kramer, I.H. Stairs, R.N. Manchester, N. Wex, A.T. Deller, et al. Strong-Field Gravity Tests with the Double Pulsar. Phys.Rev.X, 2021, 11 (4), pp.041050. 10.1103/PhysRevX.11.041050. hal-03500470

HAL Id: hal-03500470

https://hal.science/hal-03500470

Submitted on 23 Mar 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Strong-Field Gravity Tests with the Double Pulsar

M. Kramer<sup>1,2,*</sup> I. H. Stairs<sup>3</sup> R. N. Manchester<sup>4</sup> N. Wex<sup>1</sup> A. T. Deller<sup>5,6</sup> W. A. Coles<sup>7</sup> M. Ali<sup>1,8</sup>

M. Burgay<sup>9</sup> F. Camilo<sup>10</sup> I. Cognard<sup>11,12</sup> T. Damour<sup>13</sup> G. Desvignes<sup>14,1</sup> R. D. Ferdman,<sup>15</sup> P. C. C. Freire<sup>1</sup>

S. Grondin<sup>3,16</sup> L. Guillemot<sup>11,12</sup> G. B. Hobbs,<sup>4</sup> G. Janssen<sup>17,18</sup> R. Karuppusamy<sup>1</sup> D. R. Lorimer<sup>19</sup> A. G. Lyne,<sup>2</sup>

J. W. McKee<sup>1,20</sup> M. McLaughlin<sup>19</sup> L. E. Münch<sup>1</sup> B. B. P. Perera<sup>21</sup> N. Pol<sup>19,22</sup> A. Possenti<sup>9,23</sup> J. Sarkissian,<sup>4</sup>

B. W. Stappers,<sup>2</sup> and G. Theureau<sup>11,12,24</sup>

<sup>1</sup>Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>2</sup>Jodrell Bank Centre for Astrophysics, The University of Manchester, Oxford M13 9PL, United Kingdom

<sup>3</sup>Department of Physics and Astronomy, University of British Columbia,
6224 Agricultural Road, Vancouver, British Columbia V6T 1Z1, Canada

<sup>4</sup>Australia Telescope National Facility, CSIRO Space and Astronomy,
P.O. Box 76, Epping, New South Wales 1710, Australia

<sup>5</sup>Centre for Astrophysics and Supercomputing, Swinburne University of Technology,
P.O. Box 218, Hawthorn, Victoria 3122, Australia

<sup>6</sup>ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav),
P.O. Box 218, Hawthorn, Victoria 3122, Australia

<sup>7</sup>Electrical and Computer Engineering, University of California at San Diego,
La Jolla, California 92093, USA

<sup>8</sup>Perimeter Institute for Theoretical Physics,
31 Caroline Street North, Waterloo, Ontario, Canada N2L 2Y5

<sup>9</sup>INAF-Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy

<sup>10</sup>South African Radio Astronomy Observatory,
2 Fir Street, Black River Park, Observatory 7925, South Africa

<sup>11</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E
CNRS-Université d'Orléans, F-45071 Orléans, France

<sup>12</sup>Station de Radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France

<sup>13</sup>Institut des Hautes Études Scientifiques, 91440 Bures-sur-Yvette, France

<sup>14</sup>LESIA, Observatoire de Paris, Universit PSL, CNRS, Sorbonne Université, Universit de Paris,
5 Place Jules Janssen, 92195, Meudon, France

<sup>15</sup>Faculty of Science, University of East Anglia, Norwich NR4 7TJ, United Kingdom

<sup>16</sup>David A. Dunlap Department of Astronomy and Astrophysics, University of Toronto,
50 St. George Street, Toronto, Ontario M5S 3H4, Canada

<sup>17</sup>ASTRON, Netherlands Institute for Radio Astronomy,
Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, Netherlands

<sup>18</sup>Department of Astrophysics/IMAPP, Radboud University,
P.O. Box 9010, 6500 GL, Nijmegen, Netherlands

<sup>19</sup>Department of Physics and Astronomy, West Virginia University,
P.O. Box 6315, Morgantown, West Virginia 26506, USA

<sup>20</sup>Canadian Institute for Theoretical Astrophysics, University of Toronto,
60 St. George Street, Toronto, Ontario M5S 3H8, Canada

<sup>21</sup>Arecibo Observatory, University of Central Florida, HC3 Box 53995, Arecibo, Puerto Rico 00612, USA

<sup>22</sup>Department of Physics and Astronomy, Vanderbilt University,
2301 Vanderbilt Place, Nashville, Tennessee 37235, USA

<sup>23</sup>Universitá di Cagliari, Dipartimento di Fisica, S.P. Monserrato-Sestu Km 0,
700-09042 Monserrato (CA), Italy

<sup>24</sup>Laboratoire Univers et Théories LUTh, Observatoire de Paris, PSL Research University,
CNRS/INSU, Université Paris Diderot, 5 place Jules Janssen, 92190 Meudon, France

<figure><img src="image_3.png" /></figure>

(Received 3 May 2021; accepted 25 October 2021; published 13 December 2021)

<footnote>*mkramer@mpifr-bonn.mpg.de</footnote>
Continued timing observations of the double pulsar PSR J0737-3039A/B, which consists of two active radio pulsars (A and B) that orbit each other with a period of 2.45 h in a mildly eccentric ($e = 0.088$) binary system, have led to large improvements in the measurement of relativistic effects in this system. With a 16-yr data span, the results enable precision tests of theories of gravity for strongly self-gravitating bodies and also reveal new relativistic effects that have been expected but are now observed for the first time. These include effects of light propagation in strong gravitational fields which are currently not testable by any other method. In particular, we observe the effects of retardation and aberrational light bending that allow determination of the spin direction of the pulsar. In total, we detect seven post-Keplerian parameters in this system, more than for any other known binary pulsar. For some of these effects, the measurement precision is now so high that for the first time we have to take higher-order contributions into account. These include the contribution of the A pulsar's effective mass loss (due to spin-down) to the observed orbital period decay, a relativistic deformation of the orbit, and the effects of the equation of state of superdense matter on the observed post-Keplerian parameters via relativistic spin-orbit coupling. We discuss the implications of our findings, including those for the moment of inertia of neutron stars, and present the currently most precise test of general relativity's quadrupolar description of gravitational waves, validating the prediction of general relativity at a level of $1.3 \times 10^{-4}$ with 95% confidence. We demonstrate the utility of the double pulsar for tests of alternative theories of gravity by focusing on two specific examples and also discuss some implications of the observations for studies of the interstellar medium and models for the formation of the double pulsar system. Finally, we provide context to other types of related experiments and prospects for the future.

DOI: 10.1103/PhysRevX.11.041050

Subject Areas: Astrophysics, Gravitation

# I. INTRODUCTION: PULSARS AS PROBES OF GRAVITATIONAL PHYSICS

The study of gravitational physics currently benefits from a number of experimental advances which provide unprecedented opportunities for constraining the underlying theory. Many of these methods focus on the study of compact objects, namely, neutron stars (NSs) and black holes (BHs), in order to confront the theories to be studied with data obtained under strong-field conditions. A prime example is the availability of ground-based gravitational-wave detectors which provided the first detection of gravitational waves (GWs) with Earth-bound detectors [1]. Almost 40 years earlier, the first evidence for GWs was provided by observations of the first binary pulsar, PSR B1913+16 [2,3]. Pulsar observations provide unique and complementary experimental constraints, especially thanks to the unrivalled precision enabled by radio astronomical observations and the method of pulsar timing (see, e.g., Ref. [4]), which permit us to trace the orbital evolution of a system over long periods of time. Pulsar observations, furthermore, provide valuable information on the structure of NSs, on plasma physics, on asymmetries in the supernova explosion of massive stars, and on the ionized content of the interstellar medium.

## A. Pulsars

The discovery of pulsars in 1967 [5] provided astronomers with natural and extraordinarily stable fly-wheel clocks, opening up the prospect of probing the composition of matter at extremely high densities. Pulsars were soon identified as rotating NSs, exceedingly dense and compact stellar remnants formed in supernova explosions [6,7].

Most known pulsars are located within our Galaxy, typically at distances of a few kpc [8] NSs, which typically have a mass of about 1.4 times the mass of the Sun ($M_{\odot}$) but a radius of only about 12 km, can and, in fact, are observed to spin very rapidly, up to approximately 700 times every second [9]. Precise timing of pulse arrival times at Earth shows that pulsars can have a rotational stability comparable to that of the best atomic clocks. But not all pulsars have highly stable periods. Observations of the rotational instabilities known as “glitches” (in mostly young pulsars) allow a kind of stellar seismology to probe the NS’s solid crust, its liquid interior, and the coupling between them (see, e.g., Refs. [10,11]).

Importantly, many pulsars, especially those with pulse periods in the millisecond range, both are stable rotators and are in a binary orbit with another star. Their great timing stability enables tiny variations in the arrival time due to, e.g., relativistic effects, to be accurately measured, thereby allowing tests of gravitational theories and investigations of many other physical phenomena.

The timing of binary pulsars also delivers measurements of NS masses with unprecedented precision, yielding evidence for NSs with maximum masses of at least $2 M_{\odot}$ [12–15]. Such large masses are incompatible with a number of theoretical equations of state (EOSs) proposed to describe the properties of superdense matter; the discovery of even more massive NSs would restrict the families of EOSs even further [16]. Other ways of constraining the EOS via pulsar timing include the potential measurement of the moment of inertia (MOI) of a NS.

Soon after the discovery of the first binary pulsar, PSR B1913+16 in 1975 [17], it was understood that measurements of relativistic effects in binary pulsars would allow
investigation of the masses and composition of NSs, with the first attempts at describing relativistic orbital effects as part of a timing model by Ref. [18] and relativistic spin precession (due to spin-orbit coupling) by Refs. [19,20]. To first order, measured orbital perturbations are determined by the line-of-sight motion as the pulsar orbits its companion, meaning that motion in the sky plane is not measurable. Consequently, these measurements permit a range of masses and orbital inclination angles. Multiple relativistic effects must be observed in order to both constrain the masses and begin testing gravitational theories, some of which may depend on the NS composition. Many aspects of gravitational theories are best tested when the strongly self-gravitating NSs are found as pulsars in systems with short orbital periods, nonzero orbital eccentricities, and orbital planes closely aligned to our line of sight.

All of these desirable characteristics are present in the still unique double pulsar system PSR J0737-3039A/B. The 2003 discovery [21] showed it to have a very tight orbit and excellent potential for gravity tests, and shortly thereafter it was shown to be a system with two orbiting active radio pulsars [22]. The system consists of a “recycled” 23-ms pulsar (“A”) and the second-born 2.8-s pulsar (“B”). Probably after being “dead” or at least undetectable for a few million years, the A pulsar was spun up and restored to detectability as its companion star, the progenitor of pulsar B, evolved and transferred matter and angular momentum to it. Subsequently, the progenitor of B exploded, leading to two NSs in the highly relativistic, slightly eccentric 2.45-h orbit that we observe today. Further details on the double pulsar system are given in Sec. II.

## B. Pulsar timing

Pulsar timing analyses begin with measurement of precise pulse times of arrival (TOAs) at the telescope. While many pulsars emit pulsed signals at high energies (x rays and gamma rays), it is the radio band which is of most interest here. Two orthogonal polarizations of the incoming electromagnetic wave are recorded by the receiver system of the telescope. These measurements are typically made at radio frequencies of hundreds of megahertz (MHz) or several gigahertz (GHz), with signals being Nyquist sampled at twice the receiver bandwidth. This sampling is often preceded by a frequency down-conversion and followed by channelization in frequency using a digital signal-processing system that can differ depending on the telescope. Astronomical radio signals suffer a dispersive delay due to free electrons in the interstellar medium (ISM), parametrized by the dispersion measure (DM), which must be taken into account. Within a given channel bandwidth, e.g., in our observations 1 MHz, the pulsar signals are “coherently dedispersed” to remove this dispersive delay and then folded at the topocentric pulsar period. Again, the methods depend on the telescope and the receiver. The folded profiles are averaged over a “subintegration” interval, which in our observations is 30 s,

resulting in a data cube of pulse amplitude versus pulse phase for each frequency channel.

The next step is to convert the observed pulse profiles into TOAs. This is done by comparing each profile, which is time stamped by the observatory's hydrogen-maser atomic clock, with a carefully prepared template, giving both a TOA and the uncertainty in that TOA [23]. The uncertainty due to template matching is independent between TOAs but depends on the strength of the observed pulsar signal, which can vary dramatically due to interstellar scintillation caused by turbulent fluctuations in the ionized ISM density.

However, the uncertainty of the template matching is not the only source of noise. Because of relatively large-scale fluctuations in the ISM electron density, the DM is variable on a timescale of months. Pulsars also display intrinsic “timing noise” due to variations in spin frequency with a similar timescale, and this is then common to all radio frequency channels. Both of these noise sources have steep power-law spectra and become negligible on short timescales such as the orbital period of PSR J0737-3039A/B.

Finally, the observed TOAs are compared with predictions from a timing model. The model describes (and, in fact, counts) the rotations of the pulsar and accounts for physical effects that modify the TOA. The corresponding computations are done using a timing analysis program, either Tempo [24] or Tempo2 [25] in our analysis.

The (usually small) differences between the observed and predicted TOAs are known as “timing residuals.” These residuals are basic to pulsar timing analyses, since they reveal effects that are not included in the timing model. The resulting signatures typically have a particular form in the timing residuals, and some of these may result from previously unknown effects.

## C. Tests of theories of gravity

The general theory of relativity, or general relativity (GR) [26], has passed its experimental tests with flying colors, so far. Despite its successes, it may not be our final answer in describing gravity on a macroscopic scale. There is a range of parameter space, from the quasistationary weak-field regime of the Solar System to the strong-field regime of compact objects like NSs and BHs, in all of which one may encounter an experiment where the theory could be falsified [27]. It is therefore important to test different aspects of the predictions of GR and alternative theories with different methods. For instance, observations with gravitational-wave detectors are able to test the highly dynamical strong-field regime and radiative aspects of gravity, but they are not able to test aspects of light propagation in strong fields. This, on the other hand, and other aspects can be tested with binary pulsars.

The equations of GR, and indeed of alternative theories of gravity, are nonlinear and must be approximated for comparison with binary pulsar data. Damour and Deruelle [28,29] provide a leading-order pulsar timing model which includes the effects expected in GR, such as the advance of
periastron and Shapiro delay, but which is parametrized in a way that does not assume the validity of GR or any other theory of gravity. Once measured, these “post-Keplerian” (PK) parameters can be used to determine masses (based on an assumed theory) and perform self-consistency tests of theories; this is the approach subsequently taken in timing the Hulse-Taylor pulsar [30] and other systems [27,31–34]. Damour and Taylor [35] expand the formalism to include parameters based on the pulse profile changes expected in relativistic spin precession and simulate measurement timescales for several parameters.

In this framework, any given relativistic theory of gravity provides a description of the PK parameters as functions of the measured Keplerian parameters and the two *a priori* unknown masses of the binary system. Measuring *n* PK parameters, where *n* > 2, overdetermines this system of equations, providing *n* − 2 independent tests of the studied theory. In the work presented here, by measuring more PK parameters than in any other system and by measuring them also more precisely than usual in any other, we find we need to go well beyond what has been considered in earlier work. As an example, for two of the PK parameters, besides their dependence on the masses and the Keplerian parameters, we need to incorporate their dependence on the MOI of pulsar A in our analyses.

### D. A coordinated gravity experiment

In many respects, we have reached a juncture in the application of binary pulsars to tests of gravitational physics. As we demonstrate in this paper, from now on, we have to consider a number of effects that could be neglected in the past but now require attention and the application of new methods. This is true now for the double pulsar but will eventually also apply for other systems in the future.

The high precision of our measurements that we describe in Sec. III forces us to consider the relativistic mass loss of the system due to the pulsar spin-down, while considering the EOS is now essential to interpret our observational results. We show explicitly that relativistic orbital deformation needs to be accounted for in gravity tests based on time dilation (a combination of second-order Doppler effect and gravitational redshift). Higher-order light-propagation effects in strong gravitational fields are also clearly evident in our data. More specifically, we measure an (higher-order) aberration effect due to the deflection of the radio signal in the gravitational field of the companion, which gives us the rotation sense of A relative to the orbital angular momentum. The advance of periastron, meanwhile, requires two important corrections: a second post-Newtonian term and a correction due to relativistic spin-orbit coupling that changes the orientation of the orbit over time [36,37]. The latter is proportional to the pulsar spin and the (generally negligible) companion spin and, therefore, encodes information about the MOI of the pulsar. This offers an opportunity to constrain this NS property,

important for determining the EOS, for the first time via pulsar timing. Consequently, we need to develop a new timing model to jointly account for all these effects.

Before we can apply such a model, however, we also need to determine the distance to the pulsar, as this is an important parameter in many applications of pulsar timing, in particular, in tests of the effects of gravitational-wave emission. Indeed, for the Hulse-Taylor pulsar, the limited precision of the known distance (and its acceleration relative to the Solar System) has prevented any improvement as a gravity experiment for a considerable time now [38]. Since the double pulsar is relatively close on a Galactic scale, there is a tiny, but measurable, curvature in the signal wavefront. This results in a small TOA modulation with a period of six months (as we track the source during Earth’s motion around the Sun) which allows us to determine the distance to the double pulsar.

Measurement of this “timing parallax” is, unfortunately, hampered by long timescale variations in both the pulsar intrinsic spin and the intervening ISM, the latter causing variations in the observed DM. The effect of DM variations can be corrected but not the intrinsic spin noise. Hence, we cannot constrain the pulsar distance using pulsar timing as tightly as we would prefer. These effects are discussed in Secs. IV B and VIII A.

Fortunately, we have a second independent method of determining the distance to the double pulsar. This method uses high-angular-resolution imaging of the pulsar with continental-scale radio interferometers, a technique known as “very long baseline interferometry” (VLBI). It also depends on Earth’s motion around the Sun, detecting small annual modulations in the pulsar’s apparent *position* on the sky (“annual geometric parallax”; see, e.g., Ref. [39]). Since the measurements are made in the plane of the sky, as opposed to pulsar timing which is sensitive to line-of-sight changes, such VLBI measurements are subject to different systematic effects and offer complementary information. In our analysis, we adopt the weighted mean of these two independent distance measurements to obtain the distance with a sufficiently small uncertainty.

With the distance measurement in hand, we develop and present the new full-relativistic timing model in Sec. V and present its application to the analysis of our TOAs in Sec. VI. Using the measurement of PK parameters, we obtain precise mass measurements and derive finally a set of independent tests of GR, finding superb internal consistency. Our most precise test is one related to GR’s prediction for GW emission, validating this prediction at a level of $1.3 \times 10^{-4}$ with 95% confidence.

The close agreement of the double pulsar timing with GR, in turn, means that one can use these observations to place tight constraints on various alternatives to GR. In Sec. VII, we demonstrate this on the basis of two well-studied gravity theories. The first such theory is “Damour–Esposito-Farése” (DEF) gravity, a two-parameter monoscalar-tensor gravity [40], containing GR as a limit with the two parameters
$α_0 = β_0 = 0$. DEF gravity exhibits the typical effects one would expect if the strong equivalence principle (SEP; see, e.g., Ref. [41]) is violated, for instance, the existence of (scalar) dipolar GWs and a location-dependent gravitational constant, both leading to characteristic modifications of the PK parameters [42]. Moreover, in certain regions of the parameter space, DEF gravity shows genuine nonperturbative strong-field effects, present only in NSs and, therefore, not testable in the weak-field regime of the Solar System. As shown in detail in Sec. VII, with the observations presented here, we limit such deviations from GR, further constraining the $α_0 - β_0$ parameter space of DEF gravity.

The second alternative to GR which we confront with our constraints from the double pulsar is Bekenstein's tensor-vector-scalar theory (TeVES) [43], a MONDian relativistic gravity theory that evades the need for dark matter in galaxies by a modification of GR. In Sec. VII, we show that the TeVeS is practically incompatible with our observations. Despite the fact that the TeVeS is already falsified by the confirmation that tensor modes of GWs travel with the speed of light [44], the double pulsar experiment has its own merits. It tests specifically the scalar sector of the theory and shows that, depending on the details of a MONDian gravity theory, such a theory can be tested by the radiative and strong-field properties of binary pulsars.

As mentioned above, the inhomogeneous turbulent ISM in between the double pulsar and Earth leads to slow variations of the DM. Associated refractive variations affect the apparent pulsar position as measured by VLBI with a corresponding impact on the VLBI distance as discussed in Sec. IV C. However, the study of the ISM can also provide independent information on the pulsar distance. These constraints are discussed in detail in Sec. VIII, where we compare the measured pulsar distance with estimates based on established models of the Galactic ionized gas [45,46]. We also investigate the inferred structure of the ISM and bring the discussion full circle, demonstrating consistency with our preferred pulsar distance. As a last consistency

check of our results, in Sec. VIII C, we discuss the system's inferred low space velocity and its implications for formation of the binary system.

In Sec. IX, we discuss the future prospects for tests of gravitational physics with the double pulsar system, and, in Sec. X, we conclude with a summary of our results and how these complement other methods for testing gravity. Figure 1 provides an overview of both our experiment and the descriptions of its various components in this paper.

## II. THE DOUBLE PULSAR

In this section, we describe some relevant aspects of the double pulsar system in more detail, including the eclipses of the A pulsar emission by the magnetosphere of B and the resultant constraints on the orientation of the spin vectors of the two NSs relative to the orbital angular momentum vector.

As the original discovery paper [22] shows, the double pulsar system PSR J0737-3039A/B is an eclipsing dual-line binary system. Pulsar B is believed to have been formed in a low-kick supernova event, the details of which are still a matter of debate [47-49]. The relatively low eccentricity ($e = 0.088$), the small system transverse velocity ($v_{\text{trans}} \sim 10 \text{ km s}^{-1}$) [31], and the very small misalignment angle of the spin vector of A relative to the total angular momentum vector ($\delta_A < 3.2^\circ$) [50,51] are indeed all signatures of a low-kick birth event of the double pulsar. A retrograde solution ($180^\circ - \delta_A$) is possible but considered less likely, as it would require a very strong kick with a fine-tuned magnitude and direction. Furthermore, the prograde rotation of A has been confirmed by using emission properties of B [52]. In this work, we confirm this independently using the newly seen relativistic effect of aberrational light bending.

The eclipses of pulsar A by pulsar B were first detected in the B discovery paper [22] and shown to have a duration of just 30 s. The discovery paper already attributed the eclipses to absorption by the magnetosphere of B, and this idea was dramatically confirmed when it was shown that the eclipses are modulated at either the rotation rate of B or twice that rate, depending on the orbital phase [53]. This paper also discusses the likelihood that the properties of the eclipse would be affected by relativistic precession [19] of the B spin axis about the total angular momentum vector.

In contrast to A, B was born with a significant misalignment of its spin axis relative to the total angular momentum (which is dominated by the orbital angular momentum). Detailed modeling of the A eclipse modulation by the magnetosphere of B gives a measurement of B's spin misalignment angle of $\delta_B \simeq 50^\circ$ [54]. [55] The rate of the relativistic precession of B's spin vector is also derived from the eclipse modelling: The observed rate given by the PK parameter, $\Omega_B^{\text{spin}} = 4^\circ.77 \pm 0^\circ.66 \text{ yr}^{-1}$, is consistent with that predicted by GR (a 70.96-yr precessional period, corresponding to a precessional rate of $5^\circ.073 \text{ yr}^{-1}$) within

<figure><img src="image_4.png" /><figcaption>FIG. 1. An overview of the experiment described in this work. The named sections provide an orientation to descriptions of the various components.</figcaption></figure>
the uncertainty of 13% [54]. The spin precession causes changes in the observed pulse shapes of B [56,57] and ultimately resulted in B's disappearance from our view in March 2008 [58]. When the B pulsar returns to our view depends on the actual beam shape as well as the precessional rate [57,58]. The beam shape is changing with time, as B's magnetosphere is severely distorted by the wind emerging from A to form a magnetotail on the downstream side [22,53]. From the length of the eclipse, presumably caused by synchrotron self-absorption of A's radio emission in B's magnetosphere [59,60], one can estimate that the magnetosphere extends only to about 40% of B's light-cylinder radius. This is consistent with estimates of the dynamic pressure of A's wind [53,61].

The present work builds on the observations and analysis of the double pulsar system presented by Kramer *et al.* [31], together with the measurement of relativistic spin precession by Breton *et al.* [54]. Making use of observations over a 2.5-yr data span, the Kramer *et al.* work presented precise determinations of the relativistic periastron precession of the system, the combined effect of time dilation and gravitational redshift, the Shapiro delay effect due to light propagation in curved spacetime, and the decay of the orbit due to GW emission. These effects are described in more detail in Sec. V.

Combining these observations led to five strong-field tests of gravity. These make use of not only six PK parameters previously measured from relativistic effects, but also of the theory-independent [62] ratio of the two NS masses uniquely available in this system [31,54].

Since then, the system has been studied continuously using a number of radio telescopes, with improved data acquisition systems and better sensitivity, resulting in much improved timing precision over time. Here, we report on new results from the timing of A, extending the data span to 16.2 yr, and describe multiple new phenomena in the system.

### III. OBSERVATIONS AND SIGNAL PROCESSING

The vast majority of the results presented in this work are based on an extensive pulse timing experiment made over a 16.2-yr interval made at six observatories around the world: the Parkes 64-m radio telescope ("Murriyang") in New South Wales, Australia, the Green Bank 100-m radio telescope (GBT) in West Virginia, USA, the Westerbork Synthesis Radio Telescope (WSRT) in Netherlands, the Nançay Radio Telescope (NRT) in France, the Effelsberg 100-m radio telescope in Germany, and the Lovell 76-m radio telescope (LT) at Jodrell Bank Observatory in the United Kingdom. In addition, we undertook interferometric imaging observations using the Very Long Baseline Array (VLBA) to obtain complementary information on the pulsar's position in the plane of the sky as a function of time. Details of the datasets from these telescopes and the signal-processing methods are described here.

Timing observations commenced with the confirmation of the A pulsar search candidate in April, 2003. The TOAs included in our timing analysis lie between 2003 May 1, modified Julian day (MJD) 52760, and 2019 July 3, MJD 58667, a total span of 16.2 years. Details of the different datasets used in the analysis are given in Table I. This table lists in order, the observatory and/or telescope used, the common designations of the receiver and signal-processing systems, the center frequency and total bandwidth of the dataset, the number of frequency channels in the processed data and whether or not the channel data are coherently dedispersed, the basic subintegration or TOA sampling time, the data span, the number of TOAs, and the weighted root-mean-square (rms) timing residual of a basic fit to each dataset using Tempo2. As discussed in more detail below (Sec. IV), two related datasets are used for the timing analyses. For the determination of astrometric parameters and DM variations, 4-min-sampled TOAs are analyzed using Tempo2; the penultimate line of Table I summarizes this dataset. Binary and relativistic parameters are determined using a modified version of Tempo which uses a new timing model and modifications not available in Tempo2 (see Sec. V below for more details) with the 30-s-sampled TOAs. The final line in Table I gives the parameters of this combined dataset. Figure 2 illustrates the time and frequency coverage of the datasets from the different observatories.

In many respects, observations and analyses are identical or similar to those used in Ref. [31]. In the following subsections, we describe the summaries and details where required for the data-acquisition systems of the contributing observatories.

#### A. Parkes Telescope

As summarized in Table I, the Parkes observations were made in three main bands, centered around 700, 1400, and 3100 MHz, respectively, using a variety of receivers and signal processors. As illustrated in Fig. 2, in terms of time coverage, the Parkes 1400-MHz dataset is the most complete of those analyzed. All datasets have a basic time sampling or subintegration length of 30 s.

The 700-MHz and most of the 3100-MHz observations were made using the dual-coaxial 10-cm/50-cm receiver [64], while the 1400-MHz observations primarily use the center beam of the 20-cm 13-beam (MB) receiver [65]. For several months in 2016 while the MB receiver was refurbished, the H-OH receiver which covers a 512-MHz band centered at 1526 MHz was used. Following retirement of the 10-cm/50-cm receiver in 2018 October, the ultra-wideband low (UWL) receiver [66] is used to extend the 1400- and 3100-MHz datasets. All of these receivers have (or had) orthogonal linearly polarized receptors. Signal-processing systems used with these receivers comprise three groups: analog filter banks (AFB), digital filter banks (DFB), and two coherent dedispersion systems used for the
<table><caption>TABLE I. Summary of datasets. Listed are the corresponding observatory, receiver, and instrument used, together with the center frequency (Ctr freq.), bandwidth (BW), subband bandwidth, the dedispersion method (DD), the integration time per TOA, the data span, the number of TOAs, and the root-mean-square residuals (rms Res).</caption><thead><tr><th>Observatory</th><th>Receiver</th><th>Instrument</th><th>Ctr freq. (MHz)</th><th>Total BW (MHz)</th><th>Subband BW (MHz)</th><th>DD<sup>a</sup></th><th>T<sub>s</sub> (s)</th><th>Data span (MJD)</th><th>No. TOAs</th><th>rms Res. (μs)</th></tr></thead><tbody><tr><td>Parkes</td><td>50-cm</td><td>AFB</td><td>680</td><td>64</td><td>16</td><td>I</td><td>30</td><td>53054–55371</td><td>2588</td><td>186.1</td></tr><tr><td>Parkes</td><td>50-cm</td><td>CPSR2</td><td>732</td><td>64</td><td>16</td><td>C</td><td>30</td><td>53954–55100</td><td>3480</td><td>110.6</td></tr><tr><td>Parkes</td><td>50-cm</td><td>DFB3</td><td>732</td><td>64</td><td>16</td><td>I</td><td>30</td><td>55042–56748</td><td>32 595</td><td>143.1</td></tr><tr><td>Parkes</td><td>50-cm</td><td>CASPSR<sup>b</sup></td><td>720</td><td>200</td><td>12.5</td><td>C</td><td>30</td><td>55617–58396</td><td>2351</td><td>103.8</td></tr><tr><td>Parkes</td><td>MB</td><td>AFB</td><td>1390</td><td>256</td><td>64</td><td>I</td><td>30</td><td>52760–55371</td><td>111 835</td><td>96.9</td></tr><tr><td>Parkes</td><td>MB</td><td>DFB1</td><td>1369</td><td>256</td><td>64</td><td>I</td><td>30</td><td>54356–54409</td><td>3179</td><td>94.5</td></tr><tr><td>Parkes</td><td>MB</td><td>DFB2</td><td>1369</td><td>256</td><td>64</td><td>I</td><td>30</td><td>54408–54930</td><td>9201</td><td>83.2</td></tr><tr><td>Parkes</td><td>MB</td><td>DFB3</td><td>1369</td><td>256</td><td>64</td><td>I</td><td>30</td><td>54680–56769</td><td>46 853</td><td>67.9</td></tr><tr><td>Parkes</td><td>MB</td><td>DFB4<sup>b</sup></td><td>1369</td><td>256</td><td>64</td><td>I</td><td>30</td><td>56829–58667</td><td>37 217</td><td>70.1</td></tr><tr><td>Parkes</td><td>H-OH</td><td>DFB4</td><td>1526</td><td>512</td><td>64</td><td>I</td><td>30</td><td>57444–57673</td><td>11 817</td><td>79.3</td></tr><tr><td>Parkes</td><td>UWL</td><td>DFB4</td><td>1369</td><td>256</td><td>64</td><td>I</td><td>30</td><td>58431–58667</td><td>12 231</td><td>72.3</td></tr><tr><td>Parkes</td><td>10-cm</td><td>AFB</td><td>3030</td><td>768</td><td>768</td><td>I</td><td>30</td><td>52987–54619</td><td>4262</td><td>94.8</td></tr><tr><td>Parkes</td><td>10-cm</td><td>DFB2</td><td>3100</td><td>1024</td><td>512</td><td>I</td><td>30</td><td>54376–54526</td><td>690</td><td>78.6</td></tr><tr><td>Parkes</td><td>10-cm</td><td>DFB3</td><td>3100</td><td>1024</td><td>512</td><td>I</td><td>30</td><td>54618–56070</td><td>1120</td><td>90.8</td></tr><tr><td>Parkes</td><td>10-cm</td><td>DFB4<sup>b</sup></td><td>3100</td><td>1024</td><td>512</td><td>I</td><td>30</td><td>54755–58396</td><td>19 342</td><td>109.0</td></tr><tr><td>Parkes</td><td>UWL</td><td>DFB4</td><td>3100</td><td>1024</td><td>512</td><td>I</td><td>30</td><td>58264–58619</td><td>519</td><td>105.0</td></tr><tr><td>GBT</td><td>820-MHz</td><td>GASP<sup>c</sup></td><td>820</td><td>64</td><td>64</td><td>C</td><td>30</td><td>53448–55587</td><td>21 531</td><td>13.6</td></tr><tr><td>GBT</td><td>1400-MHz</td><td>GASP</td><td>1404</td><td>64</td><td>64</td><td>C</td><td>30</td><td>53266–55458</td><td>7397</td><td>22.0</td></tr><tr><td>GBT</td><td>820-MHz</td><td>GUPPI</td><td>820</td><td>200</td><td>25</td><td>I</td><td>30</td><td>55002–55431</td><td>27 539</td><td>18.8</td></tr><tr><td>GBT</td><td>820-MHz</td><td>GUPPI</td><td>820</td><td>200</td><td>25</td><td>C</td><td>30</td><td>55607–58140</td><td>255 301</td><td>18.2</td></tr><tr><td>GBT</td><td>1400-MHz</td><td>GUPPI</td><td>1500</td><td>800</td><td>100</td><td>I</td><td>30</td><td>55217–56257</td><td>46 521</td><td>26.7</td></tr><tr><td>GBT</td><td>1400-MHz</td><td>GUPPI</td><td>1500</td><td>800</td><td>100</td><td>C</td><td>30</td><td>56367–58227</td><td>76 173</td><td>21.3</td></tr><tr><td>Nançay</td><td>1.4-GHz</td><td>NUPPI</td><td>1484</td><td>512</td><td>64</td><td>C</td><td>30</td><td>55818–58211</td><td>174 173</td><td>60.8</td></tr><tr><td>Nançay</td><td>2.5-GHz</td><td>NUPPI<sup>d</sup></td><td>2520</td><td>128</td><td>128</td><td>C</td><td>30</td><td>56191–57958</td><td>2594</td><td>144.9</td></tr><tr><td>WSRT</td><td>300-cm</td><td>PuMa-II<sup>d</sup></td><td>334.6</td><td>70</td><td>8.75</td><td>C</td><td>240</td><td>54519–57119</td><td>16 110</td><td>89.2</td></tr><tr><td>Effelsberg</td><td>20-cm/7-beam</td><td>PSRIX<sup>d</sup></td><td>1360</td><td>150/250</td><td>15.3/15.9</td><td>C</td><td>160/240</td><td>55722–58328</td><td>3300</td><td>40.1</td></tr><tr><td>JBO/Lovell</td><td>1.4-GHz</td><td>PDFB<sup>d</sup></td><td>1516</td><td>384</td><td>12</td><td>I</td><td>240</td><td>55668–58219</td><td>33 780</td><td>126.9</td></tr><tr><td>All 4-min</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>240</td><td>52760–58667</td><td>199 913</td><td>13.0</td></tr><tr><td>All 30-s</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>30</td><td>52760–58640</td><td>916 648</td><td>26.0</td></tr></tbody></table>

700-MHz band, CPSR [67] and CASPSR [68]. The AFB systems had channel bandwidths of 0.5 MHz for the 700-MHz band and 1 MHz for the higher-frequency bands and utilize 1-bit digitizers. The DFB systems used 8-bit digitizers and implement polyphase filter banks with field programmable gate array processors (FPGAs). Both the AFB and DFB system utilized incoherent dedispersion to sum data across subbands. More details of these signal-processing systems can be found in Manchester *et al.* [69].

## B. Green Bank Telescope

Observations at the Green Bank Telescope (GBT) are carried out at 820 and 1400 or 1500 MHz. Typically, 5-h observations were conducted with monthly cadence and at alternating frequencies. Twice a year, we usually conducted concentrated observing campaigns in April-May and October-November using observations at 820 MHz on consecutive days to separate short- from long-term orbital

changes. The Green Bank Astronomical Signal Processor (GASP [70]) carried out 8-bit Nyquist sampling of the incoming dual-linear-polarization signal, after which it performs coherent dedispersion in software on a Linux-based cluster for each of several 4-MHz channels. Later observations were carried out with the Green Bank Ultimate Pulsar Processing Instrument (GUPPI [71]). Initially, the 8-bit digitized data streams were detected and dedispersed incoherently, while subsequently a coherent dedispersion mode is available. Coherent dedispersion was implemented later at 1500 MHz than at 820 MHz. GASP and GUPPI were operated in parallel during a transition period in order to carefully calibrate an expected small instrumental time offset between the two systems. For those overlapping epochs, in the final analysis we usually choose the GUPPI dataset, unless calibration issues with GUPPI occur. We note that earlier datasets recorded with the previous SPIGOT and BCPM processors are not

<footnote><sup>a</sup>Dedispersion technique: C = coherent, I = incoherent.</footnote>

<footnote><sup>b</sup>Band reference instrument for the DM modeling.</footnote>

<footnote><sup>c</sup>Band reference instrument for the 30-s TOA analysis.</footnote>

<footnote><sup>d</sup>Datasets not included in the 30-s TOA analysis.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 2. Time and frequency coverage of the 4-min-sampled TOAs from the contributing observatories.</figcaption></figure>

included in the current analysis (cf. Ref. [31]). A noise diode was used for flux calibration, with observations of the stable continuum source B1442+10 used to calibrate the diode noise power. Flux calibration, interference excision, and subband formation are carried out for GASP using ASPFITSReader [72,73], while for GUPPI they are carried out using PSRCHIVE [74]. The GUPPI interference excision is based on early work carried out by Grandy [75].

### C. Nançay Radio Telescope

Our dataset includes the results of high-cadence observations with the Nançay Radio Telescope (NRT), which is a transit telescope of the Kraus-type design with a collecting area equivalent to a 94-m parabolic dish. While observations of the double pulsar commenced late 2004, using an earlier back end, we make use of observations conducted since late 2011 with the Nançay Ultimate Pulsar Processing Instrument (NUPPI). NUPPI is a baseband recording system, similar to GUPPI used at the GBT; see, e.g., Ref. [76] for more details. Data were taken with the 1.4- and 2.5-GHz receivers which have frequency coverage of 1.1–1.8 and 1.7–3.5 GHz, respectively. Most of the observations with the 1.4-GHz receiver were made at a central frequency of 1484 MHz, while those with the 2.5-GHz receiver were generally centered at 2520 MHz. The effective central

frequency for a given observation depends on the excision of radio-frequency interference (RFI) which is determined via visual inspection. The data are polarization calibrated using PSRCHIVE's SingleAxis method. Observations of a reference noise diode were made prior to each pulsar observation and the data are calibrated, correcting for differential phase and amplitude between the two polarizations, before being folded into subintegration profiles.

### D. Effelsberg 100-m radio telescope

Observations with the 100-m radio telescope at Effelsberg of the Max-Planck-Institut für Radioastronomie were conducted with two different 20-cm receiver systems. For most sessions, the central beam of the seven-beam system with a bandwidth of 250 MHz was used. When this was not available, a single-pixel receiver was used, providing a bandwidth of 150 MHz. All data streams were coherently dedispersed using the ROACH-based PSRIX signal processor (see Ref. [77] for more details). As the double pulsar barely rises above the hills surrounding the telescope, full orbital tracks are not possible, and the system temperature is significantly increased because of spillover. Consequently, longer integration times of 160 and 240 s were used. Standard calibration procedures are applied as described in Ref. [77].
### E. Lovell radio telescope

Observations were also conducted with the 76-m Lovell telescope at Jodrell Bank of the University of Manchester. We used a receiver covering a frequency range 1300–1700 MHz, with a maximum usable bandwidth of 400 MHz. Data were acquired using a ROACH-based system that is essentially identical to the Effelsberg PSRIX system (see Ref. [78] for details). Gain and polarization calibration are achieved by monitoring a noise-diode signal and manual adjustments to power levels. Similar to the methods used at other observatories, RFI is excised following visual inspection of the data.

### F. Westerbork Synthesis Radio Telescope

The Westerbork Synthesis Radio Telescope (WSRT) consists of 14 25-m dishes arranged in an East-West direction. Observations were made at 334 MHz, the lowest frequency included in our dataset. Dual-polarization signals were acquired with the (nearly real-time) coherently dedispersing PuMaII instrument [79] using a total bandwidth of 70 MHz with 8.75-MHZ-wide subbands. Gain and phase differences between the two polarizations were adjusted during the phased-array calibration of the system (cf. Ref. [80]).

### G. Very Long Baseline Array

The VLBA consists of ten 25-m dishes spread across the continental USA, Hawaii, and the Virgin Islands. Dual-polarization observations were taken in standard continuum mode with 256 MHz of observing bandwidth distributed over eight 32-MHz-wide subbands centered at 1.56 GHz. Further details are given in Sec. IV C and Appendix A.

## IV. PULSE TIMING AND VLBI DATA ANALYSES

In this section, we describe how the pulse timing observations discussed in Sec. III above are processed to form pulse TOAs, how these TOAs are analyzed to determine the pulsar astrometric parameters and the long-term DM variations, and the results obtained. We also outline the methods used to analyze the VLBA interferometric data and the astrometric results obtained. A key result from these observations is our best estimate of the distance to the double pulsar, an important parameter in the gravitational tests described in Sec. V. Detailed descriptions of the observations, analysis methods, and results are deferred to Appendix A.

For the pulse timing analyses, two sets of TOAs are formed: (a) a set with 30-s spacing that are used for the relativistic-parameter analyses described in Sec. V and (b) a set with 4-min spacing that are used for the timing-based astrometry and the analysis of DM variations described in this section.

Analysis of the timing datasets is carried out in two distinct phases: The first phase uses Tempo2 to derive the

pulsar astrometric parameters and the DM variations based on the 4-min-sampled TOAs, whereas the second phase uses Tempo to derive the binary and relativistic parameters based on the 30-s-sampled TOAs. As indicated in Table I, datasets from Westerbork, Effelsberg, and Jodrell Bank observatories are not included in the Tempo analysis, since they have insufficient time resolution to correctly trace the Shapiro delay curve. The Nançay 2.5-GHz dataset is also omitted in the second analysis phase, as the TOA precision is limited and their number is relatively small. For the astrometric and DM analyses, the binary and relativistic parameters are held fixed and the 4-min sampling makes the various Monte Carlo and bootstrap analyses (described in Sec. IV B below) more tractable. The 30-s sampling interval chosen for the binary and relativistic analyses is consistent with that used for the Kramer et al. (2006) analysis (Ref. [31]) and represents a good compromise between our ability to resolve orbital effects, signal-to-noise ratio (S/N) considerations (and, hence, TOA uncertainties), and the potential impact of pulse jitter (see, e.g., Ref. [81]).

### A. Measurement of TOAs

In general, similar procedures for processing observational data and formation of TOAs are adopted for all telescopes and data acquisition systems. After removal of obvious narrow-band and broadband transient radio-frequency interference (RFI), data are calibrated, in most cases making use of a short observation of an injected noise diode signal that precedes or follows the pulsar observation. The noise diode signal itself is calibrated in flux density units by comparing it to a stable continuum source, normally Hydra A (3C 218) at Parkes and B1442+10 at Green Bank, and in polarization properties by analysis of rise-to-set tracks of the strong millisecond pulsar PSR J0437-4715 at Parkes [82] and the well-studied calibrator PSR B1929+10 at Green Bank. Following calibration, the data are formed into subbands whose width varies according to frequency band as listed in Table I, with wider bandwidths used at higher frequencies (a) because of the reduced DM smearing, (b) to maintain sufficient S/N, and (c) so that the timing program (Tempo or Tempo2) properly accounts for the fact that data at different frequencies received at a given time correspond to different orbital phases at emission because of the differential dispersion delay (see, e.g., Refs. [25,31]). The subbandwidths chosen are sufficiently small to ensure that this effect is not significant.

For the Parkes datasets, TOAs are obtained using the PSRCHIVE routine PAT with the Fourier phase gradient algorithm [23]. Separate pulse profile templates derived are used for the three main receiver bands (50-, 20-, and 10-cm). These are analytic templates aligned in phase and first used in Ref. [31], derived by fitting Gaussian components to pulse profiles that are obtained from averaging a large number (>100 000) of individual pulses.
This helps to avoid “self-mirroring,” i.e., using a template based on the same pulse profiles that are to be timed, which leads to potential biases in TOA uncertainties due to correlated noise contributions [31,83]. These same templates are used to derive TOAs for data from the Green Bank GASP instrument and from the WSRT, Effelsberg, and JBO/Lovell telescopes. The wider-bandwidth GUPPI and NUPPI data generally have higher S/N than those from GASP, and use of the analytic templates results in systematics in the residuals that worsen with increasing S/N. For GUPPI and NUPPI, we therefore construct separate templates based on multiple observations aligned and smoothed using ASPFITSReader and PSRCHIVE utilities such as psrsmooth [84].

As mentioned above, 4-min-averaged TOAs are used for the astrometric and DM analyses. The method used to obtain the 4-min TOAs differs for the different observatories. For Parkes, the original 30-s-sampled datasets are directly summed to form 4-min-sampled datasets using standard Psrchive routines and TOAs determined as for the 30-s data. For data from Green Bank and Nançay, the 30-s TOAs are averaged over 4-min intervals using the Tempo2 averageData routine. This routine determines the weighted mean residual and its uncertainty for TOAs in the 4-min sample and applies the necessary pulse phase offset to a central TOA. Data from WSRT, Effelsberg, and JBO/Lovell are used directly.

Because of the short sampling time and the effects of pulse jitter, pulse profiles often vary greatly in shape, e.g., having the interpulse stronger than the main pulse. Consequently, TOAs are sometimes very discrepant in phase. These discrepant TOAs are simply removed from the dataset to produce a “clean” set for each telescope and instrument. Next, we determine scaling factors “EFAC” and “EQUAD” for each instrument (see, e.g., Ref. [25]). These factors modify the initially derived uncertainty of each TOA, $\sigma_i$, according to $\sigma'_i = \text{EFAC} \times \sqrt{\sigma_i^2 + \text{EQUAD}^2}$ such that the reduced $\chi^2$ value obtained after fits of the initial timing model is close to unity. This ensures an appropriate weight for a given dataset in the final combined least-squares fit to our timing model. For well-behaved datasets free of systematic errors or radio interference signals, one expects the TOA uncertainties to reflect the true measurement uncertainties; hence, EFAC ~ 1 and EQUAD ~ 0. For the 30-s dataset, we determine EQUAD values for each observing epoch and exclude those where the resulting EQUADs show large deviation from zero. Then, a single EFAC value is assigned to the each of the datasets listed in Table I. For the 4-min datasets, EFACs and EQUADs are determined for each instrument and each band for each observatory using the efacEquad function of Tempo2. Because of the preliminary “cleaning,” the EFACs are generally close to unity and in most cases less than 1.2. EQUADs are in most cases 10 µs or less, although for some systems they are up to 40 µs.

Because of profile frequency evolution and the different templates used for different rf bands, there are systematic offsets between TOAs for different bands. For Parkes data, for each rf band, we determine TOA offsets for each instrument relative to a reference instrument using the Tempo2 and Tempo2 “JUMPS” facility. The reference instrument for each band is indicated in Table I. In most cases, there is significant overlap between pairs of datasets, which greatly improves the precision of the JUMP determination.

For Green Bank GUPPI data, Nançay 1.4-GHz data, and WSRT data, the fractional bandwidths are much larger, and it is necessary to determine TOA offsets between subbands. In the 4-min TOA analysis, in each case a subband close to the center of the band is chosen as reference. Since these jumps are best determined from individual datasets, they are held fixed in the subsequent analyses. In the 30-s TOA analysis, we determine phase offsets between the analytic templates and the separately derived GUPPI and Nançay templates. We account for these offsets by applying phase shifts to the corresponding TOAs prior to applying the timing model, thereby effectively aligning all templates in phase.

## B. Timing astrometry and DM variations

The 4-min-sampled TOAs from all observatories are analyzed using Tempo2 to determine the astrometric parameters for the pulsar and the variations in DM across the total data span. All analyses use the DE436 Solar-System ephemeris [85] to transfer TOAs to the Solar System barycenter (SSB) frame. Studies of the effects of using different Solar-System ephemerides for this transfer (see, e.g., Refs. [86,87]) show differences much below the TOA precision of our datasets.

Initial parameters for the Tempo2 fit are obtained from the Tempo analysis of the 30-s TOAs. Pulsar astrometric parameters (right ascension, declination, proper motions, and annual parallax), the pulsar rotational frequency ($\nu$) and its first time derivative ($\dot{v}$), and the secular orbital terms ($\dot{\omega}$ and $\dot{P}_b$) are fitted. All other binary and relativistic parameters are better determined by the 30-s TOA and, hence, are held fixed and updated in following iterations.

DM and common-mode (i.e., wavelength-independent) variations are fitted for, typically at 100-d intervals in the central part of the dataset and at somewhat longer intervals at the ends where the frequency coverage is poorer, using the methods described by Keith et al. [88]. In addition, jumps relative to the Parkes 20-cm dataset are fitted for all observatory bands listed in Table I. DM variations are relative to the reference DM, held at 48.917 208 pc cm$^{-3}$. The reference solar-wind electron density is set to zero for all fits. Because of the high ecliptic latitude ($-51.2^\circ$), the annual DM variations due to the solar wind are quite smooth, so they are mostly absorbed by the fitted DM offsets. Transient solar events, such as coronal mass ejections, are not absorbed in the fit but are smaller and
short lived. We note that the worst-case solar-wind DM contribution is negligible compared with the errors in estimating DM(t) (see Appendix B). Transient events also have a negligible effect on measured TOAs. Although TOA offsets between observatory bands resulting from DM variations are absorbed by the jumps, variations within each dataset are preserved, and these determine the derived DM variations. Higher-order rotational-frequency derivatives up to the fourth (v̄) are held fixed at values determined in the Tempo analysis of the 30-s dataset because of covariances with the common-mode variations. This whole process is iterated as the system parameters are refined.

Astrometric parameters and DM and common-mode variations are initially determined by iterated fits to the entire dataset for the default 100-d sampling of the DM and common-mode variations. Because of the known covariance between astrometric parameters and the sampling of the DM and common-mode variations, we then undertake a Monte Carlo analysis, randomly varying the spacing and phase of DM and common-mode sampling within defined limits, normally 64–256 d, for the spacing and for phase within ±0.5 of the current spacing. After 4096 iterations, the distribution for each parameter is fitted by a Gaussian function to give estimates of the parameter and its 1σ uncertainty.

To test for systematic effects related to data from a particular observatory, the Monte Carlo analysis is repeated on datasets that (a) exclude the GBT data, (b) exclude the Parkes data, and (c) contain just the GBT and Parkes data. These tests have 2048 iterations, and they are analyzed in the same way as the full dataset.

As an additional check on the reliability of the parameter uncertainties, we perform a bootstrap analysis (with replacement) of the full dataset using the 100-d DM and common-mode sampling. In a similar manner to the Monte Carlo analysis, a Gaussian function is fitted to the parameter distributions after 2048 iterations.

At timescales longer than a few months, steep-spectrum or “red” noise in the intrinsic pulsar spin rate which has a power spectrum of the form $f^α$ (where $α ≪ 0$) becomes important. In the initial analysis, much of this red noise is

removed by fitting for the first four pulse frequency derivatives. However, proper analysis of timing datasets in the presence of red noise requires the estimation of a noise covariance matrix and the use of a generalized least-squares (GLSQ) approach when fitting the timing model [89]. The DM and common-mode analysis does not require GLSQ, because the common mode absorbs all the red noise and the residuals are “white,” i.e., flat spectrum. In this case, a weighted least-squares analysis is sufficient. For the GLSQ analysis, the covariance matrix has $N_{TOA}^2$ elements, and it is not feasible to perform a full Monte Carlo analysis. We choose to estimate the effects of the GLSQ analysis compared to the approximate method described above by forming smaller averaged datasets with 16- and 32-min sampling using the Tempo2 averageData routine. Even for these smaller datasets, it is necessary to modify the indexing of the covariance matrix and eliminate unnecessary matrix inversions in Tempo2. The best-fit value of $α$ is −3.0, but changes of ±0.3 make very little difference in the parameters or their uncertainties. The results show that the initial analysis underestimates the parallax uncertainty by about 75% but does not significantly bias the value. For the other astrometric parameters, the increase in uncertainty is less than 50% (zero for the position terms), but there are some biases, in all cases less than 50% of the combined uncertainty.

Results for the astrometric parameters from the Monte Carlo analysis of the full 4-min-sampled dataset modified to reflect the biases and uncertainty changes indicated by the GLSQ analyses of the shorter datasets are listed in the second column in Table II. [90] The final three columns in Table II list Monte Carlo results from the partial datasets. All results are within or close to twice the combined uncertainty of the final result, with the parameters from the GBT-omitted dataset differing the most. This shows that the GBT dataset with its high S/N has the greatest influence on the final results. It is notable, though, that the measured parallax is not significantly biased by the dataset from any one observatory.

Results from the bootstrap analysis are given in the third column in Table II. Uncertainties from the bootstrap

<table><caption>TABLE II. Astrometric results from analysis of the 4-min timing datasets [see the text for details of the different datasets and analysis procedures, i.e., bootstrap (boot.) and Monte Carlo (MC)]. Numbers in parentheses are 1σ uncertainties relative to the last quoted digit of the parameter value. For the right ascension, we show only the seconds of $07^{h}37^{m}ss.s^s$ and for the declination only the arcseconds of −30°39′ss.s″.</caption><thead><tr><th>Dataset</th><th>All data</th><th>All data Boot.</th><th>All ex. GBT MC</th><th>All ex. Pks MC</th><th>Pks + GBT MC</th></tr></thead><tbody><tr><td>Right ascension, α (J2000) (s)</td><td>51.248 115(10)</td><td>51.248 111(6)</td><td>51.248 125(10)</td><td>51.248 105(10)</td><td>51.248 105(9)</td></tr><tr><td>Declination, δ (J2000) (′′)</td><td>40.704 85(17)</td><td>40.704 96(11)</td><td>40.705 44(21)</td><td>40.705 29(23)</td><td>40.705 25(18)</td></tr><tr><td>Proper motion R.A., μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-2.567(30)</td><td>-2.570(16)</td><td>-2.618(15)</td><td>-2.581(16)</td><td>-2.582(16)</td></tr><tr><td>Proper motion Dec., μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>2.082(38)</td><td>2.081(19)</td><td>2.178(35)</td><td>2.102(35)</td><td>2.105(25)</td></tr><tr><td>Parallax, π<sub>t</sub> (mas)</td><td>2.15(48)</td><td>2.40(11)</td><td>2.41(10)</td><td>2.01(22)</td><td>2.22(24)</td></tr><tr><td>Position epoch (MJD)</td><td>55045</td><td>55045</td><td>55045</td><td>55045</td><td>55045</td></tr><tr><td>Data span (MJD)</td><td>52760–58667</td><td>52760–58667</td><td>52760–58667</td><td>53266–58328</td><td>52760–58667</td></tr><tr><td>N<sub>TOA</sub></td><td>199 913</td><td>199 913</td><td>129 787</td><td>151 883</td><td>118 156</td></tr></tbody></table>
analysis are typically about half those from the Monte Carlo (MC) analysis. The mean values are also offset, although generally within the combined uncertainties. Both of these differences result from the covariance of the astrometric parameters and the sampling of the DM and common-mode variations.

Figure 3 shows the DM variations derived by averaging the Monte Carlo DM-offset results for the 4-min-sampled dataset in 100-d bins across the data span. As noted above, the DM results are unaffected by the GLSQ analysis. Apart from a few points with large uncertainties near the ends of the data span where there is less radio-frequency coverage, the DM shows fluctuations on timescales of a few hundred days with amplitude of up to $0.0015 \text{ pc cm}^{-3}$ and a downward trend of approximately $-1.41 \times 10^{-4} \text{ pc cm}^{-3} \text{ yr}^{-1}$ over the data span.

For the final Tempo analysis of the 30-s TOAs, the astrometric parameters and DM variations from the final Tempo2 analysis of the 4-min data are applied using a method described in Sec. VI A. The final analysis takes into account the DM variations, which we can describe by fitting a smoothly varying curve that reflects the uncertainties of the DM measurements. This is achieved using a Gaussian learning process as implemented in the PYTHON library George [91]. We explored a number of stationary kernels with very similar results. We finally adopted the Matern32Kernel leading to the modeling of the DM offsets as shown in Fig. 3. Each TOA is modified with a DM offset dithered according to the interpolated DM uncertainty at that epoch.

The observed DM variations imply the presence of transverse phase gradients which will shift the apparent position of the pulsar on a timescale of months. These shifts can be estimated from the measured intensity scintillation

<figure><img src="image_6.png" /><figcaption>FIG. 3. Variations in DM relative to the reference value, 48.917 208 pc cm$^{-3}$, derived by averaging over values from the Monte Carlo analysis at 100-d intervals (points with error bars). The solid line is the result of using a Gaussian learning process to describe the data with resulting uncertainties indicated by the gray band.</figcaption></figure>

[92]. They are of the same order as the parallax uncertainty and will affect both timing and VLBI astrometry, which we present below in Sec. IV C. We discuss these effects in Sec. VIII A. We believe that the interstellar electron density fluctuations responsible for both the DM variations and the scintillation occur in or close to the Gum Nebula, as discussed in Sec. VIII B. To avoid the poorly sampled data at the ends of the full data span, this latter discussion is based on a restricted dataset between MJDs of 53500 and 58300.

### C. VLBI astrometry

Astrometric terms in the pulsar timing model can be compared to independent estimates for these parameters obtained via radio imaging. PSR J0737–3039A/B was previously targeted by an astrometric program using the Australian Long Baseline Array [93], which measured an annual geometric parallax of $0.87 \pm 0.14$ mas. This measurement is in tension with the timing parallax presented in Table II. In an attempt to resolve this, we undertook an extended astrometric program with the VLBA. The details of these observations, data reduction, and robustness testing using in-beam reference sources are given in Appendix A, and the resultant astrometric parameters are summarized in Table III.

For this analysis, the proper motion measured by pulsar timing (and listed in Table II) is applied as a prior to the VLBI fit. While some previous studies highlight substantial disagreement between timing and VLBI proper motions (e.g., Ref. [94]), the uncertainty in the timing proper motion for PSR J0737–3039A/B is sufficiently small that even an error exceeding those seen previously for other millisecond pulsars would not materially affect the VLBI parallax uncertainty. Details of the effects of proper motion fitting are included in Appendix A. The timing and VLBI reference positions are consistent when propagated to the same reference epoch, although the precision of the VLBI reference position for PSR J0737–3039A/B is relatively low, for the reasons discussed in Appendix A.

The quantity of most interest, the parallax, lies between the previous VLBI result and the current timing value but is more precise than either. We discuss the difference between the distance implied by a straightforward inversion

<table><caption>TABLE III. VLBA astrometric results. Listed are right ascension (R.A.), declination (Dec.), proper motion in both coordinates, parallax and position epoch.</caption><tr><td>R.A., α (J2000)</td><td>07h37m51s247(1)</td></tr><tr><td>Dec, δ (J2000)</td><td>-30°39'40"68(1)</td></tr><tr><td>Parallax, π<sub>v</sub> (mas)</td><td>1.30<sup>+0.13</sup><sub>-0.11</sub></td></tr><tr><td>Proper motion in R.A., μ<sub>α</sub> (mas yr<sup>-1</sup>)<sup>a</sup></td><td>-2.567 ± 0.030</td></tr><tr><td>Proper motion in Decl., μ<sub>δ</sub> (mas yr<sup>-1</sup>)<sup>a</sup></td><td>2.082 ± 0.038</td></tr><tr><td>Position epoch (MJD)</td><td>58000</td></tr></table>

<footnote><sup>a</sup>A prior based on the timing proper motion was applied as described in Appendix A.</footnote>
of the VLBA parallax (770 ± 70 pc) and the timing parallax 465<sup>+134</sup><sub>-85</sub> pc in Sec. VIII A. For the reasons discussed in Appendix A, we adopt the weighted mean of our new VLBA parallax and the timing parallax as our best estimate of the system parallax and use the resultant distance (735 ± 60 pc) and transverse velocity (11.5 ± 1.0 km s<sup>-1</sup>) for the remainder of the paper. We stress that the precision of the distance measurement does not (yet) affect the gravitational tests presented below.

## V. BINARY AND RELATIVISTIC-PARAMETER ANALYSIS METHODS

In this section, we describe the methods used to derive the relativistic terms that form the basis of the tests of gravitational theories, the key results of this paper, from the observational data and results described in Secs. III and IV. The astrometric, pulsar spin, binary, and PK relativistic parameters, derived using a Monte Carlo analysis to allow for the uncertainties in the derived DM variations (Sec. IV), are presented in this section. The results of the gravitational tests are described in the following Sec. VI for tests of GR and in Sec. VII for tests of alternative gravity theories.

A timing model describing the rotational properties of pulsar A, its position and apparent movement on the sky, the propagation delays in the ISM, and the orbital motion and relativistic effects that cause deviations from an arrival time expected for a simple Keplerian orbit is applied to TOA measurements for pulsar A. Before the timing model can be applied, the measured topocentric arrival times need to be transferred to the SSB to account for the varying position of our telescopes during the course of the year and for relativistic effects in the Solar System [4]. Before we present the results of such a procedure, we explain the required relativistic binary timing model. We do this in some detail, as we introduce a number of new considerations, which eventually also become relevant for other systems in the future.

The barycentric arrival time $t_b$ of a pulsar signal is related to the proper time $T$ of the pulsar (modulo a constant factor) by a phenomenological timing model that accounts for all relevant contributions related to the orbital motion of the pulsar and the signal propagation in and near the binary system [29,35]:

$$t_b - t_0 = \frac{1}{D} [T + \Delta_R(T) + \Delta_E(T) + \Delta_S(T) + \Delta_A(T)]. \quad (1)$$

The time $t_0$ denotes a (chosen) reference epoch. The different $\Delta_a$ ($a = R, E, S, A$) on the right-hand side denote, respectively, the Rømer delay resulting (solely) from the orbital motion of the pulsar, the Einstein delay caused by time dilation along the pulsar's worldline, the signal propagation delay due to relativistic contributions to the propagation of the radio signal, and the aberration delay due to the fact that the radio signals originate from a

rotating moving source [18,29,95]. These contributions depend on the proper time $T$ and a set of Keplerian and PK parameters, which we explain in detail below. $D$ is the Doppler factor due to the velocity $v_b$ of the pulsar system with respect to the SSB [29]:

$$D = \frac{1 - \hat{\mathbf{K}}_0 \cdot \mathbf{v}_b/c}{\sqrt{1 - \mathbf{v}_b^2/c^2}}, \quad (2)$$

which to leading order is dominated by the unknown radial velocity $v_R = \hat{\mathbf{K}}_0 \cdot \mathbf{v}_b$ (>0 if the Earth-pulsar distance increases). $\hat{\mathbf{K}}_0$ is the line-of-sight (LOS) unit vector pointing from the SSB to the binary system. Only temporal changes of $D$ are of interest here, and, therefore, $D \equiv 1$ can be chosen at a given epoch (see discussion in Ref. [29]). The proper time of the pulsar is related to the rotational phase $\phi$ of the pulsar by

$$\frac{\phi(T)}{2\pi} = N_0 + \nu(T - t_0) + \frac{1}{2}\dot{\nu}(T - t_0)^2 + \frac{1}{6}\ddot{\nu}(T - t_0)^3 + \dots, \quad (3)$$

where $N_0$ is the pulse number at the reference epoch $t_0$. In the following, we discuss the individual contributions $\Delta_a$ and all significant effects that need to be accounted for in those contributions when analyzing the timing data presented here. In terms of theory-specific interpretation of the different $\Delta_a$ and their PK parameters, we primarily focus on GR. As it turns out, for several contributions we need to account for next-to-leading-order (NLO) terms in order to achieve a correct interpretation of our timing data. In the following sections, we therefore discuss the individual contributions to the SSB arrival times in detail.

### A. Rømer delay

In terms of the orbital motion, pulsar timing analysis of binary pulsars is based on a particularly compact solution of the first post-Newtonian (PN) two-body problem (Damour-Deruelle solution [28]). As it turns out, this quasi-Keplerian parametrization of the orbital motion is valid for a wide class of alternative theories, in particular, the class of (massless) scalar-tensor theories [29,35]. The Rømer delay, $\Delta_R$ in Eq. (1), related to the orbital motion of the pulsar, is given by [29]

$$\Delta_R = x \left[ \sin \omega (\cos u - e_r) + \sqrt{1 - e_r^2} \cos \omega \sin u \right], \quad (4)$$

where $x = a_p \sin i/c$ is the projected semimajor axis of the pulsar orbit. The angle $\omega$ is the longitude of periastron (measured from the ascending node [96]), and $u$ is the “relativistic” eccentric anomaly, related to the proper time of the pulsar via Kepler's equation
$$u - e_T \sin u = 2\pi \left[ \left( \frac{T - T_0}{P_b} \right) - \frac{\dot{P}_b}{2} \left( \frac{T - T_0}{P_b} \right)^2 \right]. \quad (5)$$

The $\dot{P}_b$ term accounts for any secular (intrinsic or apparent) change in the orbital period $P_b$. It therefore accounts for the decay in the orbital period due to GW damping, which enters the GR equations of motion at the 2.5PN level, i.e., the order of $v^5/c^5$. The parameter $T_0$ denotes a time of periastron passage. As becomes clear below, for the double pulsar it is important to distinguish between the three different eccentricities $e_r$, $e_\theta$, and $e_T$ that enter the two equations (4) and (5) and are different at their first PN correction [29]:

$$e_r \equiv e_T(1 + \delta_r), \quad e_\theta \equiv e_T(1 + \delta_\theta). \quad (6)$$

While the parameter $\delta_r$ is not observable in the double pulsar (it is absorbed by higher-order frequency derivatives), the “relativistic deformation of the orbit” $\delta_\theta$ has to be accounted for in the analysis (Sec. VI B 4). The “time eccentricity” $e_T$ corresponds to the fitting parameter called “eccentricity of orbit” in the Tempo and Tempo2 implementation of the Damour-Deruelle (DD) timing model [24,98] and, therefore, can be considered as the observed eccentricity in pulsar timing, a Keplerian parameter.

The longitude of periastron $\omega$ that enters the Rømer delay (4) changes over time, due to the relativistic precession of the apsidal line of the double pulsar. In the quasi-Keplerian parametrization of the orbital motion, the longitude of periastron is simply related to the “relativistic” true anomaly $\theta$ [28], i.e.,

$$\omega = \omega_0 + k\theta, \quad (7)$$

$$\theta = 2 \arctan \left[ \sqrt{\frac{1 + e_\theta}{1 - e_\theta}} \tan \left( \frac{u}{2} \right) \right]. \quad (8)$$

By adding or subtracting multiples of $2\pi$, the angle $\theta$ is adjusted such that it matches the number of full phase cycles of $u$, as determined by Eq. (5). The structure of the Rømer delay presented above is expected to be the same for many viable relativistic theories of gravity. What changes is how the PK parameters ($k$, $\dot{P}_b$, $\delta_r$, and $\delta_\theta$) depend on the Keplerian parameters ($P_b$, $x$, and $e_T$) and the inertial masses $m_A$ and $m_B$ of A and B, respectively. Detailed expressions for GR can be found in Refs. [4,29], where one also finds expressions for the PK parameters which we introduce below, in combination with the Einstein and Shapiro delays in the signal propagation.

Instead of the PK parameter $k$, in pulsar timing one generally quotes the value for the (secular) periastron advance parameter defined as $\dot{\omega} \equiv n_b k$, where $n_b \equiv 2\pi/P_b$ is the measured orbital frequency [cf. Eq. (5.8) in Ref. [36]]. Given the high numerical precision of $n_b$

normally obtained in pulsar timing observations, $\dot{\omega}$ and $k$ can be used synonymously. For that reason, whenever comparison with former publications is required, we use $\dot{\omega}$ for convenience.

## 1. Second post-Newtonian corrections

The high precision achieved for several PK parameters raises the question of whether higher-order corrections need to be taken into account. In terms of the Rømer delay, 2PN corrections to the quasi-Keplerian parametrization are calculated in Refs. [36,99–101]. The result is a generalized quasi-Keplerian description of the orbital motion, in which Kepler’s equation and the equation for the true anomaly are augmented by periodic terms of the order of $v^4/c^4$. These 2PN periodic terms, however, are well below the timing precision so far achieved in the double pulsar system and can be ignored. This also means that we can replace $e_\theta$ with $e_T$ in Eq. (8) when calculating the longitude of periastron $\omega$ via Eq. (7).

Apart from these additional periodic terms, parameters that are already present in the 1PN solution get extended by 2PN corrections. Of particular importance for the analysis presented here is the secular advance of periastron. The (intrinsic) advance of periastron to 2PN approximation was first calculated within GR in Refs. [36,102]. When expressed as a function of the masses and the (directly observed) Keplerian parameters, it is given by

$$k = k^{\text{1PN}} + k^{\text{2PN}} + \mathcal{O}(\beta_O^6), \quad (9)$$

where

$$k^{\text{1PN}} = \frac{3\beta_O^2}{1 - e_T^2}, \quad k^{\text{2PN}} = \frac{3f_O\beta_O^4}{1 - e_T^2}, \quad (10)$$

and

$$\beta_O = \frac{(GMn_b)^{1/3}}{c}, \quad (11)$$

$$f_O = \frac{1}{1 - e_T^2} \left( \frac{39}{4} X_A^2 + \frac{27}{4} X_B^2 + 15 X_A X_B \right) - \left( \frac{13}{4} X_A^2 + \frac{1}{4} X_B^2 + \frac{13}{3} X_A X_B \right), \quad (12)$$

with $M \equiv m_A + m_B$, $X_A \equiv m_A/M$, and $X_B \equiv m_B/M$ [cf. Eq. (5.18) in Ref. [36]]. Note that $X_A + X_B = 1$. With the numbers from the timing parameter (Table IV), one finds $\dot{\omega}^{\text{2PN}} \equiv n_b k^{\text{2PN}} \simeq 4.39 \times 10^{-4} \text{ deg yr}^{-1}$, which is about 35 times the measurement error of $\dot{\omega}$, quoted in Table IV. More importantly, other PK parameters, in particular, $s$ and $\dot{P}_b$, have now reached a precision that, for the first time in a binary pulsar, we have to account for $k^{\text{2PN}}$ in order to perform, for instance, consistent gravity tests.
The 3PN corrections to $k$ are calculated in Refs. [106,107]. They are of the order of $\beta_O^2 \approx 4 \times 10^{-6}$ times the 2PN terms and, therefore, absolutely negligible.

## 2. Spin-orbit contribution and equation of state

So far, we have ignored the influence of spin, i.e., the proper rotations of A and B, on the orbital dynamics. Several effects related to spin can lead to significant modifications of the equations of motion of a binary system, in particular, the relativistic spin-orbit and spin-spin coupling and the rotationally induced mass quadrupole moments [37]. For the analysis in this paper, only the coupling between the spin of A and the orbital motion is of any relevance (the spin of B is about 135 times smaller). It is numerically of 2PN order [33,108,109], which is typical for binary pulsars in relativistic double-NS systems [101]. Furthermore, since the spin of A is practically parallel to the orbital angular momentum (see Sec. II), there is only a contribution to the precession of the periastron that is of relevance here. We refer the reader to Refs. [33,110] for more details.

In order to incorporate spin-orbit coupling in our analysis, Eq. (9) needs to be extended by the Lense-Thirring (LT) term, i.e.,

$$k = k^{\text{1PN}} + k^{\text{2PN}} + k^{\text{LT.A}}, \quad (13)$$

where within GR the LT contribution is given by [36,37,111]

$$k^{\text{LT.A}} = -\frac{3\beta_O^2\beta_{S_A}}{1 - e_T^2}g_{S_A}^{\parallel}, \quad (14)$$

with

$$\beta_{S_A} = 2\pi\nu \frac{cI_A}{Gm_A^2}, \quad (15)$$

$$g_{S_A}^{\parallel} = \frac{X_A(1 + \frac{1}{3}X_A)}{\sqrt{1 - e_T^2}}. \quad (16)$$

Apart from the MOI $I_A$, all quantities in the above equations are known with high precision. $I_A$ depends on the EOS for NS matter, which is still afflicted by considerable uncertainty. Consequently, there is a range in the prediction for $k^{\text{LT.A}}$. From Eq. (14), one finds the numerical expression

$$\dot{\omega}^{\text{LT.A}} \equiv n_b k^{\text{LT.A}} \approx -3.77 \times 10^{-4} \times I_A^{(45)} \text{ deg yr}^{-1}, \quad (17)$$

where $I_A^{(45)} \equiv I_A/(10^{45} \text{ g cm}^2)$. Using the multimessenger constraints on the radius that can be inferred from Ref. [112] (probability distribution function of case F in their Fig. 1) [113] in combination with the radius-MOI relation for A given in Ref. [114], we find a range of $I_A^{(45)} \approx 1.15-1.48$ (95% confidence). [115] Alternatively, Eq. (17) can be used to infer limits for the MOI of A purely from the timing

observations of the double pulsar, if combined with two other suitable PK parameters [33,110,119]. A corresponding analysis is given in Sec. VI B 3.

## 3. Proper motion contributions

The proper motion of a binary pulsar leads to a change in its orientation with respect to the observer on Earth. Such a change leads to an apparent change in the longitude of periastron $\omega$ and the orbital inclination $i$ [120,121]. The change in $\omega$ leads to a proper motion related offset $\dot{\omega}^{\text{pm}}$ between the intrinsic and the observed advance of periastron. Using the proper motion and orbital inclination from Table IV, in combination with the longitude of the ascending node obtained from scintillation measurements [92], one finds $\dot{\omega}^{\text{pm}} \approx -4 \times 10^{-7} \text{ deg yr}^{-1}$ (see also Ref. [110]). This is about a factor of 30 smaller than the current measurement error for $\dot{\omega}$ (see Table IV) and can, therefore, be ignored. As a consequence, there is no need to distinguish between the observed and the intrinsic $\dot{\omega}$.

The change in the orbital inclination enters the timing model through a temporal change in the projected semimajor axis of the pulsar orbit, showing up as a $\dot{x}$ in the timing solution, if significant. However, this contribution is even smaller than the contribution to the advance of periastron, since it is greatly suppressed by the fact that $i$ is close to $90^{\circ}$ ($\dot{x}^{\text{pm}} \propto \cot i \approx 0.01$; see Table IV).

## 4. Next-to-leading-order contributions in the mass function

The inclination of the binary orbit is linked to the projected semimajor axis $x$ in Eq. (4) via the binary mass function. In Newtonian gravity, one finds (see, e.g., Refs. [4,35])

$$\sin i = \frac{n_b x}{\beta_O X_B}, \quad (18)$$

where $n_b$ and $x$ are both (observable) Keplerian parameters, generally known to very high precision for a binary pulsar. As we discuss later (Sec. V C), the measurement of the Shapiro delay in the double pulsar gives access to $\sin i$, and, therefore, Eq. (18) leads to an additional constraint for the two masses $m_A$ and $m_B$.

In the 1PN approximation, Kepler's third law, which enters the derivation of Eq. (18), gets modified by an additional term [see Eq. (3.7) in Ref. [28] and Eq. (3.7) in Ref. [35]]. Consequently, Eq. (18) gets modified as well at the 1PN level. Using the 1PN expression for Kepler's third law, one finds

$$\sin i = \frac{n_b x}{\beta_O X_B} \left[ 1 + \left( 3 - \frac{1}{3} X_A X_B \right) \beta_O^2 \right]. \quad (19)$$

We use the fact that for the Damour-Deruelle solution the Newtonian relation between the semimajor axis of the pulsar orbit and the semimajor axis of the relative orbit also holds at the 1PN level, i.e., $a_A = (m_B/M)a_R + O(v^4/c^4)$
<table><caption>TABLE IV. Timing parameters for PSR J0737–3039A in TDB units (see the text). Except for astrometry and DM, the parameters are derived using Tempo with the 30-s TOA dataset. Numbers in parentheses are $1\sigma$ uncertainties referred to the last quoted digit. The overall reduced $\chi^2$ is 0.97.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Right ascension (R.A.), α (J2000)</td><td>07h37m51s248115(10)<sup>a</sup></td></tr><tr><td>Declination (Dec), δ (J2000)</td><td>-30°39′40″/70485(17)<sup>a</sup></td></tr><tr><td>Proper motion R.A., μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-2.567(30)<sup>a</sup></td></tr><tr><td>Proper motion Dec., μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>2.082(38)<sup>a</sup></td></tr><tr><td>Parallax, π<sub>c</sub> (mas)</td><td>1.36(+0.12, -0.10)<sup>a</sup></td></tr><tr><td>Position epoch (MJD)</td><td>55045.0000</td></tr><tr><td>Rotational frequency (freq.), ν (Hz)</td><td>44.054 068 641 962 81(17)<sup>b</sup></td></tr><tr><td>First freq. derivative, ṷ (Hzs<sup>-1</sup>)</td><td>-3.4158071(11) × 10<sup>-15</sup><sup>b</sup></td></tr><tr><td>Second freq. derivative, Ṱ (Hzs<sup>-2</sup>)</td><td>-2.286(29) × 10<sup>-27</sup><sup>b</sup></td></tr><tr><td>Third freq. derivative, Ṱ̇ (Hzs<sup>-3</sup>)</td><td>1.28(26) × 10<sup>-36</sup><sup>b</sup></td></tr><tr><td>Fourth freq. derivative, Ṱ̇̇ (Hzs<sup>-4</sup>)</td><td>4.580(86) × 10<sup>-43</sup> <sup>b</sup></td></tr><tr><td>Timing epoch, t<sub>0</sub> (MJD)</td><td>55700.0</td></tr><tr><td>Profile evolution, FD parameter c<sub>1</sub></td><td>0.0000180(75)</td></tr><tr><td>Profile evolution, FD parameter c<sub>2</sub></td><td>-0.0001034(10)</td></tr><tr><td>Profile evolution, FD parameter c<sub>3</sub></td><td>0.0000474(26)</td></tr><tr><td>Dispersion measure, DM (pccm<sup>-3</sup>)</td><td>48.917 208</td></tr><tr><td>Orbital period, P<sub>b</sub> (day)</td><td>0.102 251 559 297 3(10)</td></tr><tr><td>Projected semimajor axis, x (s)</td><td>1.415 028 603(92)</td></tr><tr><td>Eccentricity (Kepler equation), e<sub>T</sub></td><td>0.087 777 023(61)</td></tr><tr><td>Epoch of periastron, T<sub>0</sub> (MJD)</td><td>55 700.233 017 540(13)</td></tr><tr><td>Longitude of periastron, ω<sub>0</sub> (deg)</td><td>204.753 686(47)</td></tr><tr><td>Periastron advance, Ṡ (deg yr<sup>-1</sup>) <sup>c</sup></td><td>16.899 323(13)</td></tr><tr><td>Change of orbital period, Ṗ<sub>b</sub></td><td>-1.247920(78) × 10<sup>-12</sup></td></tr><tr><td>Einstein delay amplitude, γ<sub>E</sub> (ms)</td><td>0.384 045(94)</td></tr><tr><td>Logarithmic Shapiro shape, z<sub>s</sub></td><td>9.65(15)</td></tr><tr><td>Range of Shapiro delay, r (μs)</td><td>6.162(21)</td></tr><tr><td>NLO factor for signal prop., q<sub>NLO</sub></td><td>1.15(13)</td></tr><tr><td>Relativistic deformation of orbit, δ<sub>θ</sub></td><td>13(13) × 10<sup>-6</sup></td></tr><tr><td>Change of proj. semimajor axis, ṡ</td><td>8(7) × 10<sup>-16</sup></td></tr><tr><td>Change of eccentricity, ṩ<sub>T</sub> (s<sup>-1</sup>)</td><td>3(6) × 10<sup>-16</sup></td></tr><tr><td colspan="2"><em>Derived parameters</em></td></tr><tr><td>sin i = 1 - exp(-z<sub>s</sub>)</td><td>0.999936(+9/-10)</td></tr><tr><td>Orbital inclination, i (deg)</td><td>89.35(5) or 90.65(5)</td></tr><tr><td>Total mass, M (M<sub>☉</sub>)<sup>d</sup></td><td>2.587052(+9/-7)</td></tr><tr><td>Mass of pulsar A, m<sub>A</sub> (M<sub>☉</sub>)<sup>d</sup></td><td>1.338 185(+12/-14)</td></tr><tr><td>Mass of pulsar B, m<sub>B</sub> (M<sub>☉</sub>)<sup>d</sup></td><td>1.248 868(+13/-11)</td></tr><tr><td>Galactic longitude, l (deg)</td><td>245.2357</td></tr><tr><td>Galactic latitude, b (deg)</td><td>-4.5049</td></tr><tr><td>Proper motion in l, μ<sub>l</sub> (mas yr<sup>-1</sup>)</td><td>-3.066(35)</td></tr><tr><td>Proper motion in b, μ<sub>b</sub> (mas yr<sup>-1</sup>)</td><td>-1.233(31)</td></tr><tr><td>Distance from π<sub>c</sub>, d (pc)</td><td>735(60)</td></tr><tr><td>Transverse velocity, v<sub>T</sub> (kms<sup>-1</sup>)</td><td>11.5(10)</td></tr></tbody></table>

<sup>a</sup>See Secs. IV B and IV C for the derivation of these values.

<sup>b</sup>See Ref. [103].

<sup>c</sup> Ṡ ≡ 2πk/P<sub>b</sub>. k is the PK timing parameter in Eq. (7).

<sup>d</sup>See Ref. [104].

(see Refs. [28,29] and the discussion in Ref. [63]). With P<sub>b</sub>, x, and the masses from Table IV, one finds for the 1PN correction in Eq. (19) approximately 1.3 × 10<sup>-5</sup>, which is only about a factor of 1.3 larger than the error for sin i in Table IV. For that reason, we use the full 1PN mass function (19). This is the first time that 1PN corrections to the mass function become relevant for any binary pulsar.

## 5. *Secular changes in orbital period*

The observed change in the orbital period is a combination of effects intrinsic to the system and apparent changes related to a temporal change in the Doppler factor D in Eq. (1) [122]. For the double pulsar, by far the dominant contribution to Ṗ<sub>b</sub> is the orbital period decay due to the emission of GWs. In GR, GW damping enters at the 2.5PN level in the equations of motion (see Ref. [123] and references therein). The explicit expression for the leading-order changes due to GW emissions for an eccentric orbit are worked out in Ref. [124] (see Ref. [4] for a more pulsar-astronomy-adapted expression). The Hulse-Taylor pulsar is the first binary system where the leading-order GW damping has been tested [2,30].

The NLO correction to the change in the orbital period corresponds to the 3.5PN terms in the equations of motion and, hence, to the 1PN corrections in the radiation reaction force [125–127]. It is calculated in Ref. [128]. For the double pulsar, this contribution amounts to about -1.75 × 10<sup>-17</sup> [110]. This is about a factor of 4.5 smaller than the error in Ṗ<sub>b</sub> (see Table IV). Although this higher-order correction is, in principle, still negligible, we include it in our analysis. This is of particular interest for the comparison with the LIGO/Virgo results in Sec. VI B 2. In the near future, however, that contribution will become of importance (see Ref. [110]).

Yet another intrinsic effect that changes the orbital period in the double pulsar is the mass loss related to the spin-down of the pulsars [122]. This mass loss is a result of Einstein's energy-mass equivalence in the sense that here one is seeing the loss of mass associated with the loss of rotational (kinetic) energy of the pulsar. In Ref. [110], these contributions are calculated based on Eq. (4.1) in Ref. [122]. While for B this is negligible, for A one has [110]

$$ \dot{P}_{\text{b}}^{\dot{m}_{\text{A}}} = 2.3 \times 10^{-17} \times I_{\text{A}}^{(45)}. \quad (20) $$

For two reasons, it is important to include this contribution in the analysis below. First, given the range for $I_{\text{A}}^{(45)}$ (see the end of Sec. VA 2), $\dot{P}_{\text{b}}^{\dot{m}_{\text{A}}}$ can be as large as 3.1 × 10<sup>-17</sup>, which is a fair fraction of the measurement error of Ṗ<sub>b</sub> (see Table IV). Second, and more importantly, when estimating a MOI constraint based solely on the double pulsar
observations, $\dot{P}_b^{\text{mA}}$ and its dependence on $I_A$ has a significant impact on the result (see discussion in Ref. [110]).

Finally, there are the external contributions to the observed $\dot{P}_b$, related to a temporal change of $D$, where $\dot{P}_b/P_b = \dot{D}/D$ [122]. There is an apparent radial acceleration due to the transverse motion of the double pulsar with respect to the SSB, which leads to the Shklovskii effect [129]:

$$ \dot{P}_b^{\text{Shk}} = \frac{(\mu_{\alpha}^2 + \mu_{\delta}^2)d}{c} P_b, \qquad (21) $$

and a physical radial acceleration due to the Galactic gravitational potential, leading to a $\dot{P}_b^{\text{Gal}}$. Both effects depend on the distance $d$, which is the main source of uncertainty in those corrections (see Sec. VIII A below). For the Galactic potential, we use the model provided in Ref. [130]. We also test the analytic correction of Ref. [122] (adapted for $b < 0$) with the latest Galactic parameters from Ref. [131]. The result is practically the same. Numerical details for these contributions are given in Sec. VI B 2.

## B. Einstein delay

The Einstein delay $\Delta_E$ in Eq. (1) links the proper time of the pulsar $T$ to the coordinate time of the binary system $t$. To leading order, it can be viewed as a combination of the second-order Doppler effect and the gravitational redshift caused by the companion. If $T$ is renormalized such that the orbital period is the same as measured in $t$, one finds to leading order the following expression [29,95]:

$$ \Delta_E = \gamma_E \sin u. \qquad (22) $$

The amplitude of the Einstein delay $\gamma_E$ is a PK parameter that in GR can be calculated according to

$$ \gamma_E = e \frac{P_b}{2\pi} X_B(1 + X_B) \beta_O^2. \qquad (23) $$

In alternative theories of gravity, $\gamma_E$ is generally more complicated, since it also includes a periodic variation of the pulsar rotation due to a variation of the pulsar's MOI along the orbit [42,132,133]. The variation of the MOI is caused by a variation of the local gravitational constant along the (eccentric) orbit, a result of a violation of the SEP [41].

NLO contributions are calculated for GR in Refs. [101,134]. Similarly to the periodic 2PN terms in the orbital motion (see Sec. VA 1), these contributions are currently absolutely negligible. They are of the order of 10 ns, which is an order of magnitude below the precision for $\gamma_E$ (see Table IV). They are not expected to be of relevance before the full Square Kilometre Array [135] is in operation.

## C. Signal propagation delay

For binary pulsars which are sufficiently edge-on, i.e., $i$ close to 90°, the curved spacetime of the companion star has a significant effect on the propagation of the pulsar's radio signal. To leading order, we have the well-known Shapiro delay [136], which for binary pulsars is expressed in the following form [29,95]:

$$ \Delta_S^{(\text{LO})} = -2r \ln \Lambda_u, \qquad (24) $$

$$ \Lambda_u = 1 - e_T \cos u - s[\sin \omega(\cos u - e_T) + (1 - e_T^2)^{1/2} \cos \omega \sin u], \qquad (25) $$

with the two PK parameters $r$ and $s$, called the *range* and the *shape* of the Shapiro delay, respectively. The Shapiro shape can quite generally be identified with the sine of the orbital inclination, i.e., $s = \sin i$ [35,137]. The range of the Shapiro delay is linked to the mass of the companion. In GR, one finds $r = Gm_B/c^3$ [29,95].

The leading-order expression for the Shapiro delay in a binary pulsar is obtained in Ref. [95], by integrating along a straight line (in an isotropic coordinate system) and by assuming a static mass distribution while the signal propagates away from the pulsar (static limit). Relaxing the first assumption and accounting for the fact that the radio signal propagates along a curved path due to the deflection in the gravitational field of pulsar B leads to a modification of Eq. (24) [138,139]. This *lensing correction* to the Shapiro delay, given the current timing precision in the double pulsar, is, however, not yet observable due to a strong covariance with $s$. We test this against the real data as well as in mock data simulations [140] and find a shift in $s$ which is yet insignificant (less than $0.5\sigma$). One has to keep in mind that these effects are strongest around conjunction, in a part of the orbit where the timing precision is significantly reduced due to the (partly intermittent) eclipses approximately $\pm 1^\circ$ around superior conjunction, caused by the rotating plasma-filled magnetosphere of B [53,54,142]. [143] In spite of that, we extend the Shapiro delay implementation in Tempo by an adapted version of the elegant approximation in Eq. (73) in Ref. [144], i.e., $\Lambda_u \rightarrow \Lambda_u + \delta\Lambda_u^{\text{len}}$ with $\delta\Lambda_u^{\text{len}} = 2rc/a_R$, which accounts to good approximation (only few percent error) for the delay related to a curved signal path and, therefore, removes nearly completely the already insignificant bias to $s$. It is interesting to note that the lensing correction leads to a reduction in the (calculated) propagation time, which is a result of Fermat's principle (see, e.g., Ref. [145]).

A higher-order signal propagation effect that becomes relevant for the double pulsar timing in the meantime, and cannot be absorbed into other timing parameters, is the *retardation effect*, i.e., specifically the velocity-dependent contributions in the Shapiro time delay (1.5PN corrections) [146,147]. It results from the fact that the “lens,”
i.e., pulsar B, moves while the signal of A propagates across the binary system toward Earth. The result derived in Ref. [146] is absolutely sufficient for modeling the TOAs around superior conjunction, [148] and one has

$$ \Delta_S = -2r \ln(\Lambda_u + \delta\Lambda_u^{\text{len}} + \delta\Lambda_u^{\text{ret}}), \quad (26) $$

where retardation correction $\delta\Lambda_u^{\text{ret}}$ is of the order of $x/P_b$ and can directly be taken from Eq. (130) in Ref. [146]. As in Ref. [146], we omit the velocity-dependent term in front of the logarithmic function, since its contribution is negligible (less than about 0.015 $\mu$s).

For systems with $i$ close to 90°, like the double pulsar, one often obtains rather asymmetric error bars for the Shapiro shape parameter $s$. For that reason, in Ref. [149] the logarithmic Shapiro shape parameter is introduced:

$$ z_s \equiv -\ln(1-s). \quad (27) $$

It not only leads to a more Gaussian distribution for the fitted parameter, i.e., $z_s$ instead of $s$, it also ensures that the a priori free parameter $s$ always fulfills $s < 1$, which is important when interpreted as the sine of the inclination angle.

## D. Aberration

The term $\Delta_A$ in Eq. (1) collectively refers to delays that are related to the fact that the pulsar is acting as a moving “lighthouse” with beamed radio emission. These effects (with one exception) are proportional to the rotational period of the pulsar, $P \equiv 1/\nu$, and would be absent if the variability of the radio signal were caused by an (intrinsic) oscillation or variability instead of the rotation of the NS. For pulsar timing, aberration means that the proper angle of emission of an observed pulse changes along the orbit [18,29]. To leading order, the aberration delay can be written as

$$ \Delta_A^{(\text{LO})} = \mathcal{A}[\sin \psi + e_T \sin \omega] + \mathcal{B}[\cos \psi + e_T \cos \omega], \quad (28) $$

where $\psi \equiv \omega + \theta$ is the longitude of the pulsar with respect to the plane of the sky. The aberration coefficients $\mathcal{A}$ and $\mathcal{B}$ are proportional to $P$ and depend on the orientation of the pulsar with respect to the orbit and the observer. Since, as discussed in Sec. II, the spin of A is (practically) parallel to the orbital angular momentum, one has $\mathcal{B} = 0$ and

$$ \mathcal{A} = \frac{Px}{P_b(1 - e_T^2)^{1/2} \sin^2 i} \approx 3.65 \ \mu\text{s}. \quad (29) $$

In practice, $\mathcal{A}$ is a nonobservable parameter. If $\mathcal{A}$ is not provided in the timing model, it gets absorbed by a shift in various timing parameters (see Refs. [29,35] for details). For that reason, we a priori add the aberration coefficient $\mathcal{A}$

to our model, as a fixed parameter with the value given in the equation above.

Classical aberration, as expressed in Eq. (28), assumes a flat spacetime in which the signal propagates toward Earth. However, this is only a first approximation, which is no longer sufficient for an adequate description of the double pulsar timing observations, in particular, for TOAs near the superior conjunction of pulsar A. It is necessary to account for the deflection of the radio signal in the gravitational field of pulsar B. This gravitational signal deflection leads to a small change in the proper angle of emission, which, in turn, leads to a lensing correction to the classical (“longitudinal”) aberration of Eq. (28). The potential importance of such a correction in close to edge-on binary systems has first been discussed in Ref. [150]. For the specific situation of A (spin aligned with orbital angular momentum), one finds

$$ \Delta_A = \mathcal{A}[\sin \psi + e_T \sin \omega] + \mathcal{D} \frac{\cos \psi}{\Lambda_u}, \quad (30) $$

where

$$ \mathcal{D} = \frac{P}{\pi x + x_B}. \quad (31) $$

The quantity $x_B$ denotes the projected semimajor axis of pulsar B, which is an observed Keplerian parameter in the double pulsar [31]. The NLO contribution in $\Delta_A$ contains the same two PK parameters as the Shapiro delay in Sec. V C, i.e., the range $r$ ($\mathcal{D} \propto r$) and the shape $s$ [in $\Lambda_u$; see Eq. (25)]. Like in the Shapiro delay (see Sec. V C), also here we account for retardation in our Tempo implementation, i.e., use the position of B when the signal reaches the minimum distance to B. This leads to small corrections in the last term of Eq. (30), i.e., a corresponding shift in $\psi$ defined by the new position of B (and the position of A at emission) and a correction of $\Lambda_u$ like in Eq. (26). Although these corrections are not yet significant, we implement them along with the retardation in the Shapiro delay, for completeness.

In Refs. [147,151], it is pointed out that the approximation for this lensing correction to the aberration delay, as given in Ref. [150], loses its validity if the impact parameter for the radio signal becomes comparable to the Einstein radius $R_E \approx \sqrt{4Gm_Ba_R/c^2} \approx 2500$ km. [152] Our tests, based on simulated and the real timing data, show, however, that Eq. (30) is absolutely sufficient for the current double pulsar timing analysis. In fact, the main difference between the (retardation-corrected) approximation (30) and the full calculations (cf. Ref. [140]) is below 50 ns (postfit) and confined to a region very close to superior conjunction where, as explained above, the timing precision is anyway significantly compromised. One has to keep in mind that for outside the eclipse region one finds the impact parameter $\gtrsim 7R_E$, which for the deflection angle means $\lesssim 0.026^\circ$.
Finally, there is the lensing correction to the “latitudinal” aberration delay, related to a potential change in the observed pulse profile, as the LOS cuts a different part of pulsar A’s magnetosphere due to the signal deflection [147,151]. This effect does not directly depend on the rotational period of pulsar A but is linked to A’s rotation by having the LOS crossing the emission region of the pulsar. At this stage, it is not clear if there is at all a significant profile evolution across the range of the LOS variation. This is a question currently under investigation, utilizing the improved data quality provided by the MeerKAT telescope [154]. Because of a strong covariance with s and the eclipses near conjunction caused by the magnetosphere of pulsar B, this effect is not (yet) expected to be of any relevance [155].

### E. Fitting for NLO contributions in the Shapiro and aberration delay

As it turns out, the NLO contributions in the Shapiro and aberration delays cannot be tested separately in the double pulsar. We have mentioned already in Sec. V C that the lensing correction to the propagation delay is covariant with the Shapiro shape s. Furthermore, the NLO contribution to $\Delta_S$ is quite similar to the NLO contribution to $\Delta_A$, which renders a separation of $\delta\Lambda_u^{\text{ret}}$ and $\mathcal{D}$, at least under current timing precision, impossible. For that reason, in order to test for the significance of NLO signal propagation contributions in our data, we collectively rescale these contributions with a common factor $q_{\text{NLO}}$ that can be fitted for in our new timing model, i.e., for the retardation effect in the Shapiro delay

$$ \delta\Lambda_u^{\text{ret}} \rightarrow \delta\Lambda_u^{\text{ret}} \times q_{\text{NLO}}, \qquad (32) $$

and for effects related to signal deflection

$$ \delta\Lambda_u^{\text{len}} \rightarrow \delta\Lambda_u^{\text{len}} \times q_{\text{NLO}}, \qquad (33) $$

$$ \mathcal{D} \rightarrow \mathcal{D} \times q_{\text{NLO}}. \qquad (34) $$

In GR, one has $q_{\text{NLO}} = 1$. The NLO corrections above depend on Keplerian parameters and the PK parameters r and s [note that $a_R = c(x + x_B)/s$]. Therefore, $q_{\text{NLO}}$ simply measures the relevance of these NLO corrections and is not considered as a new PK parameter.

## VI. BINARY AND RELATIVISTIC-PARAMETER RESULTS

In this section, we present key results from this work relating to tests of Einstein’s general relativity. These results are based on the methods described in Sec. V above and include a number of previously undetected relativistic effects. All observed effects are shown to be fully consistent with GR, a remarkable confirmation of Einstein’s

extraordinary intuition and vision. In the first timing parameters subsection, we give a full list of the parameters derived from the analysis of the 30-s sampled dataset. The next subsection presents results depending on our best determined values for the masses of the A and B pulsars, including constraints on the MOI of pulsar A ($I_A$) and, hence, on the neutron-star EOS, and the newly detected relativistic effects. Using these masses, we can constrain the rate of change of the orbital period due to GW damping $\dot{P}_b^{\text{GW}}$ and probe light propagation effects in the gravitational field of a NS, key tests of GR, to unprecedented levels. We compare limits of deviation from GR in GW emission obtained with the double pulsar with limits from the LIGO/Virgo observations of BH-BH and NS-NS coalescences, showing that, at leading PN order (quadrupole formula), the double pulsar is more sensitive by 3 orders of magnitude, whereas at high PN orders, LIGO/Virgo gives more stringent limits, demonstrating the complementarity of these experiments.

### A. Timing parameters

For the fit of the timing model to the TOAs, one does not use Eq. (1) in its form $t_b = T + \Delta(T)$ directly. What is needed is the inverted timing model, i.e., something of the form $T = t_b - \tilde{\Delta}(t_b)$. In Ref. [29], an analytic inversion of the timing model is worked out, giving an analytic expression for $\tilde{\Delta}(t_b)$. However, the approximation used in Ref. [29], in particular, with the Shapiro delay, is no longer sufficient for the double pulsar, mainly due to the higher-order signal propagation and aberration delays. For that reason, we modify the standard routine of the DDS model [157] to perform a numerical inversion of model Eq. (1) (for details, see the public Tempo distribution after November 27, 2015).

The results of applying the timing model are shown in Table IV and Fig. 4. In addition to the binary parameters introduced and discussed in detail above, we fit for spin frequency and higher-order derivatives ($\nu, \dot{\nu}, \ddot{\nu}, \dddot{\nu}$, and $\dddot{\ddot{\nu}}$), FD parameters [158,159], and time offsets of the TOA datasets relative to the reference datasets as indicated in Table I.

FD parameters mitigate the effect of possible unmodeled frequency evolution of the pulse profile by allowing for a time delay between different subbands according to $\Delta t_{\text{FD}} = \sum_{i=1}^{N} c_i (\log f)^i$ [158]. Here, $c_i$ are the FD parameters and f the subband frequency in GHz, when $\Delta t_{\text{FD}}$ is expressed in seconds. Following Ref. [159], we employ an F test to determine that three FD parameters are needed. Similarly, we employ an F test to determine that four frequency derivatives are required to model the spin evolution of pulsar A.

Nonzero FD parameters suggest that future analyses of wideband data of the double pulsar will benefit from deploying smoothly varying, frequency-dependent
<figure><img src="image_7.png" /><figcaption>FIG. 4. Residuals obtained after applying the timing model to the 30-s dataset. From bottom to top: GBT 820 MHz (black, GASP; blue, GUPPI incoherent; red, GUPPI coherent), GBT 1400 MHz (same colors as for 820 MHz), Parkes (black, 50 cm; blue, 20 cm; red, 10 cm), and Nançay. The topmost shows residuals as a function of orbital phase, averaged to 15-s-long phase bins. The large uncertainties at $\Psi = 90^{\circ}$ are the result of the presence of the eclipse. For each, we add a histogram of the distribution of residuals on the right-hand side.</figcaption></figure>

templates [160,161]. But, in general, profile evolution and dispersion measure are highly covariant, with a fit for DM being able to (partly) absorb the effects of profile evolution and vice versa. Indeed, allowing the DM reference value to float in our final fit reduces the number of significant FD parameters to one. We note that the described choice does not have an impact on any of the binary timing parameters. We choose to fix the reference DM to the value given in Table IV and determine the DM variations as described in Sec. IV B.

The parameters shown in Table IV are the result of a large number of MC runs. Before each run, we select a random realization of the astrometry as determined in Sec. IV B, i.e., position, proper motion, and parallax that are consistent with the determined distributions. While these parameters are kept fixed for the given MC run, we use a modified version of Tempo to fit the described timing model to the remainder of the 20 timing parameters (see Table IV) and

22 constant time offsets (JUMPS) between the different datasets (see Sec. IVA). In each run, the DM value of a given TOA epoch is varied according to the DM curve and the estimated (time-varying) uncertainties as determined in Sec. IV B, thereby accounting for the particular DM value and its uncertainty as determined by the Gaussian learning process. The entries in Table IV are the mean of the resulting distributions of each parameter after a few thousand MC runs. The uncertainties are the standard deviation of the distributions or the maximum Tempo error encountered in all MC runs, whichever is larger.

In order to allow direct comparisons with previous publications (especially Ref. [31]), parameters in Table IV are measured within the timescale known as “Barycentric Dynamical Time” (TDB) as implemented in Tempo. TDB runs at a slower rate than the “Barycentric Coordinate Time” (TCB), which is recommended by IAU 2006 Resolution B3 [162]. This choice does not have any consequences for the gravity tests or discussions presented below, as all (dimensionful) parameters determined from TOAs measured using TDB are multiplied by a constant factor, which either drops out or is (still) too small to be relevant for the discussion of masses or PK parameters. In order to transfer from TDB to TCB, parameters with units of time shown in Table IV need to be divided by $\kappa = (1 - 1.550519768 \times 10^{-8})$ and the values adjusted accordingly [25,162,163]. As for the astrometric timing, the transfer of the TOAs from the topocentric to the barycentric reference frame is made using the DE436 Solar-System ephemeris.

While our timing results are perfectly consistent with those based on 2.5 years of timing data presented earlier [31], the increased length and density of our dataset leads to unprecedented precision in the measured parameters. For instance, the orbital period is measured with a precision of 86 ns. Most importantly, the Keplerian and PK parameters reach a precision that leads to a very significant improvement in our ability to conduct precision tests of strong-field gravity, including radiative and light-propagation aspects, as shown in the following sections. The need to fit up to the fourth spin frequency derivative (cf. Sec. V) reflects, on one hand, the exceptional duration and density of our dataset but also indicates a certain degree of timing noise, but at a level that is consistent with other pulsars of this age [103,165].

The observed improvements in PK parameters are in line with the expectation based on our earlier measurements [33], although a detection of the relativistic deformation of the orbit, described by PK parameter $\delta_\theta$, has occurred somewhat earlier than predicted. A consequence is a slightly smaller improvement in the precision of PK parameter $\gamma_E$, due to a correlation between the two parameters (see Sec. VI B 4). The precision in the measurement of the periastron advance, PK parameter $\dot{\omega}$ (or $k$), has improved beyond the level of 2PN contributions, making it necessary for the first time to include considerations of the Lense-Thirring effect and EOS of superdense matter in our
analysis (see Sec. VIB 3). Similarly, we have to take higher-order effects into account when studying light propagation in the double pulsar system (see Sec. VIB 5). Here, we point out that the derived orbital inclination angle $i = 89^{\circ}.35 \pm 0^{\circ}.05$ (or $180^{\circ} - i$) is consistent with the value derived in our earlier timing analysis, i.e., $i_{KSM06} = 88^{\circ}.69(+0^{\circ}.50/-0^{\circ}.76)$ [31], but now also in perfect agreement with the value obtained from modeling of the eclipse pattern ($i_{eclipse} = 89^{\circ}.3 \pm 0^{\circ}.1$ [54,166]). In contrast, both values are in tension with measurements based on the changes in the intensity scintillation pattern as a function of orbital phase ($i_{scint} = 88^{\circ}.1 \pm 0^{\circ}.5$ [92]). The latter method has the potential to resolve the $i \rightarrow 180^{\circ} - i$ ambiguity caused by the Shapiro delay's dependency on $\sin i$. Ongoing MeerKAT measurements promise to improve on the scintillation results to resolve this difference [154].

## B. Masses and tests of general relativity

### 1. Mass measurement

The standard procedure to obtain the masses of a ("clean") relativistic binary pulsar system is the measurement of two PK parameters. Under the assumption of GR, and using the well-measured Keplerian parameters, one can calculate the two a priori unknown masses [4]. Alternatively, one can fit the DDGR model [30], which is based on GR and has the companion and total masses as free parameters. For the double pulsar, the situation is more complicated, since the most precisely known PK parameter, $\dot{\omega} \equiv n_b k$, has a contribution from the Lense-Thirring effect due to the rotation of pulsar A, which is more than 30 times larger than the measurement error [see Eq. (17) and $\dot{\omega}$ in Table IV]. However, a calculation of the Lense-Thirring contribution requires the knowledge of the MOI of pulsar A, $I_A$, which comes with a significant uncertainty due to our imperfect knowledge of the EOS for NS matter (cf. Refs. [119,167]). Figure 5 demonstrates how the masses, determined from $k$ and the Shapiro shape $s$ (the second-most precise PK parameter) depend on $I_A$. To constrain $I_A$, one can use constraints on the EOS obtained from the double-NS merger GW170817 [168]. As a first approximation, we use the limits on the radius of a NS obtained in Ref. [112] and convert them into a probability distribution for $I_A$ with the help of the radius-MOI relation for A given in Ref. [114]. Doing so, one obtains from Eq. (17)

$$ \dot{\omega}^{\text{LT,A}} = -4.83(+29/-35) \times 10^{-4} \text{ deg yr}^{-1}, \quad (35) $$

and in combination with the PK parameter $s$ one then gets for the (Doppler-shifted) masses [104]

$$ m_A = 1.338185(+12/-14) M_{\odot}, \quad (36) $$

$$ m_B = 1.248868(+13/-11) M_{\odot}, \quad (37) $$

$$ M = 2.587052(+9/-7) M_{\odot}. \quad (38) $$

In these calculations, for both of the PK parameters $k$ and $s$, we account for NLO contributions as given in Sec. V. While for $k$ this 2PN correction is highly significant (approximately 35$\sigma$), for $s$ the 1PN term amounts to a correction of only about $1\sigma$. It is interesting to note, while to leading order a measurement of $k$ can directly be converted into a measurement of the total mass $M$, this is no longer the case once 2PN and Lense-Thirring corrections need to be accounted for [cf. Eqs. (9) and (14)].

Concerning the precision obtained in Eqs. (36)–(38), one has to keep in mind that, strictly speaking, they need to be rescaled with the unknown Doppler factor $D$ of Eq. (1), in order get the actual (intrinsic) inertial masses [29,35]. [170] To calculate $D$, one would need the radial velocity $v_R$ of the double pulsar system with respect to the SSB, which is unknown. Given the small transverse velocity $v_T$ (see Table IV and also Sec. VIII C), one would naturally expect $v_R \sim 11 \text{ km s}^{-1}$ (cf. Fig. 20), which converts into a mass uncertainty of approximately $3 \times 10^{-5} M_{\odot}$. It is important

<figure><img src="image_8.png" /><figcaption>FIG. 5. Dependence on the masses of the double pulsar system on the MOI of pulsar A, $I_A$. Mass calculations are based on GR and the two PK parameters $k(m_A, m_B, I_A)$ [Eq. (13)] and $s(m_A, m_B)$ [Eq. (19)]. The figure shows the total mass $M = m_A + m_B$ (upper) and the masses of A (middle) and B (lower) as a function of $I_A$. The blue bands indicate the 68.3% confidence range in the masses for each given $I_A$. The gray bands indicate the likely MOI range, as suggested by the radius limits of Ref. [112] (90% confidence), when using the radius-MOI relation of Ref. [114]. The red areas are excluded by the causality condition for the EOS [169].</figcaption></figure>
to note that in gravity tests the unknown factor D drops out and is, therefore, irrelevant [29,35].

The PK parameter $\gamma_E$ (within the approximation required) is not affected by the uncertainty in the EOS. However, a mass determination based on s and $\gamma_E$ yields considerably larger errors in the masses: $m_A = 1.3393(19)$ M$_{\odot}$ and $m_B = 1.2494(9)$ M$_{\odot}$. The same is true for a mass determination based on s and R [the latter taken from Ref. [31], where $R \equiv m_A/m_B = 1.0714(11)$ is given]: $m_A = 1.3379(29)$ M$_{\odot}$ and $m_B = 1.2487(14)$ M$_{\odot}$. Note that a mass measurement using $\dot{P}_b$ is to some extent also EOS dependent, due to the spin-down mass loss of A [see Eq. (20)], although currently this is still at a negligible level. Combining s and $\dot{P}_b$ (from Sec. VI B 2) gives $m_A = 1.338176(68)$ M$_{\odot}$ and $m_B = 1.248842(34)$ M$_{\odot}$, which is obviously less precise than the masses in Eqs. (36) and (37). In summary, in spite of the uncertainty in the EOS, s and k give by far the best estimate for the double pulsar masses.

## 2. Gravitational-wave emission

With the (GR) masses (36) and (37), the binary system is fully determined, and any further PK parameter can be used for a test of GR. Of particular importance is the test of GW damping. In their backreaction on the system, GWs extract orbital energy and angular momentum from the binary motion, leading to a secular change of the Keplerian parameters $P_b$, x, and e. Temporal changes in the latter two are still well below the measurement precision, which in terms of x means that the shrinkage of the orbit due to GW damping is not yet seen directly through a change in the size of the orbit. The (observed) change of the orbital period $\dot{P}_b$, on the other hand, is highly significant, measured with a fractional precision of $6 \times 10^{-5}$ (see Table IV). This is nearly an order of magnitude better than for the Hulse-Taylor pulsar [172], and, as in that system, it is dominated by the GW damping.

For a GW test, one needs the intrinsic change of the orbital period, $\dot{P}_b^{\text{int}}$. Therefore, external contributions $\dot{P}_b^{\text{ext}}$, resulting from differential accelerations between the SSB and the double pulsar system, need to be subtracted from the observed orbital period change [122]. The Shklovskii effect [129] leads to an apparent acceleration due to the transverse motion of the double pulsar system and is given by

$$ \dot{P}_b^{\text{Shk}} = \frac{d(\mu_\alpha^2 + \mu_\delta^2)}{c} P_b = 1.72(15) \times 10^{-16}. \quad (39) $$

The distance d is calculated from the weighted mean probability distribution in Fig. 23. The differential acceleration in the Galactic gravitational potential leads to a correction of similar magnitude. If one takes the potential of Ref. [130] to calculate the Galactic gravitational accelerations of the SSB and the double pulsar, $\tilde{g}_{\text{SSB}}$ and $\tilde{g}_{\text{DP}}$, respectively, then one finds

$$ \begin{aligned} \dot{P}_b^{\text{Gal}} &= \frac{\vec{K}_0 \cdot (\vec{g}_{\text{DP}} - \vec{g}_{\text{SSB}})}{c} P_b \\ &= -3.40(+25/-24) \times 10^{-16}, \end{aligned} \quad (40) $$

where $\vec{K}_0$ is a unit vector directed from the SSB to the double pulsar. Instead of the potential in Ref. [130], one can also use a more analytic approach for the Galactic correction, like in Refs. [122,173,174], which in combination with the latest Galactic parameters based on the results in Refs. [131,175,176] gives very similar numbers. Other external contributions to $\dot{P}_b$, as, for instance, discussed in the context of the Hulse-Taylor pulsar in Ref. [122], are still negligible (see also the discussions in Refs. [110,177]). Combining Eqs. (39) and (40), while accounting for correlations, gives

$$ \dot{P}_b^{\text{ext}} = \dot{P}_b^{\text{Shk}} + \dot{P}_b^{\text{Gal}} = -1.68(+11/-10) \times 10^{-16} \quad (41) $$

and, consequently,

$$ \dot{P}_b^{\text{int}} = \dot{P}_b - \dot{P}_b^{\text{ext}} = -1.247752(79) \times 10^{-12}. \quad (42) $$

At this point, it is important to note that the error in the calculated $\dot{P}_b^{\text{int}}$ is still clearly dominated by the error in the observed orbital period change $\dot{P}_b$ (cf. Table IV). In other words, the uncertainty in our knowledge of the external contributions to the observed change in the orbital period [Eq. (41)] is still significantly smaller than the measurement error of the timing parameter $\dot{P}_b$.

The second-largest contribution to $\dot{P}_b^{\text{int}}$, besides GW damping, is related to the mass loss of A as given in Eq. (20). Using the limits of Ref. [112] and the radius-MOI relation in Ref. [114] gives

$$ \dot{P}_b^{\text{int}} = 2.9(2) \times 10^{-17}, \quad (43) $$

which is still a factor of 3 smaller than the error of $\dot{P}_b^{\text{int}}$. Nevertheless, we include the mass loss contributions in our calculations, in particular, since this is relevant in Sec. VI B 3 below, where instead of using a priori constraints on $I_A$ we use $\dot{P}_b^{\text{int}}$ to obtain a probability distribution for the MOI of A. The measured change in the orbital period due to GW emission amounts to

$$ \dot{P}_b^{\text{GW}} = \dot{P}_b^{\text{int}} - \dot{P}_b^{\text{int}_A} = -1.247782(79) \times 10^{-12}, \quad (44) $$

meaning that GW damping in the double pulsar is measured with a fractional precision of $6 \times 10^{-5}$. This improves the last published value for the double pulsar [31] by more than a factor of 200. This is also about a factor of 25 more precise than in the Hulse-Taylor pulsar [172], which is presently limited by the uncertainties in the correction for the external effects.
The high precision in Eq. (44) leads to important tests of the radiative properties of gravity theories. In Sec. VII, we use this for constraining two alternatives to GR. In this section, we provide only a test for GR. In GR, to leading order, $\dot{P}_b^{GW}$ is given by the quadrupole formula [124], which corresponds to the 2.5PN contributions in the equations of motion. Using the masses in Eqs. (36) and (37), one obtains

$$ \dot{P}_b^{GW,GR(2.5PN)} = -1.247810(+6/-7) \times 10^{-12}. \quad (45) $$

The next-higher correction (3.5PN in the equations of motion) is calculated in Ref. [128] and amounts to

$$ \dot{P}_b^{GW,GR(3.5PN)} = -1.75 \times 10^{-17}, \quad (46) $$

which is larger than the error in Eq. (45), and, therefore, we include it in our calculation for the predicted value:

$$ \begin{aligned} \dot{P}_b^{GW,GR} &= \dot{P}_b^{GW,GR(2.5PN)} + \dot{P}_b^{GW,GR(3.5PN)} \\ &= -1.247827(+6/-7) \times 10^{-12}. \end{aligned} \quad (47) $$

This value is in perfect agreement with the observed value given in Eq. (44), leading to a GR test of

$$ \dot{P}_b^{GW}/\dot{P}_b^{GW,GR} = 0.999963(63). \quad (48) $$

This is by far the most precise test of GW emission, about a factor of 25 better than in the Hulse-Taylor pulsar [172]. The agreement with GR is further demonstrated in Fig. 6, which shows the cumulative shift in periastron time due to the acceleration in the orbital phase evolution as a result of $\dot{P}_b$ [cf. Eq. (5)]. Since this is the most precise of the different tests one obtains from the double pulsar, one can say that GR is in agreement with the double pulsar at a level of $1.3 \times 10^{-4}$ with 95% confidence. We emphasize that the result (48) is calculated in a Monte Carlo run, where we simultaneously randomize the observed parameters and $I_A$ [probability distribution like for Eq. (43)] in order to calculate the masses (from k and s) and all the $\dot{P}_b$ contributions.

We also test the dependence of result (48) on our assumptions about the NS EOS. In a separate Monte Carlo run, instead of the limits in Ref. [112] we assume a uniform distribution for the MOI, in the (rather extreme) range of $I_A^{(45)} = 0.875$ (causality limit [169]) and $I_A^{(45)} = 1.981$ ($R_A = 15$ km) [114,167]. The upper limit with its large radius is in clear tension with, e.g., LIGO/Virgo observations [168]. Nevertheless, even for such a distribution for $I_A$ we obtain $\dot{P}_b^{GW}/\dot{P}_b^{GW,GR} = 0.999958(64)$, which shows that Eq. (48) depends only very weakly on the EOS uncertainty.

<figure><img src="image_9.png" /><figcaption>FIG. 6. Cumulative shift of the times of periastron passage relative to a nondissipative model. Each data point covers a time span of 60 d. To each of these subsets, we fit a Keplerian orbit optimizing only for the orbital period and the time and angle of periastron passage. We include the advance of periastron, light-propagation effects, and the Einstein delay in the orbital model but keep the values fixed to those in Table IV. We plot the difference between the measured time of periastron and a periastron time near its discovery (i.e., MJD 52759.89, or 2003.33). The red curve in the top is the GR prediction based on the masses in Sec. VI B 1. The bottom shows the deviation from this prediction, characterized by a normalized $\chi^2 = 0.76$.</figcaption></figure>

In the following, we use the accelerated phase evolution in the orbital motion as shown in Fig. 6 to update the PN parameter test introduced by the LIGO/Virgo Collaboration in Ref. [178]. In particular, it is the tests from the early-inspiral stage, where the phase evolution of the merger of compact objects is analytically described by the PN approximation [123], that can be used in such a comparison with the GW damping measured in the double pulsar. Most importantly, the measurement precision of $\dot{P}_b$ has greatly improved compared to Ref. [31], which are the data that are used in Ref. [178]. Figure 7 gives an updated version of Fig. 6 in Ref. [178] (including the new results [179] from Refs. [180,181]). The bounds are plotted for the different
PN levels, allowing for possible GR violations at different PN levels (i.e., different powers of frequency), one at a time. Note that Fig. 7 uses the “relative” PN order in the radiation reaction (i.e., PN order beyond the Einstein quadrupole formula), where the leading order, i.e., 0PN, occurs at the 2.5PN order in the binary equations of motion (see, e.g., Ref. [123] for a detailed discussion). Because of the many orbits since 2003 (approximately 60 000), which can be tracked with high precision in a phase-coherent timing solution, the double pulsar leads to considerably tighter constraints at low PN orders, whereas it becomes very quickly less constraining for higher PN orders, due to its comparatively small velocity ($v \sim 0.002c$).

While Fig. 7 certainly serves as a comparison on how much a given PN parameter of the inspiral phase evolution can (each at a time) deviate from its GR value in the different experiments, that figure has to be taken with a grain of salt when it comes to interpreting these bounds as limits on deviations from GR predicted by alternative theories of gravity. First, such a comparison mixes tests from two different types of compact objects, i.e., NSs and BHs, which might behave quite differently depending on

<figure><img src="image_10.png" /><figcaption>FIG. 7. Update of Fig. 6 in Ref. [178] (including data from Refs. [180,181]), which shows the 90% upper bounds on the absolute magnitude of the GR violation parameters $\delta\hat{q}_i$, from 0PN through 3.5PN (“relative” order) in the inspiral phase (see, e.g., Ref. [182] for the definition of the PN phase coefficients and Ref. [178] for further details on the method). As discussed in Ref. [178], the 0.5PN parameter is zero in GR and, therefore, understood not as a relative but as an absolute shift. Black circles show the combined limits from the double BH mergers, blue squares are the limits from the double-NS merger GW170817, and red triangles give the limits derived from the double pulsar GW test in this paper. The PN order on the x axis is in the GR radiation reaction, where the leading contribution (0PN) corresponds to the dissipative 2.5PN term in the equations of motion. Note that such a comparison of tests with different compact objects (BHs vs NSs) as well as different gravity regimes (mildly relativistic vs highly relativistic strong field) does come with a caveat, which is explained in more detail in the text.</figcaption></figure>

how GR is broken. Hence, constraints from experiments with material bodies might not apply to BH dynamics and vice versa. Particularly obvious cases are alternative theories where BH binaries behave like in GR (e.g., Ref. [183]) or alternative theories where NSs do not carry any scalar charge, while BHs do [184]. Second, the double pulsar tests a different gravity regime (mildly relativistic strong field) compared to the GW merger events (highly relativistic strong field). For instance, the double pulsar test would generally be insensitive to modifications of GR that lead only to short-range effects (e.g., Refs. [185,186]); see also Ref. [187]. Nevertheless, at least to some extent, such a comparison illustrates the complementarity of binary pulsar experiments and merger observations by GW detectors, as long as one keeps in mind the qualitative differences of the various experiments, which are closely linked to the details of a given theory of gravity.

### 3. Lense-Thirring effect and equation of state

In Sec. VI B 1, we use constraints on the MOI of pulsar A, $I_A$, derived from the multimessenger analysis in Ref. [112], in order to obtain the best mass estimates for the double pulsar, as given in Eqs. (36)–(38). In this section, at first, we ignore any existing constraints on the EOS of NSs and simultaneously determine $m_A$, $m_B$, and $I_A$, following the procedure outlined in Ref. [110]. As in Sec. VI B 1, we assume GR to be the correct theory of gravity and use the three best PK parameters to simultaneously calculate the individual masses of the double pulsar and the MOI of A. From the calculations in Sec. VI B 1, it is already obvious that the combination of the PK parameters $k$, $s$, and $\dot{P}_b$ is expected to give by far the best results. In a way, we use $s$ and $\dot{P}_b$ to determine the masses $m_A$ and $m_B$ and then use $m_A$ to extract $I_A$ from the observed advance of periastron $k^{\text{obs}}$ (see $\dot{\omega} \equiv n_B k$ in Table IV), a procedure already proposed for the double pulsar in Ref. [33]. In practice, the calculations are slightly more complicated, as $\dot{P}_b$ also has a contribution proportional to $I_A$ [see Eq. (20)]. Although that contribution is still smaller than the error in $\dot{P}_b$, we nevertheless account for it and follow the procedure in Ref. [110], i.e., calculate $m_A$, $m_B$, and $I_A$ by simultaneously solving the three equations $k^{\text{obs}} = k(m_A, m_B, I_A)$, $s^{\text{obs}} = s(m_A, m_B)$, and $\dot{P}_b^{\text{int}} = \dot{P}_b(m_A, m_B, I_A)$. By this, we obtain probability distributions for the double pulsar masses and the MOI of pulsar A. For the MOI, we find $I_A < 3.0 \times 10^{45} \text{ g cm}^2$ with 90% confidence. Figure 8 compares our result with those derived from the GW170817 LIGO/Virgo merger and from NICER x-ray timing. Using a universal relation, like the one in Ref. [114], one can convert the probability distribution of $I_A$ into a probability distribution for A’s radius. With 90% confidence, this gives an upper limit for A’s radius of 22 km, a value outside any physically valid EOS and clearly exceeding the range used in Ref. [114].
<figure><img src="image_11.png" /><figcaption>FIG. 8. Probability distribution for the MOI $I_A$ of A derived from the $k-s-\dot{P}_b$ test. The vertical dashed red line indicates the upper bound with 90% confidence. $I_A > 0$ is assumed as a prior. The light gray band shows the 90% credible interval one obtains with the limits from Ref. [112] using the radius-MOI relation of Ref. [114]. As a comparison, the horizontal black lines indicate different 90% ranges derived from (top to bottom) tidal-deformability constraints from GW170817 [117], Bayesian modeling of a range of EOSs [116], and two different constraints from NICER observations [118]. The red area is excluded by the causality condition for the EOS [169].</figcaption></figure>

While these numbers cannot compete with those obtained from LIGO/Virgo and NICER observations (see, e.g., Refs. [112,116–118]), they show that the double pulsar constraints are narrowing in on realistic values for NS radii. This is expected to improve considerably over the next years (see Hu et al. [110]). For the masses, we get $m_A = 1.338183(+78/-52) M_\odot$ and $m_B = 1.248869(+38/-27) M_\odot$. These masses are clearly more uncertain than in Eqs. (36) and (37), but they do not require any assumption about the EOS for matter at supranuclear densities.

Instead of using the Lense-Thirring effect to constrain the MOI, one can conversely use existing constraints on the NS EOS to test the contribution of spin-orbit coupling to the precession of the pulsar orbit (cf. Sec. VII in Ref. [110]). For this, we introduce a scaling factor $\lambda_{LT}$ for the Lense-Thirring part in Eq. (14), i.e.,

$$k = k^{\text{1PN}} + k^{\text{2PN}} + \lambda_{\text{LT}} k^{\text{LT,A}}, \quad (49)$$

For GR, one has $\lambda_{LT} = 1$. Similar to the procedure above, we can now calculate $m_A$, $m_B$, and $\lambda_{LT}$ by simultaneously solving the three equations $k^{\text{obs}} = k(m_A, m_B, \lambda_{LT}|I_A)$, $s^{\text{obs}} = s(m_A, m_B)$, and $\dot{P}_b^{\text{int}} = \dot{P}_b(m_A, m_B|I_A)$, where this time in our Monte Carlo runs $I_A$ is randomly chosen from the probability distribution based on the EOS constraints in Ref. [112] (like in Sec. VI B 1, when determining the masses from $k$ and $s$). In a sense, we determine the masses of A and B from the PK parameters $s$ and $\dot{P}_b$ to extract the

Lense-Thirring contribution form the observed periastron advance $k$, by assuming a range of values for $I_A$. We obtain

$$\lambda_{\text{LT}} = 0.7 \pm 0.9 \quad (68\% \text{ C.L.}), \quad (50)$$

This is consistent with GR, but admittedly not very constraining, at least when compared to the weak-field (test-particle-type) frame-dragging experiments in the Solar System (Lense-Thirring precession of satellite orbits [188] and Pugh-Schiff precession of gyroscopes [189]). Moreover, the above result is not generic, in the sense of the parametrized post-Newtonian (PPN) tests in the Solar System. In order to extract the Lense-Thirring contribution to $k$, we assume that all three, the advance of periastron without Lense-Thirring contribution [i.e., Eq. (9)], the GW damping, and the Shapiro shape, can be calculated from GR. This is generally not expected to be the case in alternative gravity theories. As a result, one would have to implement the equivalent analysis for any other gravity theory, in order to obtain a fully consistent Lense-Thirring test. Given the low precision of the LT test, one may ask if, in view of the other tests with the double pulsar, mostly based on very precisely measured PK parameters, this then yields any additional constraints for the gravity theories under consideration. For the alternative theories discussed in Sec. VII, at least, spin-orbit effects observed in the double pulsar (including $\Omega_B^{\text{spin}}$) do not contribute to the constraints. See Ref. [110] for a more detailed discussion, in particular, on this test in the context of near-field modifications of GR.

Finally, from the analysis outlined above, one can also extract a test of the total advance of periastron $k$, in an $s-\dot{P}_b-k$ test. From this, one obtains $k^{\text{obs}}/k^{\text{GR}} = 1.000015(26)$.

## 4. Relativistic deformation of the orbit

The relativistic deformation of the orbit discussed in context of the Rømer delay (Sec. VA) is detected in our measurements. Similar to the report of a $1.5\sigma$ detection of $\delta_\theta$ for PSR B1913+16 [172], our value is also formally detected only just above the $1\sigma$ level. However, as we demonstrate in Fig. 9, the parameter is, in fact, well constrained and in full agreement with the expectation from GR. The figure also demonstrates a correlation between $\delta_\theta$ and $\gamma_E$. Indeed, due to this correlation, the value derived for $\gamma_E$ in a timing model ignoring the relativistic deformation, i.e., assuming $\delta_\theta = 0$, would lead to a value for $\gamma_E$ that would be inconsistent with GR at the $2\sigma$ level, while both $\gamma_E$ and $\delta_\theta$ are, in fact, in perfect agreement when including $\delta_\theta$ in our timing model (see the white cross in Fig. 9). In order to illustrate the effect of a nonzero $\delta_\theta$ value, we compare a deformed orbit ($\delta_\theta > 0$) with an elliptical orbit ($\delta_\theta = 0$) in Fig. 10 for a hugely exaggerated effect.
<figure><img src="image_12.png" /><figcaption>FIG. 9. Map of $\chi^2$ contours demonstrating a correlation between the PK parameters $\gamma_E$ and $\delta_\theta$. In order to produce this map, the two parameters are held fixed at their grid position, while fitting for all remaining parameters in the timing model as described earlier. The contours indicate 68.3%, 95.4%, and 99.7% confidence levels, respectively. The expected GR value is indicated by the cross with uncertainties too small to be visible on this scale: $\gamma_E = 383.997(5)$ µs and $\delta_\theta = 12.60889(4) \times 10^{-6}$ (masses from k and s; cf. Sec. VIB 1).</figcaption></figure>

## 5. Signal propagation

Already in the discovery paper of A [21], it was suspected that the double pulsar system is seen nearly edge-on. This was confirmed by the discovery of a prominent Shapiro delay in the TOAs of A, as well as by the observation of eclipses of A near its superior conjunction, caused by the plasma-filled magnetosphere of B [22]. The measurement of the Shapiro delay gave access to two additional PK parameters, i.e., the “Shapiro range” $r$ and the “Shapiro shape” $s$ [29] [see Eq. (24)]. By 2006, these two parameters were measured with high precision, and, in particular $s$, as part of an $\dot{\omega}-R-s$ test, provided a $10^{-3}$ (95% C.L.) GR test [31]. Unfortunately, the measurement of the mass ratio $R$ has not improved since then, mostly because of the fact that B disappeared from view in early 2008 as a result of relativistic spin precession [58]. In the meantime, the measurement of the Shapiro delay, and, therefore, $s$, has improved significantly (see Fig. 11), making it the second-most constraining parameter in the mass-mass plane (see Fig. 13 below), which is also obvious from the mass determinations in Sec. VIB 1. Since $\gamma_E$ by now is more constraining than $R$, one can test $s$ as part of a $\dot{\omega}-\gamma_E-s$ test. This leads to

$$s^{\text{obs}}/s^{\text{GR}} = 1.000\,09(18), \qquad (51)$$

<figure><img src="image_13.png" /><figcaption>FIG. 10. Hugely exaggerated illustration of the effect of the relativistic deformation $\delta_\theta$ on the shape of the orbit, which modifies the usually expected elliptical shape. For the red (deformed) orbit, we choose $a = 1$, $e_T = e_r = 0.6$, and $\delta_\theta = 0.5$, while the dashed blue curve has $\delta_\theta = 0$.</figcaption></figure>

meaning a $4 \times 10^{-4}$ (95% C.L.) test in agreement with GR, which is somewhat better than the $\dot{\omega}-R-s$ test of Ref. [31], while mainly dominated by the error in $\gamma_E$. Moreover, from Eq. (3.15) in Ref. [35], it is clear that, for a large range of alternatives to GR, $s$ is linked only to the effective gravitational constant $G_{AB}$ of the orbital dynamics, as $s$ can be identified with the purely geometric quantity $\sin i$.

A test of the Shapiro range $r$ directly probes the interaction between a strongly self-gravitating NS—in our case, pulsar B—and a photon. This can be easily seen within the modified Einstein-Infeld-Hoffmann (mEIH) formalism [35,137], where $r = (1 + \gamma_{B0})G_{B0}m_B/c^3$ [see, in particular, Eq. (3.14) in Ref. [35] with $\epsilon_{B0} = 2\gamma_{B0} + 1$]. $G_{B0}$ is the effective gravitational constant and $\gamma_{B0}$ the strong-field extension of the PPN formalism parameter $\gamma^{\text{PPN}}$ (see, e.g., Ref. [41]), which enter the leading-order spacetime metric $g_{\mu\nu}$ of pulsar B. In GR, we have $G_{B0} = G$ and $\gamma_{B0} = 1$. Since $r$ is by far the most uncertain of the PK parameters, apart from $\delta_\theta$ and $\Omega_B^{\text{spin}}$, there are several combinations of PK parameters that can be used in the test. In other words, the agreement of $m_B^{(r)} \equiv rc^3/G = 1.2510(43)$ M$_\odot$ (cf. Table IV) with the companion mass determined in Sec. VIB 1 gives

$$r^{\text{obs}}/r^{\text{GR}} = 1.0016(34), \qquad (52)$$

leading to a $7 \times 10^{-3}$ (95% C.L.) test of $r$ in agreement with GR.

It is interesting to note that, near superior conjunction, the signals of the pulsar propagate through a region with a
<figure><img src="image_14.png" /><figcaption>FIG. 11. Signature of the Shapiro delay in the timing residuals, plotted against A's angular orbital position $\psi = \omega + \theta$ with respect to the ascending node. The bottom shows the aggregated residuals as a function of orbital phase when not including the Shapiro delay parameters $z_s$ (or, equivalently, $\sin i$) and $r$. The shape differs from the expected Shapiro delay curve shown in the top, as some signal power is absorbed in a fit for the orbital parameters, namely, the Rømer delay (see, e.g., Ref. [98]). The expected shape is recovered by computing the residuals when setting the Shapiro delay parameters to zero in the final timing solution. Note that the shown residuals are corrected for a slow variation in the curve due to the relativistic orbital precession. The red curves show the GR expectation using values of $\sin i$ and $r$ derived from $\dot{\omega}$ and $\gamma_E$. These are $(\sin i)^{GR} = 0.99988(18)$ and $r^{GR} = 6.1518(11) \mu s$, where the uncertainties are determined by those in $\dot{\omega}$ and $\gamma_E$. Note that the uncertainty in the prediction of $\sin i$ is correspondingly larger than that of the actually measured $\sin i$.</figcaption></figure>

spacetime curvature (as, for instance, measured by the Kretschmann scalar) that is many orders of magnitude stronger than in other experiments that test photon propagation in gravitational fields and, therefore, the coupling between gravitational and the electromagnetic fields. This is of interest, for instance, in the presence of interactions between the electromagnetic field and curvature tensor as studied in Ref. [190].

As it turns out, the leading-order expression (24) no longer provides a complete description of the signal

<figure><img src="image_15.png" /><figcaption>FIG. 12. Aggregated residuals (blue) due to NLO contributions in the Shapiro delay and the aberration, plotted against A's angular orbital position $\psi = \omega + \theta$ with respect to the ascending node. The red curve shows the theoretical prediction, and the black curve corresponds to the fitted $q_{NLO}$ (see Table IV), with a $2\sigma$ range indicated by the light gray band. Residuals are rescaled by $(1 + e_T \cos \theta)^{-1}$ to account (to leading order) for a secular variation in the amplitude due to the precession of $\omega$. A data point generally results from approximately 3000 aggregated residuals. The data point at $\psi = 90^\circ$ (i.e., superior conjunction of A) has a much larger uncertainty, as this bin coincides with the 30-s-long eclipse of pulsar A, resulting in both many fewer TOAs available for aggregation and residuals of opposite sign being averaged.</figcaption></figure>

propagation in the double pulsar. Such a model leads to significant residuals near conjunction, as can be seen in Fig. 12. These residuals are fully accounted for by the expected NLO contributions discussed in Secs. VC and VD (red curve in Fig. 12). To test the significance of the NLO corrections in $\Delta_S$ and $\Delta_A$, we scale them collectively with a factor $q_{NLO}$ [cf. Eqs. (32) and (34)], which is unity in GR. As discussed in Sec. VE, the similarity of the NLO contributions to $\Delta_S$ and $\Delta_A$ makes it impossible to test these two contributions separately. For the common factor, which also scales the NLO contribution from lensing [see Eq. (33)], we find

$$q_{NLO} = 1.15(13), \qquad (53)$$

which can be interpreted as an approximately $8\sigma$ measurement of NLO contributions in the signal propagation, in agreement with GR. This clearly shows that the timing of the double pulsar has by now reached a precision where NLO contributions to the Shapiro and aberration delay have to be taken into account.

As as consequence of its definition in Sec. VE, a fit for $q_{NLO}$ combines two aspects of gravity, 1.5PN corrections to
<figure><img src="image_16.png" /><figcaption>FIG. 13. Mass-mass diagram for the double pulsar based on GR, for the six PK parameters $\dot{\omega} \equiv n_b k$, $\gamma_E$, $\dot{P}_b$, $r$, $s$, $\Omega_B^{spin}$, and the mass ratio $R$. The width of each curve indicates the measurement uncertainty of the corresponding parameter. The seventh PK parameter $\delta_{\theta}$ is not shown here, since its limits still lie outside the mass ranges shown. For the solid black $\dot{\omega}$ line, $I_A^{(45)} = 1.32$ is used, which corresponds to a NS radius of 12 km (cf. Ref. [112]). The inset is an expanded view of the region of principal interest, where only the four best PK parameters are shown. To illustrate the influence of the LT effect, we draw in the inset a dashed black $\dot{\omega}_0$ line where the LT contribution is ignored [see Eq. (9)]. The gray band indicates the range for the $\dot{\omega}$ line under the variation of $I_A$, from the causality limit of Ref. [169] (left border) to a NS radius of 15 km (right border). The conversions between radius and MOI are based on the relation in Ref. [114]. There is a small $I_A$ dependence of the $\dot{P}_b$ curves, which is ignored. The intersection of all line pairs is consistent with a small region that corresponds to the masses of A and B.</figcaption></figure>

the Shapiro delay related to the motion of the masses [Eq. (32)] and corrections related to the deflection of the signal beam in the gravitational field of the companion [Eqs. (33) and (34)]. Although, as discussed in Sec. V E, these contributions cannot be tested separately in a simultaneous fit, in a phenomenological approach one can still test for one at a time while keeping the other one fixed. If

we apply the scaling factor $q_{NLO}$ only to the deflection-related contributions, one finds

$$q_{NLO}[\text{deflection}] = 1.26(24). \qquad (54)$$

This limit comes solely from the NLO aberration contributions [Eq. (34)], since a rescaling of $\delta\Lambda_u^{len}$ is covariant
with s. Given the large uncertainty, the above test is admittedly of little interest in all those theories or frameworks where deviations from GR in the signal deflection (caused by a strongly self-gravitating mass) are already well constrained by the measurement of the Shapiro delay PK parameters r and s, i.e., Eqs. (30) and (31) [191].

Alternatively, one can keep the deflection contributions fixed, i.e., $q_{NLO} = 1$ in Eqs. (33) and (34), and fit only for the retardation effect. In this case, one obtains

$$q_{NLO}[\text{retardation}] = 1.32(28). \qquad (55)$$

This result clearly demonstrates the significance of the motion of the companion mass while the signal of A propagates across the binary system on its way toward the observer.

Finally, when assuming that the signal deflection is sufficiently well described by GR, then Eq. (54) provides independent evidence (besides Refs. [52,57]) for the fact that the spin of A is aligned with the orbital angular momentum and not antialigned [$D \to -D$ in Eq. (30)] [193].

## 6. Mass-mass diagram

The PK parameters, once external contributions are removed and a value for $I_A$ is chosen, depend only on the a priori unknown masses of pulsars A and B and the well-measured (and, hence, determined) Keplerian parameters. Hence, each PK parameter defines a curve in a two-dimensional mass plane displaying their dependence on $m_A$ (plotted on the X axis) and $m_B$ (plotted on the Y axis) for a given theory of gravity. If the theory of gravity used to describe the dependence of the n curves on the masses is consistent with the experimental data, and if no other effects are present, then all curves will meet in a single point. As two curves each define one intersection point, n - 2 curves have still the potential to miss the intersection point and to falsify the assumed theory. This “mass-mass” diagram, as shown for GR in Fig. 13, is a graphical representation of the gravity tests using PK parameters described in Sec. II, with $\dot{P}_b$ corrected as described in Sec. VI B 2. In the case of the double pulsar, the mass ratio adds yet another (largely theory-independent) curve to this diagram in Fig. 13. The thickness of the lines represents the uncertainty in the measured PK parameters (68% C.L.). For all PK parameters but $\Omega_B^{spin}$, which is determined from eclipse modeling (see Sec. II), we have reached a point where we need to enlarge to show the intersection point in more detail to recognize the relative uncertainties. While the large panel in Fig. 13 shows most of the mass range of interest for NSs, the inset shows such an enlarged area, where we choose not to show R, r, and $\Omega_B^{spin}$ for clarity. It is obvious that all lines are consistent with a single intersection point given by the masses of Eqs. (36) and (37). The inset, however, highlights the fact that we clearly need to

take the Lense-Thirring contribution into account (see discussion in Sec. VA 2). With n = 7 constraints (i.e., the six PK parameters measured here and the mass ratio R), we can conduct n - 2 = 5 independent tests of theories of gravity, which are all passed with flying colors as shown above [194].

On a final note, in the inset in Fig. 13, one can see that the measured $\gamma_E$ nicely agrees with the intersection of the s curve with either the $\dot{\omega}$ curve or the $\dot{P}_b$ curve. One can convert this into a test of $\gamma_E$:

$$\gamma_E/\gamma_E^{GR} = 1.000\,12(25). \qquad (56)$$

The main interest in such a test lies in the fact that it tests an invariance of the local gravitational constant, which is a specific aspect of local position invariance. We discuss this in more detail in the following section, in the context of two specific alternatives to GR.

## VII. ALTERNATIVE GRAVITY THEORIES

In the previous section, we demonstrate that Einstein's GR is in perfect agreement with current observations of the double pulsar. This shows that, in contrast to Solar-System tests which are in the weak-gravity regime, GR also is an accurate description of gravitational interactions for two strongly self-gravitating bodies in a mildly relativistic strong-gravity regime. However, since GR was devised in 1915, many other theories of relativistic gravity have been proposed. In this section, we consider these alternative theories. Since it is impossible to give an even remotely complete discussion, we focus on two examples. The first, known as “Damour–Esposito-Farèse” (DEF) gravity, is a two-parameter class of theories, whereas the second, known as “Tensor-Vector-Scalar (TeVES)” gravity, represents the class of “Modified Newtonian Dynamics” or MONDian gravitational theories.

These theories are chosen since they are especially suitable for studying the deviations one typically expects when going beyond GR by breaking the SEP and how they are being constrained by the different PK parameters measured in the double pulsar. In this, these theories provide a theory-based approach for putting different GR experiments into context and allow for a better understanding which structures and symmetries of GR are actually being probed by pulsars and, more specifically, the double pulsar (see Sec. V in Ref. [63]).

Before entering the theory-specific discussion, we give a general statement about constraints on deviations from GR in the radiative regime that can come from the precise measurement of the GW emission [see Eq. (48)]. In alternative theories, one often has gravitational dipolar radiation due to a violation of the SEP as a significant contribution to GW damping. In the equations of motion, that contribution enters at the 1.5PN level [$O(\beta_O^3)$]; see, e.g., Ref. [183]]. If we parametrize a deviation from GR led by
dipolar GW damping by $\dot{P}_b^{GW} = \dot{P}_b^{GW,GR}(1 + B_D\beta_O^{-2})$, then we find $B_D \lesssim 4 \times 10^{-10}$ (95% C.L.). This is nearly 5 orders of magnitude more constraining than the limit from the double-NS merger GW170817 [195]. Note, however, that the two different limits are coming from different gravity regimes (see the discussion at the end of Sec. VI B 2). How such a limit on $B_D$ converts into actual constraints on a given gravity theory depends on the specifics of that theory. In many theories, it is expected that, due to the similarity of the two masses in the double pulsar, the limit is actually weaker than one would assume directly from the tight limit on $B_D$.

A violation of the SEP not only gives rise to dipolar GWs, it generally also results in a temporal modification of the (effective) gravitational constant due to the expansion of the Universe (see, e.g., Ref. [196]), by this having an addition effect on the orbital period. The orbital period change due to a time-varying gravitational constant is given by $\dot{P}_b^{\tilde{G}}/P_b = -2(\dot{G}/G)\mathcal{F}_{AB}$, where $\mathcal{F}_{AB}$ accounts for corrections related to the strong gravitational fields of the NSs in the double pulsar system [196–198]. [199] For weakly self-gravitating bodies, as in the Solar System, $\mathcal{F}_{AB} \approx 1$. $\mathcal{F}_{AB}$ depends on the details of the gravity theory under consideration, as well as on the EOS. While for Jordan-Fierz-Brans-Dicke gravity $\mathcal{F}_{AB} \approx 0.6$, for some theories $\mathcal{F}_{AB}$ can be significantly larger than one (cf. discussion in Ref. [27]). If one assumes that $\dot{P}_b^{\tilde{G}}$ is the only significant non-GR contribution to the orbital period change, then from Eq. (48) one finds that $|\dot{G}/G| < -0.8(14) \times 10^{-13} \text{ yr}^{-1}/\mathcal{F}_{AB}$. This limit is about a factor of 2 better than the currently best pulsar limit [201], although the exact factor is theory dependent. In cases where $\mathcal{F}_{AB}$ is large, this limit can even exceed the best limit obtained in the Solar System (i.e., $|\dot{G}/G| < 4 \times 10^{-14} \text{ yr}^{-1}$ [202]). At this stage, these considerations need to be taken with a grain of salt. As stated above, theories with non-vanishing $\dot{G}$ are expected to also predict dipolar GWs, which, in turn, leads to an additional non-GR modification of $\dot{P}_b$. As a conclusion, in order to perform a $\dot{G}$ test with the double pulsar, one has to do a fully consistent analysis within a certain class of gravity theories and combine it with other binary pulsars (similar to the method proposed in Ref. [174]). However, such an analysis goes beyond the scope of this paper and will be presented in a future publication.

## A. Damour–Esposito-Farèse Gravity

The first alternative to GR which we confront with our double pulsar results is the two-parameter monoscalar-tensor gravity $T_1(\alpha_0, \beta_0)$ introduced in Ref. [40] (“DEF gravity”). This two-parameter class of theories exhibits various effects one typically expects from alternatives to GR, including genuine strong-field effects related to the self-gravity of NSs. Hence, $T_1(\alpha_0, \beta_0)$ can be viewed as an

ideal minispace of well-studied gravity theories in which GR is embedded ($\alpha_0 = \beta_0 = 0$). They are well suited to contrasting the predictions of GR with the predictions of alternatives and allowing for the comparison of GR tests across different gravity regimes [63]. The constants $\alpha_0$ and $\beta_0$ measure, respectively, the linear and the quadratic coupling strength between matter and the scalar field. For simplicity, we assume that the potential of the scalar field can be ignored on the typical scales for the double pulsar test, while it could still be of importance on Galactic or cosmological scales. We further assume that we can neglect a temporal change in the background scalar field, which among others would lead to a long-term variation of the gravitational constant. For sufficiently negative $\beta_0$, DEF gravity shows genuine strong-field effects in binary pulsars, related to “spontaneous scalarization” [40,42]. The presence of the scalar field leads to a modification of the PK parameters in different ways [42,203].

Already at the Newtonian level, the body-dependent effective scalar couplings $\alpha_A$ and $\alpha_B$ enter via the (effective) gravitational constant $G_{AB} = G_*(1 + \alpha_A\alpha_B)$, where $G_*$ is the bare gravitational constant, linked to Newton’s gravitational constant $G$ (as measured in a Cavendish-type experiment) via $G = G_*(1 + \alpha_0^2)$. By this, the acceleration of a body depends on the gravitational binding energy, leading to a violation of the SEP [41]. As a result, one has a strong-field modification of the mass function of the binary system that enters the expression for the Shapiro shape parameter $s$ [cf. Eq. (5.8) in Ref. [42]].

At the first PN level, one has three body-dependent strong-field counterparts of Eddington’s weak-field PPN parameters $\gamma^{\text{PPN}}$ and $\beta^{\text{PPN}}$ [204] and a periodic modification of the MOI due to a violation of local position invariance. All this leads to genuine strong-field modifications of the corresponding PK parameters (see Appendix C for details).

Beyond the first PN level, there are modifications to the conservative dynamics at 2PN [183], which are not (yet) relevant for the test presented here. Hence, it is sufficient to use the corresponding GR terms. However, the presence of the dynamical scalar field $\varphi$ considerably modifies the loss of energy and angular momentum due to the emission of scalar GWs. As a result, in addition to the quadrupolar tensor GWs, one has monopolar, dipolar, and quadrupolar scalar GWs, entering the equations of motion at the 2.5PN, 1.5PN, and 2.5PN level, respectively. These dissipative contributions to the equations of motion result, most importantly, in a modification of $\dot{P}_b$ (see Appendix C for details). The leading term, i.e., scalar dipole radiation, is proportional to $(\alpha_A - \alpha_B)^2$ and, therefore, quite sensitive to the asymmetry in gravitational binding energy in a system. [205] Even though the difference in mass between A and B is only 7%, the small size of the effect is compensated for by the high precision $(1.3 \times 10^{-4})$ of the GW emission test (see Sec. VI B 2). Consequently, the GW damping test
<figure><img src="image_17.png" /></figure>

FIG. 14. Constraints on the DEF gravity from different experiments: Shapiro delay with the *Cassini* spacecraft [207], dipolar radiation (J1738+0333) [208], gravitational weak equivalence principle (J0337+1715) [76], and the double pulsar (this paper). Areas above a curve are excluded. (See Refs. [42,206] for details. We use $\chi^2 = 4$ as a conservative limit.) Solid lines are for a comparably stiff EOS (MPA1 in Ref. [167]). The dashed lines are based on a rather soft EOS (WFF1 in Ref. [167]), to illustrate the EOS dependence of the different limits. GR corresponds to $\alpha_0 = \beta_0 = 0$, and Jordan-Fierz-Brans-Dicke (JFBD) theory is along the vertical $\beta_0=0$ line with Brans-Dicke parameter $\omega_{BD} = (\alpha_0^{-2}-3)/2$. The black star indicates the parameters used in Fig. 15.

provides an important contribution to the constraints on DEF gravity from the double pulsar.

The PK parameters $r$ and $\delta_\theta$ are not used in the tests presented in this section. Because of their large uncertainty, they are not relevant for any of those tests. The same is the case for effects related to the scaling parameter $q_{\text{NLO}}$.

Finally, we note that the expression for the mass ratio remains unchanged, since $R \equiv m_A/m_B = x_B/x + O(c^{-4})$ [63].

Figure 14 shows the constraints on DEF gravity from the double pulsar, presented as a curve (boundary of allowed region) in the $\alpha_0 - \beta_0$ parameter plane. When calculating the limits, we followed the procedure outlined in Refs. [42,206]. As one can see from Fig. 14, the double pulsar provides the best limit on DEF gravity for $\beta_0 \lesssim -3$. In Fig. 15, for illustration purposes, we plot a mass-mass diagram for a point in the $\alpha_0 - \beta_0$ plane that was previous to this paper not excluded: $\alpha_0 = 5 \times 10^{-4}$ and $\beta_0 = -4$ (marked by * in Fig. 14). As one can see for this specific example, $T_1(5 \times 10^{-4}, -4)$ is falsified by the double pulsar, in particular, by the GW test.

## B. Bekenstein's TeVeS

As a second alternative to GR, we discuss Bekenstein's TeVeS [43]. In Ref. [208], TeVeS is extended to a TeVeS-like class of gravity theories with two coupling parameters $\alpha_0$ and $\beta_0$, similar to the class of standard scalar-tensor theories discussed in Sec. VII A. Here, we discuss only

<figure><img src="image_18.png" /></figure>

FIG. 15. Mass-mass diagram for DEF gravity with $\alpha_0 = 5 \times 10^{-4}$ and $\beta_0 = -4$. EOS MPA1 is used to calculate the curves. The curves fail to agree on a common region in the mass-mass plane (see, in particular, $\dot{\omega}$ and $\dot{P}_b$ curves), meaning that this specific scalar-tensor theory is excluded. The chosen point in the $\alpha_0 - \beta_0$ plane is not excluded by other experiments (see Fig. 14). This figure is merely to illustrate how a specific theory fails the double pulsar test. In this case, it is due to the additional energy loss from scalar GWs, predominantly the dipolar contribution. The $\dot{P}_b$ curve is based on Eq. (44). The 2PN and Lense-Thirring contributions to $\dot{\omega}$ are calculated assuming GR, since any corrections from the scalar field at this PN level are insignificant.

Bekenstein's original TeVeS, which corresponds to $\beta_0 = 0$, since this is the case where the double pulsar provides the most important contribution. Although Bekenstein's TeVeS has recently been falsified by the multimessenger observation of the LIGO/Virgo merger event GW170817 [44], we still use it here as an example of a MONDian gravity theory that passes Solar-System tests but can be excluded by a binary pulsar experiment. At the same time, this is a significant update on the result presented in Ref. [208].

The crucial difference between DEF gravity in Sec. VII A and TeVeS is the specific form of the physical metric, which ensures that $\gamma^{\text{PPN}} = \beta^{\text{PPN}} = 1$, making TeVeS indistinguishable from GR in the weak-field slow-motion regime. Therefore, by design, TeVeS passes Solar-System tests, like the one from the *Cassini* spacecraft [207]. Interestingly, the special form of the physical metric also leads to the fact that $\gamma_{AB} = 1$ in TeVeS-like gravity theories [208], different from Eq. (C2).

As discussed in Ref. [208], the predictions for the Einstein delay amplitude $\gamma_E$ and the GW damping $\dot{P}_b$ keep the same form as in standard scalar-tensor theories. However, $\alpha_A = \alpha_B = \alpha_0$, and, therefore, there is no dipolar radiation. As a consequence, dipolar radiation tests like in
<figure><img src="image_19.png" /><figcaption>FIG. 16. Mass-mass diagram for Bekenstein's TeVeS with a "natural" transition from the Newtonian to the MONDian regime ($\alpha_0 = 0.04$) [210]. Such a theory is obviously in contradiction with the double pulsar observations. The $\dot{P}_b$ curve is based on Eq. (44). The 2PN and Lense-Thirring contributions to $\dot{\omega}$ are calculated assuming GR, since any TeVeS-related corrections for them are negligible in this test.</figcaption></figure>

Ref. [208] do not provide any constraints on Bekenstein's TeVeS. The same is the case for the universality of free fall test with the pulsar in the stellar triple [76,209]. The double pulsar, on the other hand, is sensitive to the modification of the MOI due to the periodic variation in the local gravitational constant (modification of $\gamma_E$) related to the scalar field of TeVeS and has a sufficiently precise GW test that is sensitive to modifications of the quadrupolar GW damping due to the presence of the dynamical scalar field. This is obvious from Fig. 16, which gives the mass-mass diagram for Bekenstein's TeVeS with a coupling parameter $\alpha_0$ that guarantees a natural transition from the Newtonian to the MONDian regime (see discussion in Ref. [210]). The deviations in $\gamma_E$ and $\dot{P}_b$ clearly falsify that theory. The current double pulsar observations require $|\alpha_0| < 9 \times 10^{-3}$ (95% C.L.) which requires a highly unnatural behavior of the scalar-field contribution to the gravitational acceleration (cf. Fig. 8 in Ref. [208] and Fig. 3 in Ref. [210]).

TeVeS-like theories that modify the conformal factor between the Einstein and the Jordan frame, as introduced in Ref. [208], are similarly excluded. For $\beta_0$ values significantly different from zero, the constraints do primarily come from dipolar radiation and universality of free fall tests, e.g., from limits in Refs. [76,208].

As a final comment, we emphasize that like in Ref. [208] we neglect the contributions of the vector field to the dynamics and the gravitational radiation. We consider this as a conservative approach, since a dynamical vector field

is expected to increase the GW damping by extracting additional energy from the orbital dynamics. Furthermore, modifications of TeVeS that modify only the vector part and leave the scalar sector of the theory unchanged should, therefore, generally be excluded by this test as well (see discussion in Ref. [210] on the role of the scalar field to reproduce the MONDian dynamics).

## VIII. ASTRONOMICAL IMPLICATIONS

In this section, we complete our analysis by revisiting several important properties of the double pulsar system. We first discuss results relating to the system distance, comparing our best estimate with the DM-based estimates based on the two principal Galactic electron-density models. The observed discrepancies inform how the model parameters of ISM structures in the path to the system should be modified. We next use the observed scintillation properties and DM variations over a 13-yr data span (Sec. IV B) to investigate the structure of the ISM over 6 orders of magnitude in fluctuation scale. Finally, we consider the implications of the observed system proper motion and implied transverse velocity on formation models for the system and summarize the evidence for a prograde rotation sense for the A pulsar (Sec. VI B 5).

### A. The distance to the double pulsar

The astrometric results derived from timing (Sec. IV B) and from interferometric data (Sec. IV C) can be compared against each other and with earlier VLBI observations using the Australian Long Baseline Array [93]. A detailed comparison between the two VLBI measurements is presented in Appendix A, in which we demonstrate that our new VLBA results are more reliable than the previous Long Baseline Array results, in part due to a previously unappreciated source of systematic error (refractive image wander). Accordingly, we focus on the VLBA results presented in this work.

The estimated parallax from the VLBA measurements is $\pi_v = 1.30_{-0.11}^{+0.13}$ mas (Sec. IV C). This contrasts with the value estimated from the pulsar timing of $\pi_t = 2.15 \pm 0.48$ mas (Sec. IV B). The implied pulsar distances are, respectively, $770 \pm 70$ and $465_{-85}^{+134}$ pc. As described in Appendix A, we adopt a weighted mean of these two distances, $735 \pm 60$ pc, as our best estimate of the pulsar distance. Note that none of these distance estimates change beyond their uncertainties if one applies the Lutz-Kelker correction (cf. Refs. [211,212]). Refractive image wander also potentially has an effect on the timing-based astrometry. We show in Appendix B 2 that this effect is negligible, changing the values of the astrometric parameters by $\lesssim 0.2\sigma$.

We now consider the impact of the pulsar distance on other observables to see if they can assist in determining the most likely parallax value.
The estimated distance for the pulsar based on the NE2001 Galactic electron-density model [45] is 516 pc and, based on the YMW16 model [46], 1105 pc. We note, however, that the derivation of the YMW16 electron-density model includes the earlier VLBI-based distance [93] (corrected for the Lutz-Kelker bias) of 1100 pc in its set of independently measured distances used to calibrate the model. Removing the double pulsar from the list of independent distances and redetermining the parameters of the YMW16 model has a surprisingly large effect on the estimated distance, changing it from 1105 to 463 pc [213]. Further investigation shows that the reason for this large change is that, in the modified analysis, the radius of the model component representing the Gum Nebula changed from 125 to 128 pc. This has the effect of placing the LOS to the double pulsar just within the Gum Nebula, which is centered at about 450 pc from the Sun [214], rather than just outside it. The resulting additional electron column density is sufficient to place the pulsar within the Nebula at 463 pc. The original YMW16 distance for the double pulsar, 1105 pc, is very close to the adopted independent distance of 1100 pc. This, and the large change when this independent distance is removed from the model, highlights the fact that the number of independent distances used to build the original YMW16 model (viz., 189) is relatively small. In the intervening years, many more independent estimates of pulsar distances have been published, most notably the VLBI pulsar parallax study PSR$\pi$ [215]. Consequently, future generations of the YMW16 model will be less susceptible to this problem.

Of course, the real Gum Nebula does not have the relatively sharp and well-defined edge assumed by the YMW16 model. It is a complex region of ionized gas that is believed to be a greatly expanded remnant of a supernova explosion, a fossil H II region, a wind-blown bubble, or some combination of these (see Ref. [214] for an $H\alpha$ image of the nebula and a summary of its properties). The Vela supernova remnant is superimposed on the Nebula and probably lies within it but is not the source of the bulk of the ionization. The Gum Nebula is roughly circular in shape, with an angular radius of about 23°, corresponding to a physical radius of about 190 pc. It has large extension to the west (i.e., lower right ascension or Galactic longitude) which covers the projected direction of the double pulsar. The identification of this extended nebulosity as the location of the electron density fluctuations responsible for the DM variations and scintillation of the double pulsar is discussed in Sec. VIII B below.

## B. ISM properties

In addition to the *average* DM value used above, we can also make use of the observed *variations* of the DM shown in Fig. 3, in order to infer properties of the intervening ISM and to compare these with expectations. These variations are consistent with our LOS moving across electron-density

enhancements and deficits in the ISM (see, e.g., Ref. [216]). These density fluctuations are part of the same plasma turbulence that causes angular scattering and intensity scintillation. The diffractive intensity scintillations of pulsar A have been thoroughly analyzed [92], and our observations do not have sufficient frequency resolution to improve on this work. In the following analysis, diffractive parameters are taken from Ref. [92]. They show that the scattering screen is located at about 30% of the distance from Earth to the pulsar or about 220 pc from Earth for the adopted pulsar distance of 735 pc. If the extended nebulosity on the western side were somewhat closer to Earth than the Gum Nebula itself, we could identify the scattering screen with this nebulosity.

We also note that the rotation measure (RM) of the double pulsar is $+120.8 \pm 0.2 \text{ rad m}^{-2}$ [154], which is consistent both with the RM of other pulsars located near the Gum Nebula (in projection) and with the RM of extragalactic sources located behind the western nebulosity near the double pulsar [214]. Both of these are consistent with the double pulsar being located behind the western extension and much of the Faraday rotation occurring in this region.

Structure functions (SFs) are useful here, because the diffractive scattering can be calculated directly from the SF, and stochastic processes with a power-law spectrum have a power-law SF. The SF for DM variations is defined by

$$D_{\text{DM}}(\tau) = \langle [\text{DM}(t + \tau) - \text{DM}(t)]^2 \rangle, \quad (57)$$

where $\tau$ is the lag between DM samples (see, e.g., Ref. [217]).

To obtain the most reliable estimates of the DM variations and the corresponding SF, we restrict our analysis to the interval MJD 53500–58300, which has better sampling in both time and frequency compared to earlier and later times (see Fig. 2). We use the results of the Monte Carlo analysis described in Sec. IV B to estimate the DM variations at 24-d intervals. The resulting DM-variation curve is fully consistent with that shown in Fig. 3 for the corresponding interval.

The structure function estimated directly from Eq. (57) is shown at the upper right in Fig. 17. At lags between 100 and 1000 d, the observed SF is approximately power law with an exponent much smaller than the Kolmogorov value of 5/3. [218] A linear fit to the region between 96 and 1008 d is shown as a red line on the plot. The SF is less useful for lags greater than about 2500 d (i.e., half the data span), because there are few independent samples to be averaged.

The LOS velocity (excluding the binary motion) is about 37 km/s [92], so the timescales that are easily measurable in DM(t), 100–1000 d, correspond to spatial scales of approximately 2.0 to 20 A.U. [219] The intensity scintillations, unlike DM(t), have a characteristic spatial scale $s_0$
<figure><img src="image_20.png" /><figcaption>FIG. 17. Structure function of DM from MJD 53500 to 58300 with 24-d sampling. The point marked with a + is derived from diffractive timescale measurements at 820 MHz [92]. The blue solid line shows the extrapolation from the diffractive scale up to the measured lags assuming a Kolmogorov electron-density fluctuation spectrum. The dashed lines indicate ±0.1 in the spectral exponent. The red line shows the result of a weighted fit of a power law to the measured structure function for lags between 96 and 1008 d, as indicated by the red dots.</figcaption></figure>

which is very much smaller at approximately 3700 km. The spatial autocorrelation function (ACF) of the scintillations is given by $ACF(s) = \exp[-D_{\phi}(s)]$, where $D_{\phi}(s)$ is the structure function of the phase shift $\phi$ caused by the scattering plasma. The phase delay is the negative of the group delay, so we can obtain the phase directly from DM by computing the group delay and multiplying by $2\pi\nu$. A measurement of $s_0$ at $\nu$ provides $D_{\phi}(s_0) = 1$, which defines $D_{\text{DM}}(s_0)$. It should be noted that the timescale of the scintillations is extremely variable, because the binary velocity is much higher than the mean LOS velocity. The analysis of Ref. [92] fits the binary velocity to extract $s_0$ which is not time variable.

In this way, the observed diffractive intensity scintillations can be used to define the point in Fig. 17 marked with a + near lag 110 s. A blue line with the Kolmogorov slope of 5/3 is drawn from the diffractive point to meet the red line at a lag of approximately 30 d. The dashed lines in Fig. 17 illustrate the effect of an uncertainty in the slope of the extrapolation from the diffractive point, indicating a likely range for the intersection lag of between 10 and 100 d, or a spatial scale of 0.2 to 2.0 A.U. using the mean ISM velocity given by Ref. [92].

The location of the spectral break is important, because it can be used to estimate the thickness of the scattering region which is causing both the DM($t$) fluctuations and the intensity scintillations. Accordingly, we simulate 100 realizations of a pure power-law process with an SF exponent of +0.5 and apply the convolutions by a 100-d triangle implicit in the Tempo2 analysis and the 24-d binning

<figure><img src="image_21.png" /><figcaption>FIG. 18. Power spectral density for the DM variations used to form the SF in Fig. 17 (segmented black line). The blue and red solid lines are means of 100 simulations of the 1D PSD of DM variations given by Tempo2 for fluctuation spectral exponents of -8/3 (Kolmogorov) in blue and -1.5 in red. The dashed red lines give the 90% range of the -1.5 simulations. The vertical black line denotes the annual frequency.</figcaption></figure>

of the Monte Carlo results. Details of the simulation methods are described in Appendix B 1.

We first adjust the amplitude and spectral exponent of the simulated DM($t$) so its power spectral density (PSD) matches that of the observations. [220] The match is much clearer in the PSD than in the SF, because the PSD estimation errors are statistically independent and those of the SF are highly correlated. Figure 18 shows the PSD of the observed DM variations as a segmented black line. The mean PSD of the simulations is overplotted with the 90% confidence limits, and both the level and the spectral exponent are adjusted so the observed PSD fits within these confidence limits. A Kolmogorov PSD simulation is also shown to make it clear that the observations cannot be fitted with a Kolmogorov spectrum. The flatter spectrum with exponent -1.5 agrees well with the observations up to frequencies of approximately $7 \times 10^{-3}$ d⁻¹. At higher frequencies, the observed PSD is dominated by white noise, which is not included in the simulations.

The observed PSD also shows a statistically significant peak at a frequency close to the annual frequency. This peak arises from the annual motion of Earth moving the LOS to the pulsar across transverse A.U.-scale gradients in the integrated screen density. For further details on both the effect of the Tempo2 convolutions and the annual feature in the PSD, see Appendix B 1.

Figure 19 shows the measured SF from Fig. 17 in black on an expanded scale. The structure functions of the best-fit simulations with a spectral exponent of -1.5 are calculated, and the mean and 90% confidence limits are shown in red.
<figure><img src="image_22.png" /><figcaption>FIG. 19. Comparison of the observed structure function for DM variations with white noise removed (in black) with simulations for an input SF exponent of +0.5. The solid red line is the mean of 100 simulations including the effect of the time-domain convolutions, and the red dashed lines give the 90% range. The cyan line is the mean of the simulations before the time-domain convolutions. The blue line is the Kolmogorov extrapolation (SF exponent 5/3) from the diffractive scale.</figcaption></figure>

The theoretical SF simulated before convolutions required in the data analysis is indicated by the cyan line, and the blue line is the extrapolation from the diffractive timescale, also shown in blue in Fig. 17.

The steepening of the SF toward lags shortward of a few hundred days results from the time-domain convolutions. However, there must be a true break in the SF slope from a value close to the Kolmogorov 5/3 at short lags to a flatter slope at lags greater than a few hundred days. The intersection of the blue and cyan lines at a lag of approximately 170 d, which corresponds to a linear scale of about 3 A.U., probably gives the best estimate of the frequency of spectral break in the ISM fluctuations.

A break in the PSD from a Kolmogorov spectral exponent of -8/3 to an exponent of -5/3 occurs naturally in turbulence in a thin layer, and the transition frequency is directly related to the thickness of the layer. As discussed by Lay [221], the inverse of the transition frequency for the corresponding PSD corresponds to twice the screen thickness. To relate the SF transition lag and PSD transition frequency, we simulate our observed DM(t) 1000 times with an SF transition at 170 d and compute the average PSD, which we find has a transition scale of 1/800 d. Accordingly, the screen depth must correspond to a lag of about 400 d or a spatial scale of about 8 A.U. The mean electron density integrated through the scattering layer must exceed its rms fluctuation, since it is positive definite and the distribution is roughly symmetric. The rms of DM(t) is about 0.001 pc cm⁻³ (Fig. 3), and, therefore, a

lower limit on the mean density in the scattering region is about 25 cm⁻³.

As mentioned in Sec. VIII A above, the scattering screen can be identified with the nebulosity to the west of the Gum Nebula, which is probably associated with the Nebula but somewhat closer to Earth. From the observed Hα intensity, Purcell et al. [214] infer emission measures in the fainter parts of the Nebula that correspond to a mean electron density of about 1.5 cm⁻³. This estimate is appropriate for the nebulosity in front of J0737-3039A/B. Recognizing that a filament with a transverse scale of even 10 000 A.U. would be invisible on the Hα image (which has a resolution of 6 arcmin), an overdensity by a factor of 20 or even more in the scattering screen of depth 8 A.U. is quite plausible. We, therefore, conclude that the western nebulosity associated with the Gum Nebula is a likely location for the scattering screen. This is consistent with J0737-3039A/B being located at our preferred distance of 735 pc.

The clear flattening of the spectral exponent at a scale of 3 A.U. is very unusual in ISM observations, which usually follow the Kolmogorov exponent. Of the pulsars observed by pulsar timing arrays, which provide DM(t) with the required precision, only PSR J1713+0747 shows similar behavior. It should be noted that the diffractive scintillation of the double pulsar is entirely normal for pulsars with similar DM. Only the flattening of the exponent due to the unusually thin scattering region is abnormal.

### C. System velocity and geometry

Using the distance adopted from the weighted mean probability distribution in Fig. 23, i.e., 735 ± 60 pc, and taking the parameters of Ref. [130], in particular, the Solar distance to the Galactic midplane, $z_{SSB} = +14$ pc, [222] we find that the double pulsar is about 43 pc below the Galactic plane.

The proper motion of the double pulsar is extremely well constrained. Its total value 3.304 ± 0.033 mas yr⁻¹ converts into a transverse velocity with respect to the SSB of just 11.5 ± 1.0 km s⁻¹, where the uncertainty is dominated by that of the distance measurement. The components in Galactic longitude and latitude are $v_l = -10.7 \pm 0.9$ km s⁻¹ and $v_b = -4.3 \pm 0.4$ km s⁻¹, respectively.

The radial velocity $v_R$ with respect to the SSB is unknown, but, for reasonable values for $v_R$, the double pulsar has a peculiar velocity of ≲50 km s⁻¹ with respect to its standard of rest (see Fig. 20). In particular, the vertical component of the full Galactic velocity $v_z$ is small (≪5 km s⁻¹), meaning that the double pulsar moves practically along the Galactic plane, information that can be used to constrain possible birthplaces of the double pulsar system.

The small peculiar velocity of the double pulsar is consistent with the suggested low-kick supernova explosion for the birth of B (e.g., Refs. [48,49]), although this interpretation has been challenged [47]. A low-kick
<figure><img src="image_23.png" /><figcaption>FIG. 20. Velocity of the double pulsar with respect to its standard of rest ($\Delta v_{LSR}$; blue) and its velocity perpendicular to the Galactic plane ($v_z$; red), as a function of the unknown radial $v_R$. The gray area shows the 90% probability range for the radial velocity, assuming an (a priori) uniform probability distribution for the direction of the 3D velocity vector with respect to the SSB (cf. Ref. [48]). Calculations are based on the parameters and Galactic potential in Ref. [130], including the components for the peculiar velocity of the Sun ($U_\odot, V_\odot, W_\odot$) = (8.6, 13.9, 7.1) km s$^{-1}$.</figcaption></figure>

scenario is, however, also supported by the small misalignment angle $\delta_A$ between A's spin and the total angular momentum vector, which is based on the lack of profile changes in A [50,51]. Visible secular changes in the pulse profile of A are still lacking, even with the doubling of the observation data span, leading to even tighter constraints that will be presented elsewhere.

The previous argument that the prograde solution $\delta_A < 3.2^\circ$ is more likely than the retrograde solution ($176.5^\circ < \delta_A < 180^\circ$) is based on the justifiable notion that a large and fortuitous kick would be required for the observed alignment [51]. Here, we use our measurement of the relativistic effect of aberrational light bending to show unequivocally that A is rotating prograde (see Sec. VI B 5 for details), supporting the earlier findings of Pol et al. [52] using a fully independent method based on the emission properties of pulsar B. We note that B is also rotating prograde [57] but with a significant misalignment angle of $\delta_B \approx 50^\circ$ [54,57].

We conclude that our results are fully consistent with theoretical considerations for the formation of the double pulsar system, as well as with the observed stability of the mean pulse profile.

## IX. PROSPECTS

The results presented in the previous sections verify our past expectations on the development of the precision in the timing parameters (cf. Ref. [33]). This indicates that not only has the instrumentation improved as expected, but also that possible, as yet undetected, effects have not yet impacted on our timing observations or their interpretation.

In particular, this is true for relativistic spin precession, which could lead to a changing pulse profile, hence making the timing more difficult and potentially less accurate (see, e.g., Refs. [224,225]). In fact, spin precession rendered B temporarily undetectable [58,226], so that the precision of the (direct) mass ratio measurement [22,31], which requires timing of both A and B, cannot easily be improved. We refer to the recent work by Noutsos et al. [57] for a detailed study of possible future improvements in this additional constraint.

The results presented here are not affected by the disappearance of B. The contrasting lack of any profile changes in A further strengthens the notion that the spin direction of A is closely aligned with the total angular momentum vector of the system (see Sec. VIII C). In this case, A will continue to be visible and available for precision timing for the foreseeable future, enabling improvements on all parameters and effects studied here. This includes the new effects presented in this work, which will become even more relevant in the future. For instance, in the future, the system's acceleration relative to the SSB will become a crucial parameter, while we now also have to take into account the relativistic mass loss due to the pulsar spin down and the NS structure. This will allow us eventually to convert our current constraints on the MOI and the NS radius into concrete measurements. Corresponding detailed computations and simulations are presented by Hu et al. [110], who show that an MOI measurement with 11% accuracy by 2030 is likely. Combination with results from other sources, such as GW emission from NS-NS mergers, to determine the EOS sufficiently well, would allow for a 7% test of Lense-Thirring precession or, alternatively, provide a 3$\sigma$ measurement of the NLO GW damping in GR.

Obviously, future studies of the double pulsar will benefit greatly from more sensitive telescopes. For example, with observations such as those with MeerKAT [110,154], we can study possible profile changes of A due to signal deflection at B near superior conjunction. Moreover, the much-improved sensitivity should also result in a better timing performance during the eclipse, in general, allowing us to probe this important part of the Shapiro delay curve in more detail for improved NLO tests. Detailed eclipse studies will allow us to track the rotations of B even without seeing its radio signal, also helping to improve the mass-ratio measurement [54,154]. Finally, continued studies (timing, eclipse, and profile studies) of A will improve the precision of all relativistic effects

Future studies will also include a Bayesian-based analysis of the timing data. Currently, the number of TOAs to be studied typically overwhelms the available algorithms and computing power. Optimization work, including studies to deploy smoothly varying, frequency-dependent templates [160,161], is underway, with initial results on a subset of data being fully consistent with the results shown here.
In addition to timing observations, information obtained via different routes will further aid the exploitation of the system. First, further VLBA astrometric imaging observations can reduce the current uncertainty of the VLBI parallax and, hence, of the distance. A factor-of-2 reduction in the VLBI parallax uncertainty is within reach of an extended campaign, which would lead to a distance precision of 5%. Second, as mentioned in Sec. VI A, scintillation measurements can result in an independent measurement of the orbital inclination angle. This is one of the goals of ongoing MeerKAT observations [154].

## X. SUMMARY AND CONCLUSIONS

We have presented results from a timing campaign for the double pulsar with a duration of more than 16 years, using data from six different telescopes. In order to analyze the data adequately, we have introduced a number of new methods. To facilitate analyses of the system astrometric parameters and DM variations, we formed a dataset with 4-min integration time per TOA, which provided the widest possible coverage in terms of time and observing frequency. Analysis of the binary parameters, including the important PK parameters, requires a second dataset with TOAs based on 30-s subintegrations in order to optimally resolve the fast compact orbit of PSR J0737-3039A. The analysis of this latter dataset required also the development of a new timing model and an improved implementation in Tempo. This takes into account higher-order effects, including velocity-dependent terms in the Shapiro delay. With these efforts, our results not only improve on precision of the previously measured parameters [31], but also reveal newly measured effects.

Using a Monte Carlo analysis, we obtained a probability distribution for the weighted mean of the VLBI and pulse timing annual parallaxes to obtain our best estimate of $1.36_{-0.10}^{+0.12}$ mas for the parallax and $735 \pm 60$ pc for the distance of the double pulsar. We emphasize again that, because of the fortunate arrangement of the secular (Shklovskii) and the Galactic acceleration being of opposite sign for the double pulsar, the impact of distance uncertainties is minimal for our current GW emission test [177]. This will eventually change, but continued VLBI and timing observations should converge onto a higher-precision distance well before any limitations are reached for the GW emission test.

We then discussed in detail the contributions to the timing model that need to be considered given the much-improved precision of our measurements. Here, we considered, in particular, higher-order timing effects, extrinsic modifications to our observed timing parameters, and also spin contributions. The latter requires us, for the first time, to pay attention to the EOS of NSs when interpreting pulsar timing data.

In this work, we obtained more independent tests of GR than are possible in any other system. We summarize the

<table><caption>TABLE V. Summary of the relativistic effects measured in our analysis and list of the resulting independent strong-field tests of GR. For each test, the remaining PK parameters and the mass ratio are used to determine the masses of pulsars A and B as input for GR predictions. In addition, parameters that test the significance of specific higher-order contributions in the advance of periastron and the signal propagation are given.</caption><thead><tr><th>Relativistic effect</th><th>Parameter</th><th>Obs./GR pred.</th></tr></thead><tbody><tr><td>Shapiro delay shape</td><td>s</td><td>1.00009(18)</td></tr><tr><td>Shapiro delay range</td><td>r</td><td>1.0016(34)</td></tr><tr><td>Time dilation</td><td>γ<sub>E</sub></td><td>1.00012(25)</td></tr><tr><td>Periastron advance</td><td>ω̇ ≡ n<sub>b</sub>k</td><td>1.000015(26)</td></tr><tr><td>GW emission</td><td>Ṗ<sub>b</sub></td><td>0.999963(63)</td></tr><tr><td>Orbital deformation</td><td>δ<sub>θ</sub></td><td>1.3(13)</td></tr><tr><td>Spin precession</td><td>Ω<sup>spin</sup><sub>B</sub></td><td>0.94(13)<sup>a</sup></td></tr><tr><td colspan="3"><em>Tests of higher-order contributions</em></td></tr><tr><td>Lense-Thirring contrib. to k</td><td>λ<sub>LT</sub></td><td>0.7(9)</td></tr><tr><td>NLO signal propagation</td><td>q<sub>NLO</sub>[total]</td><td>1.15(13)</td></tr><tr><td>... from signal deflection</td><td>q<sub>NLO</sub>[deflect]</td><td>1.26(24)</td></tr><tr><td>... from signal retardation</td><td>q<sub>NLO</sub>[retard]</td><td>1.32(24)</td></tr></tbody></table>

six PK parameters measured in this work in Table V, additionally including the test of relativistic spin precession. The measurements allow us to test conservative aspects of the orbital dynamics of two strongly self-gravitating masses up to 2PN order, including an approximately $1\sigma$ constraint on the Lense-Thirring contribution, which, in turn, could be used to constrain the MOI of pulsar A under the assumption of GR. A signal propagation test resulted in a confirmation of GR at the $4 \times 10^{-4}$ level for the propagation of photons in the gravitational field of a strongly self-gravitating (material) body (see Sec. VI B 5). Moreover, NLO contributions are clearly present in the timing data, confirmed with a precision of about 10%. The most precise GR test available with our data probes the radiative aspects of GR, yielding a test at 2.5PN level in the equations of motion with a precision of $1.3 \times 10^{-4}$ (95% C.L.; see Sec. VI B 2). In terms of overall fractional precision, this is the most precise test of GR's predictions for GW emission currently available.

The new effects that we detected include a relativistic deformation of the orbit and higher-order contributions to the Shapiro and aberration delay. The latter allow us to infer the spin direction of A as prograde, confirming earlier results that require a low-kick supernova as the formation process for pulsar B. We can determine the masses of pulsars A and B with a precision of $10^{-5}$ (modulo an unknown Doppler factor, which is expected to deviate by less than $10^{-4}$ from unity), which we can combine with the first constraints on the NS MOI ever obtained from pulsar timing. With the double pulsar, we are able to greatly improve the measurement of orbital period decay resulting from GW damping compared to the current best

<footnote><sup>a</sup>Determined in Ref. [54].</footnote>