A MeerKAT view of the double pulsar eclipses Geodetic precession of pulsar B and system geometry

The double pulsar system, PSR J0737 − 3039A / B, consists of two neutron stars bound together in a highly relativistic orbit that is viewed nearly edge-on from the Earth. This alignment results in brief radio eclipses of the fast-rotating pulsar A when it passes behind the toroidal magnetosphere of the slow-rotating pulsar B. The morphology of these eclipses is strongly dependent on the geometric orientation and rotation phase of pulsar B, and their time evolution can be used to constrain the geodetic precession rate of the pulsar. We demonstrate a Bayesian inference framework for modelling high-sensitivity eclipse light curves obtained with MeerKAT between 2019 and 2023. Using a hierarchical inference approach, we obtained a precession rate of Ω BSO = 5 . 16 ◦ + 0 . 32 ◦ − 0 . 34 ◦ yr − 1 (68% confidence intervals) for pulsar B, consistent with predictions from general relativity to a relative uncertainty of 6.5%. This updated measurement provides a 6.1% test of relativistic spin-orbit coupling in the strong-field regime. We show that a simultaneous fit to all of our observed eclipses can in principle return a ∼ 1.5% test of spin-orbit coupling. However, systematic e ff ects introduced by the current geometric orientation of pulsar B along with inconsistencies between the observed and predicted eclipse light curves result in di ffi cult to quantify uncertainties when using this approach. Assuming the validity of general relativity, we definitively show that the spin axis of pulsar B is misaligned from the total angular momentum vector by 40 . 6 ◦ ± 0 . 1 ◦ and that the orbit of the system is inclined by approximately 90 . 5 ◦ from the direction of our line of sight. Our measured geometry for pulsar B suggests the largely empty emission cone contains an elongated horseshoe-shaped beam centred on the magnetic axis, and that it may not be re-detected as a radio pulsar until early 2035.


