© The Authors 2024

1. Introduction

PSR J0737−3039A/B is a highly relativistic double neutron star binary with a short, 2.45 h mildly eccentric (e = 0.088) orbit (Burgay et al. 2003). Uniquely, both neutron stars have been detected as radio pulsars (hereafter referred to as pulsars A and B) with respective spin periods of 22.7 ms and 2.77 s (Lyne et al. 2004). As a result, the system was dubbed the ‘double pulsar’. High-precision timing of the two pulsars resulted in four independent tests of Einstein’s general theory of relativity (GR) in the strong-field regime within only 2.7 yr of its initial discovery (Kramer et al. 2006). Continued timing of the faster rotating pulsar A enabled a preliminary measurement of Lense-Thirring precession in the system, which was used to place limits on the pulsar’s moment of inertia, along with several higher-order relativistic effects (Kramer et al. 2021a; Hu et al. 2022). These detections have made the double pulsar one of the most successful astrophysical laboratories for testing our theories of gravity and the behaviour of matter at super-nuclear densities.

A substantial contributor to our ability to test fundamental physics with this system comes from its remarkably edge-on orientation of the orbital plane, which is inclined at an angle of i = 89.35° ±0.05° (or equivalently, 90.65° ±0.05°) from our perspective on Earth (Kramer et al. 2021a). This chance geometric alignment results in a 30–40 s long eclipse of pulsar A by the magnetosphere of pulsar B around superior conjunction, offering a unique means to directly probe the plasma environment around an active pulsar (Lyne et al. 2004). The duration of these eclipses corresponds to a region of space that is ∼1.7 × 10⁷ m wide, spanning only ∼10% of the light-cylinder radius of an equivalent isolated pulsar with the same rotational properties as pulsar B (Lyne et al. 2004; Kaspi et al. 2004; Breton et al. 2012). This smaller than expected eclipsing region arises from the relativistic wind from pulsar A penetrating deep into the magnetosphere of pulsar B, compressing the ‘windward’ side facing pulsar A and blowing the ‘leeward’ side backwards into a cometary magnetotail (Arons et al. 2005).

