© H. Hu et al. 2022

1. Introduction

The double pulsar PSR J0737−3039A/B is a rich laboratory for strong-field gravity experiments. The system consists of a 23-ms recycled pulsar (‘A’) and a 2.8-s ‘normal’ pulsar (‘B’) in a nearly edge-on and slightly eccentric 2.45-h orbit (Burgay et al. 2003; Lyne et al. 2004). Various relativistic effects have been precisely measured in previous works (Kramer et al. 2006a, 2021a), including periastron precession, time dilation (gravitational redshift and second-order Doppler effect), Shapiro delay due to light propagation in the curved spacetime of the companion, and the orbital period decay, which currently provides the most precise test of quadrupolar gravitational waves predicted by general relativity (GR). In addition, the relativistic spin precession of B was measured by Breton et al. (2008) and the relativistic deformation of the orbit was newly detected in this system (Kramer et al. 2021a). All these make it a still unique system for gravity experiments.

Comparing with other gravity experiments, the Shapiro delay measured in the double pulsar probes the strongest spacetime curvature (∼10⁻²¹ cm⁻²) in a precision experiment with photons, that is, the interaction between gravitational and electromagnetic fields (Wex & Kramer 2020). In addition, with 16 yr of data, Kramer et al. (2021a) were able, for the first time, to measure higher-order effects of signal propagation in the strong gravitational field of a neutron star, which are currently not accessible via any other method. These include the retardation effect due to the movement of the companion (B) and aberrational light deflection by the gravitation of the companion. The latter confirms the prograde rotation of A, which is consistent with the results measured by Pol et al. (2018) using the emission properties of B and is in line with the expectations from binary evolution models.

In this work, we present observations of PSR J0737−3039A with the MeerKAT telescope, a precursor for the Square Kilometre Array (SKA) at mid-frequency range. Thanks to its location in the Southern Hemisphere, it permits a timing precision more than two times better than that of the Green Bank Telescope for this pulsar (Bailes et al. 2020; Kramer et al. 2021b). This superior precision enables an independent and improved measurement of signal propagation effects within a very short time span. We also investigate effects that have been expected but not been studied in detail before. These include potential profile variations due to latitudinal deflection, the detectability of latitudinal deflection delay, and the prospects of measuring the effect of lensing on the propagation time separately.

This paper is organised as follows. Section 2 describes the MeerKAT observations on the double pulsar and data processing. In Sect. 3, we introduce the concepts of gravitational signal propagation effects, including higher-order contributions to the Shapiro and aberration delay. The timing results and mass measurements are presented in Sect. 4. We then provide an in-depth study on the higher-order signal propagation effects in Sect. 5, with a focus on latitudinal deflection and lensing. In addition, we provide an improved analytical description for the signal propagation in the double pulsar. Finally, we discuss the results and future prospects in Sect. 6.

2. Observations and data processing

2.1. MeerKAT observations

The observations presented in this paper come from the MeerKAT telescope as part of the MeerTIME project (Bailes et al. 2020), which performs timing of known pulsars with various scientific themes. Observations on PSR J0737−3039A are conducted under the Relativistic Binary theme (RelBin, Kramer et al. 2021b), which focuses on testing the relativistic effects in binary pulsars to achieve measurements of neutron star masses and tests of theories of gravity. MeerTime observations are generally recorded using the Pulsar Timing User Supplied Equipment (PTUSE) signal processor. This processor receives channelised tied-array beamformed voltages from the correlator-beamformer engine of the MeerKAT observing system and is capable of producing coherently de-dispersed full-Stokes data in both the filterbank (search) mode and the fold (timing) mode, where the data are folded at the topocentric period of the pulsar. Details on the pulsar observing setup with MeerTime are explained by Bailes et al. (2020).

PSR J0737−3039A is regularly observed with a typical cadence of one month and duration of 3 h. As the orbital period of this pulsar is ∼2.45 h, the observations are scheduled to start shortly before an eclipse and finish after the second eclipse, in order to observe the eclipses twice in one observing session. The session is typically composed of a 30 min observation with the fold mode and search mode in parallel, followed by a 2 h fold-mode observation and another 30 min fold-search dual-mode observation. This specific arrangement is designed to maximise our sensitivity in detecting signal propagation effects, as well as in studying the magnetosphere of pulsar B (Lower et al., in prep.). Observations are performed with two receivers: the L-band receiver that covers the frequency range 856–1712 MHz, and the UHF receiver that covers the frequency range 544–1088 MHz, both with 1024 channels. The data presented here started in March 2019 and ran until May 2022. For the analysis in this paper, we use 29 full-orbit timing observations and 62 search-mode eclipse data sets, which amounts to a total of ∼87 h.

2.2. Timing data reduction

The raw timing data from the PTUSE machines are folded every eight seconds, which are then processed with the MEERPIPE data reduction pipeline. MEERPIPE carries out radio frequency interference (RFI) removal using a modified COASTGUARD algorithm (Lazarus et al. 2016), followed by flux and polarisation calibration. Details on polarisation and flux calibration are described in Serylak et al. (2021) and Spiewak et al. (2022), respectively.

After processing with MEERPIPE, the calibrated data products are reduced using the pulsar software package PSRCHIVE¹ (Hotan et al. 2004). We first correct for the rotation measure (RM) with the value measured in Kramer et al. (2021b), that is, RM = 120.84 rad m⁻². As the L-band observations between March 2019 and February 2020 were restricted to 928 frequency channels (dropping 48 channels each from the top and bottom bands), to maintain consistency throughout the analysis, we reduce the later L-band data to the same frequency channels. We treat the UHF-band data in the same way, as the roll-off adversely affects sensitivity of the top and bottom bands.

For this system, a complete timing model is only available in the pulsar timing software TEMPO² (Nice et al. 2015, more details are given in Sect. 4). Therefore, to fold the data more accurately, all data are supplied with a polyco-format ephemeris with the values measured in Kramer et al. (2021a). Since the double pulsar rapidly changes its orbital phase, the time span (TSPAN) of a predicted pulse phase solution has to be as small as possible to retain good precision³. With the PSRCHIVE version 2022-01-14, we set TSPAN to the minimum possible value, which is 3 min.

There is a known data processing issue with this pulsar, which is that the pulsar moves rapidly to a different orbital phase during the dispersion delay time. Thus, the pulses received at the same time at different frequencies correspond to different orbital phases and cannot be folded with the same phase prediction. If not properly accounted for, this folding issue will cause frequency-dependent orbital smearing. Standard pulsar software such as PSRCHIVE do not take this effect fully into consideration, even with the frequency-resolved TEMPO2 predictor⁴. To avoid this issue, we first de-disperse the total intensity data so that all frequencies correspond to the same orbital phase, then average the data, first in frequency and then in time⁵. Because of the profile frequency evolution and scintillation effects, data are sub-banded in frequency, with 32 sub-bands for the UHF-band and 16 sub-bands for the L-band.