Introduction
PSR J0737−3039A/B is a highly relativistic double neutron star binary with a short, 2.45 hr 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 7 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 (Ω B SO ) 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.Team Compas 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. 2010andJ1906+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 Section 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 Section 4, including an outline of the challenges faced given the current system geometry and model compatibility.In Section 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 Section 6 and outline potential future directions for modelling of the eclipses.

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 Meer-Guard1 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 frequencychannels 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.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.

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 2 Solar System ephemeris.The barycentred arrival times were then converted to orbital phases using the latest pul-2 https://ssd.jpl.nasa.gov/planets/eph_export.html  et al. (2008).The minimum radial extent of the synchrotron absorbing plasma (the cooling radius) is represented by R min .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).
sar A timing ephemeris from Hu et al. (2022) 3 .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 Figure 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.

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 0 .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 cos θ = cos(90 sin φ = − sin(180 • − δ B ) sin φ so sin θ . (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 Figure 2. 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 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 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 Equation 3.Both of these terms are related to the dipole unit vector magnetic polar angle (θ µ ) as where r = {x, y(t), z}, r is the distance between pulsars A and B in spherical coordinates (r = |r| = x 2 + y 2 (t) + z 2 ) and μ is the dipole unit vector, the components of which are given by μx = ( μΩ Here, α is the magnetic inclination angle of pulsar B and ϕ B is the rotation phase of pulsar B, which is related to the spinvector direction as ϕ B = Ω B t = 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 and before being passed to Equation 3 where they are integrated over x.

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-toepoch 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 where d represents the input eclipse light curves, µ is the eclipse model, Θ contains the model parameters, and σ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 pulseto-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 lightcurves 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.

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 angularmomentum 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 Figure 2).As a result, the geodetic precession of pulsar B manifests as a timevarying 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 and where θ 0 and φ 0 are the co-latitude and longitude of B's spin axis at some reference epoch, t 0 = 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.

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 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.

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 where the time dependence of Θ comes from computing φ at each eclipse epoch via Equation 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 Ω B SO and φ 0 when fitting the eclipses.

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 Figure 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 Figure 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.
In Figure 4, we show that the recovered φ 0 , precession rate and σ φ from hierarchical fits to both versions of the eclipse data vary depending on whether data is 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 −1 can be computed as (Barker & O'Connell 1975) 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 Ω B SO = 5.074005 • ±0.000003 • yr −1 , which is shown as the dashed lines in Figure 4. Comparing the different fits in Figure 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 Ω B SO 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 Figures 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 Figure 2 of McLaughlin et al. 2004 andFigure 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 lightcurve 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 Figure 5 we show the recovered geometric parameters for pulsar B from our fits to individual pairs of eclipses.As expected, α, µ, z 0 and ξ all remain at near constant values, while φ displays a significant decrease from 13.0 • ± 0.5 • to −6.2 •+1.6 • −1.4 • 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 θ = θ 0 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 Figure 5.Our recovered precession rate of Ω B SO = 5.16 •+0.32 • −0.34 • yr −1 agrees with the expected value of 5.074005 • yr −1 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 Ω B SO = 5.31 •+0.42 • −0.41 • yr −1 .In principle, the iterative method from Section 3.4.2should 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 While the 68% confidence intervals represent a ∼1.5% uncertainty on Ω B SO , 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.

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 −1 , 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 Figure 6 we present the one-and two-dimensional posterior distributions on the model parameters.For comparison, we also show the same model parameters where Ω B SO was searched over as a free parameter.
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 B = 0.05 in favour of the 'Not Assuming GR' model indicates the timeevolution 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 φ 0 (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 0 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 0 , resulting in bimodal marginalised posterior distributions for the individual parameters and four modes in the two-dimensional posteriors with each other.These Orbital phase (deg) 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.symmetries appear as ∆ϕ → ∆ϕ B ± 0.5, {α, θ, φ} → 180 • − {α, θ, φ} and −z 0 → +z 0 .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 Section 3.3 with wider prior ranges on the aforementioned parameters of interest.Our recovered posterior distributions for α, θ, φ, and z 0 were all well constrained and confined to a single mode peaking at the previously reported val- ues 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.
Table 2: Recovered model parameters from the iterative fits to the Fullband dataset using ±0.72 • cuts with and without the assumption that GR is the correct theory of gravity.Uncertainties on each parameter represent the 68% confidence intervals.The value of φ 0 is referenced to MJD 59289 (2021 March 16).In the last row we give the Bayesian evidences for the two models.

Parameter (units)
Assuming GR Not Assuming GR

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 Figure 7, we have plotted the inferred masses of pulsars A and B from Ω B SO 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 Ω B SO 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.
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 Ω B SO 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 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 G is a generalised gravitational constant for the interaction between the two pulsars.If GR is the correct theory of gravity, then we expect c 2 σ so 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 Ω B SO .
et al. (2021a) into Equation 16, we obtain c 2 σ so G = 3.54 ± 0.27.Taking the ratio of the observed and predicted values, we find = (1 + 0.017) ± 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 phasecoherent 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.

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 geom- ) 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. (2006Kramer et al. ( , 2021a) ) are shown by pairs of lines.These include the periastron precession rate ( ω), 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).etry of PSR J0737−3039A/B.As noted in Section 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), cos(∆S (t)) = cos(δ A ) cos(δ B ) + sin(δ A ) sin(δ B ) cos(φ SO (t)), (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 Equation 12. Using the Perera et al. (2014) upper-limit of δ A ≤ 2 • , our measurements of θ, φ 0 and the GRpredicted value of Ω B SO , 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 φ 0 depend on the sign of z 0 .Previously it was thought that symmetries between model parameters would make it difficult to determine the correct sign of z 0 (Breton 2009).However we demonstrated in Section 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 = −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 0 ξ = 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 0 ).The full system geometry represented in the Cartesian coordinate system utilised by the eclipse model is illustrated in Figure 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 0 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.
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 timevariations 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.

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 z O rb ita l pl an e (To observer) x Fig. 8: Diagram showing our inferred system geometry for the double pulsar using the same Cartesian coordinate system that is depicted in Figure 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 (K 0 ; 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.
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 doublepeaked 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(Perera et al. , 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), 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 Figure 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 Figure 17 of Perera et al. 2010 and Figure 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.
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 yrs 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 twodimensional 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 fanlike 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 Figure 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.

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 Ω B SO 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 Ω B SO = 5.16 •+0.32 • −0.34 • yr −1 , 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 c 2 σ so G = 3.54 ± 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 Ω B SO 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 Midfrequency 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 Ω B SO , enabling associated tests of gravity that are on-par with the measurements of other relativistic effects in this unique system.

Fig. 1 :
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 Section 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.

Fig. 2 :
Fig. 2: Diagram depicting the double pulsar geometry according to the eclipse light-curve model, adapted from Figure 1 of Bretonet al. (2008).The minimum radial extent of the synchrotron absorbing plasma (the cooling radius) is represented by R min .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).

Fig. 4 :
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 Ω B SO from general relativity.

Fig. 5 :
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.
3.6076796±0.0000021.Substituting in our best value of Ω B SO and measurements for the post-Keplerian parameters from Kramer

Fig. 7 :
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 (Ω B SO ) 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 fromKramer et al. (2006Kramer et al. ( , 2021a) ) are shown by pairs of lines.These include the periastron precession rate ( ω), 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).

Fig. 9 :
Fig. 9: Predicted change in the impact 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.