High-time resolution observations of the eclipses analysed by McLaughlin et al. (2004) with the 100-m Green Bank Telescope (GBT) revealed the light curve of pulsar A exhibits peaks and troughs in its flux density that are spaced apart by once and half the 2.77 s rotation period of pulsar B at different eclipse phases. This modulation can be almost entirely explained through a simple geometric model, in which the radio pulses undergo synchrotron absorption by the relativistic pair plasma that is confined to the toroidal, closed-field region of pulsar B’s magnetosphere (Lyutikov & Thompson 2005). The success of the model provided not only the first direct evidence for a dipole magnetic field geometry around a pulsar, but was later used to track temporal changes seen in a set of eclipses collected with the GBT over 3.9 yr, resulting in a method of detecting the geodetic precession rate of pulsar B (${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$) and an associated fifth independent test of GR (Breton et al. 2008). While the effects of geodetic precession have been detected in six other relativistic binary pulsars to date (Kramer 1998; Kirsten et al. 2014; Fonseca et al. 2014; Venkatraman Krishnan et al. 2019; van Leeuwen et al. 2015; Cameron et al. 2023), the corresponding precession rate measurements are largely indirect and primarily based on modelling the pulse profile width and polarisation properties. Interpretation of these observables can be heavily reliant on the assumed pulsar beam shape, system geometry, and the applicability of the rotating vector model (Radhakrishnan & Cooke 1969), all of which can display large deviations from simple models. Breton et al. (2008) demonstrated the precession of pulsar B has a significant impact on the observed eclipse light curve, which was fitted by adding a simple linear drift in the spin-axis longitude of pulsar B over time when computing the model templates.

In addition to enabling tests of spin-orbit coupling in the strong-field regime, the eclipses allow us to develop an improved picture of the overall geometry of the double pulsar. Such measurements provide an important input into the construction of stellar binary population synthesis models (e.g. Riley et al. 2022), improving our understanding of the formation and evolutionary history of such systems (Stairs et al. 2006; Kim et al. 2015; Tauris et al. 2017; Vigna-Gómez et al. 2018), and generating astrophysically motivated priors for analysing gravitational waves from double neutron star mergers (Zhu & Ashton 2020). Within months of pulsar B being discovered, it became obvious that the overall shape, intensity and position of the profile within the brief ‘orbital visibility windows’ were evolving with time (Burgay et al. 2005). This behaviour is linked to our changing line of sight through the emission cone as the pulsar undergoes geodetic precession. While this ultimately resulted in the disappearance of radio pulses from pulsar B sometime in late-2008 (Perera et al. 2010), it did open the rare opportunity to develop a map of the radio beam. Similar beam maps have only been obtained for three other pulsars in relativistic binaries to date (B1913+16, Weisberg & Taylor 2002; J1141−6145, Manchester et al. 2010; and J1906+0746, Desvignes et al. 2019). Accurate geometric constraints are required when creating such radio-beam maps and predicting when radio pulses from pulsar B may once again be detectable (Breton 2009; Perera et al. 2012; Noutsos et al. 2020).

In this work, we present the first results of an eclipse monitoring campaign that utilises the low system temperature and high gain of MeerKAT to capture the subtle time-evolution of the eclipse modulation pattern in the greatest detail yet. Section 2 describes the observations and data processing steps. In Sect. 3 we outline our improved approach to modelling the eclipse light curves, and the use of hierarchical inference methods to infer the geodetic precession rate of pulsar B. The results from our eclipse modelling are detailed in Sect. 4, including an outline of the challenges faced given the current system geometry and model compatibility. In Sect. 5 we discuss updates to the tests of GR and constraints on the system geometry that stem from our improved measurement of geodetic precession, implications for the radio beam shape of pulsar B and when we may again detect radio pulses from it. We summarise the results in Sect. 6 and outline potential future directions for modelling of the eclipses.

2. Observations

We have performed monthly monitoring observations of the double pulsar from July 2019 to September 2022 with the Meer Karoo Array Telescope (MeerKAT), a 64-element radio interferometer located in the Northern Cape province of South Africa (Jonas & the MeerKAT Team 2016). Our observations were collected under the Relativistic Binary theme of the MeerTime large science project (Bailes et al. 2020; Kramer et al. 2021b). Observations at MeerKAT from March 2019 to March 2020 were performed using the 1284 MHz central frequency L-band receiver system, after which the majority of observations were performed with the 816 MHz central frequency UHF receivers. These data were collected using the PTUSE instrument (Bailes et al. 2020) and are coherently dedispersed to account for the frequency-dependent delay induced by the passage of pulsar A’s radiation through the interstellar medium. PTUSE provides 1024 frequency channel filterbank data across the 856 MHz and 544 MHz bandwidths of the L-band and UHF receiver fleets respectively, along with ∼9 μs time sampling, and full Stokes information (Bailes et al. 2020).

The coherent search-mode data were folded at the predicted rotation period of pulsar A using the DSPSR software package (van Straten & Bailes 2011), where individual single pulses were saved to PSRCHIVE format archive files (Hotan et al. 2004; van Straten et al. 2012). Frequency channels that were affected by radio-frequency interference were excised using the MEERGUARD¹ pulsar data cleaning package. We then averaged the data in frequency and binned in time by four rotations of pulsar A for an effective time resolution of ∼91 ms. Only the total intensity data were analysed in this work. Analysis of the polarimetry of pulsar A throughout the eclipses is the subject of a separate work. We also saved a copy of the data set where frequency-channels that were outside a 961–1088 MHz subband that is covered by both the UHF and L-band receivers were excised. This was done to test whether or not frequency dependencies in the eclipse width and depth affect our ability to infer secular changes in the modulation pattern. The two data sets are distinguished by the FULLBAND and SUBBAND designations hereafter.

3. Methods

3.1. Eclipse light-curve extraction

To model the total intensity light curve, we first extracted the flux of pulsar A via a matched-filtering process. This was performed using the PSRFLUX tool in PSRCHIVE (Hotan et al. 2004; van Straten et al. 2012), where a high signal-to-noise template (generated from the integrated pulse profile from many hours of observations) was cross-correlated with a frequency-averaged copy of the total intensity data. The resulting flux densities were normalised by the median off-eclipse value so the points where pulsar A is unobstructed, have values of order unity. We used the PAT tool in PSRCHIVE to compute topocentric arrival times at MeerKAT for pulsar A, which were converted to equivalent arrival times at the Solar System barycentre using TEMPO and the JPL-DE436² Solar System ephemeris. The barycentred arrival times were then converted to orbital phases using the latest pulsar A timing ephemeris from Hu et al. (2022)³. As in Kramer et al. (2021a), we have taken into account retardation effects due to the motion of pulsar B as the signal from A propagates through the binary system towards Earth. From here the total-intensity data were then fit using a variation of the recipe outlined in Breton et al. (2008) in order to infer the geometry of pulsar B.

In Fig. 1 we show the complete sample of eclipse light curves corresponding to the FULLBAND version of the data. The eclipses are largely dominated by transparency windows, regions where pulsar A is not obscured, that are separated by a full rotation of pulsar B. Windows separated by only a half-rotation are restricted to the ingress and egress phases. This is in contrast to the data analysed by Breton et al. (2008), where the precession phase of pulsar B meant that the transition from the half to full rotational separation occurred closer to superior conjunction.

Fig. 1. Normalised intensity of pulsar A throughout each of the eclipses recorded by MeerKAT. The red vertical line indicates the orbital phase of superior conjunction. Regions contained within the white lines correspond to where the light curve model best matches the data (see Sect. 4.1). Cyan and orange indicators on the right-hand side denote light-curves obtained using the L-band and UHF receivers respectively. We note that the observation numbers on the y-axis do not increase linearly with time.

3.2. Light-curve model

Here we provide a summary of the light-curve model developed by Lyutikov & Thompson (2005). Pulsar B is positioned at the centre of a Cartesian coordinate system where the x-axis points towards our line of sight to the system, the z-axis is in the plane of the sky as viewed from Earth, and the y-axis runs parallel to the apparent motion of pulsar A (Breton et al. 2008). The motion of pulsar A is offset vertically from the origin by a constant value of z = z₀. The direction of Ω_B can be fully described by two angles: the spin-axis co-latitude (δ_B) with respect to the z-axis, and spin-axis longitude (φ_so) with respect to the x-y-plane (Damour & Taylor 1992). Ordinarily, δ_B and φ_so would be related to the Cartesian coordinate system via

$$\begin{array}{cc}& cos\theta =cos({90}^{°}-i)cos({180}^{°}-{\delta }_{\mathrm{B}})\hfill \\ \hfill & -sin({90}^{°}-i)sin({180}^{°}-{\delta }_{\mathrm{B}})cos{\phi }_{\mathrm{so}},\hfill \end{array}$$(1)

$$\begin{array}{cc}& sin\phi =-\frac{sin({180}^{°}-{\delta }_{\mathrm{B}})sin{\phi }_{\mathrm{so}}}{sin\theta }.\hfill \end{array}$$(2)

However, the double pulsar is nearly edge-on from our perspective (so either i = 89.35° ±0.05° or 90.65° ±0.05°; Kramer et al. 2021a), meaning the z-axis is effectively (anti-)aligned with the total angular momentum vector. As a result, we can make the simplifying assumption that θ ≈ 180° −δ and φ ≈ −φ_so. Note that in the coordinate frame of Breton et al. (2008), where A moves in the positive y direction during the eclipse, the orbital angular momentum is (almost) aligned with the negative z direction, that is, pointing downwards in Fig. 2.

Fig. 2. Diagram depicting the double pulsar geometry according to the eclipse light-curve model, adapted from Fig. 1 of Breton et al. (2008). The minimum radial extent of the synchrotron absorbing plasma (the cooling radius) is represented by Rmin. Longitude and latitude of the spin axis of pulsar B (ΩB) is given by the angles φ and θ respectively. The magnetic axis (μB) is offset from ΩB by angle α. In our model we have neglected the deflection of A’s radio signals in the gravitational field of B, which amounts to a change in the impact parameter at B of no more than about 600 km (Kramer et al. 2021a).

The total intensity models of pulsar A throughout the eclipses are computed from the synchrotron optical depth (τ) of the plasma trapped within the closed-field region of pulsar B’s magnetosphere. Following the prescription of Breton et al. (2008), the optical depth at a given point during the eclipse is computed as

$$\begin{array}{c}\hfill \tau =\frac{\mu }{{\nu }_{\mathrm{GHz}}^{-\ell }}{\displaystyle {\int }_{-R_{\mathrm{mag}}}^{+R_{\mathrm{mag}}}(\frac{Bsin\kappa }{B_{\mathrm{mag}}})\mathrm{d}(\frac{x}{R_{\mathrm{mag}}}).}\end{array}$$(3)

The optical depth depends upon radio frequency (ν) with a power-law index of ℓ, R_mag is the truncation radius of B’s magnetosphere, B is the local magnetic field strength along the line of sight in units of B_mag (magnetic-field strength at R_mag), κ is the angle between the local magnetic field and our line of sight and x is the radial position of pulsar A from our perspective in units of R_mag. Note that R_mag is not fit for directly but is inferred via the parameter ξ that scales the size of the magnetosphere to the orbital distance between pulsars A and B (Breton 2009). The scaling parameter μ combines various parameters that describe the physical properties of B’s magnetosphere as

$$\begin{array}{c}\hfill \mu =\frac{4.5\times {10}^{-6}{\lambda }_{\mathrm{mag}}}{N_{\mathrm{B}}^{1/4}}\frac{k_{\mathrm{B}}{T}_{\mathrm{e}}}{10m_{\mathrm{e}}{c}^2},\end{array}$$(4)

in which λ_mag is the pair plasma electron multiplicity, N_B alters the size of the magnetosphere based on the impact of the wind from pulsar A, k_B is Boltzmann’s constant, T_e and m_e the electron temperature in the plasma and the electron mass, and c is the vacuum speed of light. Variations in each of these parameters, which are not practicable to be fit for individually, serve to alter the depth of the eclipses as the intensity of pulsar A is computed from the optical depth as e^−τ. In Lyutikov & Thompson (2005) the frequency dependence of the optical depth is set at ℓ = 5/3. However, observations by Breton et al. (2012) showed the observed eclipse depth follows a much shallower frequency dependence than expected from the standard model. We inferred a value of ℓ ∼ 1/3 from a multi-band fit to 16 subbands across a pair of eclipses detected with the UHF and L-band receivers. We adopt this value throughout the remainder of this work.

The modulation pattern of the light curve depends strongly on the changing line-of-sight geometry of pulsar B as it rotates, which is modelled through the corresponding variations in both B and κ in Eq. (3). Both of these terms are related to the dipole unit vector magnetic polar angle (θ_μ) as

$$\begin{array}{c}\hfill cos{\theta }_{\mu }=\frac{\widehat{\mu }·\mathbf{r}}{r},\end{array}$$(5)

where r = {x, y(t),z}, r is the distance between pulsars A and B in spherical coordinates ($r=|\mathbf{r}|=\sqrt{x^2+{\mathit{y}}^2(t)+z^2}$) and $\widehat{\mu }$ is the dipole unit vector, the components of which are given by

$$\begin{array}{cc}& {\widehat{\mu }}_x=({\widehat{\mu }}_x^{\mathrm{\Omega }}cos\theta +{\widehat{\mu }}_z^{\mathrm{\Omega }}sin\theta )cos\phi -{\widehat{\mu }}_y^{\mathrm{\Omega }}sin\phi ,\hfill \\ \hfill & {\widehat{\mu }}_y=({\widehat{\mu }}_x^{\mathrm{\Omega }}cos\theta +{\widehat{\mu }}_z^{\mathrm{\Omega }}sin\theta )sin\phi +{\widehat{\mu }}_y^{\mathrm{\Omega }}cos\phi ,\hfill \\ & {\widehat{\mu }}_z={\widehat{\mu }}_z^{\mathrm{\Omega }}cos\theta -{\widehat{\mu }}_x^{\mathrm{\Omega }}sin\theta ,\hfill \end{array}$$(6)

with

$$\begin{array}{cc}& {\widehat{\mu }}_x^{\mathrm{\Omega }}=sin\alpha cos({\varphi }_{\mathrm{B}}+\mathrm{\Delta }{\varphi }_{\mathrm{B}}),\hfill \\ \hfill & {\widehat{\mu }}_y^{\mathrm{\Omega }}=sin\alpha sin({\varphi }_{\mathrm{B}}+\mathrm{\Delta }{\varphi }_{\mathrm{B}}),\hfill \\ & {\widehat{\mu }}_z^{\mathrm{\Omega }}=cos\alpha .\hfill \end{array}$$(7)

Here, α is the magnetic inclination angle of pulsar B and ϕ_B is the rotation phase of pulsar B, which is related to the spin-vector direction as ϕ_B = Ω_Bt = 2πt/P_B where P_B is the spin period of pulsar B. The parameter Δϕ_B accounts for the offset from ϕ_B = 0 at an assumed reference time. The values of B and κ are computed at each step in r as

$$\begin{array}{c}\hfill B=\frac{\sqrt{1+3{cos}^2{\theta }_{\mu }}}{r^3}{\mu }_{\mathrm{B}},\end{array}$$(8)

and

$$\begin{array}{c}\hfill cos\kappa =\frac{3cos{\theta }_{\mu }(x/r)-{\widehat{\mu }}_x}{1+3{cos}^2{\theta }_{\mu }},\end{array}$$(9)

before being passed to Eq. (3) where they are integrated over x.

3.3. Joint-fitting of eclipse pairs

Unlike in Breton et al. (2008), radio emission is not currently detected from pulsar B. Without an a priori timing solution for pulsar B, we cannot predict the rotation-phase of pulsar B throughout our observations. However, the S/N of the eclipses detected by MeerKAT is sufficiently high that we can instead perform direct fitting of individual eclipses, where the period of pulsar B is held fixed at P = 2.773 s and the rotation phase at the start of an eclipse timeseries included as a free parameter.

Our initial attempts at fitting the eclipses independently of one another were hampered by ‘eclipse weather’, epoch-to-epoch stochasticity in the eclipse envelope likely originating from fluctuations in the plasma content and radial extent of the closed magnetic field lines of pulsar B. This phenomena resulted in excess scatter in the individual geometric constraints inferred from one eclipse to the next. We were able to mitigate a substantial amount of this behaviour by performing joint fits to pairs of eclipses that were separated in time by only a single orbit. This was conducted using a Gaussian likelihood function of the form

$$\begin{array}{c}\hfill \mathcal{L}(\mathbf{d}|\mathrm{\Theta },{\sigma }_Q)={\displaystyle \prod _{i=1}^2\prod _{j=1}^N\frac{1}{\sqrt{2\pi {\widehat{\sigma }}_{i,j}^2}}exp[-\frac{{(d_{i,j}-{\mu }_{i,j}(\mathrm{\Theta }))}^2}{2{\widehat{\sigma }}_{i,j}}],}\end{array}$$(10)

where d represents the input eclipse light curves, μ is the eclipse model, Θ contains the model parameters, and ${\widehat{\sigma }}_{i,j}^2={\sigma }_{i,j}^2+{\sigma }_{Q,i}^2$ are the uncertainties on the light-curve fluxes added in quadrature with an additional error in quadrature parameter (EQUAD; σ_Q, i). This extra uncertainty parameter accounts for both pulse-to-pulse flux variations of pulsar A and unaccounted systematic errors. Our priors for the model parameters are summarised in Table 1. Posterior samples for the model parameters were generated using BILBY (Ashton et al. 2019) as a front-end to the PYMULTINEST sampler, a PYTHON-based implementation of the MULTINEST nested-sampling algorithm (Buchner et al. 2014; Feroz et al. 2009).

From our initial test fits to the data, we found the posterior distributions for the cooling radius consistently displayed equal a posteriori support for values of R_min ≲ 0.5 before running up against the edge of the prior. While this does not give us a precise radius at which the plasma becomes transparent to radio waves, we can at least constrain it to less than half the radius of pulsar B’s truncated magnetosphere. Given the template light-curves are identical for R_min ≲ 0.5, we opted to fix the value of R_min = 0.5 for the remainder of our light-curve modelling.

3.4. Measuring the spin-precession rate of pulsar B

The spin vector of a rotating body moving in the curved spacetime of a companion precesses about the total angular-momentum vector of the binary system. In the case of the double pulsar, the total angular momentum vector is almost exactly aligned with the orbit-normal vector (z-axis in Fig. 2). As a result, the geodetic precession of pulsar B manifests as a time-varying change in the spin-axis longitude while the spin-axis latitude remains unchanged. Hence, the time-evolution of these two parameters can be written as

$$\begin{array}{c}\hfill \theta ={\theta }_0,\end{array}$$(11)

and

$$\begin{array}{c}\hfill \phi (t)={\phi }_0-{\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}(t-t_0),\end{array}$$(12)

where θ₀ and φ₀ are the co-latitude and longitude of B’s spin axis at some reference epoch, t₀= MJD 59289.

As a phase-connected timing solution for pulsar B across our dataset is unavailable, we could not use the exact same simultaneous joint-fit method utilised by Breton et al. (2008) to infer the precession rate of pulsar B. Instead, we explored two alternative means for measuring the change in spin-axis longitude with time. As with the joint-eclipse fits, both approaches made use of the PYMULTINEST sampler with BILBY.

3.4.1. Hierarchical inference

Our first method involved a second-stage hierarchical fit to the posterior samples we obtained for φ at each observing epoch. We used a hyper-likelihood function of the form

$$\begin{array}{cc}& {\mathcal{L}}_{\mathrm{tot}}({\phi }_0,{\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}},{\sigma }_{\phi }|\phi (t))=\hfill \\ & {\displaystyle \prod _{i=1}^{N_{\mathrm{e}}}\frac{1}{n_i}\sum _k^{n_i}exp[-\frac{{(\phi {(t_i)}_k-{\phi }_0+{\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}(t_i-t_0))}^2}{2{\sigma }_{\phi }^2}],}\hfill \end{array}$$(13)

