Issue 
A&A
Volume 667, November 2022



Article Number  A149  
Number of page(s)  13  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202244825  
Published online  21 November 2022 
Gravitational signal propagation in the double pulsar studied with the MeerKAT telescope
^{1}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
email: huhu@mpifrbonn.mpg.de
^{2}
Jodrell Bank Centre for Astrophysics, The University of Manchester, Oxford Road, Manchester M13 9PL, UK
^{3}
Institute of Physics, Eötvös Loránd University, Pázmány P.s. 1/A, 1117 Budapest, Hungary
^{4}
Institute for Radio Astronomy & Space Research, Auckland University of Technology, Private Bag 92006, Auckland 1142, New Zealand
^{5}
INAFOsservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius, Italy
^{6}
Australia Telescope National Facility, CSIRO Space and Astronomy, PO Box 76 Epping, NSW 1710, Australia
^{7}
Dept. of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
^{8}
Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218 Hawthorn, VIC 3122, Australia
^{9}
ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Mail H29, Swinburne University of Technology, PO Box 218 Hawthorn, VIC 3122, Australia
^{10}
South African Radio Astronomy Observatory, Cape Town 7925, South Africa
^{11}
SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK
^{12}
Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa
Received:
29
August
2022
Accepted:
22
September
2022
The double pulsar PSR J0737−3039A/B has offered a wealth of gravitational experiments in the strongfield regime, all of which general relativity has passed with flying colours. In particular, among current gravity experiments that test photon propagation, the double pulsar probes the strongest spacetime curvature. Observations with MeerKAT and, in the future, the Square Kilometre Array (SKA) can greatly improve the accuracy of current tests and facilitate tests of nexttoleadingorder (NLO) contributions in both orbital motion and signal propagation. We present our timing analysis of new observations of PSR J0737−3039A, made using the MeerKAT telescope over the last three years. The increased timing precision offered by MeerKAT yields a measurement of Shapiro delay parameter s that it twice as good, and an improved mass measurements compared to previous studies. In addition, our results provide an independent confirmation of the NLO signal propagation effects and already surpass the previous measurement from 16 yr data by a factor of 1.65. These effects include the retardation effect due to the movement of the companion and the deflection of the signal by the gravitational field of the companion. We also investigate the novel effects that have been expected. For instance, we search for potential profile variations near superior conjunctions caused by shifts of the line of sight due to latitudinal signal deflection, and we find insignificant evidence with our current data. With simulations, we find that the latitudinal deflection delay is unlikely to be measured with timing because of its correlation with Shapiro delay. Furthermore, although it is currently not possible to detect the expected lensing correction to the Shapiro delay, our simulations suggest that this effect may be measured with the full SKA. Finally, we provide an improved analytical description for the signal propagation in the double pulsar system that meets the timing precision expected from future instruments such as the full SKA.
Key words: stars: neutron / pulsars: individual: J0737–3039A / gravitation / binaries: eclipsing
© H. Hu et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model.
Open access funding provided by Max Planck Society.
1. Introduction
The double pulsar PSR J0737−3039A/B is a rich laboratory for strongfield gravity experiments. The system consists of a 23ms recycled pulsar (‘A’) and a 2.8s ‘normal’ pulsar (‘B’) in a nearly edgeon and slightly eccentric 2.45h 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 secondorder 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^{−21} cm^{−2}) 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 higherorder 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 midfrequency 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 higherorder contributions to the Shapiro and aberration delay. The timing results and mass measurements are presented in Sect. 4. We then provide an indepth study on the higherorder 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 tiedarray beamformed voltages from the correlatorbeamformer engine of the MeerKAT observing system and is capable of producing coherently dedispersed fullStokes 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 foldmode observation and another 30 min foldsearch dualmode 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 Lband 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 fullorbit timing observations and 62 searchmode 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^{1} (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^{−2}. As the Lband 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 Lband data to the same frequency channels. We treat the UHFband data in the same way, as the rolloff 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^{2} (Nice et al. 2015, more details are given in Sect. 4). Therefore, to fold the data more accurately, all data are supplied with a polycoformat 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^{3}. With the PSRCHIVE version 20220114, 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 frequencydependent orbital smearing. Standard pulsar software such as PSRCHIVE do not take this effect fully into consideration, even with the frequencyresolved TEMPO2 predictor^{4}. To avoid this issue, we first dedisperse 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^{5}. Because of the profile frequency evolution and scintillation effects, data are subbanded in frequency, with 32 subbands for the UHFband and 16 subbands for the Lband.
As for the timeaveraging, the integration time needs to be short enough in order to properly resolve the Shapiro delay and the nexttoleadingorder (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 32s integration time, consistent with the analysis by Kramer et al. (2021a). After frequency and time averaging, data are redispersed to allow measurements of dispersion measure (DM) in the timing analysis.
2.3. Wideband templates and TOA extraction
Wideband 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 frequencydependent 2D templates. With this technique, DM measurements are dependent on the DM value used to align (dedisperse) 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 Lband and UHFband data. This could potentially be solved with a simultaneous observation with Lband and UHFband 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 dedisperse the corresponding fullbandwidth data. This minimises the DM offset between Lband and UHFband data, which can be seen in Fig. 2. These data are then subbanded and averaged in time. Finally, by smoothing the profiles with psrsmooth/PSRCHIVE, we obtain 2D templates, with 16 subbands at Lband and 32 subbands at UHFband. These templates are then used to measure frequencyresolved TOAs by crosscorrelating 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 Lband receivers. 
Information on the MeerKAT observations and data set for PSR J0737−3039A.
2.4. DM variation
The wideband observation and high precision of MeerKAT telescope make it possible to obtain an accurate DM measurement on a perepoch 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 4min 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^{6}, and this will be revisited by Hu, Porayko et al. (in prep.). A more thorough investigation of this effect, as well as the frequencydependent 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 edgeon 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 leadingorder this is the wellknown Shapiro delay (Shapiro 1964), which is expressed in the following form for binary pulsars (Blandford & Teukolsky 1976; Damour & Deruelle 1986):
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 DamourDeruelle (DD) timing model (Damour & Deruelle 1986) that can be fitted in the pulsar timing software TEMPO or TEMPO2 (Hobbs et al. 2006). The two postKeplerian (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 , where c is the speed of light in vacuum and cm^{3} s^{−2} 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 , where G is the gravitational constant.
The leadingorder 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)): with
where the semimajor axis of the relative orbit a_{R} = (x + x_{B})/s, with x and x_{B}^{7} being the projected semimajor 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
where the retardation correction can be taken directly from Kopeikin & Schäfer (1999, Eq. (130)) as
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 theoryindependent 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 lowkick 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
The first term on the righthand 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
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 in Eq. (6) is the higherorder correction originating from the ‘longitudinal deflection delay’, and can be written as (Doroshenko & Kopeikin 1995)
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:
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 xaxis. 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 () 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
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 dropout 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 16yr 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^{8}. 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. Postfit 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 Lband data are plotted in blue, whereas the UHFband 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} ≈ 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.
Timing parameters for PSR J0737−3039A using the TEMPO DDS binary model.
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)^{9} 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^{10} (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^{11} is used.
All binary parameters in Table 2 are consistent with the 16yr 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 χ^{2} of the individual observation is close to one, suggesting that our result is not affected by jitter noise. We also perform simulations and singlepulse 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 and the Shapiro shape parameter s.
For the advance of periastron, in addition to the 1PN contribution, we also need to account for higherorder corrections due to 2PN effects and LenseThirring (LT) precession caused by spinorbit coupling of pulsar A, as they are much larger than the measurement error of (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)
The first and second postNewtonian (PN) terms and are functions of masses and observed Keplerian parameters. The situation is more complicated for the LT contribution , 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 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^{45} g cm^{2} 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 from Eq. (35) in Kramer et al. (2021a), which uses the constraints on the EOS from Dietrich et al. (2020):
The Shapiro shape parameter s is the sine of the orbital inclination i. In Newtonian gravity, the orbital inclination is linked to the projected semimajor axis x via the binary mass function (e.g. Lorimer & Kramer 2004):
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):
Taking the measurements of P_{b}, x, and masses from Table 2, one can calculate that the 1PN correction is approximately 1.27 × 10^{−5}. 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 and s, one obtains the (Dopplershifted) 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 , NLO gravitational wave damping and massloss 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 and the mass loss contribution to Ṗ_{b}. For periastron advance , 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^{45} g cm^{2} 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.
Mass measurements with a new modified DDGR model that accounts for NLO contributions in the orbital motion and signal propagation in this system.
Following Kramer et al. (2021a), one could test the agreement of r in GR by comparing (cf. Table 2) with the companion mass determined here, which gives
This leads to a 5.3 × 10^{−3} (95% confidence) test of GR.
5. Studying NLO signal propagation effects
Because the double pulsar system is nearly edgeon 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 leadingorder 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,
Fig. 5. Aggregated residuals (blue) due to NLO contributions in Shapiro and aberration delay, shown in the orbital phase ψ. Residuals are rescaled by (1 + e_{T} cos θ)^{−1} to account for secular variations in amplitude due to the precession of periastron. The black curve indicates the fitted q_{NLO} (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 , and corrections related to the signal deflection in the gravitational field of the companion and . 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
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 colatitude 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 signaltonoise 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 subintegration. 
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 Lband data and in Fig. 9 for the UHFband 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 noneclipse 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 Lband 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. 
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
where χ_{0} 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 χ_{0} = 45° ^{12} (i.e. a relatively large latitudinal deflection delay) and scale it with a factor q^{latdef}. We simulate highprecision TOAs using this test model and fit for the scaling factor. The prefit 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 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 1ns 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: prefit 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: postfit 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 q^{len} 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 , where the quantity α_{E} is the angle corresponding to the Einstein radius, and is given by
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 plasmafilled 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′}^{13}. Consequently, one obtains
Since ΔΘ ≤ α_{E} ≪ 1 for all angles of Θ, we can expand the equation above in ΔΘ while keeping terms only up to order . This leads to a quadratic equation:
which has, under the assumptions made, one solution:
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 retardationcorrected orbital phase via cosΘ = sin isin(ψ + δψ^{ret}). The longitudinal and latitudinal deflection delay are given by
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 . 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 3yr timing observations of the double pulsar using the MeerKAT telescope, with a specific focus on studying higherorder signal propagation effects in the gravitational field of the companion. In order to minimise the effects from profile evolution and DM variation, we used frequencydependent 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^{−3} (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 , 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 subbanded 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 Lband data from MeerKAT (see Table 4). This makes their predictions conservative and we are likely to achieve even better measurements with future observations.
Comparison of the MeerKAT timing precision σ_{rms} assumed in Hu et al. (2020) and from real observations with Lband and UHFband receivers, scaled to a 5 min integration time over the full bandwidth.
x_{B} had been observed in Kramer et al. (2006a).
It should be noted that χ_{0} is not a constant, but our purpose here is to get a feeling for the measurability of the effect in timing rather than a proper account of the effect, which requires knowledge of the latitudinal variation in the emission pattern of the pulsar. Furthermore, as discussed at the beginning of Sect. 5.1, the beam geometry adopted by Rafikov & Lai (2006b) is not the correct one anyway.
We note that the D_{d} in Eq. (24) of Wucknitz (2008) is a typo and should not be there since it is already part of the definition of m.
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 MaxPlanck Society as part of the ‘LEGACY’ collaboration with the Chinese Academy of Sciences on lowfrequency 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 PrinSKA/CTA approved with the Presidential Decree 70/2016.
References
 Bailes, M., Jameson, A., Abbate, F., et al. 2020, PASA, 37, e028 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R., & Teukolsky, S. A. 1976, ApJ, 205, 580 [Google Scholar]
 Breton, R. P., Kaspi, V. M., Kramer, M., et al. 2008, Science, 321, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Burgay, M., D’Amico, N., Possenti, A., et al. 2003, Nature, 426, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T. 2007, arXiv eprints [arXiv:0704.0749] [Google Scholar]
 Damour, T., & Deruelle, N. 1985, Ann. Inst. Henri Poincaré Phys. Théor., 43, 107 [Google Scholar]
 Damour, T., & Deruelle, N. 1986, Ann. Inst. Henri Poincaré Phys. Théor., 44, 263 [Google Scholar]
 Damour, T., & Schäfer, G. 1988, Nuovo Cimento B Serie, 101B, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840 [NASA ADS] [CrossRef] [Google Scholar]
 Demorest, P., Ramachandran, R., Backer, D. C., et al. 2004, ApJ, 615, L137 [NASA ADS] [CrossRef] [Google Scholar]
 Dietrich, T., Coughlin, M. W., Pang, P. T. H., et al. 2020, Science, 370, 1450 [Google Scholar]
 Doroshenko, O. V., & Kopeikin, S. M. 1995, MNRAS, 274, 1029 [NASA ADS] [Google Scholar]
 Ferdman, R. D., Stairs, I. H., Kramer, M., et al. 2008, AIP Conf. Ser., 983, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Ferdman, R. D., Stairs, I. H., Kramer, M., et al. 2013, ApJ, 767, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Guillemot, L., Kramer, M., Johnson, T. J., et al. 2013, ApJ, 768, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655 [Google Scholar]
 Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302 [Google Scholar]
 Hu, H., Kramer, M., Wex, N., Champion, D. J., & Kehl, M. S. 2020, MNRAS, 497, 3118 [NASA ADS] [CrossRef] [Google Scholar]
 Klioner, S. A., & Kopeikin, S. M. 1992, AJ, 104, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Klioner, S. A., & Zschocke, S. 2010, Class. Quant. Grav., 27, 075015 [NASA ADS] [CrossRef] [Google Scholar]
 Kopeikin, S. M., & Schäfer, G. 1999, Phys. Rev. D, 60, 124002 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2006a, Science, 314, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2006b, Ann. Phys., 15, 34 [CrossRef] [Google Scholar]
 Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2021a, Phys. Rev. X, 11, 041050 [NASA ADS] [Google Scholar]
 Kramer, M., Stairs, I. H., Venkatraman Krishnan, V., et al. 2021b, MNRAS, 504, 2094 [CrossRef] [Google Scholar]
 Lai, D., & Rafikov, R. R. 2005, ApJ, 621, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Lazarus, P., Karuppusamy, R., Graikou, E., et al. 2016, MNRAS, 458, 868 [Google Scholar]
 Lorimer, D. R., & Kramer, M. 2004, Handbook of Pulsar Astronomy (Cambridge: Cambridge University Press), 4 [Google Scholar]
 Lyne, A. G., Burgay, M., Kramer, M., et al. 2004, Science, 303, 1153 [CrossRef] [PubMed] [Google Scholar]
 Nice, D., Demorest, P., Stairs, I., et al. 2015, Astrophysics Source Code Library [record ascl:1509.002] [Google Scholar]
 Parthasarathy, A., Bailes, M., Shannon, R. M., et al. 2021, MNRAS, 502, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Perlick, V. 2004, Liv. Rev. Relativ., 7, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Piran, T., & Shaviv, N. J. 2004, arXiv eprints [arXiv:astroph/0401553] [Google Scholar]
 Pol, N., McLaughlin, M., Kramer, M., et al. 2018, AJ, 853, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Prša, A., Harmanec, P., Torres, G., et al. 2016, AJ, 152, 41 [Google Scholar]
 Rafikov, R. R., & Lai, D. 2006a, Phys. Rev. D, 73, 063003 [Google Scholar]
 Rafikov, R. R., & Lai, D. 2006b, ApJ, 641, 438 [NASA ADS] [CrossRef] [Google Scholar]
 Ransom, S. M., Backer, D. C., Demorest, P., et al. 2004, arXiv eprints [arXiv:astroph/0406321] [Google Scholar]
 Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses (New York: SpringerVerlag) [Google Scholar]
 Serylak, M., Johnston, S., Kramer, M., et al. 2021, MNRAS, 505, 4483 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789 [Google Scholar]
 Shapiro, I. I. 1967, Science, 157, 806 [NASA ADS] [CrossRef] [Google Scholar]
 Smarr, L. L., & Blandford, R. 1976, ApJ, 207, 574 [Google Scholar]
 Soffel, M., Klioner, S. A., Petit, G., et al. 2003, AJ, 126, 2687 [Google Scholar]
 Spiewak, R., Bailes, M., Miles, M. T., et al. 2022, PASA, 39, e027 [NASA ADS] [CrossRef] [Google Scholar]
 Stairs, I. H., Thorsett, S. E., Dewey, R. J., Kramer, M., & McPhee, C. A. 2006, MNRAS, 373, L50 [CrossRef] [Google Scholar]
 Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170 [Google Scholar]
 Taylor, J. H., & Weisberg, J. M. 1989, ApJ, 345, 434 [NASA ADS] [CrossRef] [Google Scholar]
 Ward, W. R. 1970, ApJ, 162, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Wex, N., & Kramer, M. 2020, Universe, 6, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Willems, B., & Kalogera, V. 2004, ApJ, 603, L101 [NASA ADS] [CrossRef] [Google Scholar]
 Willems, B., Kaplan, J., Fragos, T., Kalogera, V., & Belczynski, K. 2006, Phys. Rev. D, 74, 043003 [NASA ADS] [CrossRef] [Google Scholar]
 Wucknitz, O. 2008, MNRAS, 386, 230 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Mass measurements with a new modified DDGR model that accounts for NLO contributions in the orbital motion and signal propagation in this system.
Comparison of the MeerKAT timing precision σ_{rms} assumed in Hu et al. (2020) and from real observations with Lband and UHFband receivers, scaled to a 5 min integration time over the full bandwidth.
All Figures
Fig. 1. Pulse profile of PSR J0737−3039A observed at multiple frequencies with the MeerKAT UHF and Lband 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 xaxis. 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 () 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. Postfit 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 Lband data are plotted in blue, whereas the UHFband 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 rescaled by (1 + e_{T} cos θ)^{−1} to account for secular variations in amplitude due to the precession of periastron. The black curve indicates the fitted q_{NLO} (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 subintegration. 

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 Lband 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 UHFband data summed from 37 eclipses. 

In the text 
Fig. 10. Lensing simulation: prefit 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: postfit residual plotted against orbital phase ψ. 

In the text 
Fig. 12. Uncertainty of factor q^{len} 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 plasmafilled magnetosphere of B (cf. Lai & Rafikov 2005; Rafikov & Lai 2006b). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.