As for the time-averaging, the integration time needs to be short enough in order to properly resolve the Shapiro delay and the next-to-leading-order (NLO) signal propagation contributions (q_NLO, see Sect. 3), which are largest at superior conjunction. We perform a simulation to test the measurability of these NLO contributions with different integration times. The results show that a good measurement of Shapiro delay and q_NLO can be achieved if the integration time is ⪅32 s, but it becomes significantly worse if the integration time is longer than 1 min for q_NLO and 2 min for Shapiro delay. Therefore, we average all data with a 32-s integration time, consistent with the analysis by Kramer et al. (2021a). After frequency and time averaging, data are re-dispersed to allow measurements of dispersion measure (DM) in the timing analysis.

2.3. Wide-band templates and TOA extraction

Wide-band observations such as MeerTIME can suffer significant profile evolution in frequency, and thus the traditional 1D template is not favoured. To best determine the pulse time of arrivals (TOAs) at multiple frequencies, we employ frequency-dependent 2D templates. With this technique, DM measurements are dependent on the DM value used to align (de-disperse) the 2D template. Due to the correlation between the DM and profile evolution, DM measurements are to some extent frequency dependent, which can lead to a DM offset between L-band and UHF-band data. This could potentially be solved with a simultaneous observation with L-band and UHF-band receivers, which is missing in our case. Therefore, to avoid this problem, we choose a bright observation from each band for making 2D templates, and measure DM using data from their overlapping frequencies. Then, we use these DM values to de-disperse the corresponding full-bandwidth data. This minimises the DM offset between L-band and UHF-band data, which can be seen in Fig. 2. These data are then sub-banded and averaged in time. Finally, by smoothing the profiles with psrsmooth/PSRCHIVE, we obtain 2D templates, with 16 sub-bands at L-band and 32 sub-bands at UHF-band. These templates are then used to measure frequency-resolved TOAs by cross-correlating with the reduced data using pat/PSRCHIVE. The pulse profile of PSR J0737−3039A at multiple frequencies is shown in Fig. 1. More information on the observing systems and data sets is given in Table 1.

Fig. 1. Pulse profile of PSR J0737−3039A observed at multiple frequencies with the MeerKAT UHF- and L-band receivers.

2.4. DM variation

The wide-band observation and high precision of MeerKAT telescope make it possible to obtain an accurate DM measurement on a per-epoch basis so as to minimise the influence of DM noise in the data. To do so, we fit for only DM and spin frequency ν for each observing epoch using 4-min TOAs, and keep the other parameters fixed. The DM measurements are shown in Fig. 2. Following Kramer et al. (2021a), we use a modified version of TEMPO for our timing analysis, which corrects dispersive delays for each TOA based on the exact DM measurement of that epoch.

Fig. 2. DM measurement per observing epoch relative to the reference value 48.913 pc cm−3 (see Sect. 2.4).

One should note that our data set does not show the apparent DM variation as a function of orbital phase, as was seen in Ransom et al. (2004). It had been demonstrated that this effect occurs due to an unaccounted Doppler shift of the observational frequency as the pulsar moves in a binary system⁶, and this will be revisited by Hu, Porayko et al. (in prep.). A more thorough investigation of this effect, as well as the frequency-dependent orbital smearing (see Sect. 2.2), is ongoing and will be presented in detail in the future publication.

3. Signal propagation effects at superior conjunction

In this section, we recapitulate the necessary concepts of signal propagation effects in the double pulsar, including the NLO contributions in the Shapiro delay and aberration delay. These concepts were described in greater detail in Kramer et al. (2021a).

Being a nearly edge-on binary system (i.e. i ∼ 90°), the curved spacetime of the companion star (pulsar B) has a significant effect on the propagation of the pulsar’s signal. To leading-order this is the well-known Shapiro delay (Shapiro 1964), which is expressed in the following form for binary pulsars (Blandford & Teukolsky 1976; Damour & Deruelle 1986):

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{S}}^{(\mathrm{LO})}& =-2rln{\mathrm{\Lambda }}_u,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_u& =1-e_{\mathrm{T}}cosu-s[sin\omega (cosu-e_{\mathrm{T}})\hfill \\ \hfill & +{(1-e_{\mathrm{T}}^2)}^{1/2}cos\omega sinu].\hfill \end{array}$$(2)

Here, u denotes the eccentric anomaly (from Kepler’s equation with eccentricity e_T), and ω denotes the longitude of periastron measured from the ascending node. The time eccentricity e_T corresponds to the eccentricity parameter in the Damour-Deruelle (DD) timing model (Damour & Deruelle 1986) that can be fitted in the pulsar timing software TEMPO or TEMPO2 (Hobbs et al. 2006). The two post-Keplerian (PK) parameters r and s represent the ‘range’ and ‘shape’ of Shapiro delay, respectively. The shape parameter is generally identified with the sine of the orbital inclination i as s ≡ sin i, whereas the range parameter is linked to the mass of the companion m_B, which in GR follows r = T_⊙ m_B. The constant $T_{\odot }\equiv {(\mathcal{GM})}_{\odot }^{\mathrm{N}}/c^3$, where c is the speed of light in vacuum and ${(\mathcal{GM})}_{\odot }^{\mathrm{N}}\equiv 1.3271244\times {10}^{26}$ cm³ s⁻² is the nominal solar mass parameter defined by the IAU 2015 Resolution B3 (Prša et al. 2016). Throughout the paper, all masses expressed in solar mass M_⊙ are referred to the nominal solar mass by taking the ratio $G{m}_i/{(\mathcal{GM})}_{\odot }^{\mathrm{N}}(i=A,B)$, where G is the gravitational constant.

The leading-order expression Eqs. (1) and (2) were obtained by integrating along a straight line (in harmonic coordinates) and assuming a static mass distribution when the pulsar’s signal propagates through the system (Blandford & Teukolsky 1976). In reality, the pulsar’s signal propagates along a curved path due to the deflection in the gravitational field of the companion and leads to a ‘lensing correction’ to the Shapiro delay. This actually results in a reduced propagation time as a consequence of Fermat’s principle (Perlick 2004). The effect of lensing is not yet observable in any pulsar systems, but for completeness, one can extend Eq. (2) by an adapted version of the approximation in Klioner & Zschocke (2010, Eq. (73)): ${\mathrm{\Lambda }}_u\to {\mathrm{\Lambda }}_u+\delta {\mathrm{\Lambda }}_u^{\text{len}}$ with

$$\begin{array}{c}\hfill \delta {\mathrm{\Lambda }}_u^{\mathrm{len}}=2rc/a_{\mathrm{R}},\end{array}$$(3)

where the semi-major axis of the relative orbit a_R = (x + x_B)/s, with x and x_B⁷ being the projected semi-major axes of pulsar A and pulsar B, respectively. For the double pulsar, one needs to account for the fact that the companion star moves while the pulsar’s signal propagates across the system. This effect is known as the ‘retardation effect’ or 1.5PN correction to the Shapiro delay (Kopeikin & Schäfer 1999; Rafikov & Lai 2006a). To sufficient approximation, the signal propagation delay can be extended to

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{S}}=-2rln({\mathrm{\Lambda }}_u+\delta {\mathrm{\Lambda }}_u^{\mathrm{len}}+\delta {\mathrm{\Lambda }}_u^{\mathrm{ret}}),\end{array}$$(4)

