Issue 
A&A
Volume 643, November 2020



Article Number  A143  
Number of page(s)  38  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202038566  
Published online  17 November 2020 
Understanding and improving the timing of PSR J0737−3039B
^{1}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
email: anoutsos@mpifrbonn.mpg.de
^{2}
Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
^{3}
Department of Physics, West Virginia University, Morgantown, WV 26506, USA
^{4}
Australia Telescope National Facility, Commonwealth Scientific and Industrial Research Organisation (CSIRO), PO Box 76, Epping, NSW 1710, Australia
^{5}
INAFOsservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy
^{6}
University of Manchester, Jodrell Bank Observatory, Macclesfield SK11 9DL, UK
^{7}
University of Manchester, School of Physics and Astronomy, Jodrell Bank Center for Astrophysics, Alan Turing Building, Manchester M13 9PL, UK
^{8}
Alfred P. Sloan Research Fellow
^{9}
National Radio Astronomy Observatory, Green Bank, WV 24944, USA
^{10}
LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France
^{11}
Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA
^{12}
Università di Cagliari, Dip di Fisica, S.P. MonserratoSestu Km 0,700, 09042 Monserrato, Italy
^{13}
School of Chemistry, University of East Anglia, Norwich NR4 7TJ, UK
Received:
2
June
2020
Accepted:
31
August
2020
The double pulsar (PSR J0737−3039A/B) provides some of the most stringent tests of general relativity (GR) and its alternatives. The success of this system in tests of GR is largely due to the highprecision, longterm timing of its recycledpulsar member, pulsar A. On the other hand, pulsar B is a young pulsar that exhibits significant shortterm and longterm timing variations due to the electromagneticwind interaction with its companion and geodetic precession. Improving pulsar B’s timing precision is a key step towards improving the precision in a number of GR tests with PSR J0737−3039A/B. In this paper, red noise signatures in the timing of pulsar B are investigated using roughly a fouryear time span, from 2004 to 2008, beyond which time the pulsar’s radio beam precessed out of view. In particular, we discuss the profile variations seen on timescales ranging from minutes – during the socalled “bright” orbital phases – to hours – during its full 2.5 h orbit – to years, as geodetic precession displaces the pulsar’s beam with respect to our line of sight. Also, we present our efforts to model the orbitwide, harmonic modulation that has been previously seen in the timing residuals of pulsar B, using simple geometry and the impact of a radial electromagnetic wind originating from pulsar A. Our model successfully accounts for the longterm precessional changes in the amplitude of the timing residuals but does not attempt to describe the fast profile changes observed during each of the bright phases, nor is it able to reproduce the lack of observable emission between phases. Using a nested sampling analysis, our simple analytical model allowed us to extract information about the general properties of pulsar B’s emission beam, such as its approximate shape and intensity, as well as the magnitude of the deflection of that beam, caused by pulsar A’s wind. We also determined for the first time that the most likely sense of rotation of pulsar B, consistent with our model, is prograde with respect to its orbital motion. Finally, we discuss the potential of combining our model with future timing of pulsar B, when it becomes visible again, towards improving the precision of tests of GR with the double pulsar. The timing of pulsar B presented in this paper depends on the size of the pulsar’s orbit, which was calculated from GR, in order to precisely account for orbital timing delays. Consequently, our timing cannot directly be used to test theories of gravity. However, our modelling of the beam shape and radial wind of pulsar B can indirectly aid future efforts to time this pulsar by constraining part of the additional red noise observed on top of the orbital delays. As such, we conclude that, in the idealised case of zero covariance between our model’s parameters and those of the timing model, our model can bring about a factor 2.6 improvement on the measurement precision of the mass ratio, R = m_{A}/m_{B}, between the two pulsars: a theoryindependent parameter, which is pivotal in tests of GR.
Key words: pulsars: individual: PSR J07373039A/B
© A. Noutsos et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
The PSR J0737−3039 system is commonly known as the double pulsar and is composed of an old, recycled pulsar (pulsar A), with a spin period of P_{A} ≈ 22.7 ms, and a young pulsar (pulsar B), with a spin period of P_{B} ≈ 2.77 s (Burgay et al. 2003; Lyne et al. 2004). The two pulsars orbit each other in tight, loweccentricity orbits: the orbital period and eccentricity of the orbits is P_{b} ≈ 2.45 h and e ≈ 0.088, respectively. The double pulsar exhibits a plethora of physical effects, many of them due to the intense gravitational interaction between the two neutron stars (NS), which are never separated by more than ≈3 lights. Moreover, the system is viewed almost perfectly edgeon. These properties have allowed for very precise tests of general relativity (GR; Kramer et al. 2006).
2. Previous work
2.1. Tests of GR with the double pulsar
Amongst its many properties, the double pulsar is unique in that it is the only doubleNS system we know to date, where we detect pulsed emission from both its members, thus allowing us to directly measure the length of its semimajor axes through pulsar timing. Pulsar A is a very stable rotator, exhibiting little timing noise. Solely from timing pulsar A, several orbital (Keplerian and postKeplerian) parameters have been measured with high precision for this system; and this pulsar’s timing precision has served in testing theories of gravity (Kramer et al. 2006, and in prep., Kramer & Wex 2009). The timing of pulsar B, on the other hand, exhibits significant amounts of systematic noise, over a variety of timescales (see section below). To date, the precision of a number of tests of GR with the Double Pulsar has been limited by the significantly less precise timing of pulsar B, which dominates the uncertainty of the ratio between the two pulsars’ semimajor axes, x_{B}/x_{A}: in a very large set of gravity theories, this ratio is equal to the mass ratio, R = m_{A}/m_{B}, an important theoryindependent parameter in those tests.
2.2. Timing pulsar B
Within an orbital period, pulsar B exhibits a periodic modulation of its brightness. More specifically, during two orbital phase ranges, called bright phase 1 (hereafter BP1) and bright phase 2 (hereafter BP2), the pulsar appears much brighter than anywhere else (Lyne et al. 2004). These two phases correspond to orbital phase (from the ascending node) of ≈210° and ≈280°, although it has been observed that their location and extent gradually change with time, at the rate of a few degrees per year (Burgay et al. 2005; Perera et al. 2010; hereafter PMK+10). During the rest of its orbit, for roughly 85% of the orbital period, pulsar B is barely detectable, which has granted this phase the title, weak phase (hereafter WP). Furthermore, during each of BP1 and BP2, the profile of pulsar B evolves dramatically as a function of orbital phase, in contrast to the WP, during which the profile evolution is significantly less. This significant profile evolution during each of the bright phases (BPs) was already noticed soon after its discovery (Lyne et al. 2004). The strong profile evolution during the BPs, coupled with the intermittent visibility during the orbit, has limited the amount of timing precision that can be achieved for this pulsar.
Kramer et al. (2006; hereafter KSM+06) presented a detailed timing analysis of pulsar B from observations between MJD52760 and MJD53736. The authors divided the orbit into five intervals (hereafter also “orbitalphase windows”) and generated a timing template per interval to time the pulsar. A different set of five templates was used for every subsequent threemonth period, thus accounting for secular pulseshape changes due to geodetic precession (see Sect. 2.4). They published pulsar B’s spin parameters, based only on timing data during the WP, having a much more stable profile than the BPs, while excluding BP1 and BP2 from their timing analysis on the grounds of significant profile variation during those phases.
2.3. The geometry of pulsar B
The inclination of the double pulsar’s orbit is i ≈ 89°, resulting in periodic, 30slong eclipses of pulsar A’s emission by the magnetosphere of pulsar B, as the latter pulsar moves in front of the former, at conjunction, during their orbit. Breton et al. (2008; hereafter BKK+08) successfully modelled the fluxdensity modulation of pulsar A during those eclipses, with a simple geometric model describing the relative orientation of the two pulsars. In that work, the authors were able to determine a number of parameters of pulsar B’s geometry with high precision: for example, the magnetic inclination, which is the angle between the spin and the magnetic axis, , and the angle between the spin axis and the orbital angular momentum, . It must be noted that in the geometry of BKK+08, the angle that was constrained was the “colatitude” of the spin axis, θ = 180° −δ, which is supplementary to the angle δ defined in Damour & Taylor (1992). Interestingly, the work of BKK+08 could not distinguish between δ (prograde spin) and 180° −δ (retrograde spin), both of these solutions being degenerate, and hence the sense of the pulsar’s rotation could not be uniquely determined.
2.4. Geodetic precession
The problem of timing pulsar B is exacerbated by the longterm modulation of the effects mentioned above, due to the large degree of geodetic precession that this pulsar exhibits. The misalignment of pulsar B’s spin axis with the total angular momentum of the binary system, the latter being approximately equal to the orbital angular momentum, L, results in the geodetic precession of pulsar B’s spin around L. BKK+08, through their modelling of pulsar A’s eclipses, measured the rate of geodetic precession to be yr^{−1}, which is consistent with the prediction of GR, meaning yr^{−1} (Barker & O’Connell 1975). Also, the precession phase, which is the angle between the plane containing L and the spin axis and that containing L and the line of sight (LOS), was determined to be . Geodetic precession causes the gradual shift of the trace of our LOS across the emission region of pulsar B, as a function of time, which in turn causes the observed pulse profile to change significantly over timescales of years (Burgay et al. 2005; PMK+10). Ultimately, the emission beam of pulsar B precessed out of view in March 2008 (MJD54552), hence limiting the total span of available pulsar B data, from its discovery to its disappearance, to ≈4.3 years (i.e. MJD52997–MJD54552). It should be noted that pulsar A does not exhibit geodetic precession, as its spin axis is closely aligned with L (Ferdman et al. 2013).
2.5. Theoretical modelling
The nature of the profileshape changes is poorly understood; however, it was proposed early on that the distortion of pulsar B’s emission is caused by the pressure of pulsar A’s relativistic wind (McLaughlin et al. 2004; Arons et al. 2004). In that context, there have been a number of attempts for a physical description: for example, Lyutikov (2005) proposed a model based on distorted Euler potentials of pulsar B’s magnetic dipole in order to explain the observed profile and intensity modulation during the orbit. Although the model predicts that the phase of the pulsed emission drifts by ∼15 ms across each of the BPs, those drifts progress in opposite directions for BP1 and BP2, according to the model; in contrast, observations show that the pulse drifts towards later pulse phases during both BP1 and BP2 (e.g. compare Fig. 3 of Ransom et al. 2005 to Fig. 5 of Lyutikov 2005).
More recently, notable efforts were made by Perera et al. (2012) and Lomiashvili & Lyutikov (2014; hereafter LL14) to model the shape and location of the emission region of Pulsar B. The former authors used all available pulse profiles at 820 MHz, from the pulsar’s discovery to its disappearance, to calculate the geometry of the bow shock that is formed at the equilibrium points between the wind pressure of pulsar A and the magnetic pressure of pulsar B. By tracing the last closed field lines of pulsar B’s confined magnetosphere, those authors were able to estimate emission heights, assuming they are produced at the locations where the tangent to the field lines coincide with the observer’s LOS. The estimated heights ranged from ≈100 to 400 km. Furthermore, that work tried to model the beam shape of pulsar B’s emission, assuming it follows the topology of concentric, elliptic hollow cones of diminishing intensity towards the edge of the beam. Using solely the observed pulseprofile widths during BP1, the authors performed a maximum likelihood analysis to determine the beam ellipticity: the beam’s semimajor to semiminor axis ratio was determined to be . In addition, their analysis constrained the values of α and δ to be and , which the authors claimed were both consistent within 2σ with those derived by BKK+08. We note, however, that in the geometry of Perera et al. (2012), the angle δ (the angle between the orbital momentum and the spin) is identified as θ, which is in fact supplementary to the angle θ defined in BKK+08. In that sense, the constraint placed on δ is only consistent with the degenerate solution of BKK+08 where the pulsar spin is retrograde with respect to the orbital motion. Finally, Perera et al. used the magnetohydrodynamic confinement model of Tsyganenko (2002a,b) to set an upper limit on the emission height, based on the maximum deflection angle of the polar magneticfield lines. The latter was , which is equal to the full width, along the semimajor axis, at 10% of the maximum intensity of their calculated elliptical beam. The upper limit on the emission height was ∼25 000 km, which is compatible with the values derived by the same authors from fieldline tracing, although less constraining.
Additionally, the work of LL14 modelled the orbital and secular variations of pulsar B, using the MHD model of Tsyganenko (2002a,b) and a Dungeytype model (Dungey 1961) to describe the deformation of pulsar B’s magnetosphere under the strong wind pressure of pulsar A. Motivated by the work of PMK+10, LL14 assumed a smooth, analytic shape for the emission beam (at rest) of pulsar B, whose intensity distribution resembles a “horseshoe” convolved with a 3D Gaussian profile. In particular, this shape is consistent with the secular evolution of the average BP1 and BP2 profiles of pulsar B (i.e. averaged over each of the BPs), which change from predominantly single peaked to double peaked, over the span of ≈4 years (PMK+10). LL14 determined the best beamshape and axialorientation parameters of pulsar B, as well as the height of radio emission, by fitting the model’s predictions of the peak pulse intensity in a pulse period to the observed intensity variations, as a function of orbital and precession phase. It is worth noting that their model did not consider the pulsephase delays with respect to an inertial system, and hence pulsar timing was not attempted in that work.
2.6. Paper layout
The present paper is an attempt to model the periodic and secular variations observed in the pulse profiles of pulsar B, using a geometric model of the interaction with its companion (pulsar A) and a synthetic model of the structure of its emission beam. The ultimate goal of this paper is to provide a semianalytic description of those variations, as a function of time and orbital phase, which can then be used to improve the timing of pulsar B and hence the precision of tests of gravity with the doublepulsar system. Towards that purpose, this work is organised as follows. In Sect. 3, we describe the data we have used and present an initial timing analysis of pulsar B in order to obtain a qualitative characterisation of the systematic timing variations, which we attempt to model in the following section. In addition, we present our measurements of physical properties such as the dependence of pulsar B’s average flux density on orbital phase and observation epoch. We also provide an estimate of the spectral index of pulsar B’s fluxdensity spectrum. In Sect. 4, we first employ an empirical, toymodel description of the observed systematic variations, based only on the initial timing analysis: this contributes to quantitative estimates of the magnitude and functional dependence of the systematic variations. Following that, we lay out the full 3D geometry of our main model, including the parametrisation of pulsar B’s emission beam. We then provide a description of our method of parameter estimation, based on a nested sampling algorithm and the observed fluxdensity profiles of pulsar B. Finally, we present the resulting most likely configuration of the pulsar’s emission geometry, as well as the most likely magnitude of the deflection of pulsar B’s beam by pulsar A’s wind. In Sect. 5, we discuss the magnitude of the timing improvements resulting from the use of our model in timing pulsar B. We also discuss the prospects of future modelling improvements that would give rise to measurements of postKeplerian parameters, such as beam aberration. At the end of the section, we present tantalising evidence for a systematic displacement of pulsar B’s BPs, likely caused by geodetic precession, using a simple model of the secular changes in their orbital locations. The paper concludes with Sect. 6, wherein we summarise and discuss the main results.
3. Data reduction
The data available for this work came from archival observations with the Parkes and GBT telescopes, at 685 and 1400 MHz, for Parkes, and 820 MHz, for GBT. A description of the observations can be found in KSM+06 and Kramer et al. (in prep.). The total span of the data was MJD53004 to MJD54496. The original data products of the GBT observations comprised subintegrations of 20 s or 180 s in length, which contained dedispersed and averaged profiles of pulsar B across a 48 MHz band centred at 820 MHz. The Parkes observations during MJD53004–53035 and MJD53102–53192 comprised subintegrations of 60 s in length, across a 64 MHz band centred at 685 MHz, and across a 256 MHz band centred at 1375 MHz, respectively. In total, the original data contained roughly 29 500 subintegrations. The data contained no polarisation information. In a complete reprocessing of the data, the GBT observations were converted from the proprietary format of the original BCPM data to PSRFITS archives, for further processing with PSRCHIVE (Hotan et al. 2004). To maintain a uniform temporal resolution across all our profiles, we set the pulsephase resolution to N_{bin} = 512 bins, downsampling where necessary: the corresponding temporal resolution is ≈5.5 ms. Following visual inspection of the data, we excised those subintegrations and frequency channels that were plagued by RFI. Subsequently, the remaining frequency channels in each archive were averaged together, leaving only the subintegrations in the data.
In order to determine the longterm temporal evolution of the emission properties of pulsar B, we uniformly grouped our data into 14 100day intervals (hereafter MJD bins). Table 1 shows the instruments and observation properties corresponding to the data products within each MJD bin. Apart from the poorly sampled interval MJD54300–54400, and the interval MJD54200–54300, where no data were available, all other bins contain at least 20 subintegrations. Our data set provides a nearly continuous coverage of pulsar B’s profile evolution, with a 100day resolution, of approximately four years.
Properties of the observations and data products used in the analysis of this paper, listed per 100day MJD bin.
The chosen length of the MJD bins was optimised to retain a useful amount of data in each MJD bin (see last column of Table 1), while at the same time limiting the profile variations across an MJD bin due to geodetic precession. Theoretically, across the MJD range covered by our data, geodetic precession causes a drift of pulsar B’s magnetic axis with respect to our LOS at the rate of dβ/dt ∼ sin ϕ_{SO}sin(δ − α)Ω_{SO} ≲ 1° per MJD bin, where β is the impact angle, that is, the angle between the LOS and the magnetic axis at the spin phase of the closest approach. Observationally, PMK+10 determined the rate of broadening pulsar B’s profile to be yr^{−1}, for ca. MJD53900–54500. Given that the profile’s FWHM in BP1 and ca. MJD54000 is ≈6°, the expected intrabin smearing across an MJD bin is less than 15%.
As has been reported in previous studies, pulsar B is very bright during BP1 and BP2, while it is barely detectable elsewhere, particularly during the WP. Moreover, for at least 800 days, from MJD53100 to MJD53900, pulsar B appears to have an orbitalphase window of intermediate brightness, roughly an order of magnitude less than that of BP1, which we call the intermediate phase (hereafter IP). To determine the locations and extents of BP1, BP2, and the IP across our data, first we calculated the signaltonoise ratio (hereafter S/N) of the profiles contained in each MJD bin. Then, we determined the orbital phase corresponding to each profile, from its time stamp and the orbital parameters of pulsar A, appropriately adding 180° to the longitude of periastron (Kramer et al., in prep.). More specifically, we converted the site arrival times (SATs) to barycentric arrival times (BATs) using TEMPO2 and calculated for each SAT the phase from the ascending node (ϕ_{asc}) by solving Kepler’s equations, at each step correcting for the periastron advance, . The distribution of S/N as a function of ϕ_{asc}, for each MJD bin, is shown in Fig. B.1, where it can be seen that in discrete segments of the orbit the distribution exhibits a roughly symmetric rise and fall of the S/N, with relatively sharp maxima. We approximated those features with a normal distribution, and fitted a Gaussian function to the S/N distribution of each of the BPs, and, depending on whether it was visible, that of the IP. The centroids and 3σ widths of BP1, BP2, and the IP for each MJD bin, corresponding to the maxima and standard deviations of the fitted Gaussians, are shown in Table 2. In this work, the WP was defined as the orbitalphase interval between the upper 3σ bound of the IP and the lower 3σ bound of BP1, which corresponds to roughly ϕ_{asc}/2π ∈ [0.15, 0.53]. As can be seen in Fig. 1, the location of BP1 and BP2 changes significantly with time, whereas the orbitalphase evolution of the IP is more uncertain.
Fig. 1. BP1,BP2 and IP: Orbitalphase evolution of the centroids of BP1, BP2, and the IP as a function of MJD epoch. The abscissa values correspond to the mean MJD across each bin (see Table 2). The dashed red lines are the best linear fits to the data; the slope and reduced chisquared of each fit is shown in the legend. ⟨BP1+BP2⟩: The average orbitalphase shift of the centroids of BP1 and BP2 as a function of MJD epoch; the best linear fit is shown with a dashed red line. 
In addition to the phase shifts of BP1 and BP2, there is significant brightness evolution of BP1, BP2, and the IP. We used the radiometer equation to calculate the mean flux density of the observed profiles, based on their S/N, integration length, bandwidth, and pulse width at 10% of the maximum. In particular, for the 820 MHz observations with the GBT we used T_{sys} = 35 K (based on PMK+10) and G = 2 K Jy^{−1}. In the MJD intervals 53000–53100 and 53100–53200, there were no available observations at 820 MHz, so we decided to use the 685 MHz and 1400 MHz data, respectively, from Parkes observations. For those observations, we used T_{sys} = 45 K and G = 0.66 K Jy^{−1} at 685 MHz, and T_{sys} = 22 K and G = 0.74 K Jy^{−1} at 1400 MHz. Figure 2 shows the pulseaveraged flux density of pulsar B, averaged over BP1, BP2, the IP, and the total orbit, as a function of MJD. In Fig. 3, the pulseaveraged flux density as a function of orbital phase and MJD is presented in a greyscale plot.
Fig. 2. Pulseaveraged flux density of pulsar B at 820 MHz (filled circles), as a function of MJD, averaged over BP1 (red), BP2 (blue), the IP (green), and over the entire orbit (black). The fluxdensity calculation of the individual pulse profiles is based on the radiometer equation, assuming a duty cycle equal to the full pulse width at 10% of the pulse maximum. For the intervals MJD53000–53100 and MJD53100–53200, the flux densities (empty circles) were scaled to 820 MHz, using the spectral indices reported in Sect. 3.3, from observations at 685 and 1400 MHz, respectively. In the interval MJD54200–54300, there were no available data. The dashed lines have been used as guides only. 
Fig. 3. Greyscale map of pulseaveraged flux density of pulsar B, as in Fig. 2, as a function of orbital phase and MJD. Bilinear interpolation between bins has been applied. The red, blue and green dashed lines delineate the orbitalphase regions of BP1, BP2, and IP as a function of MJD, assuming the corresponding 3σ widths shown in Table 2. We note that the IP is undetectable before MJD53100 and after MJD53900. The WP is defined as the orbitalphase region between the upper bound of IP and lower bound of BP1. 
3.1. Timing analysis
In KSM+06, the BPs were excluded from the timing analysis of pulsar B, as the pulse profile evolves dramatically as a function of orbital phase, during each of the orbitalphase intervals corresponding to BP1 and BP2. On the other hand, although much weaker, the profile during the WP appeared to be more stable. Based on those findings, we created an analytic template by fitting a threeGaussiancomponent model to the average WP profile of pulsar B (see Fig. 4): the latter was generated by averaging all available WP data in our data set, over the range of epochs where the WP emission was detectable. More specifically, we averaged all 820 MHz BCPM data that were taken between MJD53200 and MJD53700, noting that although the WP was also detectable during the 1400 MHz Parkes observations, between MJD53100 and MJD53200, we wanted to avoid mixing profiles at different frequencies, causing possible profile smearing due to frequency evolution. Subsequently, we generated the pulse times of arrival (TOAs) by crosscorrelating the WP template profile of Fig. 4 with the data profiles of all subintegrations. The crosscorrelation was performed within PSRCHIVE using the Fourier domain fitting algorithm of Taylor (1992) and Markov chain Monte Carlo (MCMC) sampling to estimate the uncertainties. We note here that in most observations the S/N of the subintegrations during the IP and the WP is very low. In order to obtain useful fits of the WP template to the profiles of those phases, we increased the S/N by averaging over a number of subintegrations depending on the amount of signal in the observation. At that point, many of the TOAs still came from random fits to noise and are unwanted. To further excise data that are unrelated to the pulsar’s emission, we rejected all TOAs with uncertainties (σ_{TOA}) exceeding the upper 68% confidence limit of the log(σ_{TOA}) distribution: σ_{TOA} > 9.3 ms (see Fig. 5). Finally, we visually inspected all profiles corresponding to the remaining TOAs and removed those that were produced by impulsive RFI. The final data set contained 4115 TOAs.
Fig. 4. Analytic template constructed from averaging pulsar B observations during the WP (ϕ_{asc}/2π ≈ 0.15–0.53), between MJD53200 and MJD53700. 
Fig. 5. Probability density function (PDF) of the logarithm of TOA uncertainties, σ_{TOA} (in μs), from the crosscorrelation of all of the original data profiles (after RFI excision) with the WP template profile of Fig. 4. The solid, dashed, and dasheddotted vertical lines correspond to the median, the 68%, and 95% confidence intervals of the distribution, respectively. The unshaded area of the histogram corresponds to the range of uncertainties (σ_{TOA} > 9.3 ms) corresponding to TOAs that were excluded from further analysis. 
Using the spin parameters for pulsar B, published by KSM+06, and the orbital parameters, calculated from the ephemeris of pulsar A (Kramer et al., in prep.) and assuming GR, we used TEMPO2 (Hobbs et al. 2006) and the DDGR model (Taylor & Weisberg 1989) to calculate the residuals of pulsar B for each MJD bin. We would like to clarify at this point that although it was also possible to use the orbital parameters from KSM+06, which were derived from directly timing pulsar B, this would have introduced systematics due to the noisy timing behaviour of this pulsar – which is the subject of this paper’s investigation. Using the DDGR model and pulsar A’s ephemeris provides much higher precision and avoids the systematics of pulsar B. In particular, the DDGR model requires only the Keplerian parameters and the two pulsar masses (determined from pulsar A’s timing to a high precision) to calculate all the relativistic effects affecting pulsar B’s timing. Equations (1)–(5), below, show the definitions of the ephemeris parameters we used to calculate the orbital parameters of pulsar B from those of pulsar A.
where is the solarmass parameter, expressed in units of time^{1}; n = 2π/P_{b} is the orbital frequency; M_{A}, M_{B}, and M are the inertial masses of pulsar A, pulsar B, and the binary system, respectively; x_{A} is the projected semimajor axis of pulsar A’s orbit; and ω_{A} and ω_{B} are the longitudes of periastron of pulsar A and pulsar B, respectively.
We must clarify here that, although we have assumed GR to calculate x_{B} from the precession of periastron and the Shapiro shape, the result actually covers all fully conservative gravity theories where the generalised Eddington parameters are sufficiently close to the values assumed by GR (cf. Sect. III B in Damour & Taylor 1992, where we actually only require ϵ − ξ/2 to be sufficiently close to the GR value). In particular, the modelling of the beam shape of pulsar B and the wind of pulsar A, which we present in Sect. 4.2, is fairly independent of the choice of a gravity theory, as the precision we require for x_{B} in order to calculate orbital phases would be sufficiently good, even with the measured value of this parameter. Later on in this paper (Sect. 5.1), where we discuss the improvements on the precision measurements of R that our modelling could bring about, the precisely calculated value of x_{B} becomes more important. Therefore, we must recognise that, because the timing performed in this paper relies upon the GR value of x_{B}, our results cannot directly be used to test theories of gravity. However, since x_{B} is the only theorydependent timing parameter in our analysis, and since as was mentioned above, the assumed values in our timing analysis are consistent across a range of gravity theories, so we are confident that our modelling and the conclusions that stem from it can contribute to future timing observations of pulsar B upon its reappearance (see discussion in Sect. 6.2).
Besides the timing corrections of the above ephemeris, no further fits for any of the pulsar parameters were attempted. The residuals in each of the 14 100day bins are shown in Fig. A.1. A visual inspection of the residuals reveals two types of delay as a function of orbital phase: (a) a slow, harmonic delay across the entire orbit, with an amplitude of ∼10 ms, and (b) a fast and, to firstorder, approximately linear delay across each of BP1, BP2, and the IP, of the same order of magnitude as (a). In addition, a secular increase of the amplitude of the above delays can be seen on a timescale of years, which is most likely caused by geodetic precession (e.g. PMK+10).
One possible source of harmonic variations in pulsar timing is the beam aberration due to the pulsar rotation, which modifies the intrinsic direction of the emitted radiation, as seen by the inertial observer. The amplitude of this effect for pulsar B is ∼(P_{B}/P_{b})x_{B} ≲ 0.5 ms, where x_{B} is the projected semimajor axis of pulsar B’s orbit (Damour & Taylor 1992); therefore, it cannot account for the observed delays, which have larger amplitudes of 1–2 orders of magnitude. In Sect. 4, we propose that the presence of the harmonic delay is mainly due to the external action of pulsar A’s wind, deflecting pulsar B’s emission beam relative to our LOS, and thus modulating the observable part of the pulsar’s emission. In this work, we mainly attempt to model that slow, harmonic modulation across the orbit. The fast profile changes across the BPs and the IP are briefly discussed in Sect. 5.
3.2. Pulse profile analysis
3.2.1. Average profiles
The strong profile evolution of pulsar B across its orbit was already noted by Lyne et al. (2004), for example. In Sect. 4, we try to draw conclusions about the structure of pulsar B’s emission by modelling this evolution. However, the large number of profiles in the original data would render our modelling approach computationally very expensive. In order to reduce the complexity of the problem, we have decided to generate average profiles across each of the four orbitalphase windows, for every MJD bin.
The significant profileshape modulation of pulsar B’s emission as a function of orbital phase and MJD limits the integration time that can be used to obtain average profiles before averaging smears out any intrinsic features of the pulsed emission. The modelling performed in this paper relies on mapping the evolution of those features as a function of orbital phase and epoch. Therefore, we have tried to limit the amount of smearing caused by averaging over those two parameters. The original data profiles used in this work come from subintegrations with a typical length of t_{sub} ≈ 7P_{B} (≈0.002P_{b}), while only a small subset of profiles (before MJD53400) have t_{sub} ≈ 65P_{B} (≈0.02P_{b}: see column 8 of Table 1). Such short integrations are not enough to obtain a stable profile, as typically a few hundred to a few thousand pulses are required for this (Lorimer & Kramer 2005). In order to obtain more stable average profiles in each MJD bin, we have further averaged the original data in orbitalphase bins of Δϕ_{asc}/2π = 0.02. The number of profiles in the original data, contained in a given combination of MJD and orbitalphase bins, varied between 1 and 210.
In Fig. 6, we show the evolution of the average profiles across BP1, BP2, and the IP, for a number of orbitalphase bins and for two MJD bins: MJD53400–53500 (top panels) and MJD53700–53800 (bottom panels). All profiles shown come from folding and averaging the original data with the timing model presented in the previous section. Since the main purpose of this figure is to highlight the relative average profile changes across each of BP1, BP2, and the IP, each group of profiles, corresponding to a combination of MJD bin and an orbitalphase window, has been equally rotated in phase, such that the peak flux density of the centroid profile – that is, the profile corresponding to the orbitalphase bin that includes the centroid phase, as is tabulated in Table 2 – occurs at pulse phase 0.5 (vertical dashed lines). In that figure, the reader can also compare the analytic WP template of Fig. 4 (red profiles), which was constructed by fitting a threeGaussiancomponent model to the average WP profile of pulsar B, with the observed average profiles. The corresponding number of single pulses averaged to obtain those profiles is also shown. In the original data, after the RFI excision and the application of our selection criteria, the fraction of subintegrations, during the WP, which contained a significant signal (> 5σ) amounts to less than 1% of the total data set. In addition, during the WP, pulsar B shows little profile evolution. Therefore, for each MJD bin we averaged all available WP data to produce a single, average WP profile per MJD bin.
Fig. 6. Profile evolution of pulsar B at 820 MHz, during BP1, BP2, and the IP, corresponding to two MJD bins: (top half) MJD53400–53500 and (bottom half) MJD53700−53800. Each profile comes from averaging the original data in the corresponding MJD bin and an orbitalphase interval of Δϕ_{asc}/2π = 0.02, centred around the orbital phase, ϕ_{asc}, shown in each panel. The number of single pulses averaged to obtain each profile are also shown at the topleft corner of each panel. The grey, vertical dashed lines at pulse phase 0.5 are used here to guide the eye and to highlight the relative changes between profiles observed during most BPs. We note that the profile alignment in this figure can only be used to compare the profiles of a single MJD bin and orbitalphase window (see main text for details). The blue, vertical error bar next to each profile corresponds to the quadrature sum of the offpulse RMS of the profile shown, σ_{stat}, and the average RMS due to systematic pulsetopulse variations of the profiles that were summed to obtain the profile shown (see Sect. 3.2.2). For comparison, we have overlaid the analytic WP template profile of Fig. 4 with the observed average profiles corresponding roughly to the middle of BP1, BP2, and the IP. 
3.2.2. Profile uncertainties
The modelling we perform later in this paper is based on the minimisation of the χ^{2} between the model and the observed profiles. As such, our method depends on both the magnitude and the uncertainty of each profile’s fluxdensity values. In particular, it is important to consider both the statistical uncertainties (σ_{stat}) arising from radiometer noise, as well as the systematic uncertainties (σ_{sys}) borne from pulsetopulse variations within the integration length of our average profiles.
Although the offpulse RMS in the average profiles (=σ_{stat}) is a good indicator of the radiometer noise present in our observations, it does not contain any information about the pulsetopulse variability across the set of profiles that were averaged together to produce the final profile. Hence, we have decided to make an estimate of the systematic uncertainties corresponding to the flux density of each phase bin, I_{i}, in each average profile, by calculating the RMS of I_{i} across the n_{prof} profiles that were averaged together. The final RMS of I_{i}, including both statistical (σ_{stat}) and systematic (σ_{sys(i)}) uncertainties is then given by
where
and ⟨I_{i}⟩ is the average flux density of the ith phase bin, across n_{prof} profiles.
Figure B.2 shows a grid of the 41 observed average profiles (black lines), corresponding to the centroids of BP1, BP2, and the IP, and to the WP. The average, centroid profiles of BP1, BP2, and the IP were generated from all originaldata profiles contained within the orbitalphase bins that included , , and , respectively, according to Table 2. The orbitalphase bins from which the above average profiles were calculated are indicated with blue horizontal error bars in Fig. A.1. Each column in the grid corresponds to BP1, BP2, the IP, or the WP, and each row corresponds to a 100day MJD bin. For each profile, the figure also shows the 1σ_{sys} fluxdensity envelope via grey lines. The 1σ_{stat} value for each profile is shown with a blue vertical error bar.
3.2.3. Profile templates
As a further step towards characterising the average centroidprofile evolution, as a function of orbital and precession phase, we fitted one or twocomponent Gaussian templates, depending on the complexity of the profile, to the observed average profiles of Fig. B.2. We hence represented the noisy observed profiles with the smooth analytical versions shown in Fig. B.3 (red lines), thus filtering out highfrequency noise and offpulse artefacts. The description of pulsar B’s profile evolution via noiseless templates is particularly advantageous in our modelling of Sect. 4, wherein we achieve modelparameter convergence by using such templates instead of the observed average profiles. The fitting of the templates was performed using PSRCHIVE, and, at this stage, considered only the statistical uncertainties (σ_{stat}). Across all 41 profiles, the best fit Gaussian templates were a good match to the observed average profiles, with the total of the summed differences between them corresponding to . Indeed, the construction of those templates, and the decision, for example, of how many Gaussian components are needed to describe the data, may still have been influenced by the presence of nonGaussian noise. However, based on the of the individual fits, we estimate that for most profiles the contribution of such systematics is of the order of , which is at least an order of magnitude smaller than the average σ_{sys} of the profiles. The exception is the IP profiles at MJD53600–53700 and MJD53700–53800 and the BP1 profile at MJD54300–54400, which come from a single profile in the original data and thus have zero systematic noise: these are profiles with very low S/N and it is difficult to be confident about the reliability of the fit.
3.2.4. Profile evolution
Using the Gaussian templates derived in the previous subsection (Fig. B.3), we calculated the peak flux density, the pulse width at 10% maximum, W_{10}, and the peak separation between the leading and trailing components (Δ_{P2P}) as a function of MJD, for each orbitalphase window. The evolution of these parameters across our data set is shown in Fig. 7. It is interesting to note that only the average centroidprofiles of BP1 show a clear fluxdensity decrease with increasing MJD (see Fig. 7a). The W_{10} values for all orbitalphase windows do not exhibit any clear trends as a function of time over four years. More precisely, over that time interval, we have σ_{W10}/⟨W_{10}⟩ ≈ 0.13 − 0.36, where σ_{W10} is the RMS of the width over the data span and ⟨W_{10}⟩ is the mean value. Locally, the exception is the evolution of the BP1 profiles, during the interval MJD53300–53900. During that time, we observe an evolution from a profile with a weak leading component and bright trailing one (ca. MJD53300–53600) – and with a significant decrease of Δ_{P2P} to roughly zero, during ca. MJD53400–53600 – to a profile with a bright leading component and weak trailing one, with Δ_{P2P} monotonically increasing thereafter, for the remainder of the data span. This exchange of the relative position of the brightest component is represented in Fig. 7c with a change in the sign of Δ_{P2P}.
Fig. 7. (a) Peak flux density of pulsar B at 820 MHz (filled circles), as a function of MJD, calculated from the average pulse profile corresponding to the centre of BP1 (red), BP2 (blue), IP (green) and to the WP (grey). For the intervals MJD53000–53100 and MJD53100–53200, the flux densities (empty circles) were scaled to 820 MHz from observations at 685 MHz and 1400 MHz, respectively, using the spectral indices reported in Sect. 3.3. In the interval MJD54200–54300, there were no available data. The dashed lines show the best fit light curves of the peak flux density of the model. (b) Full pulse width at 10% of the maximum (W_{10}) as a function of MJD, for the average profiles corresponding roughly to the centre of the orbitalphase range of BP1 (red), BP2 (blue), the IP (green), and the WP (grey). (c) Peak separation (Δ_{P2P}) between the leading and trailing component of the average profiles. Open circles indicate that the brightest component is trailing. The dashdotted grey line is the best linear fit to the BP1 data from this work (see Sect. 3.2); the solid and dotted black lines are the best linear fits to the BP1 and BP2 data, respectively, across the corresponding ranges shown by Perera et al. (2010). All values shown have been calculated from the best fit Gaussian templates of Fig. B.3. 
Moreover, PMK+10 measured the average rate of change of Δ_{P2P} for BP1 and BP2 during the interval MJD53900–54500. Their measured rates for BP1 and BP2 were yr^{−1} and yr^{−1}, respectively. We can compare those rates with our measurements based on the average profiles at the centroids of BP1 and BP2 for a similar range of epochs. Our calculations considered only the absolute component separation (Δ_{P2P}) thus ignoring their relative intensity. We find that for BP1, for the above range of epochs, yr^{−1}. For BP2, we could not detect a significant leading component beyond MJD54400; restricting the epoch range to MJD53800–54400, we estimated that yr^{−1}. Furthermore, over the entire data span, we find that BP1 shows the most systematic change of Δ_{P2P}. In particular, if we exclude MJD53000–53100 and MJD53100–53200, where Δ_{P2P} = 0, and MJD53500–53600, where the second component is marginally detected, a linear fit yields yr^{−1}. The fits by PMK+10 and our fit to the BP1 values are shown in Fig. 7c. The difference between our value of the global gradient of Δ_{P2P}, for BP1, being roughly half of that published by PMK+10, is most likely the result of the different MJD ranges considered.
3.2.5. Beam morphology
According to the geometry of BKK+08, the direction of geodetic precession leads to an increasing separation between our LOS and the visible magnetic pole of pulsar B over the span of our data. In other words, the absolute value of the impact angle (β) increases with time (β itself becomes more negative). The approximate amount by which this happens is Δβ=β(MJD 54500)−β(MJD 53000) ≈  − 15° −(−2° ) = 13° (see also Fig. 17b). This fact, together with the evolution of the average BP1 and BP2 profiles, from mostly singlepeaked profiles (before MJD53700) to predominantly doublepeaked profiles (after MJD53800), suggests a convex emission beam, where the separation of the active regions increases with the distance from the magnetic axis. Therefore, at least based on these simple observations, the emission beam does not appear to be consistent with a concentric ring or a wedge, centred on the magnetic axis. In particular, it contradicts the proposed horseshoe model by LL14, that is, a concave wedge centred on the magnetic axis, which would result in a profile evolution from a double to a singlepeaked profile, based on the geometry of BKK+08. Such a geometry would also lead to a decrease of the overall profile width with time.
A simple calculation can be made in support of the above statements. Assuming that at the start of our data the LOS trace was tangent to an emission region with a circularring geometry of radius R – where the circular ring has a negligible thickness compared to R – then, according to Fig. 7c, at the end of the data the peak separation would reach a maximum of Δ_{P2P} ≈ 8°. The radius R is then simply given by , where Δβ is the absolute value of the change of β across our data, as was defined above. For such a circular beam geometry, the rate of change of the pulse width with respect to time would be d(Δ_{P2P})/dt = 4(R − Δβ)(dβ/dt)/Δ_{P2P}. The average rate is ⟨d(Δ_{P2P})/dt⟩ = 4(R − Δβ)(Δβ/T)/Δ_{P2P}, where T ≈ 4 yr is the length of our data. Using the values corresponding to the end of our data set, the rate is yr^{−1}. Compared to the observed evolution of the peaktopeak separation of BP1 across our data set, we find that this rate is over three times higher, suggesting that the emission region must indeed be elongated rather than circular. The exact shape of pulsar B’s beam is investigated in Sect. 4.2.
3.3. Spectral index
The modelling we perform in Sect. 4 partly depends on the fluxdensity evolution of pulsar B’s profile. Therefore, in order to have comparable flux densities across our data set, we must scale the flux density of our multifrequency data set to a common reference frequency. To achieve that, we need to estimate the spectral index of pulsar B’s fluxdensity spectrum. For this estimate, we have used the peak flux density of multifrequency, fluxcalibrated, average profiles at 685, 1400, and 3100 MHz. We note that the 3100 MHz data were only available during MJD53000–53100 and came from observations with Parkes and the WBCORR backend, using 768 MHz of bandwidth; these data were calibrated using T_{sys} = 28 K and G = 0.62 K Jy^{−1} (see KSM+06). As is shown in Fig. 3, the radio flux is significantly modulated both across the orbit and during the BPs, and, due to precession, it diminishes with increasing MJD. Hence, for the spectralindex fitting, we required an overlap in orbital phase and MJD. The optimal data combinations, for which we deemed this calculation to be reliable, contained two frequencies, had an orbital phase overlap with Δϕ_{asc}/2π < 0.01, and their epochs were separated by 30 days. We calculated the spectral index, p, assuming a single power law. The results were as follows: (a) from the combination of the 685 MHz data at MJD53005 and ϕ_{asc}/2π = 0.6058 with the 1400 MHz data at MJD53033 and ϕ_{asc}/2π = 0.6068, the spectral index was p = −1.65(18); (b) from the combination of the 685 MHz data at MJD53005 and ϕ_{asc}/2π = 0.5658 with the 3100 MHz data at MJD53034 and ϕ_{asc}/2π = 0.5721, the spectral index was p = −1.39(6). The probability distributions of p from the individual data combinations and the joint probability distribution are shown in Fig. 8. The median value of the joint distribution was and was used to scale the average profiles at 685 and 1400 MHz to their corresponding values at 820 MHz. All profiles at those two frequencies presented herein have thus been scaled.
Fig. 8. Probability distributions of the spectral index of pulsar B’s radio emission, calculated assuming a single powerlaw between 685 and 1400 MHz (light grey shade) and 685 and 3100 MHz (dark grey shade), from multifrequency profiles corresponding to roughly the same orbital phase and separated in time by roughly 30 days. The joint probability distribution is shown with a solid black line, and the corresponding median spectral index and 68% confidence interval, with solid and dashed vertical lines, respectively. 
4. Simulating the effects of pulsar A’s wind
4.1. Toy model description of the harmonic delays
We explored simple linear relations between the observed harmonic delays in the timing residuals and orbital quantities that vary harmonically as a function of orbital phase. In particular, we tested functions such as the separation between the pulsars (∝[1 + ecos(ϕ_{asc} − ω_{B})]^{−1}) and the radiation pressure of the dipole (∝[1 + ecos(ϕ_{asc} − ω_{B})]^{2}). Such relations, although varying harmonically across the orbit like our data, are out of phase with the variations in the data. In Fig. A.1, it can be seen that for the epochs where the BPs, the IP, and the WP are detected, the maximum delay occurs near ϕ_{asc} = nπ, where n = 0, 1, 2… In contrast, the pulsar separation and the radiation pressure have maxima at ϕ_{asc} = ω_{B} + nπ, with ω_{B} increasing monotonically from 260° to 326° during the fouryear span of our data due to the precession of periastron. This observation suggests that the harmonic delays are ∝cos ϕ_{asc}. Such a dependence possibly implies a relation between the timing delays and the magnitude of the skyprojected component of an orbital property.
In this work, we assume that a radial wind, directed from pulsar A to pulsar B, deflects the emission beam via its relativistic wind pressure. This assumption motivated previous studies of this system based on plasma physics (e.g. Lyutikov 2005; Lomiashvili & Lyutikov 2014). Moreover, the harmonic delays in the residuals are caused by the component of the wind that is perpendicular to the LOS. That is to say,
where we define as the righthanded orthogonal Cartesian reference system, with being coplanar with the orbital plane, being coincident with the orbital angular momentum vector, L, and pointing towards the observer. The geometry assumed in this toy model is shown in Fig. 9. In the above expression, w_{0} is a dimensionless coefficient, equal to the orbitaveraged displacement perpendicular to the LOS (for e ≪ 1), as a fraction of the radius of the emission site, r_{em}. To first order, we expect that such a beam deflection introduces a phase delay (δϕ_{s}), where ϕ_{s} = Ω_{B}t is the spin phase of pulsar B, with Ω_{B} being the instantaneous spin frequency accounting for spindown. The amplitude of this delay varies harmonically across the orbit, thus reflecting what we observe in the residuals of pulsar B when a fixed template is used. The corresponding timing delay due to this deflection is
Fig. 9. Geometry of the toy model used to fit the harmonic variation of the timing residuals in Fig. A.1. In this figure, pulsar B spins clockwise with a spin frequency of Ω_{B}, and it has a velocity of v_{B}, directed away from the LOS. The configuration shown corresponds to ϕ_{asc} = 0, where the wind of pulsar A maximally deflects the emission beam of pulsar B, perpendicularly to the LOS (direction ) by angle δϕ_{s}. The LOS direction is defined by . In this model, the emission is assumed to be generated at a distance (r_{em}) from the centre of the star. The wind of pulsar A, with magnitude equal to w_{⊥} (in units of r_{em}), deflects the emission site from position E to E′. 
We investigated the dependence of w_{0}, as a function of MJD, by fitting Eq. (9) to the residuals in each of the 14 MJD bins. Figure A.1 shows the best fit functions and their 1σ confidence interval for each MJD bin, overlaid with the timing residuals. It can be seen that the amplitude of the function increases monotonically with time. Further analysis shows that the evolution of w_{0} with time can be approximated well, within the considered MJD range, with a linear regression of the form w_{0} = a(t − t_{0}) (see Fig. 11). The slope of the regression was determined to be equal to 0.016(1) yr^{−1} and the epoch when w_{0} = 0 was determined to be MJD52852(65). The MJD at which the amplitude of the harmonic variations of the residuals becomes zero, according to the best linear fit, corresponds to an impact angle of ; furthermore, β < 0 corresponds to the configuration where the LOS at its minimum approach to the currently visible magnetic pole lies between said magnetic pole and the pulsar’s north pole. In addition, we mapped the probability distribution of a and t_{0}, given the entire data set of available residuals, using the Bayesian inference tool, MULTINEST (Feroz et al. 2009), with uniform parameter priors. For this step, we used the function of Eq. (9), but with a linear dependence of w_{0} on time, as shown above, to determine the parameter values corresponding to the maximum likelihood. The 2D joint probabilitydensity map for a and t_{0} is shown as an inset in Fig. 11. The most likely values from that analysis were a = 0.015 yr^{−1} and t_{0} = MJD 52820; the latter epoch corresponds to . The small value of β at ≈MJD 52850 implies that the discovery of pulsar B occurred when our LOS was tracing emission very near the magnetic pole.
To first order, the above results imply that the impact of the wind displaces the emission region of pulsar B relative to our LOS, with an increasing orbitaveraged magnitude, across the span of our data. This conclusion is most likely the consequence of the simplified 2D model that we used here to fit the data, leading to Eq. (9), rather than that of a wind with a timedependent orbitaveraged magnitude, w_{0}(t). In reality, in contrast to the cartoon representation of Fig. 9, the spin axis of pulsar B is nonorthogonal to the orbital plane, and geodetic precession changes the angles between the spin axis and the wind direction and the spin axis and our LOS, as a function of time. This results in a changing amount of spinphase delay (δϕ_{s}) as a function of time, as is depicted in the cartoon representation of Fig 10. This effect is calculated more accurately, as part of the 3D modelling that we perform in the following sections.
Fig. 10. Cartoon representation of spinphase delay (δϕ_{s}) produced by the impact of the wind of pulsar A, with a perpendicular component to the LOS (w_{⊥}) for two different phases of geodetic precession: at time (t_{0}) and at a later time (t_{0} + δt). In this representation, the observer’s LOS is perpendicular to the plane of the figure, when looking down onto it; the plane of the orbit is perpendicular to the orbital angular momentum vector, L; and pulsar A is located to the left of pulsar B. The precession of the spin axis, s, about L, proceeds along the path shown with a dashed, grey line, in the direction indicated by the arrow. The NS’s equator is shown with a solid, grey line, with the grey arrow indicating the spin direction; the magnetic axis, μ, rotates about the spin axis in the same direction. In both instances of time, the action of w_{⊥} causes the deflection of the emission (assumed here to be at the magnetic pole) from position E (solid black point) to E′ (open circle). The equivalent spinphase delay, δϕ_{s} at (t_{0}) and at (t_{0} + δt), is measured as the polar angle between E and E′. For the same wind magnitude (w_{⊥}), it can be seen that . 
Fig. 11. Orbitaveraged wind magnitude (w_{0}) as a function of MJD (data points) from fits of Eq. (9) to the residuals in each of the MJD bins shown in Fig. A.1. The solid grey line and the shaded area correspond to the linear fit to the data and the 1σ confidence interval of the best fit, respectively. The inset figure shows the 2D joint probabilitydensity function of the rate of change of w_{0} with time (a) and the reference MJD, t_{0} (at which w_{0} = 0), given the entire data set of residuals. The red star symbol corresponds to the location of the most likely values of a and t_{0}. The dashed line corresponds to the most probable function that describes the evolution of w_{0} with time. 
4.2. Simulation of the emission geometry
As a further step, we attempted to describe the profile evolution of pulsar B, as a function of orbital phase and epoch, with a geometric 3D model of the emission region. The intensity distribution of the emission of pulsar B, relative to the magneticaxis direction () was parametrised with a twocomponent 2D Gaussian function:
where
and
The motivation for using such a general Gaussianbeam model, instead of a more specific “horseshoe” or fanbeam shape, came primarily from the high number of degrees of freedom that this model has, even allowing for beam shapes that are similar to those mentioned above. In addition, the LOS traces of such a beam produce one or twocomponent Gaussian fluxdensity profiles, which is a close approximation of the observed profiles of pulsar B, as can be seen in Fig. B.3.
The definitions of the all the angles used in our parametrisation are depicted in Figs. C.1a,b and 12. The polar coordinates (hereafter beam coordinates) in this parametrisation are γ, which is the beam colatitude, equal to the angle between the magnetic axis and the direction of emission; and ϕ, which is the beam longitude, equal to the angle between the plane containing the pulsar’s spin and magnetic axes (fiducial plane) and the plane containing the magnetic axis and the direction of the emission. The angle ϕ – where ϕ ∈ [0, 2π) – is measured counterclockwise on the plane of the sky, from the pulsar north through to the pulsar south. The parameter I_{0ℓ} corresponds to the peak intensity of each of the two Gaussian components; Γ_{0ℓ} and Φ_{0ℓ} are the beam colatitude and longitude of the peakintensity location of the Gaussian components, respectively; σ_{Γℓ} and σ_{Φℓ} = (1 − f_{ℓ})σ_{Γℓ} are the halfbeam widths of the Gaussian components, along the direction of increasing Γ_{0ℓ} and Φ_{0ℓ}, respectively (represented by the unit vectors, and ) – where f_{ℓ} is the flatness parameter, with f_{ℓ} = 0 corresponding to a circular Gaussian. Finally, ζ_{ℓ} ∈ [0, π) is the angle of rotation about the direction of peak intensity of the Gaussian, measured anticlockwise on the plane of the sky.
Fig. 12. Schematic illustration of the emission beam of pulsar B, modelled in this work as a pair of Gaussian components, of which the intensity profile is represented here as greyscale contours. This figure shows a 2D projection of the emission beam on the plane of the sky and the definitions of the beam parameters that were used in our model to describe the location and shape of the emission (Eq. (10); compare with Fig. C.1). The white circles at the centres of the elliptical contours show the locations of maximum beam intensity. In this representation, the observer’s LOS is perpendicular to the plane of the paper (as denoted by the circled cross). The circled dot at the top of the figure shows the direction of the magnetic axis () in this projection pointing out of the plane of the paper, towards the observer. The black and grey dotted lines show two examples of the trace of the LOS across the beam, at epochs t_{1} and t_{2}(> t_{1}), respectively. 
In addition, as can be seen in Fig. C.1b, d_{ℓ} is the angle between the location of maximum intensity of each of the Gaussians, with beam coordinates (Γ_{0ℓ}, Φ_{0ℓ}), and a given location on the sky, with coordinates (γ, ϕ):
ψ_{ℓ} is the angle between and the direction of increasing d_{ℓ} (represented by the unit vector, ) measured anticlockwise on the plane of the sky.
In order to parametrise the interaction of the wind with the pulsar beam, we have defined two orthogonal Cartesian reference systems, and . C is defined such that coincides with the direction of the orbital angular momentum (L) and points towards the observer. Here, we have approximated . In reality, the orbital inclination of the double pulsar system is (Breton 2009)^{2}. C^{Ω} is defined such that coincides with the spin axis of pulsar B (Ω) and always lies in the plane defined by Ω and (hereafter fiducial plane), while . Hence, C^{Ω} is inclined with respect to C by an angle δ, meaning L ∧ Ω = δ. Figure C.1 shows the full 3D geometry assumed in our model.
The calculation of the beam intensity towards the LOS (direction ) is defined with respect to the magnetic axis (γ = ϕ = 0); it requires the transformation of , from C^{Ω} to C. The transformation from C to C^{Ω} can be described as a rotation by ϕ_{SO} around the zaxis, then a rotation by δ around the y axis, and, finally, a rotation by π − ϕ_{0} about the z axis, where ϕ_{0} = tan^{−1}(tanϕ_{SO}/cos δ): the last rotation ensures that Ω, and are coplanar, so that ϕ_{s} = 0 at the closest approach of the observer’s LOS to the magnetic pole. Hence, the magnetic axis in C is expressed as
where R_{k}(ξ) is the 3 × 3 rotation matrix that rotates a 3D vector counterclockwise by angle ξ, around direction k. The above operation is noncommutative. The precession phase, ϕ_{SO}, is a timedependent quantity, meaning ϕ_{SO} = Ω_{SO}(t − t_{0}), where Ω_{SO} is the precession rate and t_{0} is the reference epoch corresponding to ϕ_{SO} = 0.
At any instant of time, the beam coordinates, (γ_{x}, ϕ_{x}), of the trace of the LOS on the unit sphere centred on pulsar B are given by
Here, we have defined the orthogonal Cartesian system , having its origin at the location of the magnetic pole; is the unit tangent vector pointing in the direction of the star’s rotation (i.e. eastwards), and is the unit tangent vector pointing in the direction of the north pole (see Fig. C.1):
Similarly, the unit tangent vector directed towards is given by
The normal unit vector,
at the position of the magnetic pole ensures the correct sense of ϕ_{x}, as was defined at the beginning of this section.
The radial wind vector, with components
and
deflects the beam in the direction of , such that the deflected emission forms an angle with , which is given by
where
defines the direction of the deflection as anticlockwise – when viewed from the celestial north – when ϕ_{asc} ∈ [0, π/2)∪[3π/2, 2π), and clockwise elsewhere. To calculate the pulse profile after the effect of the wind, it is sufficient to only rotate by −δϕ_{x} (i.e. in the opposite direction to the deflection) and recalculate the pulse profile from Eq. (10), using the coordinates of Eq. (15).
Ultimately, for a given set of beamshape parameters, wind magnitude, and epoch, our simulation produces a fluxdensity profile: F(t; I_{0}, Γ_{0}, Φ_{0}, σ_{Γ}, f, ζ, w_{0}). We stress here that our model does not account for the rapid profile evolution that is observed during each of the BPs, and, to a certain degree, during the IP.
4.3. Simulation setup
We combined Eqs. (10) and (15) to calculate the pulse profile of pulsar B, as function of ϕ_{asc} and MJD, under the influence of pulsar A’s wind, of average magnitude w_{0}. The orientation of pulsar B’s spin and magnetic axes relative to the orbital plane, and its rate and phase of precession have been estimated with high precision in the work of BKK+08, using a model of the eclipses of pulsar A. We used those values as fixed parameters in our simulation. BKK+08 calculated the above values in a global selfconsistent fit. Although we are aware that the GR prediction for Ω_{SO} (i.e. ) yields a much more precise value, to remain consistent with the global solution across all model parameters we decided to use the fitted value for Ω_{SO}. Specifically, we fixed the values^{3} for , , and yr^{−1}. We note that we define , which is in the opposite direction to that in BKK+08; therefore, in our model, ϕ_{SO} increases with time, and, accordingly, ϕ_{SO} has the opposite value at the reference epoch.
The free parameters of our simulation were I_{{01, 02}}, Γ_{{01, 02}}, Φ_{{01, 02}}, σ_{Γ{1, 2}}, f_{{1, 2}}, ζ_{{1, 2}}, w_{0}. In order to account for unmodelled physical effects that amplify or suppress the beam intensity during the four different orbitalphase windows, BP1, BP2, IP, and WP, we assigned separate free parameters for I_{{01, 02}}, to each of those orbitalphase intervals: for example, , , , . For a beam model and wind magnitude that exactly describes the profile evolution of pulsar B, these 19 parameters should be constant and independent of the orbital phase or epoch of observation.
In LL14, the authors constrained the beamshape parameters by fitting the simulated intensity variations of pulsar B, as a function of epoch, to the observed peakintensity maps of pulsar B, published by PMK+10. The best fit parameters were determined via the computation of the correlation coefficient between the simulated and the observed intensity maps, across a multidimensional parameter grid with iteratively increasing resolution around the highest coefficient. Apart from the precession rate, which was fixed to yr^{−1}, LL14 assumed no prior knowledge of the pulsar’s orientation and also fitted for α and θ ( = 180° −δ), which they determined to be 56° and 122°, respectively. An additional parameter that LL14 determined was the height of the emission, r_{em} = 3750 R_{NS}, where R_{NS} ∼ 10 km is the NS radius. In this work, we did not consider an emission height or any associated phase shift due to the lag between the rotating frame of the NS and that of an inertial observer. Indeed, for the slowrotating pulsar B, such a phase shift only becomes significant if the emission is generated close to the light cylinder (Perera et al. 2012). We would like to note that the high value of r_{em} in the work of LL14 resulted in perpetual visibility for pulsar B, which contradicts observations. To mitigate this problem, the authors forced a cutoff on the intensity by means of a Gaussian filter function, allowing emission only from near the nullcharge surface, meaning regions for which Ω ⋅ B ∼ 0, where B is the local magnetic field at the region of interest. It is interesting to note that the horseshoe beam of LL14, when adapted to the geometry of BKK+08, leads to no emission during the period from the pulsar’s discovery to its disappearance. The possibility of emission with the beam of LL14 still exists, of course, but it requires the beam to be rotated by ≈180° about the magnetic axis. In that case, precession would evolve the pulse profile from doublepeaked to singlepeaked, which is the opposite of what is observed. Notably, in LL14, the direction of precession is reversed (see Sect. 8.3 in that paper), leading to the desired pulseshape evolution.
4.4. Parameter estimation
The significant profileshape evolution of pulsar B across its orbit means that timing with a fixed profile template, as was done to derive the timing residuals shown in Fig. A.1, will introduce systematic offsets induced by the crosscorrelation procedure: the latter tries to find the best phase offset that minimises the difference between the shape of the template and that of the data profile. If the shape of the observed profile changes significantly (see e.g. Fig. 6), the phase offset does not only represent the relative shift between the template and the data but also a systematic offset that reflects the difference in pulse width, the number and the relative amplitude of the components, etc. (see e.g. Hassall et al. 2012). As was mentioned in the introduction, to avoid this problem in timing pulsar B, KSM+06 generated a set of different templates for different epochs and orbital phases. Generated templates that are good approximations of the observed profile shapes would yield smaller systematic offsets between the model and data profiles, and the corresponding model parameters would be weighted higher. However, the templates of KSM+06 were not motivated by an underlying physical model but were constructed from the data. Our geometric model allows us to derive analytic templates for a given epoch and orbital phase, based on a coherent evolution of the viewing geometry, a parametrisation of the beam shape of pulsar B, and the impact of a radial wind. As is also discussed in Sect. 5.1, the systematic uncertainties of the crosscorrelation procedure represent only a fraction of the “true” uncertainties, since they do not account for a potential plethora of alternative geometries that could provide an equally good or better fit to the observations.
We attempted to constrain our model’s parameters, given the observed profiles of pulsar B, by mapping their probability distribution via the following likelihood function:
where and represent the flux density of the ith pulsephase bin of the observed and model profile, respectively, corresponding to orbital phase ϕ_{asc(j)} and epoch t_{k}. The model profiles are calculated from Eqs. (10) and (15), by calculating I(γ_{x}, ϕ_{x}) for every pulsephase bin. The model parameters are represented here as a vector , where J ∈ {BP1, BP2, IP, WP}, for each orbitalphase window, and ℓ ∈ {1, 2}, for each Gaussian component of the beam. Finally, σ_{ijk} are the standard deviations of the observed profiles (calculated according to Eq. (6)).
Our aim was to determine the most likely parameters that describe the profile evolution both across the orbit, due to pulsar A’s wind, and over the years, due to geodetic precession. For that reason, our model must be able to account for the profile changes that are observed across the entire span of our data. However, it is clear that the effect of the radial wind in our model can only produce a harmonic modulation of the profile shape (at the orbital period) via the periodic displacement of the beam as a function of orbital phase. As such, the fast evolution across each of the BPs and the IP cannot be accounted for, in our simple model. Therefore, we decided to only use the average profiles nearest to the centre of BP1, BP2, and the IP (see Fig. B.3). As was mentioned earlier, when the pulsar is detectable during the WP, the average profile was generated from all the available data in the corresponding orbitalphase window. Furthermore, as was explained in Sect. 3.2, to avoid biasing the parameter estimation due to offpulse artefacts and other nonGaussian noise unrelated to the pulsar’s emission, we used the noiseless Gaussian templates of Fig. B.3 as instead of the average data profiles. Our choice to use noiseless templates instead of the observed profiles leads to overestimated values of ℒ, because the magnitudes of the statistical and systematic noise, σ_{ijk}, which are considered in Eq. (23), are not reflected in the differences, , leading to smaller values of the arguments being summed. In terms of the statistical noise, this decision has little impact, as σ_{stat} is roughly the same across our data set, with mJy across the 41 profiles. More significant is the systematic noise, σ_{sys}, which is both higher in the majority of cases and varies significantly between profiles: for example, by excluding all pulse phases outside the FWHM of the pulse profiles, and ignoring phase bins with σ_{sys} < ⟨σ_{stat}⟩, where ⟨σ_{stat}⟩ is the median of the statistical noise, as shown above, we calculate mJy. However, the contribution of the systematic noise could only be properly accounted for in Eq. (23), if all of the 1467 original data profiles that were averaged to produce the final set of 41 averaged profiles were used in the calculation of ℒ. Unfortunately, this would have been prohibitively expensive in terms of computation.
The sampling of the parameter space was done using POLYCHORD (Handley et al. 2015), which is a nested sampling algorithm that is tailored for problems with high dimensionality. For a given set of modelparameter values, the code calculates the likelihood of Eq. (23); it then converges towards the most likely region of the parameter space by means of sequentially sampling that space with 500 of socalled live points. These are updated sequentially towards increasingly constraining regions around the global maximum likelihood.
In addition to the 19 model parameters, we also constrained the global constant phase offset (δϕ_{s0}) between the model profiles and the corresponding observed profiles: this was necessary because the reference phase of the observed profiles, which is typically defined by the start of the observations and the timing ephemeris, does not necessarily match the reference phase of the model. Moreover, we used a more conservative approach and introduced separate δϕ_{s0} parameters for the 685, 820, and 1400 MHz profiles: , , and . This decision was motivated by the frequencydependent phase offsets that can arise from differing instrumentation and/or observatories.
Finally, to account for the unknown amount of covariance between σ_{sys} and σ_{stat}, for each of the 41 average profiles we introduced an additional error coefficient, q_{kj}, where . In total, the number of free parameters that were simultaneously constrained given the 41 observed profiles was . For all parameters, we chose uniform priors: the corresponding ranges are shown in the last column of Table 3. For q_{kj}, the range of the priors was [0,10].
Mean and most likely values from our analysis, for the beam intensity, beamshape, and average wind magnitude.
4.5. Results
4.5.1. Parameter distributions
In Fig. 13, we show a small subset of the PDFs that were generated by our POLYCHORD run, focusing on the shape and location parameters of the brightest Gaussian component of the BP1 beam, Γ_{02}, Φ_{02}, f_{2}, and ζ_{2}, as well as the wind parameter, w_{0}. The complete set of PDFs can be found in Figs. C.2 and C.3. The values corresponding to the maximum global likelihood are indicated with a star symbol, in each of the panels. It is clear from these plots that certain parameter combinations, like ζ_{2} and Φ_{02}, are highly covariant, whereas others, like w_{0} and Φ_{02}, are close to orthogonal. It can also be seen, most clearly in Fig. 13a, that for certain distributions the most likely values lie outside the 1σ contour. The reason for this is that the nonGaussianity of some of the distributions – most notably, the one shown in Fig. 13d, for example – biases the joint likelihood of certain parameter combinations. Nevertheless, these values are still the most likely, when the likelihood over the entire multidimensional space is considered. The mean values, 1σ uncertainties, and the values corresponding to the solution with the highest global likelihood are shown in Table 3. It is worth noting that the constant phase offset, δϕ_{s0}, was found to be identical for all three frequencies and equal to . Apart from the three profiles corresponding to the IP at MJD53600–53700 and MJD53700–53800 and the BP1 at MJD54300–54400, which as mentioned earlier have σ_{sys} = 0 and for which q_{kj} was unconstrained, the majority of the rest of the profiles have q_{kj} < 0.5.
Fig. 13. Probability density plots of four pairs of model parameters, corresponding to the second Gaussian beam component, from our analysis: (a) Γ_{02}–w_{0}, the beam colatitude of the peak intensity of the Gaussian component against the wind magnitude; (b) f_{2}–σ_{Γ2}, the flatness parameter against the halfbeam width along the direction of increasing Γ_{02}; (c) Φ_{02}–w_{0}, the beam longitude of the peak intensity of the Gaussian component against the wind magnitude; and (d) Φ_{02}–ζ_{2}, the beam longitude of the peak intensity of the Gaussian component against the angle of rotation about the direction of peak intensity. The star symbols indicate the values corresponding to the solution with the maximum global likelihood. 
4.5.2. Beam structure
We used the values corresponding to the most likely solution for the beam’s parameters, to provide a visual representation of pulsar B’s emission beam. Since our model parametrised the intensities of each component individually, for each orbitalphase window, the beam shape due to the different and values, for each J ∈ {BP1, BP2, IP, WP}, also differs. Using the most likely values for the beam intensities, in Fig. 14 we show the 2D projections of the beam of BP1, BP2, the IP, and the WP, where for clarity we have normalised the intensity to 1. Due to this normalisation, we stress that these plots show only the relative intensity between the Gaussian beam components in each orbitalphase window and should not be used to compare intensities between orbital phases. The trace of the observer’s LOS at MJD53500, MJD53750, and MJD54500, for an unperturbed beam (w_{0} = 0) and for w_{0} ≈ 0.017, is shown with solid and dashed white lines, respectively. The beam intensity maps reveal that in BP1 and the IP, the most likely beam consists of a dominant primary component that, according to the values of I_{0ℓ}, is roughly 20 and 100 times brighter, respectively, than the secondary; the absolute intensity of the dominant component in BP1 is roughly ten times brighter than in the IP. In contrast, in BP2 and the WP, the intensity of the primary and secondary components is roughly equal. Except for the IP, where the secondary component is practically invisible, in all other orbitalphase windows the two components are offset with respect to each other by roughly 15° (combine Γ_{0ℓ} and Φ_{0ℓ} with Eq. (13)); also, the two components form an angle of ≈20° (combine parameters Φ_{0ℓ} and ζ_{ℓ}). Interestingly, although both the BP1 and IP beams comprise a clearly dominant component, for BP1 this corresponds to ℓ = 1, while for the IP it corresponds to ℓ = 2. In both cases, the secondary component is practically eclipsed by the brightness of the primary. For illustration purposes, Fig. 15 shows a 3D representation of the most likely beam during BP1, BP2, the IP, and the WP, rendered on a sphere, where the viewing geometry corresponds to MJD53000, MJD53750, and MJD54500.
Fig. 14. Twodimensional projection on plane of the sky of most likely emission beams of pulsar B, corresponding to BP1, BP2, the IP, and the WP, as they were determined from our analysis. The Cartesian coordinates shown (X, Y) are defined in Sect. 4.2. In each panel, the intensity has been colourcoded and normalised to a maximum value of 1 (see colour bar) to provide a clear representation of the beam shape; the relative intensity between the orbital phases is not shown. The solid white lines indicate the traces of the unperturbed LOS for MJD53000, MJD53750, and MJD54500; the dashed white lines show the same traces after applying the deflection caused by the most likely magnitude of the wind (i.e. w_{0} = 0.017). The circleddot symbol indicates the position of the magnetic axis relative to the beam. A 3D rendering of these beams is provided in Fig. 15. 
Fig. 15. (a) 3D representation of most likely emission beams during BP1, BP2, the IP, and the WP, derived from this work, shown at the minimum approach of the magnetic axis () to the LOS (). The figure shows three viewing configurations, corresponding to three geodeticprecession phases, at MJD53000, MJD53750, and MJD54500. The trace of the LOS is shown with a white line. The colour scale represents the beam intensity in arbitrary units, normalised between 0 (blue) and 1(red), in steps of 0.1. 
Fig. 15. continued. 
4.5.3. Profile evolution
An important aspect of any beam model is how closely it reproduces the observed fluxdensity profiles as a function of orbital and precessional phase. Using again the most likely beam and wind parameters, in Fig. B.4 we show the model profiles (red lines) overlaid with the observed profiles (black lines). Similarly to Fig. B.2, the 1σ confidence interval of the systematic uncertainties, multiplied by the most likely error coefficient (q_{kj}), such as , is shown for each profile, with grey lines.
Qualitatively, the model tracks the precessional evolution of the profile well, beginning with mainly singlepeaked profiles at the earliest epochs (< MJD53500), and evolving to almost exclusively doublepeaked profiles after MJD53900. Equally successful is the model’s ability to track the phase drift of the profile as a function of orbital phase, although there are a few cases, almost exclusively in the IP (e.g. MJD53400–53500 & MJD53800–53900), where the peak of the model profile is shifted towards earlier phases, by a few phase bins, relative to the observed profile. The reason behind this discrepancy could be related to the following characteristics of the IP. Firstly, its orbitalphase extent is significantly larger than that of the BP1 and BP2 (see Table 2). As a result, it is less certain that a single, average profile at the centroid of the IP would reflect the magnitude of the harmonic profile evolution at that orbital phase. Furthermore, the profile evolution during the IP, during the epochs where the largest discrepancies between the model and the data are seen, is more complex than that of BP1 and BP2. Indeed, the timing residuals during the IP, at MJD53400–53500 and MJD53500–53600, exhibit a more irregular drift pattern than the monotonic drifts observed during the BPs (see Fig. A.1). Moreover, at ϕ_{asc}/2π ≈ 0.02, a discontinuity can be seen, accompanied with a decrease in S/N and the number of detections. These timing irregularities imply a more stochastic profile evolution across the IP compared to the other orbitalphase windows, a conclusion that is partly supported by the profiles shown in Fig. 6. Lastly, as can be seen in Fig. B.1, the S/N during the IP is significantly lower than that of BP1 and BP2. As a result, fewer profiles were averaged together to create the average, centroid profiles for that orbital window, lending to less stable profiles that are possibly less representative of the average profile shape at that orbital phase. Despite those shortcomings, our work did not try to further characterise the complex structure of the IP and was restricted to considering it as a single interval, due to the limited statistics. Overall, the reduced chisquared between the model and observed profiles was (calculated from 20 992 data points and a model with 63 free parameters, i.e. using 20 929 degrees of freedom).
4.5.4. Flux evolution
As was mentioned earlier, our analysis has determined the most likely amplitudes of the Gaussian components, for each orbitalphase window. However, as well as changing the profile shape, geodetic precession also causes a change in the profiles’ amplitudes, as our LOS traces different parts of the beam. It is therefore interesting to compare the peak fluxdensity evolution of the model profile to that observed. In Fig. 7, together with the measured flux densities of the observed profiles (data points) we have plotted the peak fluxdensity of the model profiles with dashed lines, for each orbitalphase window. The roughly constant flux density of the profiles of BP2, the IP, and the WP is well reproduced by the model. More interestingly, the fading of the BP1 profiles with increasing MJD is, overall, tracked by the model. However, at around MJD53500 and after MJD53800 the model appears to systematically underestimate the BP1 flux density by as much as 2–3σ. It must be stressed, however, that this comparison does not consider the complex uncertainty of the model flux density, arising from the combination of the uncertainties of the model parameters (see Table 3). To first order, we can estimate that uncertainty, based on the 1σ uncertainty of , as being roughly 10%. Although this alone does not eliminate the inconsistencies between the observed and model flux density, in particular around MJD53100, 53400, and 53800, it reduces them to within 2σ.
4.5.5. Sense of rotation
Another important property that was determined by our parameter estimation is the sense of pulsar B’s rotation. Such a determination in binary systems could provide useful information for studies of binarystar evolution and formation of such doubleNS systems. Very recently, the sense of pulsar A’s rotation was determined by Pol et al. (2018), using the drifting features in the subpulse structure of pulsar B, caused by the electromagnetic radiation of pulsar A. The authors concluded that pulsar A spins in a prograde fashion relative to its orbital motion.
As was mentioned in the introduction, the analysis of Breton (2009) constrained the angle, δ, between the orbital angular momentum and the spin axis of pulsar B. However, the solution was not unique and there was a degeneracy between δ and 180° −δ, hence the sense of rotation was not uniquely determined. In our analysis, we tested both the retrograde and prograde model, while constraining the model parameters. In our model, reversing the spin of the pulsar affects the direction of the phase delay caused by the wind vector at a given orbital phase. For example, as is shown in Fig. 9, at ϕ_{asc} = 0 the radial wind deflects the pulsar beam counter to the pulsar’s retrograde rotation, causing a phase delay in the pulse arrival time. This is of course what we observe in the residuals if a fixed template is used. However, in reality, the wind deflects the entire emission beam relative to the LOS, causing not only a translation but also a profile evolution, as our LOS traces different parts of the beam. Hence, depending on the complexity of the emission region, both retrograde and prograde configurations could result in the observed profile evolution and must be therefore tested.
A direct comparison could be made between the logevidence values, log Z, corresponding to the retrograde and prograde configurations, where Z is the integral of the likelihood over the entire parameter space. For the retrograde case, POLYCHORD reported log Z_{↺} = 19 345.0(8), while for the prograde case the value was log Z_{↻} = 19 923.7(8). Converting the logevidence values to a probability of the retrograde solution compared to the prograde, we obtain P = R/(1 + R) ≈ 10^{−578} ≈ 0, where log R = log Z_{↺} − log Z_{↻} = −578 is the logarithm of the Bayes factor. In addition, the most likely profiles in the retrograde case were exclusively singlecomponent Gaussians, across the entire MJD range, and did not reproduce the observed profile evolution due to geodetic precession (see Fig. B.4). In particular, the total reduced chisquared between the observed profiles and the retrograde model was (cf. , for the prograde case). Based on those facts, we have concluded that the most likely sense of rotation for pulsar B, relative to its orbital motion and relative to pulsar A, is prograde.
5. Discussion
5.1. Timing improvements
The model profiles of Fig. B.4 can be used as analytic templates for timing pulsar B. We have generated TOAs for each of the original 4115 observed profiles, using the model profiles corresponding to the orbitalphase and MJD range as timing templates. Here, we recall that the model templates we used for timing pulsar B represent only a small fraction of the pulsar’s orbit and therefore we may not expect that they track the profile changes across each orbitalphase window. Nevertheless, we are interested in the model’s performance with regards to the harmonic delays that are observed across the entire orbit. Ideally, a beam model must be able to match the observed profile shapes as well as the Gaussian templates of Fig. B.3. A visual representation of such an ideal case can be seen in Fig. A.2, where we show the timing residuals derived from timing pulsar B with those Gaussian templates: compared to the fixed template, the total RMS of the residuals improves from RMS = 9.33 ms to RMS = 4 ms. This can be considered as the nearoptimal case of correcting for the harmonic delays. However, such a heuristic method is not informative as to the physical processes behind those delays: our model attempts to interpret those delays in terms of the geometry of pulsar B’s emission beam and how it is affected by pulsar A’s wind pressure.
The performance of our model can also be quantified in terms of the improvement on the RMS of the timing residuals, relative to the case where a fixed template was used to time pulsar B (see Fig. A.1). In Fig. A.3, we show the residuals per MJD bin, generated from the original observed profiles and the most likely model profiles used as templates. The corresponding RMS was 6.6 ms, this value being approximately midway between the naïve case of a fixed template and the ideal case of heuristic Gaussian templates.
An interesting comparison can be made between the timing residuals using the fixed WP template of Fig. 4, the Gaussian templates of Fig. B.3, and those using the model templates, for the MJD interval, MJD53900–54200. It is clear that the residuals corresponding to the BPs, based on the first two sets of templates, are clustered in two distinct sets separated by ≈0.01 in pulse phase. This is roughly the separation of the two distinct peaks that the BP profiles develop during that interval. Since the fixed WP template cannot account for the complexity of the observed profile during those epochs, the crosscorrelation procedure stochastically determines the phase delay based on the brightest of the two peaks, leading to two clusters. The heuristic Gaussian templates track the complexity of the BP profiles much better, leading to a less pronounced clustering (although some still remains, possibly owing to a small subset of singlepeaked profiles, for which the doublepeaked template is a bad fit). Finally, the model templates seem to track the profile shape well during BP1, resulting in practically no clustering of the residuals, in contrast to the residuals during BP2, where as expected from the fact that our model fails to reproduce the observed doublepeaked profiles the clustering is more pronounced.
An important aspect of our efforts to improve the timing of pulsar B is the corresponding improvement on the precision of the orbital parameters that play a central role in tests of GR and alternative theories of gravity. One of these parameters is the ratio of the intrinsic semimajor axes, , which is qualitatively different from the rest, as it provides theoryindependent information, which can be used to constrain a large family of theories of gravity in a generic way. Although the precision of most PK parameters used in those tests improves significantly with time, solely by the continuing timing of pulsar A (Kramer et al., in prep.), R does not: its precision is dominated by the uncertainty of the observed semimajor axis of pulsar B, , where ϵ_{aberr} is the fractional change of the semimajor axis caused by beam aberration (Damour & Taylor 1992).
Crucially, as has been mentioned earlier, the timing of pulsar B presented here depends on the precise value of x_{B}, which was calculated from the orbit of pulsar A, assuming GR. Therefore, the residuals between the GR timing model and the TOAs, whether the latter were derived using a fixed template or the templates of our model, are dependent on that assumption. It follows, then, that those residuals cannot directly be used in tests of GR and other theories of gravity. However, we would like to reiterate that our results are valid for a range of fully conservative gravity theories.
It should also be emphasised that there are a number of sources of systematic uncertainty, which we have not accounted for. Even in an ideal case where our model had been able to perfectly predict the observed profile shape as a function of orbital phase and MJD, it would still contain systematic uncertainties related to the assumption of theory and modeldependent parameters used in this work, such as x_{B}, Ω_{SO}, as well as all the parameters of pulsar B’s geometry that were adopted from BKK+08. In that hypothetical case, the RMS of the residuals would merely represent the statistical uncertainties between the observed data and the noiseless templates; the mean values of the timing parameters that one would derive from fitting a timing model to the TOAs would be consistent with those originally calculated from GR. Any additional red noise in the residuals, such as that seen using our imperfect model, would then be indication of the inability of our model to predict the exact phase and shape of the observed profiles. In this work, we have tried to account for part of that red noise using an analytical model based on geometry. As such, we cannot exclude the possibility that there are a number of alternative models that would perform equally well or better. This uncertainty mainly arises from the nature of such type of mathematical modelling, that is, one not based on physical principles derived from a deep knowledge of pulsar magnetospheric processes and how plasma winds interact with pulsar emission. Therefore, this systematic uncertainty associated with our model will always remain, even if our parametrisation is successfully applied in combination with future, independent timing of pulsar B to provide better constraints on R.
Apart from all the aforementioned shortcomings, it must also be stressed that the timing improvements that resulted from the application of our model do not automatically reflect corresponding improvements on the precision of the orbital parameters, as they do not account for the covariances between the parameters of the timing model and those of our model. To account for those covariances, a global fit including all parameters is necessary, which is beyond the scope of this paper. However, in the idealised case where there are no such covariances, we can estimate the maximal improvement expected on the precision of relative to its previous estimate by KSM+06 using all the TOAs from our analysis in TEMPO2. More specifically, we performed a fit for , while keeping the rest of the parameters fixed. The precision from the analysis of KSM+06, based on a smaller (507 TOAs) and shorter (≈2.5 yr) data set than ours, was σ_{xB} = 1.6 lightms; using our wind model and beam, we were able to reduce this value by a factor of ≈2.6, improving it to σ_{xB} = 0.61 lightms. This improvement also points towards an even more interesting prospect, the measurement of ϵ_{aberr} for pulsar B for the first time. At the moment, however, the above precision falls slightly short of its predicted value, which is lightms, ≈1.5× larger by comparison.
Finally, we briefly investigated the timing improvement that a future model could bring, if it is also able to account for the quasilinear drifts observed across BP1, BP2, and the IP. After calculating the residuals based on the Gaussian templates, we estimated the average drift across each of BP1, BP2, and the IP, assuming a linear regression, via the fits shown in Fig. A.2. The evolution of the slopes of the average drifts as a function of observation epoch is shown for each orbitalphase window in Fig. 16. We also tested a quadratic model for the drifts, but found that on the whole it was not significantly better than the linear model, given the uncertainties of our timing residuals. After correcting for the linear drifts, using the best fit functions, we once again estimated the RMS of the timing residuals. The complete heuristic correction, based on the Gaussian templates and the bestfit linear drifts, resulted in RMS = 2.5 ms. The timing residuals after those corrections can be seen in Fig. A.4.
Fig. 16. Average drift rate of timing residuals across BP1 (red), BP2 (blue), and the IP (green), as a function of MJD. The drift rates have been estimated from the linear fits shown in Fig. A.2 and are shown both in units of ms per orbit and as a fraction of pulsar B’s spin phase per orbit. 
5.2. The migration of the bright phases
An effect that is likely caused by geodetic precession is the migration of the locations of BP1 and BP2 as a function of time. However, a simple linear fit to the centroids of BP1 and BP2 shows that the average migration rate is significantly less than the precession rate (see Fig. 1). We note here that the pulseprofile evolution due to geodetic precession can only alter the phase of emission within a pulse period, that is, cause timing delays of ≲P_{B} ∼ 10^{−4}P_{b}, and it cannot explain the shifts of the BPs by ∼10P_{B} yr^{−1}. Although our data span is short compared to the period of precession, it can be seen that the migration of the centroids of BP1 and BP2 deviates from being linear with time (see values in Fig. 1) and that this deviation appears to progress in opposite directions. Indeed, towards the later epochs, the rate of change of the phase of BP1 increases, whereas at the same time that of BP2 decreases. To confirm this apparent complementarity between the migration rates, we plotted the average orbitalphase shift of BP1 and BP2 as a function of epoch: this was calculated using Table 2 as , where is the numerical average of the locations of BP1 and BP2, at the earliest epoch. As can be seen in the bottom plot of Fig. 1, the average shift can be fitted fairly well () with a linear regression with a slope of yr^{−1}, which confirms that the rates of change of the locations of BP1 and BP2 are complementary across the investigated range of epochs.
The above findings motivated us to further explore the connection between the migration rates of the BPs and the rate of geodetic precession. We would like to emphasise that the simple model employed here is entirely independent of the detailed 3D modelling of Sect. 4.
The evolution of ϕ_{SO} is linear with time (to a high level of precision). However, as we saw above, the migration rates of BP1 and BP2 are not constant across the investigated range of MJDs. It is therefore possible that the locations of BP1 and BP2, hereafter and , are not strictly proportional to ϕ_{SO}, but vary according to other defining angles of the system’s geometry. Two such angles that can possibly influence the position of the BPs are the angle between pulsar A’s wind and pulsar B’s precessing spin axis, , and the angle between pulsar B’s spin axis and the direction to the observer, . These angles are related to the phase of geodetic precession, as follows:
As can be seen in Fig. 17a, where we plot λ and ρ for ϕ_{asc}/2π = 0.58 (BP1) and ϕ_{asc}/2π = 0.8 (BP2), for a wide range of epochs, both angles are harmonic functions of time. It is thus reasonable to assume that if and are functions of λ and ρ, they will also be harmonic functions of time. Based on that assumption, we have parametrised the orbital phases of BP1 and BP2 as
Fig. 17. (a) Time evolution of λ, which is the angle between the spin axis of pulsar B and the LOS, and ρ, which is the angle between the spin axis of pulsar B and the radial wind direction, calculated at the centroids of BP1 and BP2, ca. MJD53600. Time, along the bottom horizontal axis, is expressed as days since the epoch of the earliest MJD bin (MJD53018.1). Along the top horizontal axis, we show the corresponding phase of geodetic precession during the interval considered. (b) Time evolution of the impact angle, β, for the magnetic pole that was visible until ca. 2008 (solid curve), and the opposite magnetic pole (dashed curve). (c) Best fit sinusoidal functions of the location of BP1 (red curve and data points) and BP2 (blue curve and data points) as a function of time. The light red curve is the alias of the best fit function for BP1, modulo the orbital period. The solid black curve and data points show the average of the functions and data, respectively. Finally, the vertical dotted grey lines indicate the earliest observation epoch and the subsequent epoch at which the impact angle has the same value as that at the earliest epoch; in (a) and (b), the horizontal dotted grey lines indicate the values of λ and β at those two epochs. 
Fig. 18. Joint 2D probabilitydensity map of the period (in 10^{3} days) of the sinusoidal functions that were fit to the locations of BP1 and BP2 as a function of time (see Fig. 17). The median and 1σ uncertainties of the period of geodetic precession period of pulsar B, published by Breton et al. (2008), is shown with a filled red circle with error bars. The star symbol corresponds to the most likely values of the periods from our analysis. Finally, the empty circle corresponds to the period of geodetic precession predicted by GR. 
where P_{BP1} and P_{BP2} are the periods of the harmonic movement of BP1 and BP2, respectively; t_{01} and t_{02} are the epochs when and , respectively; finally, and are the respective amplitudes of the harmonic oscillation of and .
We explored the above parameter space using nested sampling, as before, and determined the most likely values of the parameters by simultaneously maximising the likelihood of all eight parameters in Eq. (26). In Fig. 18, we show the joint probability density maps of P_{BP1} and P_{BP2} and of t_{01} and t_{02}. The most likely values of P_{BP1} and P_{BP2} are shown with a star in that plot, and they are consistent within 1σ with the value for the period of geodetic precession, published by BKK+08; red circle. It is noteworthy that the predicted value for the period of geodetic precession by GR (empty circle) is also consistent within 1σ with both the aforementioned values.
In Fig. 17c, we plot Eq. (26) overlaid with the data, over 30 000 days, centred at MJD53018.1 (the centre epoch of the earliest MJD bin), using the most likely parameter values (Table 4). Our simple model of the movement of the BPs, as can be seen in that figure, has the following implications. Firstly, as shown with a solid black line and black data points, the average shift of BP1 and BP2 has a much flatter orbitalphase evolution with precession phase compared to those of BP1 and BP2. Consequently, over short data spans compared to the period of geodetic precession, such us ours, it can be approximated with a linear function. Secondly, according to Fig. 17c, our model of the movement of the BPs suggests that when the pulsar returns to the same viewingangle configuration as that of the earliest observations, the locations of BP1 and BP2 will be at orbital phase ≈0.17 and ≈0.27, respectively. These locations are roughly half an orbit away from where they where when the pulsar was originally observed. Also, the angles ρ_{BP1} and ρ_{BP2} will also be significantly different at that time. However, without a complete description of the effect of pulsar A’s wind on pulsar B’s emission region, we cannot make any predictions as to whether these orbitalphase intervals will still be where the pulsar will appear brightest. Certainly, irrespectively of where this happens in the orbit, our LOS will intersect the beam of pulsar B at the same latitude as when it was first seen, so during that time the pulsar will again be visible. The evolution of the impact angle, , is shown in Fig. 17b, where we use grey dotted lines to indicate the epoch of the MJD bin corresponding to the earliest observations in our data and the subsequent epoch when the viewing geometry of the pulsar’s beam will return to the same configuration (i.e. to the same β; ca. MJD62530).
5.3. Emission height
According to the most likely parameter values from our analysis, the magnitude of the wind corresponds to a maximum displacement of the emission region by ≈1.7 − 2% of the emission height (Eq. (8) with w_{0} = 0.017). Equivalently, this corresponds to a maximum deflection by . In our model, we did not consider the emission height, r_{em}, as a free parameter, but rather parametrised the wind magnitude as a fractional displacement. Nevertheless, based on the above amount of deflection and the simulations of Perera et al. (2012), we can estimate r_{em}. The aforementioned work provides an analytical expression for the angular deflection of the emission, α_{defl}, as a function of the orbital phase and the emission height, r_{em}, expressed in units of the standoff radius, r_{s} (see their Eq. (34)). The latter quantity corresponds to the distance from pulsar B, where its magnetic pressure, B^{2}(r_{s})/8π, balances the dynamic pressure of pulsar A’s wind, Ė_{A}/4πc(a − r_{s})^{2}, where B(r_{s}) is the magnetic field of pulsar B at r_{s}, Ė_{A} is the spindown luminosity of pulsar A, and a ≈ x_{A} + x_{B} is the average distance between the two pulsars. In Perera et al. (2012), the standoff radius was constrained to be r_{s} ∼ 40 000 km, which is roughly equal to 30% of the undistorted lightcylinder radius, R_{LC} = cP_{B}/2π. In that work, the authors placed upper limits on r_{em} by postulating that the maximum deflection cannot exceed the width of the pulsar’s beam, meaning . Our analysis provided an estimate for α_{defl}, which can be directly converted to a value of r_{em}. Using the centroid locations of BP1 and IP, which roughly correspond to the orbital phases where the maximum deflection occurs, and the maximum deflection of , we calculate r_{em} ≈ 12 300 − 14 500 km and r_{em} ≈ 15 100 − 16 200 km, respectively, across the span of our data. In other words, the emission site is located somewhere from 30 to 40% of the standoff radius above the star.
6. Summary and future perspectives
6.1. Summary
The work presented here has made use of roughly four years of archival pulse profiles of pulsar B (PSR J0737−3039B), from Parkes and GBT observations, covering nearly the entire range of epochs from the pulsar’s discovery, to its disappearance in 2008. The main objective of this study was to make use of all available observations of this pulsar to construct a model that describes the systematic profile variations, as a function of orbital phase and observation epoch, and to investigate the level of improvement that such a model could bring about, in future tests of GR and alternative theories of gravity, with the double pulsar.
As has been noted in previous studies of the double pulsar, pulsar B is mainly detected during two narrow orbitalphase windows, “bright phase 1” (BP1; at ≈210°) and “bright phase 2” (BP2; at ≈280° ), whereas it is barely detectable elsewhere, during the so called weak phase (WP). Using available fluxcalibrated multifrequency profiles of pulsar B, observed in BP1 at nearly identical orbital phases, we estimated the spectral index of the radio emission to be . Furthermore, we determined the locations of BP1 and BP2, as a function of time, over the fouryear data span. We find that the BP1 emission is confined within ≈180° −235°, whereas the BP2 emission is within ≈260° −308°. In addition, for half of our data span, we detected an orbitalphase window of intermediate brightness, between the BPs and the WP, which we call the intermediate phase (IP): we determined its bounds to be ≈318° −360° and ≈0° −70°. Analogously to previous work, we also determined that BP1 and BP2 shift to later orbital phases as a function of time, over the span of the available observations. We measured the average rate of migration of BP1 and BP2 to be yr^{−1} and yr^{−1}, respectively, but noted that (a) the migration rate deviates significantly from being constant across the data span, and that (b) this deviation progresses in opposite directions, for BP1 and BP2. Intriguingly, we found that the migration rates of BP1 and BP2 are complementary across the data span, which we showed by calculating the numerical average of the locations of BP1 and BP2, as a function of time, and modelling it with a linear regression. The linear fit to the average locations is a better description of the data, compared to the fits to each of the BP1 and BP2 locations separately; the slope of the best fit is equal to yr^{−1}. In a further investigation, we explored the possibility that the migration rates of BP1 and BP2 are harmonic functions of time and found that the most likely periods of such functions were 74(25) and 66(16) yr, respectively, which are consistent with the published period of geodetic precession, yr, and the period of geodetic precession predicted by GR, yr.
Pulsar B exhibits dramatic profile variations on different time scales due to the interaction with the electromagnetic wind of its companion and due to geodetic precession. To characterise these variations, we generated over 4000 TOAs from the observed pulse profiles, using a synthetic template based on the pulsar’s WP emission. After correcting for all known orbital and spin delays, assuming a model based on pulsar A timing and the validity of GR, we examined the systematic delays in the timing residuals, as a function of orbital phase and epoch of observation. The timing residuals exhibit quasilinear drifts during each of BP1 and BP2, of the order of 10–20 ms. In addition, across the entire orbit, the residuals exhibit on average a harmonic variation at the period of the orbit. Also, the amplitude of this variation slowly increases with time, over the four years of data, from ≈5 ms (ca. MJD53000) to ≈30 ms (ca. MJD54500).
The present work focused on modelling solely the harmonic profile changes associated with the mean variation of the magnitude of the effect of pulsar A’s wind on pulsar B, perpendicular to the observer’s LOS. To achieve that, we employed a simple geometric model of a radial wind, directed from pulsar A to pulsar B, which deflects the emission beam of pulsar B to a varying degree, as a function of orbital phase. The assumed wind harmonically displaces pulsar B’s beam direction, such that our LOS intersects pulsar B’s emission at different magnetic latitudes, thus introducing profile variations across the orbit.
In order to reproduce the observed profile variations, in conjunction with the effect of the wind we assumed a parametric 2D beam model, consisting of two 2D Gaussian components. The peak intensity of the Gaussian components was allowed to vary between the orbitalphase windows, to account for the significant brightness changes across the orbit, which we speculated are the result of the (herein unmodelled) interaction of pulsar A’s wind with the magnetospheric plasma of pulsar B.
On the whole, our model can be defined using 19 parameters, describing the beam shape and intensity, and the wind magnitude. It can provide an analytical approximation of the observed fluxdensity profiles, as a function of orbital phase and epoch of observation. It is important to stress the main inadequacies of the model. Since it assumes a harmonically varying wind as the source of profile variations, it cannot describe the fast profile variations across the BPs and the IP; for that reason, it was deemed sufficient to represent the profile changes across the orbit using only the average profiles corresponding to the centroids of BP1, BP2, the IP, and the WP (the orbital phases corresponding to maximum S/N). Consequently, our model is informed by only a small fraction of pulsar B’s detectable emission, across its orbit. Moreover, our model cannot reproduce the lack of observed emission between the BPs, the IP and the WP. We can speculate that those gaps in the range of detectable emission are related to the interaction responsible for the fast profile variations during the BPs and the IP. For example, the timing drifts across each of these phases and the gradually diminishing S/N towards the phase bounds perhaps suggests that the beam also drifts across our LOS, becoming undetectable beyond the bounds of each phase. However, we cannot rule out that some other interaction between the wind of pulsar A and the magnetospheric plasma of pulsar B is responsible.
Given the observed average profiles corresponding to the centroids of BP1, BP2, the IP, and the WP between MJD53000 and MJD54500, we determined the most likely parameters of our model, using a nestedsampling Bayesian algorithm. Overall, the most likely template profiles from our model track the evolution of the average profile of BP1 well, from singlepeaked to doublepeaked ones, although towards later epochs the profile intensity is somewhat underestimated. Of course, it is important to take into consideration the original systematic uncertainties of the flux density (see Fig. B.2), that is, before the application of the error coefficients. In contrast, the evolution of the average BP2 profile, which clearly develops two distinct components during MJD53800–54200, is poorly tracked by our model – although the profileintensity evolution seems more consistent with the observations, compared to that of BP1. It is interesting to note, however, that during the above MJD range, the model appears to produce an additional very weak leading component, at approximately the pulse phase of the much brighter observed leading component of BP2. Lastly, the observed average profiles of the IP and the WP are much more weak and erratic, and our model was only able to approximate them with a single Gaussian – although again an additional very weak component is generated by the model whenever the observed profiles are distinctly double peaked (e.g. the WP profiles during MJD53600–53800). The inability of our model to fit the IP and WP profiles is also evident in the residuals of Fig. A.3, where timing of those orbital phases with our model templates results in systematic offsets from zero.
An interesting property that was constrained by our modelling is the shape and intensity of pulsar B’s emission beam. Our beam model was limited in that it only comprised a combination of two surface Gaussian distributions. However, the location, rotation, ellipticity, and peak intensity of those Gaussians was allowed to cover nearly all possible values of the respective parameters. We found that the most likely beam of BP1 and the IP comprises a primary component that is at least an order of magnitude brighter than the secondary Gaussian component. Interestingly, the brightness dominance of the primary and secondary Gaussian components alternates between BP1 and the IP, essentially shifting the phase of the pulse profile’s peak. As those two orbitalphase windows represent the largest amount of timing delay between any two parts of the orbit where the pulsar is detected, the corresponding phase shift between the primary components is possibly, to a certain degree, covariant with the effect of the wind. On the other hand, in BP2 and the WP, the two components seen alternating between BP1 and IP are now both present with roughly equal intensity. In hindsight, our preliminary reasoning with regards to the beam shape (see Sect. 3.2) seems justified: we indeed find that the emission beam is consistent with an elongated, elliptical region, with components that diverge as we move away from the magnetic axis.
Finally, our modelling provided tantalising evidence for a prograde rotation for pulsar B. Assuming that pulsar B spins in a prograde fashion with respect to its orbital motion yielded model profiles that were both a closer match to the observed profile evolution as a function of orbital phase and observation epoch, and at the same time demonstrated a prograde solution with significantly higher likelihood than its retrograde counterpart. If in future studies our conclusions are confirmed, then the prograde rotation of both pulsars in the double pulsar system will have interesting implications for the formation and evolutionary chain of doubleNS systems. Moreover, in population syntheses of such systems, such information may also serve as an additional constraint when modelling the gravitational waves produced by doubleNS mergers.
Our modelling is a first step towards improving the timing of pulsar B. We have estimated the improvement we can expect, if such a model is used to time pulsar B, by generating TOAs from the original 4115 observed profiles and the model profiles as timing templates. The magnitude of the improvement is roughly halfway between the best case scenario, where the harmonic delays across the orbit are completely eliminated, and the uncorrected case. More importantly, the model offers a factor 2.6 improvement on the precision of the observed size of pulsar B’s orbit, a parameter that is central in tests of GR with this system. Additionally, it brings us closer to measuring the amount of beam aberration for pulsar B for the first time, only falling short by a factor of 1.5. The timing of pulsar B presented in this work relied to some extent on the size of pulsar B’s orbit, which was calculated assuming GR, from the orbit of pulsar A, in order to account for the orbital delays. Hence, the aforementioned improvement on the measurement precision of x_{B} is not independent of this assumption. However, the value of x_{B} is the only timing parameter in our work that significantly depends on the assumed theory of gravity and, moreover, its value is consistent across a range of fully conservative theories. Nevertheless, because of this dependence on a particular set of gravity theories, our timing results cannot directly be used to test GR and other, alternative theories of gravity. We must also caution that the above timing improvements account for neither the systematic uncertainties borne from the assumption that our simple model is the only one that correctly describes the observed profile evolution, nor the covariances between the model’s parameters and those of the timing model used. As such, the above improvements can be considered as an idealised case, specific to our model. Ultimately, a more physically motivated model combined with a global fit over all the model and timing parameters is needed to estimate the true magnitude of the improvement.
Although our modelling presents a significant improvement over previous work, there are significant, unmodelled components in our timing, reflecting our lack of knowledge of pulsar B’s emission geometry and its interaction with pulsar A. However, we are confident that, within its limitations, our model reflects real physical effects, and that it is not a phenomenological exercise in absorbing the observed timing systematics. This confidence is derived from the model’s ability to reproduce longterm precessional effects, such as the amplitude of the observed harmonic delays, without introducing superfluous parametrisation, but with simple geometry.
6.2. Future perspectives
Future modelling of the profile variations observed in archival data of pulsar B, including the timing drifts observed during BP1, BP2, and the IP, will result in even higher timing precision. As an exercise of what can be expected, we combined the heuristic Gaussian templates with a linear model of the drifts during the orbitalphase windows, which resulted in roughly factor 4 improvement in the timing precision across four years. If made through physical modelling, such an improvement would reduce the uncertainty on the projected semimajor axis of pulsar B to levels comparable to or even lower than the expected value of beam aberration for this pulsar. The improved precision of x_{B} will yield an equally significant improvement on the precision of , which will enable us to place stringent theoryindependent constraints on the strongfield parameters of binary motion, as was done in Kramer & Wex (2009).
Lastly, when pulsar B inevitably becomes visible again, which based on geometric arguments^{4} should be ca. 2024 at the latest (Breton 2009), it will be possible to perform joint timing between both pulsars, A and B, and further increase the precision of the observed timing parameters. For pulsar B, increased timing precision can come from techniques that exploit the drifting subpulse features (Freire et al. 2009; Liang et al. 2014; Pol et al. 2018). Upon continuing the timing of pulsar B, after its reappearance, having useful constraints on the phase of the pulsed emission will help to bridge the gap between the periods of visibility. To that effect, our model can be used to make predictions of the phase range we expect the emission to occur, at a given epoch. If coherent timing, either side of the pulsar’s disappearance is achieved, it will provide the eclipse model of BKK+08 with a long timing baseline, which can be used in precise tests of GR via the Ω_{SO} parameter. As discussed in Sect. 5.2, when pulsar B returns, the relative orientation between the orbital angular momentum and pulsar B’s spin, as well as the orientation of the orbit with respect to our LOS, will have changed due to geodetic and periastron precession, respectively. At that time, we expect the drift pattern in pulsar B’s pulsed emission, which is caused by the impact of pulsar A’s wind, to have also changed as a result. This future perspective will provide additional information with which we can further constrain our model’s parameters.
Here, we are assuming that reappearance will occur when the impact angle (β) returns to the same value as that of the last epoch when the pulsar was still detectable (i.e. ca. 2008). Of course, an earlier reappearance is still possible if additional beam components have become active in the meantime.
Acknowledgments
This paper is partly based on data analysis that was performed on the HERCULES cluster of the Max Planck Computing & Data Facility in Garching, Germany. AN and GD acknowledge financial support by the European Research Council (ERC) for the ERC Synergy Grant BlackHoleCam. The authors would like to thank Drs. David Champion, Ralph Eatough, Ewan Barr and Nataliya Porayko for invaluable discussions and advice, during the preparation of this manuscript. Additional thanks are given by the authors to Dr. Alessandro Ridolfi for providing specialised software that greatly assisted the dataquality inspection and selection.
References
 Arons, J., Backer, D. C., Spitkovsky, A., & Kaspi, V. M. 2004, in Probing Relativistic Winds: The Case of PSR J0737–3039 A and B, eds. F. Rasio, & I. Stairs, ASP Conf. Ser., 95 [Google Scholar]
 Barker, B. M., & O’Connell, R. F. 1975, Phys. Rev. D, 12, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Breton, R. P. 2009, Ph.D. thesis, McGill University, Canada [Google Scholar]
 Breton, R. P., Kaspi, V. M., Kramer, M., et al. 2008, Science, 321, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Burgay, M., D’Amico, N., Possenti, A., et al. 2003, Nature, 426, 531 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Burgay, M., Possenti, A., Manchester, R. N., et al. 2005, ApJ, 624, L113 [CrossRef] [Google Scholar]
 Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840 [NASA ADS] [CrossRef] [Google Scholar]
 Dungey, J. W. 1961, J. Geophys. Res., 66, 1043 [CrossRef] [Google Scholar]
 Ferdman, R. D., Stairs, I. H., Kramer, M., et al. 2013, ApJ, 767, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Freire, P. C. C., Wex, N., Kramer, M., et al. 2009, MNRAS, 396, 1764 [CrossRef] [Google Scholar]
 Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015, MNRAS, 450, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Hassall, T. E., Stappers, B. W., Hessels, J. W. T., et al. 2012, A&A, 543, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655 [NASA ADS] [CrossRef] [Google Scholar]
 Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, M., & Wex, N. 2009, Classical Quantum Gravity, 26, 073001 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2006, Science, 314, 97 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Liang, Z.X., Liang, Y., & Weisberg, J. M. 2014, MNRAS, 439, 3712 [CrossRef] [Google Scholar]
 Lomiashvili, D., & Lyutikov, M. 2014, MNRAS, 441, 690 [NASA ADS] [CrossRef] [Google Scholar]
 Lorimer, D. R., & Kramer, M. 2005, Handbook of Pulsar Astronomy (Cambridge: Cambridge University Press) [Google Scholar]
 Lyne, A. G., Burgay, M., Kramer, M., et al. 2004, Science, 303, 1153 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lyutikov, M. 2005, MNRAS, 362, 1078 [CrossRef] [Google Scholar]
 McLaughlin, M. A., Kramer, M., Lyne, A. G., et al. 2004, ApJ, 613, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Perera, B. B. P., McLaughlin, M. A., Kramer, M., et al. 2010, ApJ, 721, 1193 [NASA ADS] [CrossRef] [Google Scholar]
 Perera, B. B. P., Lomiashvili, D., Gourgouliatos, K. N., McLaughlin, M. A., & Lyutikov, M. 2012, ApJ, 750, 130 [CrossRef] [Google Scholar]
 Pol, N., McLaughlin, M., Kramer, M., et al. 2018, ApJ, 853, 73 [CrossRef] [Google Scholar]
 Ransom, S., Demorest, P., Kaspi, V., Ramachandran, R., & Backer, D. 2005, in Binary Radio Pulsars, eds. F. A. Rasio, & I. H. Stairs, ASP Conf. Ser., 328, 73 [Google Scholar]
 Taylor, J. H. 1992, Philos. Trans. Royal Soc. A, 341, 117 [Google Scholar]
 Taylor, J. H., & Weisberg, J. M. 1989, ApJ, 345, 434 [NASA ADS] [CrossRef] [Google Scholar]
 Tsyganenko, N. A. 2002a, J. Geophys. Res., 107, 1179 [Google Scholar]
 Tsyganenko, N. A. 2002b, J. Geophys. Res., 107, 1176 [Google Scholar]