where N_e is the number of epochs for which we have measured φ, n_i is the total number of posterior samples from the i-th eclipse, and σ_φ is a normalising factor that accounts for the variance in our measurements of φ. Uniform priors were assumed for each of the hyper-parameters.

3.4.2. Iterative inference

Our second approach to measuring the precession rate involved performing a global ‘iterative’ fit to every eclipse light-curve simultaneously. Instead of using a timing model for pulsar B to predict its rotation phase throughout our observations, as was done by Breton et al. (2008), we used the median value for the rotation phase obtained from the initial eclipse-pair fits as our input for the phase of pulsar B at the beginning of each eclipse. This allowed us to avoid fitting for the rotation phase of the pulsar, thereby reducing the number of model parameters by N_obs. As with the eclipse-pair fits, we used a joint Gaussian likelihood function of the form

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{tot}}(\mathbf{d}|\mathrm{\Theta }(t))={\displaystyle \prod _{i=1}^{N_{\mathrm{obs}}}}\prod _{j=1}^{N_i}\frac{1}{\sqrt{2\pi {\widehat{\sigma }}_{i,j}^2}}exp[-\frac{{(d_{i,j}-{\mu }_{i,j}(\mathrm{\Theta }(t)))}^2}{2{\widehat{\sigma }}_{i,j}}],\end{array}$$(14)

