One parsec (pc) is approximately $3.086\times {10}^{16}\mathrm{m}$, One $\mathrm{kpc}=1000\mathrm{pc}$.

Note that the angle $\theta$ in Ref. [54] is (to a good approximation) $180°-{\delta }_{\mathrm{B}}$, since in the coordinate system of Ref. [54] the orbital angular momentum is closely aligned with the negative $z$ direction.

At least up to 1PN order for boost-invariant gravity theories [63].

We note that we have also explored alternative methods, e.g., using the “DMX” method implemented in tempo, in order to determine the astrometry and its uncertainties. Using simulated timing data with parallaxes of known sizes injected into our datasets, we determined that the method implemented here appears to be the most reliable in recovering the correct astrometry with appropriate uncertainties.

Depending on the literature, this is called the argument of periastron. We follow the naming that is more common in pulsar astronomy and, for instance, also used in Ref. [97] (see their Fig. 3.2).

The rotational spin frequency parameters are estimated with a standard tempo analysis. As discussed in Ref. [89], this procedure underestimates the true uncertainty in the presence of unmodeled red timing noise. Comparison of the standard analysis with the results of a full generalized least-squares analysis shows that the true uncertainty in the spin frequency parameters is about an order of magnitude larger than the values quoted in Table 4. We emphasize that the unmodeled red noise has no effect on the orbital parameters as expected because of their much shorter timescale.

Mass determination from PK parameters yields the product $G{m}_i$ ($i=\mathrm{A},\mathrm{B}$). Consequently, whenever we express the masses in terms of solar masses ${\mathrm{M}}_{\odot }$, the numerical value is calculated via $G{m}_i/(\mathcal{G}\mathcal{M}{)}_{\odot }^{\mathrm{N}}$, where $(\mathcal{G}\mathcal{M}{)}_{\odot }^{\mathrm{N}}\equiv 1.3271244\times {10}^{26}{\mathrm{cm}}^3{\mathrm{s}}^{-2}$ is the nominal solar mass parameter as defined by the IAU 2015 Resolution B3 [105].

For the precision needed, we can safely use $g_{S_{\mathrm{A}}}^{\parallel }$ instead of $g_{S_{\mathrm{A}}}$ in Ref. [36], since ${\delta }_{\mathrm{A}}<3.2°$, while $g_{S_{\mathrm{A}}}\simeq g_{S_{\mathrm{A}}}^{\parallel }\cos {\delta }_{\mathrm{A}}$ for the double pulsar (cf. Ref. [33]).

Although the mass of A ($1.34{\mathrm{M}}_{\odot }$) is slightly below $1.4{\mathrm{M}}_{\odot }$, within the accuracy needed here, the radius constraints in Ref. [112] for a $1.4-{\mathrm{M}}_{\odot }$ NS can also be applied to pulsar A.

There are other distributions that have been derived for the MOI of pulsar A; see, e.g., Refs. [116, 117, 118]. However, for the results of this paper, this does not make any difference.

Simulations are based on the integration of null geodesics in a Schwarzschild spacetime, comoving with pulsar B. Since additive constant terms can be ignored, the use of Eddington-Finkelstein coordinates turns out to be of great advantage in the numerical integration for the propagation time [141].

The (intermittent) eclipses of A by the rotating magnetosphere of B lead to a significant reduction in the number of pulses available for the integrated pulse profile. This, in turn, reduces its S/N, leading to significantly fewer viable TOAs, which at the same time have considerably larger errors.

Lensing corrections to the retardation effect as discussed in Ref. [147] are totally negligible. The corresponding postfit residuals are below $0.02\mu \mathrm{s}$, restricted to a narrow range around conjunction where, because of the eclipse, the timing precision is anyway greatly compromised.

Note that the expressions in Refs. [147, 151] are based on the lensing equation, which makes the assumption of small angles (cf., e.g., Ref. [153]), and, therefore, they are strictly speaking valid only near conjunction, which also limits their usability in a tempo implementation.

Periodic variations in the pulsar profile due to a latitudinal shift in the LOS caused by aberration are observed in PSR $\mathrm{B}1534+12$ [156].