Appendix A: Timing residuals
The figures presented in this appendix come from the timing analysis performed in Sects. 3.1 and 5.1, as well as the modelling of Sect. 4.1. All figures show the timing residuals after subtracting a timing model of the spin and orbital delays. For figures showing the residuals after the subtraction of supplementary timing models, describing the pulsar’s profile evolution, please refer to the respective captions for details. The orbital phase is shown as a fraction of the orbit (measured from the ascending node) and it is aliased over two orbital periods for clarity. The error bars on the residuals correspond to 1σ uncertainty. The values of the weighted RMS and the corresponding reduced chisquared (), shown for each data set, were calculated with TEMPO2 and refer to the total data span. For more details, please see the respective figure captions and the aforementioned sections.
Fig. A.1. Residuals of pulsar B for 14 100day intervals, calculated using the analytic WP template of Fig. 4, constructed from fitting a threeGaussiancomponent model to the average WP profile of pulsar B. The weighted RMS of the residuals, calculated over the entire data set shown in this figure, is 9.33 ms; the corresponding value is 75. The red lines correspond to the function of Eq. (9), using the best fit values for the parameter w_{0}. The shaded envelope corresponds to the 1σ confidence interval of the above function. The horizontal blue error bars indicate the orbitalphase intervals, at the location of peak brightness of BP1, BP2, and IP (centroids), that were used to calculate the average profiles of Fig. B.2; for the WP, the entire width of the orbitalphase window was used. Also, the panel corresponding to the interval MJD53400–54500 shows the extents of the orbitalphase regions, BP1, BP2, the IP, and the WP, for that interval, as red, blue, green, and grey shaded areas. 
Fig. A.2. Timing residuals using fixed WP template profile of Fig. 4 (black points) and using the Gaussian templates of Fig. B.3 (red open circles). For each data set, the weighted RMS and the corresponding , calculated over the entire data set shown in this figure, is shown in the topright corner of the topleft subplot. The blue lines show the best linear fits to the residuals of BP1, BP2, and IP, which were used to correct for the linear drifts across those orbitalphase windows (see Fig. A.4). The orbitalphase extent of the lines is equal to the W_{3σ} value of the corresponding orbitalphase window. 
Fig. A.3. Timing residuals using fixed WP template profile of Fig. 4 (black points) and the model templates of Fig. B.4 (red open circles). For each data set, the weighted RMS and the corresponding , calculated over the entire data set shown in this figure, is shown in the topright corner of the topleft subplot. 
Fig. A.4. (black points) Timing residuals using Gaussian templates of Fig. B.3 and after correcting for the linear drifts across each orbitalphase window, via the best fits shown in Fig. A.2. For comparison, light grey circles show the timing residuals using the fixed WP template profile of Fig. 4. For each data set, the weighted RMS and the corresponding , calculated over the entire data set shown in this figure, is shown in the topright corner of the topleft subplot. 
Appendix B: Profile evolution
The figures presented in this appendix show the positional and the average profile evolution of the BPs, the IP, and the WP of pulsar B, corresponding to the centroid of each orbital window (the orbital phase at peak brightness), as a function of time, binned in 14 100day bins. The details of the corresponding analyses can be found in Sects. 3, 4.2, and 4.4.
Fig. B.1. Distributions of S/N across the orbit of pulsar B, for each MJD bin (data points). For this figure, the values of the S/N have been normalised by the square root of the integration time. The locations of the two BPs, BP1 and BP2, and where visible that of the IP have been determined using Gaussian fits to the data (black curves). The vertical, solid, dashed and dot–dashed blue lines show the 3σ confidence intervals of the phase locations of BP1, BP2 and the IP, respectively. To emphasise the shift of the locations of BP1 and BP2, as a function of MJD, the dashed and dotted, vertical, grey lines mark the positions of the maxima of the fitted functions for BP1 and BP2, respectively, during MJD53000–53100. 
Fig. B.2. Observed average profiles of pulsar B (black lines), corresponding to the centroids of BP1, BP2, and the IP; the WP profiles were calculated using all the data in that orbitalphase window (see blue error bars in Fig. A.1). Each column corresponds to a different orbitalphase window (labelled along the top edge of the figure); each row corresponds to a different MJD bin (labelled along the right edge of the figure). The grey lines show the lower and upper 1σ_{sys} confidence limits (Eq. (7)). The blue vertical error bar next to each profile corresponds to the offpulse RMS, σ_{stat}, of the average profile. 
Fig. B.2. continued. 
Fig. B.3. Comparison between observed average profiles of pulsar B (black) and the best fit Gaussian templates (red), as described in Sect. 3.2.3. The profiles shown have been tabulated according to orbitalphase window (columns): BP1, BP2, IP, and WP, and the MJD range (rows). The reduced chisquared value shown for each plot has been calculated using only the offpulse RMS (σ_{stat}) as weights. 
Fig. B.3. continued. 
Fig. B.4. Comparison between observed profiles (black) and those produced by our simulation, as described in Sect. 4, using the model parameters with the highest likelihood and a prograde (PG) pulsar spin (red). Green lines show the most likely profiles of an alternative model with a retrograde (RG) pulsar spin. The profiles shown have been tabulated according to orbital phase window, that is, BP1, BP2, IP, and WP, and the MJD range to which they belong. The grey lines show the lower and upper 1σ_{sys} confidence limits of the observed flux density, multiplied by the most likely error coefficients (q_{kj}) that were derived from our analysis. 
Fig. B.4. continued. 
Appendix C: 3D geometry
The following figures accompany the description of our model’s geometry, presented in Sect. 4.2, as well as the description of the most likely beam shape of pulsar B, presented in Sect. 4.5.
Fig. C.1. (a) Geometry of the model used in this work to describe the effect of the radial wind of pulsar A (w) on the emission of pulsar B. In this figure, Ω and are the spin and magnetic moment of pulsar B, respectively, and L is the orbital angular momentum of the binary system. For clarity, we separated the elements corresponding to the frame of the orbit () and those corresponding to the frame of the pulsar () in black and red, respectively. The sense of pulsar B’s orbit around pulsar A and that of μ about Ω are shown with black and red arrows, respectively. Also shown with a black arrow is the sense of the geodetic precession of Ω about L (defined by an increasing ϕ_{SO}). We note that, compared to Breton et al. (2008), we used the opposite direction for L. The direction to the observer coincides with the unit vector . In this work, we have assumed that the double pulsar system is viewed exactly edgeon: so, the orbital inclination, i, is exactly equal to 90°. At the position of , we show an example location in the beam, with beam coordinates (γ, ϕ), corresponding to intensity I(γ, ϕ) (see Eq. (10)). For the definition of the reference frames and angles shown, the reader is directed to Sect. 4.2. (b) The definitions of the angles that are used in our parametrisation of pulsar B’s emission beam location and intensity. 
Fig. C.2. Matrix of probability density plots between the model parameters corresponding to the first Gaussian component in our beam model (indexed with “01” in Table 3), plus the wind magnitude. 
Fig. C.3. Same as in Fig. C.2, for the second Gaussian component in our beam model (indexed with “02” in Table 3). 
All Tables
Properties of the observations and data products used in the analysis of this paper, listed per 100day MJD bin.
Mean and most likely values from our analysis, for the beam intensity, beamshape, and average wind magnitude.
All Figures
Fig. 1. BP1,BP2 and IP: Orbitalphase evolution of the centroids of BP1, BP2, and the IP as a function of MJD epoch. The abscissa values correspond to the mean MJD across each bin (see Table 2). The dashed red lines are the best linear fits to the data; the slope and reduced chisquared of each fit is shown in the legend. ⟨BP1+BP2⟩: The average orbitalphase shift of the centroids of BP1 and BP2 as a function of MJD epoch; the best linear fit is shown with a dashed red line. 