where the time dependence of Θ comes from computing φ at each eclipse epoch via Eq. (12). Aside from the rotation phase of pulsar B and the EQUAD used for each eclipse (which were also set to the median values from the eclipse-pair fits), we assumed uniform priors on all of the geometric parameters. This technique allows us to directly measure for the values of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ and φ₀ when fitting the eclipses.

4. Results and analysis

4.1. Geodetic precession and model limitations

Our joint modelling of eclipse pairs were able to reproduce most of the observed peaks and troughs detected by MeerKAT at both L-band and UHF frequencies (see Fig. 3). However, the eclipse model often failed to fully capture the morphology of the ingress and egress phases, where the effects of small variations in magnetosphere size or plasma density have the largest impact on the observed light curve. Indeed, the median a posteriori models displayed in Fig. 3 show clear deviations from the data in both the ingress and egress phases. This mismatch between the observed eclipse edges and the model was noted by Lyutikov & Thompson (2005), who found the best match occurred towards the centre of the eclipses. The same effect was also seen by Breton et al. (2008), which they compensated for by restricting their fits to only data between −1.0° / + 0.75° from superior conjunction. In our case, the effects of these local distortions are amplified by the increased sensitivity of MeerKAT over the GBT.

Fig. 3. Two eclipses of pulsar A detected with MeerKAT using the UHF (left) and L-band (right) receiver systems. These light curves were generated after averaging over the full receiver bandpasses. We have inflated the flux uncertainties by the median recovered EQUAD for each eclipse. The orange traces correspond to the median a posteriori light-curve models recovered from fitting the observations. Bottom panels show the residuals after subtracting the median model.

In Fig. 4, we show that the recovered φ₀, precession rate and σ_φ from hierarchical fits to both versions of the eclipse data vary depending on whether data are cut from the ingress and egress phases at different distances from superior conjunction. Assuming the validity of general relativity, the predicted precession rate of pulsar B in rad s⁻¹ can be computed as (Barker & O’Connell 1975)

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}=\frac{1}{2}(\frac{G}{c^3}{)}^{2/3}(\frac{P_b}{2\pi }{)}^{-5/3}\frac{m_{\mathrm{A}}(4m_{\mathrm{B}}+3m_{\mathrm{A}})}{(1-e^2){(m_{\mathrm{A}}+m_{\mathrm{B}})}^{4/3}},\end{array}$$(15)

Fig. 4. Comparison between hierarchical fits to different realisations of the light curves using the FULLBAND and SUBBAND datasets. Contours in the two-dimensional posteriors represent the 68% and 95% confidence regions. Each colour represents measurements from eclipses that had different amounts of the ingress or egress phases removed. The dashed lines correspond to the predicted values of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ from general relativity.

where G is the gravitational constant, c the vacuum speed of light, P_b and e are the orbital period eccentricity of the binary, m_A is the mass of pulsar A and m_B is the mass of pulsar B. Using the pulsar masses and timing measurements from Kramer et al. (2021b), we obtain a predicted precession rate of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}=5.074005°\pm 0.000003°$ yr⁻¹, which is shown as the dashed lines in Fig. 4. Comparing the different fits in Fig. 4, it is clear that using the same restricted range as Breton et al. (2008) had little to no impact on the recovered model parameters when compared to runs that utilised the full eclipse light curves. The vast majority of the restricted orbital phase ranges appear to result in measurements of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ that are biased towards larger values, appearing to be either marginally inconsistent with the GR prediction or consistent to within the 95–99.7% confidence intervals.

The cause of this behaviour is apparent when comparing the light-curve shapes detected by MeerKAT between 2019 and 2023 to that detected by the GBT from 2003 and 2008. From Figs. 1 and 3, the light-curves appear relatively symmetric with the ‘partial’ transparency window where some of the radiation from pulsar A leaks through the magnetosphere of pulsar B appearing around superior conjunction. Notably, the transition points where the spacing between the transparency windows goes from half to a full rotation of pulsar B are restricted to the early-ingress and late-egress phases in our MeerKAT data. The eclipses seen in the earlier GBT data (see Fig. 2 of McLaughlin et al. 2004 and Fig. 3 of Breton et al. 2008) are substantially more asymmetric, with the half-to-full rotation transition taking place approximately one-third of the way through the eclipses, while the partial transparency windows appeared just prior to egress. Accurately capturing the location of these transition points in our model fits is critical for recovering the correct spin-axis longitude over time. The current locations of these features within the edges of the eclipse can account for the difficulties in measuring the pulsar geometry, as these regions are subject to the aforementioned local distortions in the magnetosphere of pulsar B, resulting in mismatches between the light-curve model and the data. Both the ±0.70° and ±0.72° cuts result in hierarchical fits that are consistent with the GR-predicted value at the 68% confidence interval with the FULLBAND uncut dataset.

With these limitations in mind, for the remainder of this work we refer to results obtained using the FULLBAND version of the eclipse data with the ±0.72° cuts unless otherwise specified. In Fig. 5 we show the recovered geometric parameters for pulsar B from our fits to individual pairs of eclipses. As expected, α, μ, z₀ and ξ all remain at near constant values, while φ displays a significant decrease from 13.0° ±0.5° to $-{6.2}_{-{1.4}^{°}}^{°+{1.6}^{°}}$ as the spin axis of pulsar B precesses throughout the 3.5 yr data set. Notably, θ appears to show a slight downward trend which is not expected under our earlier assumption that θ = θ₀ in GR. However, given the current unfavourable geometry and mismatch of the model to the data, this may simply be a result of time-dependent degeneracy between the evolving values of φ and θ. We also over-plot the predicted precession rate of pulsar B in addition to that inferred from our hierarchical analysis of φ(t) in Fig. 5. Our recovered precession rate of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}={5.16°}_{-0.34°}^{+0.32°}$ yr⁻¹ agrees with the expected value of 5.074005° yr⁻¹ from GR to within ∼6.5%. For comparison, the corresponding precession rate obtained from the SUBBAND dataset with the same restricted fits to data within the ingress and egress returned ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}={5.31°}_{-0.41°}^{+0.42°}$ yr⁻¹.