where the retardation correction $\delta {\mathrm{\Lambda }}_u^{\text{ret}}$ can be taken directly from Kopeikin & Schäfer (1999, Eq. (130)) as

$$\begin{array}{cc}\hfill \delta {\mathrm{\Lambda }}_u^{\mathrm{ret}}=& \frac{2\pi }{s}\frac{x}{P_{\mathrm{b}}}\frac{m_{\mathrm{A}}}{m_{\mathrm{B}}}{e}_{\mathrm{T}}sinu-\frac{2\pi s}{{(1-e_{\mathrm{T}}^2)}^{1/2}}\frac{x}{P_{\mathrm{b}}}\frac{m_{\mathrm{A}}}{m_{\mathrm{B}}}\hfill \\ \hfill & \times [sin\omega (cosu-e_{\mathrm{T}})+{(1-e_{\mathrm{T}}^2)}^{1/2}cos\omega sinu]\hfill \\ \hfill & \times [e_{\mathrm{T}}cos\omega +\frac{(cosu-e_{\mathrm{T}})cos\omega -{(1-e_{\mathrm{T}}^2)}^{1/2}sin\omega sinu}{1-e_{\mathrm{T}}cosu}].\hfill \end{array}$$(5)

The quantity P_b denotes the orbital period and m_A denotes the mass of pulsar A. It should be noted that in the double pulsar, the mass ratio m_A/m_B can be obtained in a theory-independent way (Kramer et al. 2006a; Damour 2007). Hence, apart from the Shapiro shape parameter s, Eq. (5) contains only Keplerian parameters.

Moreover, the classical aberration expression (Smarr & Blandford 1976) assumes a flat spacetime for the propagation of the pulsar signals, which is no longer sufficient for describing the observations of the double pulsar, particularly near the superior conjunction of pulsar A. One needs to account for the gravitational deflection of the pulsar’s signal caused by its companion (Doroshenko & Kopeikin 1995; Rafikov & Lai 2006b), which adds a lensing correction to the classical aberration. For pulsar A the misalignment angle between its spin vector and the orbital angular momentum is very small (< 3.2°, Ferdman et al. 2008, 2013), which is in line with a low-kick birth event (cf. Piran & Shaviv 2004; Willems & Kalogera 2004; Willems et al. 2006; Stairs et al. 2006; Tauris et al. 2017). Since the spin of A is practically parallel to the orbital angular momentum, the aberration delay can be simplified as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{A}}=\mathcal{A}(sin\psi +e_{\mathrm{T}}sin\omega )+{\delta }_{\mathrm{A}}^{\mathrm{londef}}.\end{array}$$(6)

The first term on the right-hand side of Eq. (6) is the classical aberration delay, where ψ = ω + θ is the longitude of the pulsar with respect to the ascending node (θ is the true anomaly, which defines the angle between the direction of the pulsar and the periastron), and the aberration coefficient

$$\begin{array}{c}\hfill \mathcal{A}=\frac{x}{\nu P_{\mathrm{b}}{(1-e_{\mathrm{T}}^2)}^{1/2}{sin}^{2}i}\simeq 3.65\mathrm{\mu }\mathrm{s}.\end{array}$$(7)

As 𝒜 is practically not observable and can be absorbed by a shift in various timing parameters (see discussions in Damour & Deruelle 1986; Damour & Taylor 1992), we a priori add the aberration coefficient 𝒜 as a fixed parameter in our timing model with the value given in Eq. (7).

The second term ${\delta }_{\mathrm{A}}^{\mathrm{londef}}$ in Eq. (6) is the higher-order correction originating from the ‘longitudinal deflection delay’, and can be written as (Doroshenko & Kopeikin 1995)

$$\begin{array}{c}\hfill {\delta }_{\mathrm{A}}^{\mathrm{londef}}=\mathcal{D}\frac{cos(\psi +\delta {\psi }^{\mathrm{ret}})}{{\mathrm{\Lambda }}_u+\delta {\mathrm{\Lambda }}_u^{\mathrm{ret}}},\mathrm{with}\mathcal{D}=\frac{1}{\pi \nu }\frac{r}{x+x_{\mathrm{B}}}·\end{array}$$(8)

Just as in the Shapiro delay, retardation correction is also accounted for here. As a sufficiently good approximation, the position of B when the signal reaches its minimum distance from B can be used (retardation corrected position; cf. Kopeikin & Schäfer 1999; Rafikov & Lai 2006a). The angle δψ^ret denotes the retardation related correction for the angle between the (coordinate) vector from B to A and the ascending node.

As already discussed in Kramer et al. (2021a), the NLO contributions in the Shapiro and aberration delays cannot be tested separately in the double pulsar due to the similarity of their effects on signal propagation. In addition, the lensing correction to the propagation delay (Eq. (3)) is challenging to measure as it can be absorbed in the fit of Shapiro shape s (see Sect. 5.3 and discussions in Kramer et al. 2021a). Therefore, to test the significance of the NLO contributions and to obtain an unbiased timing result, a common factor q_NLO is multiplied by these contributions and can be fitted for in our timing model:

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_u^{\mathrm{ret}}& ={\mathrm{\Lambda }}_u^{\mathrm{ret}}\times q_{\mathrm{NLO}},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_u^{\mathrm{len}}& ={\mathrm{\Lambda }}_u^{\mathrm{len}}\times q_{\mathrm{NLO}},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\delta }_{\mathrm{A}}^{\mathrm{londef}}& ={\delta }_{\mathrm{A}}^{\mathrm{londef}}\times q_{\mathrm{NLO}}.\hfill \end{array}$$(11)

In GR, the scaling factor q_NLO = 1. Figure 3 illustrates the different effects related to signal deflection in the double pulsar system.

Fig. 3. Simplified illustration of effects related to the deflection of A’s radio signals (solid red) in the gravitational field of B (top down and side perspective). The observer is located at a large distance along the x-axis. Apart from modifications in the propagation time due to a curved path in the gravitational field of B (lensing), one has a longitudinal deflection delay (${\delta }_{\mathrm{A}}^{\mathrm{londef}}$) due to the fact that the pulsar has to rotate by more than 360° between two pulses while approaching the conjunction. After conjunction, it is less than 360°, which makes pulsar signals arrive earlier at the observer. In addition, there is a latitudinal effect, due to a latitudinal shift in the emission direction towards the observer. This can lead to changes in the pulse profile since the line of sight cuts a different part of the emission region, which can also be accompanied by changes in the pulse arrival times (more details in Sects. 5.1 and 5.2).

4. Timing results

For the timing analysis, we use the timing model in TEMPO known as DDS, which is a modification of the DD model (Damour & Deruelle 1986) that uses a different parameterisation of the Shapiro delay. In DDS, the Shapiro shape parameter s is replaced by the logarithmic Shapiro shape parameter z_s via

$$\begin{array}{c}\hfill z_{\mathrm{s}}\equiv -ln(1-s),\end{array}$$(12)