In the text 
Fig. 2. Pulseaveraged flux density of pulsar B at 820 MHz (filled circles), as a function of MJD, averaged over BP1 (red), BP2 (blue), the IP (green), and over the entire orbit (black). The fluxdensity calculation of the individual pulse profiles is based on the radiometer equation, assuming a duty cycle equal to the full pulse width at 10% of the pulse maximum. For the intervals MJD53000–53100 and MJD53100–53200, the flux densities (empty circles) were scaled to 820 MHz, using the spectral indices reported in Sect. 3.3, from observations at 685 and 1400 MHz, respectively. In the interval MJD54200–54300, there were no available data. The dashed lines have been used as guides only. 

In the text 
Fig. 3. Greyscale map of pulseaveraged flux density of pulsar B, as in Fig. 2, as a function of orbital phase and MJD. Bilinear interpolation between bins has been applied. The red, blue and green dashed lines delineate the orbitalphase regions of BP1, BP2, and IP as a function of MJD, assuming the corresponding 3σ widths shown in Table 2. We note that the IP is undetectable before MJD53100 and after MJD53900. The WP is defined as the orbitalphase region between the upper bound of IP and lower bound of BP1. 

In the text 
Fig. 4. Analytic template constructed from averaging pulsar B observations during the WP (ϕ_{asc}/2π ≈ 0.15–0.53), between MJD53200 and MJD53700. 