Fig. 5. Evolution of our recovered geometric parameters for pulsar B over time measured using the FULLBAND (±0.72°) data set. The GR-predicted evolution of φ due to Geodetic precession is indicated by the black line. Dark orange shading indicates 68% range covered by the precession rate inferred from our hierarchical fit. Light orange shading includes scatter in the values of φ over time.

In principle, the iterative method from Sect. 3.4.2 should provide a higher precision measurement of the precession rate as we average over more of the stochastic eclipse weather. Performing this fit with the FULLBAND dataset with ±0.72° cuts, we obtain ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}={5.27°}_{-0.08°}^{+0.07°}$ yr⁻¹. While the 68% confidence intervals represent a ∼1.5% uncertainty on ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$, the marginalised posterior distribution is only consistent with the GR-predicted value to within the 99.7% confidence interval. This is insufficient to use as a precision test of gravity and instead likely reflects the current difficulties in measuring the time-evolution of the spin-axis longitude with the existing light-curve model.

4.2. Iterative eclipse fits

While the current alignment of pulsar B is not conducive to performing rigorous, high-precision tests of spin-obit coupling, we can still obtain an accurate measurement of the overall system geometry under the assumption that GR is the correct theory of gravity. Fixing the precession rate to the GR-predicted value of 5.074005° yr⁻¹, we used the iterative fitting approach with the FULLBAND data restricted to fitting the ±0.72° region of the light curves to obtain high-precision measurements of the remaining geometric and magnetospheric parameters. We present the recovered geometric and magnetospheric model parameters in Table 2. In Fig. 6 we present the one- and two-dimensional posterior distributions on the model parameters. For comparison, we also show the same model parameters where ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ was searched over as a free parameter.

Fig. 6. One- and two-dimensional posterior distributions of the eclipse light-curve parameters from iteratively fitting the FULLBAND data with the ±0.72° cuts where the precession rate was either fixed to the GR predicted value (green) or allowed to be fit as a free parameter (orange). Contours indicate the 68, 95 and 99.7% confidence intervals. Dashed black line is the predicted value of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$.

We find only marginal differences in the recovered model parameters when we do or do not assume the GR-predicted precession rate. Indeed the negligible log-Bayes factor of lnℬ = 0.05 in favour of the ‘Not Assuming GR’ model indicates the time-evolution of the eclipses is well described by pulsar B precessing at the GR-predicted rate. However, our geometric and magnetospheric constraints are substantially different to those previously reported by Breton et al. (2008) Notably, our α and θ differ by close to 10°, while φ₀ (referenced to MJD 59289) is significantly different to the value predicted by Breton et al. (2008) for the same date. These differences between eclipse modelling results likely originates from our framework fully sampling the model parameter space (except for R_min), rather than setting μ, z₀ and ξ to fixed values.

We also tested whether evolution of the light-curve due to precession of pulsar B can break various covariances that exist between model parameters (Breton 2009). Symmetries exist between Δϕ_B, α, θ, φ, and z₀, resulting in bimodal marginalised posterior distributions for the individual parameters and four modes in the two-dimensional posteriors with each other. These symmetries appear as Δϕ → Δϕ_B ± 0.5, {α, θ, φ}→180° −{α, θ, φ} and −z₀ → +z₀. Since φ varies with time due to pulsar B precessing, it is possible to break these symmetries by carrying out a simultaneous fit to several detected eclipses that are drawn from different points in our 3.5 yr campaign. We performed this joint fit to a sample of five eclipses using an expanded version of the approach detailed in Sect. 3.3 with wider prior ranges on the aforementioned parameters of interest. Our recovered posterior distributions for α, θ, φ, and z₀ were all well constrained and confined to a single mode peaking at the previously reported values in Table 2. The only symmetries that were not broken via this approach were those between Δϕ_B and α, where our recovered posterior distributions displayed two peaks separated by half a rotation of pulsar B for each Δϕ_B and α ∼ 60° or 120°. Both come down to our lack of a priori knowledge on the rotation phase of pulsar B and could be resolved through independent timing of the pulsar via its radio pulses.

5. Discussion

5.1. Testing theories of gravity

The geodetic precession rate can be used in combination with measurements of other relativistic effects to test GR and alternate theories of gravity. In Fig. 7, we have plotted the inferred masses of pulsars A and B from ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ alongside those from post-Keplerian parameters measured through precision timing of pulsar A (Kramer et al. 2021a) and the mass-ratio (Kramer et al. 2006) assuming the validity of GR. A deviation of one or more pairs of lines away from a common range of pulsar masses would indicate an inconsistency with GR, assuming the impact of higher-order effects to measurements of the post-Keplerian parameters are either negligible or have been corrected for. The pair of lines associated with the 68% confidence interval for ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ intersects the same common-point as the other post-Keplerian parameters, indicating that it is consistent with GR to our current measurement uncertainty. This provides one of the five independent tests that are possible with the measurements displayed in this diagram.

Fig. 7. Mass-mass diagram illustrating current tests of general relativity with the double pulsar. Our improved measurement of pulsar B’s spin-precession rate (${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$) is given by the (mostly) vertical blue lines where the separation between them indicates the 68% uncertainty. Other individual constraints on the masses of pulsar A and B from Kramer et al. (2006, 2021a) are shown by pairs of lines. These include the periastron precession rate ($\dot{\omega }$), shrinking of the orbital period due to gravitational radiation (Ṗb), gravitational redshift (γ), Shapiro delay range (r) and shape (s), and the mass-ratio (R). The shaded region is forbidden by the individual mass functions of the two pulsars (i.e. sin i ≤ 1).

In addition to the mass-mass comparison, the double pulsar is the only relativistic binary system in which direct constraints can be placed on the strong-field spin-orbit precession, as the measurement of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ can be combined with the orbital parameters inferred from timing both pulsars independently (Breton et al. 2008; Kramer & Wex 2009). Under the set of generic Lorentz-invariant relativistic theories introduced by Damour & Taylor (1992), the geodetic precession rate can be reformulated as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}=(\frac{P_b}{2\pi }{)}^{-3}\frac{x_{\mathrm{A}}{x}_{\mathrm{B}}}{s^2(1-e^2)}\frac{c^2{\sigma }_{\mathrm{so}}}{\mathcal{G}},\end{array}$$(16)

where x_A and x_B are the projected semi-major axes of pulsars A and B, s = sin i is the Shapiro-delay shape parameter, σ_so is the spin-orbit coupling constant and 𝒢 is a generalised gravitational constant for the interaction between the two pulsars. If GR is the correct theory of gravity, then we expect $(\frac{c^2{\sigma }_{\mathrm{so}}}{G}{)}_{\mathrm{GR}}=2+\frac{3}{2}\frac{m_{\mathrm{A}}}{m_{\mathrm{B}}}=3.6076796\pm 0.0000021$. Substituting in our best value of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ and measurements for the post-Keplerian parameters from Kramer et al. (2021a) into Eq. (16), we obtain $(\frac{c^2{\sigma }_{\mathrm{so}}}{G})=3.54\pm 0.27$. Taking the ratio of the observed and predicted values, we find $(\frac{c^2{\sigma }_{\mathrm{so}}}{G}{)}_{\mathrm{obs}}/(\frac{c^2{\sigma }_{\mathrm{so}}}{G}{)}_{\mathrm{GR}}=(1+0.017)\pm 0.061$. Hence, our measured spin-orbit coupling constant is consistent with the expectation from GR to within an uncertainty of 6.1%. The level of consistency between this test and the comparison between the predicted and observed precession rates from earlier is expected, as both are reliant on the same measured value.