which is more suitable when s is very close to one (see Kramer et al. 2006b, 2021a). The NLO contributions in the Shapiro and aberration delays are also implemented in the latest DDS model, which can be measured through a common factor q_NLO. Because the analytic inversion of the timing model developed in Damour & Deruelle (1986) is no longer sufficient for the double pulsar, primarily due to the NLO contributions, a numerical inversion of the timing model was also implemented in the latest DDS model in TEMPO (see Kramer et al. 2021a).

4.1. Timing parameters

In our analysis, the full MeerKAT data set shows a large deviation in the Shapiro range parameter r compared to the 16 yr result (Kramer et al. 2021a). We perform a drop-out analysis by removing each observing epoch and fitting the parameters. We find that r is dependent on specific observing epochs, where six epochs affect r by a significant amount while the rest of the epochs do not. These six epochs are marked as black crosses in Fig. 4. After excluding all these six epochs, r is consistent with the 16-yr result and the mass measurement in GR (Kramer et al. 2021a). Even though a number of tests and simulations have been made, we are still unclear about the cause of this problem. Possible reasons could be systematic errors in the observations or folding techniques. It should be noted that all data were folded with the TEMPO2 phase predictor during observations, which has shown outliers in this pulsar and has been doubly confirmed by our simulations. These outliers disappear after reinstalling a TEMPO polyco ephemeris in data processing (see Sect. 2.2), but we cannot rule out underlying problems due to the folding technique. The results shown here are based on data processed with the polyco scheme⁸. In any case, this issue should not affect the measurement of NLO signal propagation effects, which is the main focus of this paper. Therefore, we leave this question to future studies. In the following analysis, we use a trimmed data set that excludes these six epochs.

Fig. 4. Post-fit residuals of PSR J0737−3039A using the DDS binary model as a function of time (top panel) and the orbital phase of pulsar A with respect to the ascending node ψ (bottom panel). The MeerKAT L-band data are plotted in blue, whereas the UHF-band data are in red. The epochs of the excluded six observations are marked as black crosses.

Table 2 and Fig. 4 present the results obtained from fitting the TEMPO DDS model to the trimmed MeerKAT data set. We fix the proper motion (μ_α, μ_δ) and parallax π_x to the more precise values determined from the 16 yr timing and the very long baseline interferometry measurements (see Kramer et al. 2021a). As the time span of our data is not sufficient to obtain a reliable measurement of the orbital period derivative Ṗ_b and the orientation of the orbit (ω₀ ≈ 0°) is not at a favourable position for a precise measurement of the Einstein delay amplitude γ_E, we choose to fix these parameters to the more precise measurements from 16 yr data (Kramer et al. 2021a). Fixing the above parameters has no impact on the measurements of signal propagation effects and masses. The two PK parameters that describe the relativistic deformation of the orbit, δ_r and δ_θ (Damour & Deruelle 1986), are also held fixed at the GR value in our analysis, as δ_r cannot be measured (see Kramer et al. 2021a) and δ_θ is not yet measurable with the current MeerKAT data.

The values shown in Table 2 are the result of 1000 Monte Carlo (MC) runs, where in each run, a random realisation of proper motion, parallax, DM, γ_E, and Ṗ_b is selected. The DM value is selected according to the DM measurements and uncertainties shown in Fig. 2. We use the aforementioned modified version of TEMPO to correct the DM for each TOA and fit for all other timing parameters in each run. The numbers shown in Table 2 are the mean values of the distribution of each parameter after 1000 MC runs, whereas the uncertainties are taken from the larger one among the standard deviation of the distributions and the maximum error from TEMPO in all MC runs.

In order to allow direct comparisons with previous publications, parameters shown in Table 2 are measured with respect to the same epochs and the same terrestrial time standard UTC (NIST)⁹ within the ‘Barycentric Dynamical Time’ (TDB) timescale as implemented in TEMPO. Even though TDB runs at a slower rate than ‘Barycentric Coordinate Time (TCB)’, which was recommended by IAU 2006 Recolution B3¹⁰ (see also Soffel et al. 2003), this choice does not have any impact on the results presented in this paper (see discussions in Kramer et al. 2021a). To transfer the topocentric TOAs to the Solar System Barycentre (SSB), the DE436 Solar System ephemeris published by the Jet Propulsion Laboratory¹¹ is used.

All binary parameters in Table 2 are consistent with the 16-yr data, except for x being different by ∼3σ. This is because x is highly correlated with γ_E, which is kept fixed in our fit. This should be improved in the future, once we have enough MeerKAT data to fit for x and γ_E simultaneously. In our fit, the root mean square (rms) is very close to the mean TOA uncertainty, and the reduced χ² of the individual observation is close to one, suggesting that our result is not affected by jitter noise. We also perform simulations and single-pulse analysis following the methods in Parthasarathy et al. (2021) and find little evidence of jitter noise.

The rms of the MeerKAT data shown in Table 4 is more than two times better than that of the Green Bank Telescope (see Kramer et al. 2021b). Thanks to this much improved precision, the measurements of the Shapiro parameters improve quickly. Compared to Kramer et al. (2021a), the shape parameter s improves by a factor of two and the range parameter r improves by a factor of 1.3 (see Table 2).

4.2. Mass measurements

The standard approach for measuring the masses of a binary pulsar system is using two PK parameters. Assuming GR, one can calculate the two a priori unknown masses based on the measurements of Keplerian parameters. For the double pulsar, the two most precisely measured PK parameters are periastron advance $\dot{\omega }$ and the Shapiro shape parameter s.

For the advance of periastron, in addition to the 1PN contribution, we also need to account for higher-order corrections due to 2PN effects and Lense-Thirring (LT) precession caused by spin-orbit coupling of pulsar A, as they are much larger than the measurement error of $\dot{\omega }$ (see Hu et al. 2020; Kramer et al. 2021a). For the analysis of this paper, the total intrinsic contribution to the periastron advance can be expressed, with sufficient precision, as (Damour & Schäfer 1988)

$$\begin{array}{c}\hfill \dot{\omega }={\dot{\omega }}^{1\mathrm{PN}}+{\dot{\omega }}^{2\mathrm{PN}}+{\dot{\omega }}^{\mathrm{LT},\mathrm{A}}.\end{array}$$(13)

The first and second post-Newtonian (PN) terms ${\dot{\omega }}^{1\mathrm{PN}}$ and ${\dot{\omega }}^{2\mathrm{PN}}$ are functions of masses and observed Keplerian parameters. The situation is more complicated for the LT contribution ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}$, as it requires the knowledge of the moment of inertia (MOI) of pulsar A, which is still not very constrained because of our limited knowledge of the equation of state (EOS) of neutron stars. As discussed in detail in Hu et al. (2020), one could measure the masses and the MOI simultaneously by introducing a third PK parameter Ṗ_b into the $\dot{\omega }-s$ test. Such a test has already been carried out by Kramer et al. (2021a) using the 16 yr data with an upper limit obtained: I_A < 3.0 × 10⁴⁵ g cm² with 90% confidence. This measurement is expected to improve with the combination of the 16 yr data with MeerKAT data in a forthcoming paper, and should improve considerably in the coming years as more data becomes available. This promises an important complementary constraint on the EOS (Hu et al. 2020). For the calculations here, we take the value of ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}$ from Eq. (35) in Kramer et al. (2021a), which uses the constraints on the EOS from Dietrich et al. (2020):