In the text 
Fig. 5. Probability density function (PDF) of the logarithm of TOA uncertainties, σ_{TOA} (in μs), from the crosscorrelation of all of the original data profiles (after RFI excision) with the WP template profile of Fig. 4. The solid, dashed, and dasheddotted vertical lines correspond to the median, the 68%, and 95% confidence intervals of the distribution, respectively. The unshaded area of the histogram corresponds to the range of uncertainties (σ_{TOA} > 9.3 ms) corresponding to TOAs that were excluded from further analysis. 

In the text 
Fig. 6. Profile evolution of pulsar B at 820 MHz, during BP1, BP2, and the IP, corresponding to two MJD bins: (top half) MJD53400–53500 and (bottom half) MJD53700−53800. Each profile comes from averaging the original data in the corresponding MJD bin and an orbitalphase interval of Δϕ_{asc}/2π = 0.02, centred around the orbital phase, ϕ_{asc}, shown in each panel. The number of single pulses averaged to obtain each profile are also shown at the topleft corner of each panel. The grey, vertical dashed lines at pulse phase 0.5 are used here to guide the eye and to highlight the relative changes between profiles observed during most BPs. We note that the profile alignment in this figure can only be used to compare the profiles of a single MJD bin and orbitalphase window (see main text for details). The blue, vertical error bar next to each profile corresponds to the quadrature sum of the offpulse RMS of the profile shown, σ_{stat}, and the average RMS due to systematic pulsetopulse variations of the profiles that were summed to obtain the profile shown (see Sect. 3.2.2). For comparison, we have overlaid the analytic WP template profile of Fig. 4 with the observed average profiles corresponding roughly to the middle of BP1, BP2, and the IP. 