These measurements of the geodetic precession rate of pulsar B and the associated test of spin-orbit coupling are encouraging, though the current unfavourable geometric alignment of pulsar B and systematic effects that arise from the model inconsistencies may introduce additional uncertainties on top of our reported statistical values. Another source of systematic uncertainties we have not explored here is the impact the current lack of phase-coherent timing of pulsar B has on our measurements. Breton et al. (2008) noted that variations in the input period of pulsar B alter both the expected rotational phase during the eclipses and duration of the transparency windows, which can lead to slightly faster or slower precession rates being inferred. As a result of these effects, the values presented here may have somewhat underestimated uncertainties.

5.2. System geometry

Combined with precision timing of pulsar A over the past two decades (Kramer et al. 2021a), our eclipse modelling results provide the clearest picture to date of the three-dimensional geometry of PSR J0737−3039A/B. As noted in Sect. 4.2, the temporal evolution of the eclipse modulation pattern due to geodetic precession breaks several symmetries in the light-curve model. From this we can unambiguously state that the spins of both pulsars are pointed out of the orbital plane of the system, with the spin axis of pulsar B being misaligned from the total angular momentum vector by δ_B = 180° −θ = 40.6° ±0.1°. This suggests pulsar B is rotating prograde in its orbit around pulsar A, consistent with beam-modelling results performed by Noutsos et al. (2020). Note that previous works incorrectly assumed that δ_B ≈ θ, where the seemingly large spin-axis offset was interpreted as arising from an off-centre supernova kick causing the newly born pulsar to tumble (Farr et al. 2011). The corrected value of δ_B suggests the supernova kick resulted in only a small offset in the final spin-axis away from the total angular momentum vector. Combined with the precession rate, the angular offset between the spin axes of pulsars A and B can be computed as (Perera et al. 2014),

$$\begin{array}{c}\hfill cos(\mathrm{\Delta }S(t))=cos({\delta }_{\mathrm{A}})cos({\delta }_{\mathrm{B}})+sin({\delta }_{\mathrm{A}})sin({\delta }_{\mathrm{B}})cos({\phi }_{\mathrm{SO}}(t)),\end{array}$$(17)

where δ_A is the angle between the spin axis of pulsar A and the total angular momentum vector, δ_B = 180° −θ and φ_SO(t) is computed from Eq. (12). Using the Perera et al. (2014) upper-limit of δ_A ≤ 2°, our measurements of θ, φ₀ and the GR-predicted value of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$, we obtain ΔS = 40° ±2° at our reference epoch of MJD 59289. While our value is inconsistent with the ΔS = 138° ±5° obtained by Perera et al. (2014), we note they had assumed that δ_B ≈ θ. Correcting their measurement as 180° −ΔS returns a value of 42° ±5° which is in-line with what we obtained.

Measurements of the Shapiro delay shape parameter from timing pulsar A return two equally-likely solutions, where i = 89.35° ±0.05° or 90.65° ±0.05° (Kramer et al. 2021a). Attempts have been made to determine the ‘sense’ of i through measurements of the varying scintillation timescale of pulsar A throughout its orbit, with works by Ransom et al. (2004) and Rickett et al. (2014) finding i = 88.7° ±0.9° and 88.1° ±0.5° respectively. A separate study of the scintillation properties of both pulsars by Coles et al. (2005) found |i − 90° | = 0.29° ±0.14°. Modelling of the double pulsar eclipses provides an independent means of determining the orbital inclination angle as the recovered values of θ and φ₀ depend on the sign of z₀. Previously it was thought that symmetries between model parameters would make it difficult to determine the correct sign of z₀ (Breton 2009). However we demonstrated in Sect. 4.2 that the precession of pulsar B effectively breaks this degeneracy, and only a single peak consistent with a negative value of z₀ = −0.5408 ± 0.0009 R_mag was recovered. Assuming the small-angle approximation we can compute the orbital inclination by multiplying the z-axis offset by the magnetospheric extent, obtaining i = 90° −z₀ ξ = 90.501° ±0.001°. Note that i in this case is measured using the ‘DT92/RVM’ convention (Damour & Taylor 1992), where it corresponds to the angular offset of the orbital angular momentum (L) from a vector pointing along our line from the Earth through the centre of mass (K₀). The full system geometry represented in the Cartesian coordinate system utilised by the eclipse model is illustrated in Fig. 8. At the apparent level of precision for which we measure i, it appears to be in a 3-σ tension with the i > 90° value obtained from the Shapiro delay of pulsar A. But, as with the tests of GR, it is likely that the uncertainties on z₀ and ξ are underestimated by a substantial factor due to unmodelled systematic effects. In particular, both of these parameters should display some level of time-dependence as the projected distance between the pulsars vary due to their elliptical orbits undergoing periastron advance. Implementation of these time-dependencies will be explored in future work.

Fig. 8. Diagram showing our inferred system geometry for the double pulsar using the same Cartesian coordinate system that is depicted in Fig. 2. The inclination angle (i) is defined as the angle between the orbital angular momentum vector (L) and a vector pointing away from the observer (K0; Damour & Taylor 1992). We note that the offset of i from 90° is exaggerated by a factor of ∼20 in this image. All other angles are approximately correct.

Our recovered inclination sense disagrees with the aforementioned results from scintillation studies, where i < 90° is preferred. Determining the orbital properties of the double pulsar via scintillometry is challenging due to substantial time-variations in the observed scintillation timescale, which Rickett et al. (2014) attempted to overcome by modelling the anisotropy of the turbulent screen along the line of sight. However, it is possible that more complex geometric models are required to fully capture rapidly changing temporal variations in the scintillation pattern. These improved models are currently under development using a combination of wide-bandwidth data collected with the MeerKAT UHF, L-band and S-band receivers (Askew et al., in prep.). Preliminary results indicate i > 90° is preferred. We also note that an independent measure of the inclination angle from performing a rotating vector model fit to the linear polarisation of pulsar A under the assumption δ_A ≈ 0° returned i = 91.6° ±0.1° (Kramer et al. 2021a), albeit with uncertainties that only reflect the statistical errors in the measurement and do not take into account potential systematic effects.

5.3. Implications for beam modelling and B’s return

Geodetic precession of pulsars in relativistic binaries results in our line of sight passing through different regions of their emission cone over time. This affords us the rare opportunity to map the active radio-emitting field lines above the polar cap. The extent to which we can map a given pulsars beam depends on several factors: the shape of the active region, the beam filling factor, emission cone structure and pulsar geometry. For pulsar B this is further complicated by the intense wind from pulsar A which distorts the active magnetic field lines around the polar cap of pulsar B and resulted in the pulsar only being detected in discrete ‘visibility windows’ spread around different parts of the orbit (Lyne et al. 2004; Lyutikov 2005; Lomiashvili & Lyutikov 2014). While the profiles detected in each window displayed slight differences in overall pulse shape, they did share the same secular evolution from having a clear single-peak to a double-peaked profile (Burgay et al. 2005). Attempts to model this evolution has resulted in two competing interpretations. Early works suggested that pulsar B possesses a partially filled hollow cone, in which the radio-emitting region resembles an elongated horseshoe centred on the magnetic axis (Perera et al. 2010, 2012; Lomiashvili & Lyutikov 2014). A more recent re-analyses of the profile evolution suggested pulsar B has a wedge-shaped beam consisting of two elongated Gaussian components that diverge towards the outskirts of the emission cone (Noutsos et al. 2020). Both interpretations are strongly dependent on the assumed geometry of the pulsar, in particular the evolution of the angular offset between our line of sight and the magnetic axis (β) as the pulsar precesses and whether our line of sight was moving towards or away from the magnetic axis when pulsar B was visible.