$$\begin{array}{c}\hfill {\dot{\omega }}^{\mathrm{LT},\mathrm{A}}=-4.83(+29,-35)\times {10}^{-4}{\mathrm{deg}\mathrm{yr}}^{-1}.\end{array}$$(14)

The Shapiro shape parameter s is the sine of the orbital inclination i. In Newtonian gravity, the orbital inclination is linked to the projected semi-major axis x via the binary mass function (e.g. Lorimer & Kramer 2004):

$$\begin{array}{c}\hfill sini=\frac{{(n_{\mathrm{b}}M)}^{2/3}x}{T_{\odot }^{1/3}{m}_{\mathrm{B}}},\end{array}$$(15)

where x and the orbital frequency n_b ≡ 2π/P_b are both measured Keplerian parameters. Equation (15) gets modified by a 1PN term in the 1PN approximation for Kepler’s third law (see Eq. (3.7) in Damour & Deruelle 1985 and Eq. (3.9) in Damour & Taylor 1992):

$$\begin{array}{c}\hfill sini=\frac{{(n_{\mathrm{b}}M)}^{2/3}x}{T_{\odot }^{1/3}{m}_{\mathrm{B}}}[1+(3-\frac{m_{\mathrm{A}}{m}_{\mathrm{B}}}{3M^2}){\left(T_{\odot }M{n}_{\mathrm{b}}\right)}^{2/3}].\end{array}$$(16)

Taking the measurements of P_b, x, and masses from Table 2, one can calculate that the 1PN correction is approximately 1.27 × 10⁻⁵. This correction was considered for the first time in pulsar analysis by Kramer et al. (2021a), where the significance was about 1.3σ. Now with MeerKAT data, this 1PN correction is 2.5σ significant and cannot be ignored in the analysis. We hereby use the full 1PN mass function Eq. (16) to measure the masses.

Combining the PK parameters $\dot{\omega }$ and s, one obtains the (Doppler-shifted) masses, which are listed in Table 2. These measurements are fully consistent with those obtained with 16 yr data (Kramer et al. 2021a), and the precision of m_A and m_B are also better.

Alternatively, one can fit for masses using the timing model known as DDGR (Taylor & Weisberg 1989), which is based on the DD model where the PK parameters are explicitly calculated from the masses and the Keplerian parameters, assuming GR. Beside the Keplerian parameters, it fits explicitly for the total mass M and the companion mass m_B. To make the measurements, we modify the DDGR model so that it incorporates all NLO contributions that need to be accounted for in this system, including NLO signal propagation effects, LT contribution ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}$, NLO gravitational wave damping and mass-loss contribution to Ṗ_b (see Hu et al. 2020; Kramer et al. 2021a). An MOI needs to be provided to the model for the calculation of ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}$ and the mass loss contribution to Ṗ_b. For periastron advance $\dot{\omega }$, the uncertainty from the MOI is still smaller than that from MeerKAT observations (see Eq. (14) and Table 2). Therefore, based on the EOS constraint from Dietrich et al. (2020), we fix the MOI to I_A = 1.28 × 10⁴⁵ g cm² in our fit. Table 3 shows the mass measurements obtained using the DDGR model. The results are fully consistent with the measurements derived from the DDS model, with smaller uncertainties in m_B and M.

Following Kramer et al. (2021a), one could test the agreement of r in GR by comparing $m_{\mathrm{B}}^{(r)}=r/T_{\odot }=1.2512(33)M_{\odot }$ (cf. Table 2) with the companion mass determined here, which gives

$$\begin{array}{c}\hfill r^{\mathrm{obs}}/r^{\mathrm{GR}}=1.0019(26).\end{array}$$(17)

This leads to a 5.3 × 10⁻³ (95% confidence) test of GR.

5. Studying NLO signal propagation effects

Because the double pulsar system is nearly edge-on to our line of sight (LOS, see i in Table 2), it is ideal for measuring signal propagation effects caused by the gravitational field of the companion near superior conjunction. The leading-order expression Eq. (1) is no longer sufficient for describing the signal propagation in the double pulsar. Such a model would result in significant residuals near superior conjunction when aggregating residuals in the orbital phase, as shown in Fig. 5. These residuals agree very well with the expected NLO contributions discussed in Sect. 3, which are shown by the red curve. The significance of the NLO corrections can be tested by scaling these corrections collectively by a common factor q_NLO (cf. Eqs. (9)–(11); q_NLO = 1 in GR) and fitting for it. We find, with the trimmed data set,

$$\begin{array}{c}\hfill q_{\mathrm{NLO}}=0.999(79),\end{array}$$(18)

Fig. 5. Aggregated residuals (blue) due to NLO contributions in Shapiro and aberration delay, shown in the orbital phase ψ. Residuals are re-scaled by (1 + eT cos θ)−1 to account for secular variations in amplitude due to the precession of periastron. The black curve indicates the fitted qNLO (see Table 2) with the 2σ range shown by the grey shaded areas, which agrees very well with the theoretical prediction indicated by the red dotted line.

which has surpassed the 16 yr result by 1.65 times with only ∼3 yr of data, thanks to the much improved precision offered by MeerKAT.

Following the definition of q_NLO in Sect. 3, a fit for this parameter involves two aspects of gravity: a 1.5PN correction of the Shapiro delay due to the movement of the companion $\delta {\mathrm{\Lambda }}_u^{\text{ret}}$, and corrections related to the signal deflection in the gravitational field of the companion $\delta {\mathrm{\Lambda }}_u^{\text{len}}$ and ${\delta }_{\mathrm{A}}^{\mathrm{londef}}$. Even though these contributions cannot be tested individually in a simultaneous fit because of their similarity, one can still test one at a time while keeping the other one fixed (cf. Kramer et al. 2021a). We find

$$\begin{array}{cc}\hfill q_{\mathrm{NLO}}[\mathrm{deflection}]& =1.00(15),\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill q_{\mathrm{NLO}}[\mathrm{retardation}]& =1.00(17).\hfill \end{array}$$(20)

5.1. Searching for profile variation at eclipse

The lensing correction to the aberration delay may not only lead to a shift in time in the longitudinal aspect but can also result in a change of the co-latitude of the emission direction towards Earth, namely the ‘latitudinal deflection delay’. This would cause profile variations as the LOS cuts a different region of the pulsar beam (Rafikov & Lai 2006a,b). An illustration of the latitudinal deflection effect is shown in the right panel of Fig. 3. Various analyses have confirmed that pulsar A is an orthogonal rotator (Guillemot et al. 2013; Ferdman et al. 2013; Kramer et al. 2021b), meaning that the main pulse and the interpulse come from opposite magnetic poles. Therefore, we do not expect shifts of pulse components in phase, as discussed in Rafikov & Lai (2006b), based on the (incorrect) assumption of an aligned rotator suggested by Demorest et al. (2004).