In the text 
Fig. 7. (a) Peak flux density of pulsar B at 820 MHz (filled circles), as a function of MJD, calculated from the average pulse profile corresponding to the centre of BP1 (red), BP2 (blue), IP (green) and to the WP (grey). For the intervals MJD53000–53100 and MJD53100–53200, the flux densities (empty circles) were scaled to 820 MHz from observations at 685 MHz and 1400 MHz, respectively, using the spectral indices reported in Sect. 3.3. In the interval MJD54200–54300, there were no available data. The dashed lines show the best fit light curves of the peak flux density of the model. (b) Full pulse width at 10% of the maximum (W_{10}) as a function of MJD, for the average profiles corresponding roughly to the centre of the orbitalphase range of BP1 (red), BP2 (blue), the IP (green), and the WP (grey). (c) Peak separation (Δ_{P2P}) between the leading and trailing component of the average profiles. Open circles indicate that the brightest component is trailing. The dashdotted grey line is the best linear fit to the BP1 data from this work (see Sect. 3.2); the solid and dotted black lines are the best linear fits to the BP1 and BP2 data, respectively, across the corresponding ranges shown by Perera et al. (2010). All values shown have been calculated from the best fit Gaussian templates of Fig. B.3. 

In the text 
Fig. 8. Probability distributions of the spectral index of pulsar B’s radio emission, calculated assuming a single powerlaw between 685 and 1400 MHz (light grey shade) and 685 and 3100 MHz (dark grey shade), from multifrequency profiles corresponding to roughly the same orbital phase and separated in time by roughly 30 days. The joint probability distribution is shown with a solid black line, and the corresponding median spectral index and 68% confidence interval, with solid and dashed vertical lines, respectively. 