The geometry and precession rate of pulsar B inferred from our eclipse modelling can be used to determine β as per (Breton 2009),

$$\begin{array}{c}\hfill cos(\pi -\zeta )=sin(\theta )cos(\phi )sin(i)+cos(\theta )cos(i),\end{array}$$(18)

where β = ζ − α. Using our measured geometry in Table 2 and that reported by Breton et al. (2008), we plot the resulting β(t) curves from both geometries in Fig. 9. Tracks of β(t) are shown for both the ‘active’ magnetic pole that was visible from 2003–2008 and the antipodal pole from which we may detect radio pulses in the future. Both sets of curves appear quite similar in terms of amplitude. There is however a substantial offset in the ‘phase’ between the two assumed geometries, which comes from the ∼28° difference in φ that are predicted for a given epoch. Consequently, these two sets of results provide conflicting interpretations for the time-evolution of the pulse profile of pulsar B. For the Breton et al. (2008) geometry, which was assumed to be correct in the beam-modelling of Noutsos et al. (2020), the magnetic axis of pulsar B crossed our line of sight just prior to its discovery and the disappearance occurred when β ∼ −20°. This was interpreted as emission cone of the pulsar missing our line of sight, with the corresponding value of β setting a limit on the opening angle of the pulsar beam and motivated the aforementioned wedge-shaped beam map. In contrast, the β evolution corresponding to our inferred geometry is consistent with the elongated horseshoe beam model of Perera et al. (2010), where the magnetic axis was initially at ∼20° when the pulsar was discovered and only crossed our line of sight in late-2009, after the radio pulses disappeared. Despite the values of α and θ they recovered for the two orbital visibility windows varying in consistency with our measurements (Table 2), the φ differs by only ∼4° when computed at the same reference epoch. The values of α and θ obtained via an improved magnetospheric model for pulsar B in Perera et al. (2012) match our measurements extremely well. However, the updated φ is much closer to that of Breton et al. (2008), differing from ours (at the same reference epoch) by ∼30°. This results in a large shift in the predicted pulse-shape evolution with time (compare Fig. 17 of Perera et al. 2010 and Fig. 6 of Perera et al. 2012) and substantial differences in the predicted date that radio pulses will again be detected from pulsar B. Given our improvements in modelling the eclipses over Breton et al. (2008), we suggest that the beam shape of pulsar B most likely resembles the elliptical, partially filled hollow cone described by Perera et al. (2010), or an inverted variant of the two-component wedge model where the separation between the two components increases as our line of sight cuts closer to the magnetic axis.

Fig. 9. Predicted change in the impact parameter between the magnetic axis of pulsar B and our line of sight over time. Solid lines indicate the ‘active’ beam that was detected in from 2003–2008, dashed lines correspond to the opposite magnetic pole. The vertical shaded region indicates the time-span over which radio pulses were detected from pulsar B.

If the radio beam of pulsar B were symmetric about the magnetic axis then our geometry predicts we should have re-detected the pulsar mid-way through 2011, which did not occur. Searches for radio pulses at the GBT and Parkes (‘Murriyang’) radio telescope over the past 15 yr have so far returned non-detections (R. N. Manchester, private communication). This is interesting because little is known with certainty about the two-dimensional shape of a radio pulsar beam map (but see Desvignes et al. 2019). As about half of all pulsars show multiple pulse components, the emission intensity certainly varies along our line of sight (i.e. in a roughly longitudinal direction over the neutron star). Yet, how the emission varies with latitude is mostly unexplored. The two-dimensional emission can be modelled as roughly circular; and shine completely, or only from either cones or patches (cf. Lyne & Manchester 1988). Some models, on the other hand, predict fan-like beams (Dyks et al. 2010).

The continued absence of pulsar B suggests that the lower half of the emission cone must be entirely devoid of active field lines. This is unusual for pulsars in general as we know emission arises on either side of the magnetic pole (Johnston et al. 2023) and is not the case for any of the other precessing pulsars. The number of other pulsar systems with beam maps to compare against is limited, with considerable variation in the underlying sources. PSR B1913+16, on the one hand, is an old pulsar where recycling has diminished the magnetic field strength, and, arguably, influenced its configuration. Its beam shape is consistent with either a hollow-cone or nested-hollow-cone depending on the chosen beam structure and viewing geometry (Weisberg & Taylor 2002; Clifton & Weisberg 2008). In both cases the emission can be fit by a circularly symmetric conical beam. On the other hand, PSRs J1141−6145 and J1906+0746 are young and unrecycled. The radio beam of PSR J1141−6145 appears to be largely filled (Manchester et al. 2010); in some contrast to the behaviour we infer for pulsar B. PSR J1906+0746, however, combines both types of beam shape in a single source. This nearly orthogonal rotator allowed for the creation of beam maps from both poles (Desvignes et al. 2019). The evolution of its ‘interpulse’ suggests that emission from this pole, the only one currently visible, is generally circular, and filled; similar to PSR J1141−6145. Its ‘main pulse’, however, which has since disappeared due to geodetic precession, consists of a narrow fan-like strip of emission that extends radially from the magnetic axis. Irregular beam shapes are therefore possible and can help explain the continued absence of pulsar B.

The limited understanding of the shape of the pulsar B emission cone introduces substantial uncertainties when attempting to predict the return of its radio pulses. Breton (2009) and later Perera et al. (2012) predicted that the pulsar would re-enter the same range of β values as 2003–2008 as soon as 2024 (see the solid blue trace in Fig. 9). However, our geometry predicts the same range will not be reached until early 2035, meaning the return of detectable radio pulses from pulsar B may be a decade later than initially predicted. This is consistent with the earlier beam modelling of Perera et al. (2010). Similarly, radio emission from the opposite magnetic pole may not appear until sometime in the 2040’s to 2050’s as opposed to ∼2034 as predicted in Breton (2009). In either case, the re-appearance of pulsar B, while requiring patience, will be an important input for understanding the beam shapes of non-recycled pulsars.

6. Summary and conclusions

We presented the first results of our ongoing MeerKAT campaign to monitor the eclipses detected in the double pulsar system PSR J0737−3039A/B. Our Bayesian inference framework is capable of modelling the eclipse light-curve in the absence of a phase-coherent timing model of pulsar B, thereby allowing us to measure the pulsar geometry without a priori knowledge of the pulsar rotation phase. However, the high sensitivity of MeerKAT, combined with the current geometry of pulsar B as viewed from Earth and limitations of the Lyutikov & Thompson (2005) light-curve model presented significant challenges to obtaining a robust measurement of the geodetic precession rate of pulsar B. Using a hierarchical Bayesian inference technique, we showed the recovered precession rate is strongly dependent on how much of the ingress and egress phases of the eclipses are included in the fits, thereby limiting the robustness of our current approach to measuring ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ and the associated tests of general relativity. With this in mind, we demonstrate that only fitting for the light-curve data within ±0.72° of superior conjunction returned a precession rate of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}={5.16°}_{-0.34°}^{+0.32°}$ yr⁻¹, consistent with the expected value from GR to within 6.5% uncertainty. This in-turn provided an update to a theory-independent test of strong-field spin-orbit coupling, where our measurement of $(\frac{c^2{\sigma }_{\mathrm{so}}}{G})=3.54\pm 0.27$ is consistent with GR at the 6.1% level. We showed the precession rate measurement could be improved via the use of an iterative framework to fit all of the observed eclipses simultaneously. Yet the ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ that we obtained via this approach is only consistent with GR to within the 99.7% confidence interval. We suggest this is likely a result of the aforementioned unfavourable pulsar geometry and limitations of the existing eclipse model.