The profile variation is expected to be maximum at the superior conjunction and symmetric around ψ = 90° (retardation corrected). This study requires high time resolution, for which we use the search mode data. We select the data that are near the eclipses and fold them into single pulses using TEMPO polyco (with TSPAN = 1 min). Data are then combined, cleaned, and polarisation calibrated before integrated into total intensity and averaged in frequency. As the single pulses are still very weak, we average over every eight pulses to increase the signal-to-noise ratio (S/N). An example of eclipse data is shown in Fig. 6.

Fig. 6. Example of one eclipse observation at UHF band, plotted in intensity against orbital phase ψ and pulse phase. The intensity modulation occurs when pulsar A is eclipsed by the magnetosphere of pulsar B. Each integration is a sum of eight pulses. When plotting, the discontinuities between recordings are patched with the previous sub-integration.

In order to get a high S/N profile, we first mask the regions where pulsar A’s emission is blocked by pulsar B (see Fig. 7) and split the data into orbital phases with a step of Δψ = 0.25° for 89° < ψ < 91°. Then for each phase interval, we integrate pulses from all observations of a given band (L or UHF) together to increase the S/N. The resulting profiles are shown in Fig. 8 for the L-band data and in Fig. 9 for the UHF-band data, which are summed from 25 and 37 eclipses respectively. The difference between the added profiles at the eclipse and a 2 h integrated profile at the non-eclipse part of the orbit is insignificant. The subtle residual structures in these figures can result from interstellar medium effects (DM variation and scintillation) based on our simulations. Therefore, we conclude that the current data are not (yet) sensitive to profile variations caused by the latitudinal aberration delay, or are not significant in the region that is seen by our LOS. These profiles will be used in a subsequent study on the geometry of the system.

Fig. 7. Same as for Fig. 6, but masking out the regions without pulses (blocked by the magnetosphere of pulse B). The dashed red lines indicate the orbital phase bins used in Figs. 8 and 9.

Fig. 8. Profile variation with respect to orbital phase. Left: black lines indicate the integrated profile of L-band data summed from 25 eclipses for orbital phases 89° < ψ < 91° with an interval of Δψ = 0.25°. The baseline of the profile is placed at the centre of each interval, and the numbers on the right side of the profiles (np) indicate the estimated upper limit for the number of pulses in that interval. The red curves indicate the reference (standard) profile integrated over a 2 h observation excluding the eclipse part. Right: residuals of the added profile with respect to the reference profile.

Fig. 9. Same as for Fig. 8, but for the UHF-band data summed from 37 eclipses.

5.2. Simulation for latitudinal deflection delay

To investigate whether the deflection delay caused by latitudinal deflection is measurable from pulsar timing, we perform a simulation based on a simple emission model, which consists of a set of circular cones. Following Doroshenko & Kopeikin (1995) and Rafikov & Lai (2006b), the latitudinal deflection delay for pulsar A can be written as

$$\begin{array}{c}\hfill {\delta }_{\mathrm{A}}^{\mathrm{latdef}}=-\mathcal{D}\frac{cosisin(\psi +\delta {\psi }^{\mathrm{ret}})}{({\mathrm{\Lambda }}_u+\delta {\mathrm{\Lambda }}_u^{\mathrm{ret}})tan{\chi }_0},\end{array}$$(21)

where χ₀ is the angle between the arc connecting the LOS and spin axis and the arc connecting LOS and magnetic axis. It should be noted that Eq. (21) is based on the approximation of Doroshenko & Kopeikin (1995) for the deflection angle, and therefore assumes that the impact parameter is (sufficiently) large compared to the Einstein radius. While this is sufficient, at least until full SKA becomes operational, we present an improved description further below in Sect. 5.4.

We include this deflection time delay in our test model assuming χ₀ = 45°  ¹² (i.e. a relatively large latitudinal deflection delay) and scale it with a factor q^latdef. We simulate high-precision TOAs using this test model and fit for the scaling factor. The pre-fit residuals show an advance signature with an amplitude of −2.8 μs and symmetric to ψ = 90°, which is in a similar shape to Shapiro delay but with the opposite sign and a smaller amplitude. However, after fitting for Shapiro parameters, the signature gets mostly absorbed and leaves residuals below 42 ns at superior conjunction.

5.3. Prospects of lensing measurement

Even though the retardation and deflection effects can be tested separately while keeping the other one fixed, as shown in Eqs. (19) and (20), measuring the lensing correction to Shapiro delay $\delta {\mathrm{\Lambda }}_u^{\text{len}}$ independently is challenging. As already pointed out by Kramer et al. (2021a), this effect is difficult to observe because of its strong covariance with s, or equivalently z_s. Our simulation also confirms that the lensing signature can be mostly absorbed by z_s in timing due to its symmetry with respect to conjunction. For demonstration purposes, we simulate 1-ns TOAs, which model all NLO signal propagation contributions, then keep the retardation and deflection delay fixed in the model and fit for the scaling factor attached to the lensing correction, q^len (corresponding to q_NLO in Eq. (10)). Figure 10 shows the residuals when lensing correction is not taken into account, leading to a reduced propagation time of about 850 ns as a result of Fermat’s principle. However, after fitting for z_s, the lensing signature gets absorbed and leaves the residuals to be below 70 ns (Fig. 11), making a detection with the current timing precision certainly impossible.

Fig. 10. Lensing simulation: pre-fit residual plotted against orbital phase ψ. Data displayed here are centred on ω = 180° and span a decade. The scattering at the lower end of the curve is due to the precession of periastron.

Fig. 11. Lensing simulation: post-fit residual plotted against orbital phase ψ.

To investigate whether lensing can be measured separately in the near future, we simulate TOAs for MeerKAT, MeerKAT extension, and the first phase of the SKA (SKA 1) until 2030 based on similar assumptions made in Hu et al. (2020). In addition, as the TOA precision reduces significantly due to the intermittent signals during the eclipse (see Fig. 6), we account for this in our simulations by increasing the uncertainty of these TOAs based on MeerKAT observations. As a simple estimate, we assume GR and perform the simulation using the modified DDGR model with a grid fit to q^len. If lensing is measurable, the value of q^len should be close to one. However, it turns out that with the observed and simulated data from 2019 to 2030, the uncertainty of q^len is still larger than one.

To further push the precision, we assume that an instrument will be available in the future that is capable of providing a timing accuracy one order of magnitude better than that of the SKA 1 (i.e. 100 ns for an integration time of 30 s), for example, a future full SKA. As a rough estimate, here we only consider radiometer noise and ignore any other noise sources, such as jitter noise or scintillation noise. The uncertainty of q^len against the time span of the simulated data is shown in Fig. 12. With such precision, one would expect to get a 5-σ test of lensing with ∼4 yr of data.

Fig. 12. Uncertainty of factor qlen as a function of time span for the simulated data assumed to be ten times better than the SKA 1.

From the simulation, we also obtain an estimated uncertainty for the common factor of NLO contributions q_NLO in the near future. Assuming no jitter noise, with MeerKAT and the SKA 1, we can expect a 50-σ detection by 2030.

5.4. Improvements in the timing model for ≲50 ns precision