In the text 
Fig. 9. Geometry of the toy model used to fit the harmonic variation of the timing residuals in Fig. A.1. In this figure, pulsar B spins clockwise with a spin frequency of Ω_{B}, and it has a velocity of v_{B}, directed away from the LOS. The configuration shown corresponds to ϕ_{asc} = 0, where the wind of pulsar A maximally deflects the emission beam of pulsar B, perpendicularly to the LOS (direction ) by angle δϕ_{s}. The LOS direction is defined by . In this model, the emission is assumed to be generated at a distance (r_{em}) from the centre of the star. The wind of pulsar A, with magnitude equal to w_{⊥} (in units of r_{em}), deflects the emission site from position E to E′. 

In the text 
Fig. 10. Cartoon representation of spinphase delay (δϕ_{s}) produced by the impact of the wind of pulsar A, with a perpendicular component to the LOS (w_{⊥}) for two different phases of geodetic precession: at time (t_{0}) and at a later time (t_{0} + δt). In this representation, the observer’s LOS is perpendicular to the plane of the figure, when looking down onto it; the plane of the orbit is perpendicular to the orbital angular momentum vector, L; and pulsar A is located to the left of pulsar B. The precession of the spin axis, s, about L, proceeds along the path shown with a dashed, grey line, in the direction indicated by the arrow. The NS’s equator is shown with a solid, grey line, with the grey arrow indicating the spin direction; the magnetic axis, μ, rotates about the spin axis in the same direction. In both instances of time, the action of w_{⊥} causes the deflection of the emission (assumed here to be at the magnetic pole) from position E (solid black point) to E′ (open circle). The equivalent spinphase delay, δϕ_{s} at (t_{0}) and at (t_{0} + δt), is measured as the polar angle between E and E′. For the same wind magnitude (w_{⊥}), it can be seen that . 