The DDS timing model is a modification of the DD model [29], where the Shapiro shape parameter $s$ is replaced by the logarithmic Shapiro shape parameter $z_s$ via $s\equiv 1-\exp (-z_s)$ (see Ref. [149] for details). It is part of the tempo standard distribution.

Even though the factor ($1-\kappa$) is small, such correction is required before using the listed timing results for predictions of the folding parameters with tempo2, which implements TCB. According to IAU Resolution 2006 B3, the conversion takes place as $\mathrm{TDB}=\mathrm{TCB}-(1-\kappa )\times ({\mathrm{JD}}_{\mathrm{TCB}}-2443144.5003725)\times 86400-6.55\times {10}^{-5})$, where ${\mathrm{JD}}_{\mathrm{TCB}}$ is the relevant Julian date [164]. Note that this described conversion needs to be applied as given to the quoted epoch of periastron passage ($T_0$) value before it is used with tempo2.

This Doppler-shift rescaling, $m^{\text{measured}}=m^{\text{intrinsic}}/D$, between the pulsar-timing measured masses and the intrinsic ones [29] is due to the same effect which causes the masses directly deduced from gravitational-wave-template fitting in LIGO/Virgo data analysis to differ from the intrinsic, source-frame masses by a redshift factor: $m^{\mathrm{obs}}=(1+z)m^{\text{source}}$ [171].

We take the combined results of these references, obtained from the IMRPhenomPv2 waveform family.

Short-range modifications of GR that significantly affect the MOI could, in principle, show up in the Lense-Thirring test of Eq. (50); see discussion in Ref. [110].

This is, for instance, the case in the mEIH formalism, where (to leading order) the deflection angle is proportional to $(1+{\gamma }_{\mathrm{B}0})G_{\mathrm{B}0}$ (see e.g., Appendix C in Ref. [192]), which is already tightly constrained by the test of the Shapiro range parameter $r$. More generally, the NLO contribution to ${\mathrm{\Delta }}_{\mathrm{A}}$ is parametrized by $r$ and $s$ according to Eqs. (30) and (31).

It follows from the analysis in Ref. [51] that the rotational axis of A is (practically) perpendicular to the orbital plane. See also Sec. 8c.

An additional test is provided by the PK parameter ${\mathrm{\Omega }}_{\mathrm{B}}^{\mathrm{spin}}$; see Table 5.

There are special cases, like Barker’s constant-$G$ theory [200], where ${\dot{P}}_{\mathrm{b}}^{\dot{G}}$ cannot be written in this form.

To leading order, one finds in DEF gravity $B_D\simeq (5/96)({\alpha }_{\mathrm{A}}-{\alpha }_{\mathrm{B}}{)}^2[1+\mathcal{O}({\alpha }_{\mathrm{A}},{\alpha }_{\mathrm{B}})][1+\mathcal{O}(e^2)]$ (cf. Refs. [203, 206]).

Since the DM is the column integral of the electron density along the line-of-sight ($z$) axis, its 2D spatial spectrum is the same as that of the 3D density with the wave number in $z$ set to zero. Therefore, the 2D spectral exponent of DM ($\beta$) is the same as that of the 3D electron-density fluctuations, i.e., $-11/3$ for a Kolmogorov process. The corresponding 2D structure function has a spectral exponent $\alpha =-(\beta +2)=+5/3$ for Kolmogorov.

One astronomical unit (A.U.) is approximately $1.496\times {10}^8\mathrm{km}$.

The fluctuations in $\mathrm{DM}(t)$ are those of a 1D cut across the 2D spatial distribution. Therefore, the 1D PSD of $\mathrm{DM}(t)$ is the 2D spatial spectrum of DM (scaled by the velocity) with the perpendicular wave number integrated out. This reduction to 1D reduces the Kolmogorov spectral exponent from $-11/3$ to $-8/3$.

This value agrees well with a pulsar-based estimate [223]. There is some spread in measured values for $z_{\mathrm{SSB}}$ (see, e.g., Table 1 in Ref. [223]), but this has little impact on our conclusions.