Equations (8) and (21) are based on the approximation for the signal deflection used by Doroshenko & Kopeikin (1995). As discussed in detail in Kramer et al. (2021a), this is still sufficient for the analysis of current timing data. For that reason, the analysis in this paper is still based on Doroshenko & Kopeikin (1995), which (including corrections for retardation) is already part of the TEMPO distribution. However, in the future we can expect to obtain a timing precision of better than ∼50 ns near conjunction (±1°), so that an improved treatment of the deflection is required. In a series of papers, Rafikov and Lai used the standard lensing equation to treat the signal propagation in the double pulsar near conjunction (Lai & Rafikov 2005; Rafikov & Lai 2006a,b). This allowed them to drop the assumption that the impact parameter is much larger than the Einstein radius. The standard lensing equation, however, is based on the assumption of small angles (see e.g. Schneider et al. 1992). Therefore, strictly speaking, Lai and Rafikov’s calculations are only valid near conjunction. Similar to the calculations of Shapiro (1967) and Ward (1970), Wucknitz (2008) studied the deflection of photons in the gravitational field of a ‘point mass’ for general lensing scenarios not limited to regions close to the optical axis. Based on these results, we give an analytic expression for the signal deflection that is valid for the whole orbit, and recovers the calculations by Lai and Rafikov near conjunction, and those of Doroshenko & Kopeikin (1995) if the impact parameter is large compared to the Einstein radius.

In the following paragraphs, Θ is the angle between r_AB′ and the direction towards the observer, where r_AB′ denotes the vector from the position of pulsar A at emission to the (retardation corrected) position of pulsar B (the underlying geometry for our calculations is illustrated in Fig. 13). The deflection ΔΘ of A’s radio signal by pulsar B can be obtained from Eq. (24) in Wucknitz (2008), with the replacements α → ΔΘ, θ → Θ + ΔΘ, and $m\to {\alpha }_{\text{E}}^2$, where the quantity α_E is the angle corresponding to the Einstein radius, and is given by

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{E}}=\frac{2}{c}\sqrt{\frac{G{m}_{\mathrm{B}}}{|{\mathit{r}}_{\mathrm{AB}\prime }|}}\ll 1.\end{array}$$(22)

Fig. 13. Schematic picture of the lensing geometry as used in Sect. 5.4. B′ denotes the retardation corrected position of B (cf. Klioner & Kopeikin 1992; Kopeikin & Schäfer 1999). In principle there is a second photon path towards the observer (below B′). However, for the double pulsar, this signal is not only significantly weaker, but the path also comes so close to pulsar B that the photons are absorbed by the plasma-filled magnetosphere of B (cf. Lai & Rafikov 2005; Rafikov & Lai 2006b).

This is the maximum value ΔΘ can assume. The distance D_d in Wucknitz (2008) corresponds to |r_AB′|¹³. Consequently, one obtains

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Theta }sin(\mathrm{\Theta }+\mathrm{\Delta }\mathrm{\Theta })-\frac{{\alpha }_{\mathrm{E}}^2}{2}[1+cos(\mathrm{\Theta }+\mathrm{\Delta }\mathrm{\Theta })]=0.\end{array}$$(23)

Since ΔΘ ≤ α_E ≪ 1 for all angles of Θ, we can expand the equation above in ΔΘ while keeping terms only up to order ${\alpha }_{\text{E}}^2$. This leads to a quadratic equation:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Theta }}^2+\mathrm{\Delta }\mathrm{\Theta }sin\mathrm{\Theta }-\frac{{\alpha }_{\mathrm{E}}^2}{2}(1+cos\mathrm{\Theta })\simeq 0,\end{array}$$(24)

which has, under the assumptions made, one solution:

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Theta }\simeq \frac{1}{2}(\sqrt{{sin}^2\mathrm{\Theta }+2{\alpha }_{\mathrm{E}}^2(1+cos\mathrm{\Theta })}-sin\mathrm{\Theta }).\end{array}$$(25)

For Θ ≪ 1 this agrees with the standard lensing equation (see e.g. Schneider et al. 1992).

The angle Θ ∈ [0, π] needs to be computed from the retardation-corrected orbital phase via cosΘ = sin isin(ψ + δψ^ret). The longitudinal and latitudinal deflection delay are given by

$$\begin{array}{cc}\hfill {\delta }_{\mathrm{A}}^{\mathrm{londef}}& =\frac{\mathrm{\Delta }\mathrm{\Theta }}{2\pi \nu }\frac{cos(\psi +\delta {\psi }^{\mathrm{ret}})}{sin\mathrm{\Theta }sini},\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill {\delta }_{\mathrm{A}}^{\mathrm{latdef}}& =-\frac{\mathrm{\Delta }\mathrm{\Theta }}{2\pi \nu }\frac{sin(\psi +\delta {\psi }^{\mathrm{ret}})}{sin\mathrm{\Theta }tanitan{\chi }_0},\hfill \end{array}$$(27)

respectively (cf. Eqs. (10) and (24) in Rafikov & Lai 2006b), with ζ = π − i and η = −π/2 spin of A aligned with orbital angular momentum). If Θ is much larger than α_E, one has $\mathrm{\Delta }\mathrm{\Theta }\simeq {\alpha }_{\text{E}}^2(1+cos\mathrm{\Theta })/(2sin\mathrm{\Theta })$. This corresponds to the approximation of Doroshenko & Kopeikin (1995) for the deflection angle, and one recovers Eqs. (8) and (21).

6. Discussion

In this paper, we have presented results from 3-yr timing observations of the double pulsar using the MeerKAT telescope, with a specific focus on studying higher-order signal propagation effects in the gravitational field of the companion. In order to minimise the effects from profile evolution and DM variation, we used frequency-dependent 2D templates to generate TOAs and a DM model to correct dispersive delay in TOAs.

Thanks to its high inclination and orbital compactness, the double pulsar is a unique pulsar system for testing NLO signal propagation effects in strong fields. The significantly increased precision offered by MeerKAT permits an independent verification of NLO signal propagation effects and has already surpassed the 16 yr result with only ∼3 yr of data. In our analysis, the Shapiro shape parameter s has been improved by a factor of two compared to the previous result (Kramer et al. 2021a), which also leads to a better mass measurement. The Shapiro range parameter r agrees with GR at 5.3 × 10⁻³ (95% confidence). The precision of the measurement of NLO signal propagation effects q_NLO has been improved by a factor of 1.65. In this work, we investigated the potential profile variation due to latitudinal deflection delay and the possibility of measuring lensing correction to the Shapiro delay, which has never been studied in detail before in pulsar analysis. With the current MeerKAT data, we found little evidence of profile variation at superior conjunction. It could be that the profile variation is not significant at the region we are looking at, or our current data are not sensitive enough to identify it. We also performed simulations on latitudinal deflection delay based on a simple emission model and found it unlikely to be detected because of its correlation with Shapiro delay. As for the lensing correction $\delta {\mathrm{\Lambda }}_u^{\text{len}}$, we found it can be mostly absorbed by the Shapiro shape parameter. Our simulation showed that lensing is unlikely to be measured separately from timing before the full SKA or similarly powerful instruments, and may then be measurable with a few years of timing observations if noises such as phase jitter and scintillation do not limit our precision.