In the text 
Fig. 11. Orbitaveraged wind magnitude (w_{0}) as a function of MJD (data points) from fits of Eq. (9) to the residuals in each of the MJD bins shown in Fig. A.1. The solid grey line and the shaded area correspond to the linear fit to the data and the 1σ confidence interval of the best fit, respectively. The inset figure shows the 2D joint probabilitydensity function of the rate of change of w_{0} with time (a) and the reference MJD, t_{0} (at which w_{0} = 0), given the entire data set of residuals. The red star symbol corresponds to the location of the most likely values of a and t_{0}. The dashed line corresponds to the most probable function that describes the evolution of w_{0} with time. 

In the text 
Fig. 12. Schematic illustration of the emission beam of pulsar B, modelled in this work as a pair of Gaussian components, of which the intensity profile is represented here as greyscale contours. This figure shows a 2D projection of the emission beam on the plane of the sky and the definitions of the beam parameters that were used in our model to describe the location and shape of the emission (Eq. (10); compare with Fig. C.1). The white circles at the centres of the elliptical contours show the locations of maximum beam intensity. In this representation, the observer’s LOS is perpendicular to the plane of the paper (as denoted by the circled cross). The circled dot at the top of the figure shows the direction of the magnetic axis () in this projection pointing out of the plane of the paper, towards the observer. The black and grey dotted lines show two examples of the trace of the LOS across the beam, at epochs t_{1} and t_{2}(> t_{1}), respectively. 