Our measurements of the system geometry are consistent with the spin axis of pulsar B being offset from the total angular momentum vector by 40.6 ± 0.1°. This confirms that both pulsars A and B rotate in a prograde direction with respect to their orbital motion. We also showed that our precession measurement broke several symmetries between model parameters, allowing us to determine the sense of the inclination angle i = 90.501° ±0.001° > 90°. Upcoming scintillation measurements will confirm or rule out the presence of a tension with the i < 90° that is favoured in past studies. Finally, we discussed our improved geometry in the context of past attempts to map the radio beam of pulsar B and provided an updated prediction for when radio pulses may again be detected. Our inferred β(t) suggest the radio beam likely resembles the elongated horseshoe shape presented in Perera et al. (2010), and that the pulsar will not return to the same viewing geometry as 2008 until early 2035.

Ultimately, many of the challenges faced in this work will be addressed through developing models of the pulsar B magnetosphere that better reflect the observed eclipse phenomenology. Improvements to account for frequency evolution across the wide fractional bandwidth of MeerKAT may also help reduce systematic uncertainties such that effective ToAs could be produced for pulsar B. This would provide updated timing measurements of the double pulsar independent to pulsar A, as well as enable phase-coherent searches for radio pulses from pulsar B. Construction of these models is particularly important given the near-future MeerKAT extension project and the eventual integration of MeerKAT into the Square Kilometer Array Mid-frequency telescope. These will provide substantial increases in telescope gain, allowing us to resolve the eclipses in even greater detail than we present here. Extending the baseline over which the precession rate is measured through continued monitoring with MeerKAT and linking with the long-running GBT monitoring campaign will help overcome issues pertaining to the unfavourable viewing geometry. A combined MeerKAT+GBT dataset would return a precision measurement of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$, enabling associated tests of gravity that are on-par with the measurements of other relativistic effects in this unique system.

Acknowledgments

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory (SARAO), which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time space reference systems department of the Paris Observatory. PTUSE was developed with support from the Australian SKA Office and Swinburne University of Technology. This work made use of the OzSTAR national HPC facility at Swinburne University of Technology. MeerTime data are housed on the OzSTAR supercomputer. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government. MK acknowledges significant support from the Max-Planck Society (MPG) and the MPIfR contribution to the PTUSE hardware. 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. RMS acknowledges support through Australian Research Council (ARC) Future Fellowship FT190100155. Part of the work was undertaken with support from the ARC Centre of Excellence for Gravitational Wave Discovery (CE1700004). RPB acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 715051; Spiders) VVK acknowledges financial support from the ERC starting grant ‘COMPACT’ (grant agreement number: 101078094). JvL acknowledges funding from NWO Vici research programme ‘ARGO’ (639.043.815). MAM is supported by NSF Physcis Frontiers Center award number 2020265. Pulsar research at UBC is funded by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. MEL thanks the MPIfR fundamental physics research group for their hospitality during a visit to work on this project. This work has made use of NASA’s Astrophysics Data System.

References

All Tables

All Figures

Fig. 1. Normalised intensity of pulsar A throughout each of the eclipses recorded by MeerKAT. The red vertical line indicates the orbital phase of superior conjunction. Regions contained within the white lines correspond to where the light curve model best matches the data (see Sect. 4.1). Cyan and orange indicators on the right-hand side denote light-curves obtained using the L-band and UHF receivers respectively. We note that the observation numbers on the y-axis do not increase linearly with time.

In the text

Fig. 2. Diagram depicting the double pulsar geometry according to the eclipse light-curve model, adapted from Fig. 1 of Breton et al. (2008). The minimum radial extent of the synchrotron absorbing plasma (the cooling radius) is represented by Rmin. Longitude and latitude of the spin axis of pulsar B (ΩB) is given by the angles φ and θ respectively. The magnetic axis (μB) is offset from ΩB by angle α. In our model we have neglected the deflection of A’s radio signals in the gravitational field of B, which amounts to a change in the impact parameter at B of no more than about 600 km (Kramer et al. 2021a).

In the text

Fig. 3. Two eclipses of pulsar A detected with MeerKAT using the UHF (left) and L-band (right) receiver systems. These light curves were generated after averaging over the full receiver bandpasses. We have inflated the flux uncertainties by the median recovered EQUAD for each eclipse. The orange traces correspond to the median a posteriori light-curve models recovered from fitting the observations. Bottom panels show the residuals after subtracting the median model.

In the text

Fig. 4. Comparison between hierarchical fits to different realisations of the light curves using the FULLBAND and SUBBAND datasets. Contours in the two-dimensional posteriors represent the 68% and 95% confidence regions. Each colour represents measurements from eclipses that had different amounts of the ingress or egress phases removed. The dashed lines correspond to the predicted values of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$ from general relativity.

In the text

	Fig. 5. Evolution of our recovered geometric parameters for pulsar B over time measured using the FULLBAND (±0.72°) data set. The GR-predicted evolution of φ due to Geodetic precession is indicated by the black line. Dark orange shading indicates 68% range covered by the precession rate inferred from our hierarchical fit. Light orange shading includes scatter in the values of φ over time.
In the text

Fig. 6. One- and two-dimensional posterior distributions of the eclipse light-curve parameters from iteratively fitting the FULLBAND data with the ±0.72° cuts where the precession rate was either fixed to the GR predicted value (green) or allowed to be fit as a free parameter (orange). Contours indicate the 68, 95 and 99.7% confidence intervals. Dashed black line is the predicted value of ${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$.

In the text

Fig. 7. Mass-mass diagram illustrating current tests of general relativity with the double pulsar. Our improved measurement of pulsar B’s spin-precession rate (${\mathrm{\Omega }}_{\mathrm{SO}}^{\mathrm{B}}$) is given by the (mostly) vertical blue lines where the separation between them indicates the 68% uncertainty. Other individual constraints on the masses of pulsar A and B from Kramer et al. (2006, 2021a) are shown by pairs of lines. These include the periastron precession rate ($\dot{\omega }$), shrinking of the orbital period due to gravitational radiation (Ṗb), gravitational redshift (γ), Shapiro delay range (r) and shape (s), and the mass-ratio (R). The shaded region is forbidden by the individual mass functions of the two pulsars (i.e. sin i ≤ 1).

In the text

Fig. 8. Diagram showing our inferred system geometry for the double pulsar using the same Cartesian coordinate system that is depicted in Fig. 2. The inclination angle (i) is defined as the angle between the orbital angular momentum vector (L) and a vector pointing away from the observer (K0; Damour & Taylor 1992). We note that the offset of i from 90° is exaggerated by a factor of ∼20 in this image. All other angles are approximately correct.

In the text

	Fig. 9. Predicted change in the impact parameter between the magnetic axis of pulsar B and our line of sight over time. Solid lines indicate the ‘active’ beam that was detected in from 2003–2008, dashed lines correspond to the opposite magnetic pole. The vertical shaded region indicates the time-span over which radio pulses were detected from pulsar B.
In the text