However, our analysis also showed that adding certain epochs has a significant impact on the measured Shapiro parameters, but not on q_NLO. This could be due to the fact that the phase predicted using the polyco scheme is particularly worse at the superior conjunction, which caused the discrepancies in the Shapiro parameters. Comparison of polyco with different TSPAN values showed residuals oscillating near the superior conjunction, and we may already be limited by the precision of polyco scheme. Of course, there may exist other unknown systematic errors in the data.

To support our timing analysis and study of profile variations at eclipse due to latitudinal deflection, we also checked profiles from all observations. We found variations in the total profile from epoch to epoch. The differences in the profiles are more prominent at lower frequencies and broadband. Our simulation suggested that these profile variations are likely to be associated with DM variation and scintillation. Even though we have sub-banded data into 16/32 frequency bands and used 2D templates, profile variations may still have an impact on timing. The study of profile variations will be continued in a subsequent work to improve the constraint on the geometry of the system.

Moreover, although not discussed in this paper, we found evidence of red noise in the spectrum of timing residuals with an amplitude two orders of magnitude larger than for typical millisecond pulsars. If not taken into account, it may strongly affect astrometric parameters, as well as influence binary parameters, according to our simulations. This makes it more difficult to combine the 3 yr MeerKAT data set and the 16 yr data set. Given that the timing precision of the former significantly outperforms the latter, the weighting of MeerKAT data already exceeds the 16 yr data and can dominate noise modelling. For the purpose of this paper, we did not include 16 yr data because of their minor contribution to the q_NLO measurement (∼10% improvement). However, for studying secular relativistic effects, an appropriate noise modelling may be required to combine these data. We will investigate this in further ongoing studies.

In the future, continuing observations with MeerKAT and the SKA will further improve the precision of tests on signal propagation effects, and a 50-σ detection of q_NLO can be expected by 2030. For that, we have also provided an improved analytical description of the signal propagation in the double pulsar. Furthermore, as demonstrated by Hu et al. (2020), the precision of secular relativistic effects will also be greatly improved and will eventually enable the measurement of the MOI of pulsar A and the NLO gravitational wave damping in the near future. The timing precision of the MeerKAT data used in this work is even better than that assumed in Hu et al. (2020), which is based on early L-band data from MeerKAT (see Table 4). This makes their predictions conservative and we are likely to achieve even better measurements with future observations.

Acknowledgments

We acknowledge Kuo Liu for helpful discussions on data processing and analysis, and Olaf Wucknitz for carefully reading the manuscript and discussions on gravitational lensing which were particularly helpful in Sect. 5.4. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. MeerTime data is housed on the OzSTAR supercomputer at Swinburne University of Technology. H.H. is a member of the International Max Planck Research School for Astronomy and Astrophysics at the Universities of Bonn and Cologne. This work is supported by the Max-Planck Society as part of the ‘LEGACY’ collaboration with the Chinese Academy of Sciences on low-frequency gravitational wave astronomy. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. Part of this work has been funded using resources from the research grant ‘iPeska’ (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016.

References

All Tables

All Figures

	Fig. 1. Pulse profile of PSR J0737−3039A observed at multiple frequencies with the MeerKAT UHF- and L-band receivers.
In the text

	Fig. 2. DM measurement per observing epoch relative to the reference value 48.913 pc cm−3 (see Sect. 2.4).
In the text

Fig. 3. Simplified illustration of effects related to the deflection of A’s radio signals (solid red) in the gravitational field of B (top down and side perspective). The observer is located at a large distance along the x-axis. Apart from modifications in the propagation time due to a curved path in the gravitational field of B (lensing), one has a longitudinal deflection delay (${\delta }_{\mathrm{A}}^{\mathrm{londef}}$) due to the fact that the pulsar has to rotate by more than 360° between two pulses while approaching the conjunction. After conjunction, it is less than 360°, which makes pulsar signals arrive earlier at the observer. In addition, there is a latitudinal effect, due to a latitudinal shift in the emission direction towards the observer. This can lead to changes in the pulse profile since the line of sight cuts a different part of the emission region, which can also be accompanied by changes in the pulse arrival times (more details in Sects. 5.1 and 5.2).

In the text

	Fig. 4. Post-fit residuals of PSR J0737−3039A using the DDS binary model as a function of time (top panel) and the orbital phase of pulsar A with respect to the ascending node ψ (bottom panel). The MeerKAT L-band data are plotted in blue, whereas the UHF-band data are in red. The epochs of the excluded six observations are marked as black crosses.
In the text

Fig. 5. Aggregated residuals (blue) due to NLO contributions in Shapiro and aberration delay, shown in the orbital phase ψ. Residuals are re-scaled by (1 + eT cos θ)−1 to account for secular variations in amplitude due to the precession of periastron. The black curve indicates the fitted qNLO (see Table 2) with the 2σ range shown by the grey shaded areas, which agrees very well with the theoretical prediction indicated by the red dotted line.

In the text

	Fig. 6. Example of one eclipse observation at UHF band, plotted in intensity against orbital phase ψ and pulse phase. The intensity modulation occurs when pulsar A is eclipsed by the magnetosphere of pulsar B. Each integration is a sum of eight pulses. When plotting, the discontinuities between recordings are patched with the previous sub-integration.
In the text

	Fig. 7. Same as for Fig. 6, but masking out the regions without pulses (blocked by the magnetosphere of pulse B). The dashed red lines indicate the orbital phase bins used in Figs. 8 and 9.
In the text

Fig. 8. Profile variation with respect to orbital phase. Left: black lines indicate the integrated profile of L-band data summed from 25 eclipses for orbital phases 89° < ψ < 91° with an interval of Δψ = 0.25°. The baseline of the profile is placed at the centre of each interval, and the numbers on the right side of the profiles (np) indicate the estimated upper limit for the number of pulses in that interval. The red curves indicate the reference (standard) profile integrated over a 2 h observation excluding the eclipse part. Right: residuals of the added profile with respect to the reference profile.

In the text

	Fig. 9. Same as for Fig. 8, but for the UHF-band data summed from 37 eclipses.
In the text

	Fig. 10. Lensing simulation: pre-fit residual plotted against orbital phase ψ. Data displayed here are centred on ω = 180° and span a decade. The scattering at the lower end of the curve is due to the precession of periastron.
In the text

	Fig. 11. Lensing simulation: post-fit residual plotted against orbital phase ψ.
In the text

	Fig. 12. Uncertainty of factor qlen as a function of time span for the simulated data assumed to be ten times better than the SKA 1.
In the text

Fig. 13. Schematic picture of the lensing geometry as used in Sect. 5.4. B′ denotes the retardation corrected position of B (cf. Klioner & Kopeikin 1992; Kopeikin & Schäfer 1999). In principle there is a second photon path towards the observer (below B′). However, for the double pulsar, this signal is not only significantly weaker, but the path also comes so close to pulsar B that the photons are absorbed by the plasma-filled magnetosphere of B (cf. Lai & Rafikov 2005; Rafikov & Lai 2006b).

In the text