In the text 
Fig. 13. Probability density plots of four pairs of model parameters, corresponding to the second Gaussian beam component, from our analysis: (a) Γ_{02}–w_{0}, the beam colatitude of the peak intensity of the Gaussian component against the wind magnitude; (b) f_{2}–σ_{Γ2}, the flatness parameter against the halfbeam width along the direction of increasing Γ_{02}; (c) Φ_{02}–w_{0}, the beam longitude of the peak intensity of the Gaussian component against the wind magnitude; and (d) Φ_{02}–ζ_{2}, the beam longitude of the peak intensity of the Gaussian component against the angle of rotation about the direction of peak intensity. The star symbols indicate the values corresponding to the solution with the maximum global likelihood. 

In the text 
Fig. 14. Twodimensional projection on plane of the sky of most likely emission beams of pulsar B, corresponding to BP1, BP2, the IP, and the WP, as they were determined from our analysis. The Cartesian coordinates shown (X, Y) are defined in Sect. 4.2. In each panel, the intensity has been colourcoded and normalised to a maximum value of 1 (see colour bar) to provide a clear representation of the beam shape; the relative intensity between the orbital phases is not shown. The solid white lines indicate the traces of the unperturbed LOS for MJD53000, MJD53750, and MJD54500; the dashed white lines show the same traces after applying the deflection caused by the most likely magnitude of the wind (i.e. w_{0} = 0.017). The circleddot symbol indicates the position of the magnetic axis relative to the beam. A 3D rendering of these beams is provided in Fig. 15. 

In the text 
Fig. 15. (a) 3D representation of most likely emission beams during BP1, BP2, the IP, and the WP, derived from this work, shown at the minimum approach of the magnetic axis () to the LOS (). The figure shows three viewing configurations, corresponding to three geodeticprecession phases, at MJD53000, MJD53750, and MJD54500. The trace of the LOS is shown with a white line. The colour scale represents the beam intensity in arbitrary units, normalised between 0 (blue) and 1(red), in steps of 0.1. 

In the text 
Fig. 15. continued. 

In the text 
Fig. 16. Average drift rate of timing residuals across BP1 (red), BP2 (blue), and the IP (green), as a function of MJD. The drift rates have been estimated from the linear fits shown in Fig. A.2 and are shown both in units of ms per orbit and as a fraction of pulsar B’s spin phase per orbit. 

In the text 
Fig. 17. (a) Time evolution of λ, which is the angle between the spin axis of pulsar B and the LOS, and ρ, which is the angle between the spin axis of pulsar B and the radial wind direction, calculated at the centroids of BP1 and BP2, ca. MJD53600. Time, along the bottom horizontal axis, is expressed as days since the epoch of the earliest MJD bin (MJD53018.1). Along the top horizontal axis, we show the corresponding phase of geodetic precession during the interval considered. (b) Time evolution of the impact angle, β, for the magnetic pole that was visible until ca. 2008 (solid curve), and the opposite magnetic pole (dashed curve). (c) Best fit sinusoidal functions of the location of BP1 (red curve and data points) and BP2 (blue curve and data points) as a function of time. The light red curve is the alias of the best fit function for BP1, modulo the orbital period. The solid black curve and data points show the average of the functions and data, respectively. Finally, the vertical dotted grey lines indicate the earliest observation epoch and the subsequent epoch at which the impact angle has the same value as that at the earliest epoch; in (a) and (b), the horizontal dotted grey lines indicate the values of λ and β at those two epochs. 

In the text 
Fig. 18. Joint 2D probabilitydensity map of the period (in 10^{3} days) of the sinusoidal functions that were fit to the locations of BP1 and BP2 as a function of time (see Fig. 17). The median and 1σ uncertainties of the period of geodetic precession period of pulsar B, published by Breton et al. (2008), is shown with a filled red circle with error bars. The star symbol corresponds to the most likely values of the periods from our analysis. Finally, the empty circle corresponds to the period of geodetic precession predicted by GR. 

In the text 
Fig. A.1. Residuals of pulsar B for 14 100day intervals, calculated using the analytic WP template of Fig. 4, constructed from fitting a threeGaussiancomponent model to the average WP profile of pulsar B. The weighted RMS of the residuals, calculated over the entire data set shown in this figure, is 9.33 ms; the corresponding value is 75. The red lines correspond to the function of Eq. (9), using the best fit values for the parameter w_{0}. The shaded envelope corresponds to the 1σ confidence interval of the above function. The horizontal blue error bars indicate the orbitalphase intervals, at the location of peak brightness of BP1, BP2, and IP (centroids), that were used to calculate the average profiles of Fig. B.2; for the WP, the entire width of the orbitalphase window was used. Also, the panel corresponding to the interval MJD53400–54500 shows the extents of the orbitalphase regions, BP1, BP2, the IP, and the WP, for that interval, as red, blue, green, and grey shaded areas. 

In the text 
Fig. A.2. Timing residuals using fixed WP template profile of Fig. 4 (black points) and using the Gaussian templates of Fig. B.3 (red open circles). For each data set, the weighted RMS and the corresponding , calculated over the entire data set shown in this figure, is shown in the topright corner of the topleft subplot. The blue lines show the best linear fits to the residuals of BP1, BP2, and IP, which were used to correct for the linear drifts across those orbitalphase windows (see Fig. A.4). The orbitalphase extent of the lines is equal to the W_{3σ} value of the corresponding orbitalphase window. 

In the text 
Fig. A.3. Timing residuals using fixed WP template profile of Fig. 4 (black points) and the model templates of Fig. B.4 (red open circles). For each data set, the weighted RMS and the corresponding , calculated over the entire data set shown in this figure, is shown in the topright corner of the topleft subplot. 

In the text 
Fig. A.4. (black points) Timing residuals using Gaussian templates of Fig. B.3 and after correcting for the linear drifts across each orbitalphase window, via the best fits shown in Fig. A.2. For comparison, light grey circles show the timing residuals using the fixed WP template profile of Fig. 4. For each data set, the weighted RMS and the corresponding , calculated over the entire data set shown in this figure, is shown in the topright corner of the topleft subplot. 

In the text 
Fig. B.1. Distributions of S/N across the orbit of pulsar B, for each MJD bin (data points). For this figure, the values of the S/N have been normalised by the square root of the integration time. The locations of the two BPs, BP1 and BP2, and where visible that of the IP have been determined using Gaussian fits to the data (black curves). The vertical, solid, dashed and dot–dashed blue lines show the 3σ confidence intervals of the phase locations of BP1, BP2 and the IP, respectively. To emphasise the shift of the locations of BP1 and BP2, as a function of MJD, the dashed and dotted, vertical, grey lines mark the positions of the maxima of the fitted functions for BP1 and BP2, respectively, during MJD53000–53100. 

In the text 
Fig. B.2. Observed average profiles of pulsar B (black lines), corresponding to the centroids of BP1, BP2, and the IP; the WP profiles were calculated using all the data in that orbitalphase window (see blue error bars in Fig. A.1). Each column corresponds to a different orbitalphase window (labelled along the top edge of the figure); each row corresponds to a different MJD bin (labelled along the right edge of the figure). The grey lines show the lower and upper 1σ_{sys} confidence limits (Eq. (7)). The blue vertical error bar next to each profile corresponds to the offpulse RMS, σ_{stat}, of the average profile. 

In the text 
Fig. B.2. continued. 

In the text 
Fig. B.3. Comparison between observed average profiles of pulsar B (black) and the best fit Gaussian templates (red), as described in Sect. 3.2.3. The profiles shown have been tabulated according to orbitalphase window (columns): BP1, BP2, IP, and WP, and the MJD range (rows). The reduced chisquared value shown for each plot has been calculated using only the offpulse RMS (σ_{stat}) as weights. 

In the text 
Fig. B.3. continued. 

In the text 
Fig. B.4. Comparison between observed profiles (black) and those produced by our simulation, as described in Sect. 4, using the model parameters with the highest likelihood and a prograde (PG) pulsar spin (red). Green lines show the most likely profiles of an alternative model with a retrograde (RG) pulsar spin. The profiles shown have been tabulated according to orbital phase window, that is, BP1, BP2, IP, and WP, and the MJD range to which they belong. The grey lines show the lower and upper 1σ_{sys} confidence limits of the observed flux density, multiplied by the most likely error coefficients (q_{kj}) that were derived from our analysis. 

In the text 
Fig. B.4. continued. 

In the text 
Fig. C.1. (a) Geometry of the model used in this work to describe the effect of the radial wind of pulsar A (w) on the emission of pulsar B. In this figure, Ω and are the spin and magnetic moment of pulsar B, respectively, and L is the orbital angular momentum of the binary system. For clarity, we separated the elements corresponding to the frame of the orbit () and those corresponding to the frame of the pulsar () in black and red, respectively. The sense of pulsar B’s orbit around pulsar A and that of μ about Ω are shown with black and red arrows, respectively. Also shown with a black arrow is the sense of the geodetic precession of Ω about L (defined by an increasing ϕ_{SO}). We note that, compared to Breton et al. (2008), we used the opposite direction for L. The direction to the observer coincides with the unit vector . In this work, we have assumed that the double pulsar system is viewed exactly edgeon: so, the orbital inclination, i, is exactly equal to 90°. At the position of , we show an example location in the beam, with beam coordinates (γ, ϕ), corresponding to intensity I(γ, ϕ) (see Eq. (10)). For the definition of the reference frames and angles shown, the reader is directed to Sect. 4.2. (b) The definitions of the angles that are used in our parametrisation of pulsar B’s emission beam location and intensity. 

In the text 
Fig. C.2. Matrix of probability density plots between the model parameters corresponding to the first Gaussian component in our beam model (indexed with “01” in Table 3), plus the wind magnitude. 

In the text 
Fig. C.3. Same as in Fig. C.2, for the second Gaussian component in our beam model (indexed with “02” in Table 3). 

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