© The Authors 2023

1. Introduction

PSR J1933−6211 was discovered as part of the Parkes High Galactic Latitude Survey (Jacoby et al. 2007), which used the 64 m CSIRO Parkes Murriyang radio telescope in Parkes, New South Wales, Australia (henceforth the Parkes telescope), to search for radio pulsars at Galactic latitudes between 15 and 30 deg. The fully recycled nature of the pulsar, combined with a very low eccentricity (e ∼ 1.2 × 10⁻⁶), indicates that the companion very likely is a white dwarf (WD) star, whose progenitor recycled the pulsar. Jacoby et al. (2007) used timing observations of the pulsar with the Parkes telescope and the CPSR2 backend to derive the binary mass function, and estimated a minimum mass for the companion (M_c) of 0.32 M_⊙ by assuming that the pulsar mass (M_p) is 1.4 M_⊙.

The fast spin of PSR J1933−6211 (P₀ = 3.5 ms) is typical of millisecond pulsars (MSPs). Most MSPs formed in low-mass X-ray binaries (LMXB), which allow for the long accretion times required to spin up neutron stars (NSs) to these short spin periods. In these systems, the companions are helium WDs (He WDs), in which the binary orbital period (P_b) and the WD mass are thought to be correlated (Tauris & Savonije 1999). For the orbital period of PSR J1933−6211 (12.8 d), the correlation predicts a companion mass between 0.25 M_⊙ and 0.28 M_⊙, depending on the properties of the progenitor of the WD. This predicted value, using the He WD correlation, is lower than the originally estimated minimum companion mass (0.32 M_⊙), suggesting either an unusually light pulsar or that the companion is not a He WD, but a more massive type of WD, such as a carbon-oxygen (CO) WD, formed instead in an intermediate-mass X-ray binary (IMXB; Tauris et al. 2011).

These features make this system a relative rarity; there are only 4 other pulsars with spin rates below 6 ms with established or likely CO-WD companions (PSRs J1101−6424; Ng et al. 2015; J1614−2230; Alam et al. 2020; J1618−4624; Cameron et al. 2020; and J1943+2210; Scholz et al. 2015), compared to the 101 pulsar binary systems in this spin period range with He WD companions (Manchester et al. 2005)¹. Because the more massive companions evolve faster, MSPs resulting from IMXBs tend to be significantly slower (P > 9 ms) than those that result from LMXBs. Thus, the spin periods of these pulsars are somewhat anomalous if the companion is a CO WD.

The orbital inclination of one of these systems, PSR J1614−2230, is close to 90 deg, which allowed for a precise measurement of the Shapiro delay (Shapiro 1964) and showed the pulsar to have a mass that is very likely above 1.9 M_⊙ (Demorest et al. 2010; Arzoumanian et al. 2018), thereby introducing strong constraints on the equation of state of dense nuclear matter (Özel & Freire 2016). This raises the question of how the pulsar becomes this massive. A detailed study of this system suggested an alternative evolutionary pathway, in which the NS was spun up via case A Roche-lobe overflow (RLO; see Tauris et al. 2011). Case A RLO takes place when the companion is still a main-sequence star. Under these conditions, the accretion timescale is related to the hydrogen-burning timescale of the donor. This long accretion episode can therefore in principle allow the NS to gain much angular momentum, explaining the fast rotations observed in these systems. In contrast, in terms of mass gain, Tauris et al. (2011) estimated that PSR J1614−2230 gained 𝒪 (0.2) M_⊙ at most during accretion, that is, they concluded that PSR J1614−2230 is massive because it was born that way.

Motivated by the possibility that a CO WD companion might result in PSR J1933−6211 being as massive as PSR J1614−2230 (or by the possibility that if the companion is a He WD, the pulsar would have an unusually low mass), Graikou et al. (2017) attempted to measure the masses of the components of the PSR J1933−6211 system via the Shapiro delay. They timed the system using coherently dedispersed Parkes data, which significantly improved the timing precision. A timing residual root-mean-square (rms) of 1.23 μs was obtained using only a subset of bright observations in which the pulsar signal was boosted by interstellar scintillation. However, owing to unfavourable scintillation during many observations and the now-estimated low orbital inclination of the system, obtaining mass measurements was not possible. Based on their available data, the authors placed an upper limit of 0.44 M_⊙ on M_c.

In this work, we again attempt to measure the masses of the components of the PSR J1933−6211 system; this time successfully. The superior sensitivity of the MeerKAT telescope (Jonas 2009) was crucial for this. It yields a sensitivity for pulsars that is an order of magnitude better than that of the Parkes telescope (Bailes et al. 2020; the exact number depends on the spectral index of the pulsar and interstellar scintillation). This is especially important for pulsars like PSR J1933−6211, which are located so far south that they cannot be observed with any of the sensitive Northern Hemisphere telescopes. MeerKAT timing observations are carried out under the MeerTime Large Survey Project (LSP), pursuing a broad range of scientific topics (Bailes et al. 2020). The PSR J1933−6211 observations were made under two distinct research sub-themes within the MeerTime LSP: (1) the relativistic binary timing programme (Kramer et al. 2021, RelBin), which performs dedicated observations of binary pulsar systems to measure relativistic effects in their timing, with the aim of testing gravity theories and measuring NS masses, and (2) the pulsar timing array programme (Spiewak et al. 2022; PTA), which observes an array of southern millisecond pulsars to search for nanohertz gravitational waves. This paper reports the results of these timing measurements. Our results are aided by the aforementioned data from the Parkes telescope, and include additional measurements made with that telescope using the new Ultra-Wideband Low receiver (Hobbs et al. 2020).

The structure of this paper is as follows. In Sect. 2 we discuss the observations, the resulting data, and how these were analysed. In Sect. 3 we present the polarimetric profile of the pulsar, together with a rotating vector model (RVM) of the polarimetry. In Sect. 4 we present our timing results, in particular, a discussion of the most important timing parameters, the component masses, and the orbital orientation for this system. In Sect. 5 we discuss the main results and their implications for the nature of the system and for the evolution of MSP-CO WD systems in general. We also provide conclusions and future prospects.

2. Observations and data reduction

In this work, we used the 2003/2004 time-of-arrival values (ToAs) from the Parkes observations described by Jacoby et al. (2007), as well as 2011/2012 Parkes ToAs associated with the data in Graikou et al. (2017) and provided to us as 5 min time-averaged ToA values. These represent a curated ToA set from which outliers with a low signal-to-noise ratio (S/N) were removed. These ToAs were originally obtained using the Parkes 20 cm multibeam receiver (Staveley-Smith et al. 1996) with the Caltech Swinburne Parkes Recorder 2 (CPSR2) and the CASPER Parkes Swinburne Recorder (Venkatraman Krishnan 2019, CASPSR) backends. We describe the new observations we obtained below. Combined, the full timing baseline reported in this work is 19 years. An overview of all the data and their characteristics used in this work is presented in Table 1.

2.1. Parkes observations

More recent Parkes data of PSR J1933−6211 were collected through the P965 and P1032 Parkes observing programmes. This includes coherently dedispersed fold-mode observations using the ultra-wide bandwidth low-frequency (UWL) receiver with its Medusa backend (Hobbs et al. 2020); as well as a few fold-mode observations using the 20 cm multibeam receiver and the CASPSR backend. The latter setup is identical to the one used in the 2011/2012 data of Graikou et al. (2017) and therefore provides an overlap between the MeerKAT/PTUSE and Parkes/CASPSR data sets, which are otherwise separated by a large gap in observations from 2016 to 2019 that could hamper accurate phase connection for PSR J1933−6211.

We have a total of 17 UWL observations that vary in duration from 890 s to 1 h 4 min, taken between August 6, 2018 and March 18, 2022. The Parkes UWL receiver operates at a centre frequency of 2368 MHz and has a total bandwidth of 3328 MHz. In the fold-mode setup used here, it produces 1024 phase bins across the rotational phase of the pulsar.

2.2. MeerKAT observations

Data from the MeerKAT telescope were obtained between March 8, 2019 and May 16, 2022. The observations made by the PTA programme were 256 s each and were regularly spaced, with a mean cadence of two weeks, while the RelBin observations were longer (≥2048 s) and were aimed at obtaining good orbital coverage. In particular, the RelBin data set contains one 4 h observation (MJD 58746.80) and two 90 min observations (MJDs 58836.50 and 58823.69) taken close to and across superior conjunction to optimise for Shapiro delay measurements.

The MeerKAT observations were recorded using the L-band receiver (856−1712 MHz) in its 1K (1024) channelisation mode, using the Pulsar Timing User Supplied Equipment (PTUSE) backend (Bailes et al. 2020), which provided coherently dedispersed folded pulsar archives with 1024 phase bins across the pulse profile of 3.54 ms, or with a phase-bin resolution of 3.46 μs.

Prior to the observations, standard array calibration is applied via the MeerKAT science data processing (SDP) pipeline, as described in Serylak et al. (2021). This includes online polarisation calibration, such that (since April 9, 2020) the Tied Array Beam data stream ingested to PTUSE produces polarisation-calibrated L-band pulsar data products. Data recorded before access to the online polarisation calibration pipeline were calibrated offline according to the steps outlined in Serylak et al. (2021).

2.3. Data reduction

The data reduction and analysis in this section rely on a combination of well-established pulsar software suites, including PSRCHIVE (Hotan et al. 2004) and TEMPO2 (Hobbs et al. 2006; Edwards et al. 2006), as well as observatory or research programme-specific pipelines (e.g. MEERPIPE, PSRPYPE). We denote particular tools within software packages as tool/SOFTWARE.

2.3.1. Parkes: Multibeam/CASPSR

The CASPSR data taken in October 2020 were reduced in a similar manner as reported in Graikou et al. (2017). Band edges were removed and radio frequency interference (RFI) manually excised using pazi/PSRCHIVE before creating frequency integrated, intensity-only (Stokes I) profiles with 512 phase bins using pam.

Testament to the scintillating nature of PSR J1933−6211, of the four observations obtained (with observing lengths ranging from ∼1 − 1.5 h), only the two observations taken on 17 October produced profiles with S/N > 10. The brightest of these were reduced to two time intervals, and the second to a single averaged profile only.

2.3.2. Parkes: UWL/Medusa

Data from the UWL receiver were reduced using the PSRPYPE processing pipeline². The pipeline performs flux and polarimetric calibration, along with automated RFI excision using CLFD³. This works in a similar way to MEERPIPE and produces RFI excised, calibrated, and decimated to a number of time, frequency, and polarisation resolutions. To increase the profile S/N values leading up to computed ToA measurements, we further reduced the data products to four frequency channels, single time integrations, and full intensity only.

2.3.3. MeerKAT: L-band/1K PTUSE

The MeerTime observations were reduced using the MEERPIPE pipeline, which produces archive files cleaned from RFI (based on a modified version of COASTGUARD; Lazarus et al. 2016) of varying decimation using standard pam/PSRCHIVE commands. We started our customised data reduction from the output products containing 16 frequency channels across the inner 775.75 MHz of MeerKAT L-band, an eight-fold reduction in subintegration time, and calibrated Stokes information.

Based on our findings that the PSR J1933−6211 scintillation cycles last approximately 20 to 30 min on average, all longer-duration observations were decimated to have a minimum integration length of 500 s. To increase the S/N per ToA, we reduced the channelisation to eight frequency channels for all data. A rotation measure (RM) correction of 9.2 rad m⁻² was applied using pam, based on the measurement presented in Kramer et al. (2021).

2.4. Estimating pulse times of arrival

Three additional CASPSR ToAs were added to the data set using the same standard profile as was used to generate the CASPSR ToAs in Graikou et al. (2017), together with the data described in Sect. 2.3.1. This provided an overlap between the Parkes/CASPSR and MeerKAT/PTUSE ToAs.

To create ToA values from the UWL/Medusa data sets, a high S/N timing standard was created through the addition of the available observations (psradd/PSRCHIVE) following their reduction and cleaning by PSRPYPE as well as additional RFI removal by hand using pazi. Based on the obtained S/N values, we chose to create a total intensity standard with four frequency channels (providing a per channel profile with an S/N 70 to 350) and turned them into DM-corrected analytical templates using psrsmooth/PSRCHIVE. This template was then used in pat/PSRCHIVE to obtain ToA values at the telescope for the reduced UWL data described in Sect. 2.3.2.

A MeerKAT multi-frequency timing standard was created using all PSR J1933−6211 observations with an estimated S/N > 200. These were added using psradd and reduced to create a template with a single subintegration, eight frequency channels with four Stokes polarisations, with DM and RM corrections applied. Polarisation-resolved standards were motivated by the timing improvements Graikou et al. (2017) reported using matrix template matching (MTM; van Straten 2004, 2013) for ToA generation.

Analytical polarisation and frequency-resolved standards were generated from these high S/N templates by again applying wavelet smoothing using psrsmooth. These were subsequently used to apply MTM using pat on the MeerKAT data products described in Sect. 2.3.3, providing measurements of the ToAs. As shown by Graikou et al. (2017), the timing of PSR J1933−6211 benefits especially from using the MTM because of its sharp polarisation features, as shown in Fig. 1. Use of MTM strongly relies on an accurate polarisation calibration of the pulsar data. The well-calibrated MeerTime data products will therefore benefit from the use of MTM.

Fig. 1. Averaged flux-calibrated and RM-corrected PSR J1933−6211 pulse profile (bottom) and associated polarisation position angles (middle and top). Bottom: this averaged flux-calibrated and RM-corrected profile was obtained from adding all 24.5 h of MeerKAT L-band data, providing an S/N ≈ 3000, a peak flux density of 8.7 mJy, and a mean flux density of 1.1 mJy. The profile shows significant linear (red, L) and circular (blue, V) polarisation, the sharped-edged features of which make it an ideal candidate for MTM timing techniques when detecting it brightly. Middle and top: the linear polarisation PA are plotted for all phase bins with significant (> 5σ) linear polarisation. To visualise PA swings, we present the PA values, together with PA+180 deg and PA−180 deg (grey data points) in both panels. The solid (brown) line shows the best-fit RVM model when using a prior of ζ = 125 deg and obtaining α = 41.66 (4) deg (top), or a prior of ζ = 55 deg and obtaining α = 121.99 (4) deg (middle) as described in Sect. 3, with their ±180 deg equivalent as dashed (brown) lines. This fit was obtained using only the data points in black and removing the fainter data points of the centrally plotted PA values. The dotted lines (blue) in both panels are PA ± 90 deg to show the potential PA values of orthogonally polarised modes.

Finally, to account for the varying S/N values of the observations that are due to the high occurrence of scintillation in PSR J1933−6211, we manually removed individual ToAs with uncertainties larger than 20 μs from all data sets. This was done after visual inspections that confirmed that the large ToA uncertainties were indeed due to low S/N detections.

2.5. Timing analysis and orbital models

The analysis of the ToAs was made using TEMPO2⁴. The telescope-specific ToAs computed above were transformed into TT(BIPM2021)⁵, which is a realisation of terrestrial time as defined by the International Astronomical Union (IAU), and thereafter converted into time of arrivals at the Solar System barycentre using the most recent DE440 Solar System ephemeris of the Jet Propulsion Laboratory (JPL; Park et al. 2021).

Initial orbital and pulsar parameter estimates were found using the DDH orbital model description as implemented by the TEMPO2 software. This is an extension of the DD model (Damour & Deruelle 1986), and describes the Keplerian orbit via the parameters orbital period (P_b), length of the projected semi-major axis (x_p), orbital eccentricity (e), longitude of the periastron (ω), and the time of passage through the ascending node (T₀) along with several relativistic corrections, which are quantified by a set of phenomenological post-Keplerian (PK) parameters. In particular, DDH uses the orthometric amplitude (h₃) and the orthometric ratio (ς) to model the Shapiro delay, whereas the standard DD model describes it with the range (r) and shape (s) parameters (Freire & Wex 2010; Weisberg & Huang 2016). These parameters have the advantage of being far less strongly correlated than r and s, especially for low orbital inclinations, as is the case for PSR J1933−6211.

However, the DDH model fails to account for the full set of kinematic contributions described in Sect. 4.2; in particular, it does not describe the annual orbital parallax (AOP; Kopeikin 1995), but can only model the secular variation of x caused by the proper motion (Kopeikin 1996), ẋ. Consequently, it cannot discriminate between the multiple solutions for the orbital orientation of the system given by a measured ẋ and ς. Furthermore, unmodeled residual trends caused by the AOP pollute the very weak Shapiro delay signal whose higher harmonics are of the same order of magnitude as the AOP for PSR J1933−6211.

For these reasons, we refined our parameter estimations by using the T2 binary model, which is based on the DD model, but self-consistently accounts for all kinematic contributions to orbital and post-Keplerian parameters described in Sect. 4.2 (Edwards et al. 2006). Within the description of the T2 model, all kinematic effects caused by the proper motion are calculated internally from the orbital orientation of the system, given by the position angle of the ascending node (Ω, KOM) and orbital inclination (i, KIN) parameters. If astrometric dynamics is the only cause of the variation of the semi-major axis, then there is no need for an additional ẋ parameter under this paradigm.

We note that for systems with very low orbital eccentricities, such as PSR J1933−6211, ω and T₀ estimated through the DD or T2 model, for example, can be highly correlated. The ELL1-type orbital models (Lange et al. 2001) are a popular alternative to replace these with the time of ascending node (T_asc = T₀ − ω P_b) and the Laplace-Lagrange parameters, ϵ₁ ≡ e sin ω and ϵ₂ ≡ e cos ω. Similarly to the DD models, however, the ELL1-type orbital models fail to include the relevant kinematic contributions included in the T2 model. Consequently, following our T2 analysis, we derived T_asc, ϵ₁ and ϵ₂ to produce a full set of accurate timing parameters. We note that the TEMPO2 implementation of the T2 model can also work with the ELL1 parametrisation, which we also performed as a check and obtained consistent results.

In order to calculate reliable error bars and parameter correlations within the T2 model, we employed the TEMPONEST plugin to TEMPO2. TEMPONEST is a Bayesian parameter estimation tool that allows for physically motivated prior distributions on timing parameter values while also fitting for additional noise models to the data, including red noise and DM noise (Lentati et al. 2014). TEMPONEST internally uses the MULTINEST (Feroz et al. 2019) sampler. We set the multi-modal flag ON, as we expected multiple modes to be present for some of our parameters a priori.

2.6. Noise model selection

We settled on a best noise model to describe the PSR J1933−6211 timing data by performing Bayesian non-linear fits of timing models with varying noise characteristics to the data using TEMPONEST. The tested noise models included (1) a white-noise only model, where we fit for the noise parameters EFAC and EQUAD that add to or scale the uncertainties of the ToA measurements (as described in Lentati et al. 2014); (2) white noise plus a DM noise model characterised through a chromatic power-law model and; (3) white noise plus a stochastic achromatic (red) timing-noise model similarly described by a power law; as well as (4) white noise plus DM noise plus a red-noise model. To each model, we provided uniform priors centred on the initial best-fit tempo2 parameter value and ranging across ±40σ, where σ is the associated TEMPO2 uncertainty (i.e. we set FitSig to 40 in TEMPONEST). For a select set of parameters, we provided physically motivated uniform priors; (Ω, KOM) and (i, KIN) were set to cover their range of possible values: [0,360] deg and [0,180] deg, respectively; (ϖ, PX) was set to range from 0.1 to 2.2 and (M_c, M2) from to 0.1 to 1.5 M_⊙.

We performed Bayes factor (BF) comparisons between these models, and we find the strongest evidence for a red- and white-noise model, which compared to white-noise only has a BF of 16.6. Comparisons of the red- and white-noise model to models that include DM-noise yield a BF of 5.4 against DM and white noise and 1.8 against DM, white, and red noise. We conclude that all DM effects are well modelled through the inclusion of the TEMPO2 timing parameters, DM1 and DM2 (which describe the coefficients to the first- and second-order DM derivatives expressed as a DM Taylor series expansion). For the remainder of the results section, we therefore focus on the outcomes of the TEMPONEST posterior distributions, which include red- and white-noise parameters.

3. Results: Pulsar geometry

Throughout this paper, we use the observer’s convention to define angles and vectors, unless explicitly stated otherwise. In this framework, the position angle and the longitude of the ascending node (Ω) increase counter-clockwise on the plane of the sky, starting from north. Furthermore, the orbital inclination i is defined as the angle between the orbital angular momentum and the line from the pulsar to the Earth. This can vary between 0 and 180 deg. Figure 2 shows these angular definitions. The observer’s convention is also used by the T2 orbital model (Edwards et al. 2006). Any angle without a subscript follows this convention. We recall that this convention is different from the conventions used in Damour & Taylor (1992) and Kopeikin (1995), where Ω was measured clockwise from east and i is the angle between the angular momentum of the orbit and a vector pointing from the Earth to the pulsar. Angles in this alternate convention are explicitly denoted with the subscript DT92.

Fig. 2. Schematic of the 3D geometry for the PSR J1933−6211 binary in the observer’s convention. The orbital plane is defined by the unit vectors i and j, with the orbital angular momentum perpendicular to the plane, along the direction of k. The spin of the pulsar S is assumed to be aligned with the orbital angular momentum. The plane of the sky is defined by the east and north unit vectors, with n being the unit vector of our line of sight that is defined from the pulsar to the observer. The angle between n and k by definition is the orbital inclination angle, i. In the observer’s convention, this is equal to the polar angle of the pulsar ζ. The projection of S on the sky plane forms. The angle between north and i is the longitude of the ascending node. ρ denotes the opening angle of the emission cone, with β the impact angle of n on the cone, as measured from the magnetic axis μ. αp is the inclination of μ with respect to S.

Figure 1 provides our highest S/N profile for PSR J1933−6211 as obtained when adding 88 079 s (24.5 h) MeerKAT L-band data, cleaned from RFI. The flux-calibrated and RM-corrected profile has a mean flux density of 1.1 mJy and an estimated S/N = 2833.

The variation in position angle of the linear polarisation (PA; ψ) of the pulse profile across the pulsar longitude changes due to the viewing geometry, and under ideal assumptions, it results in an S-shaped swing. This is often described by the rotating vector model (Radhakrishnan & Cooke 1969, RVM), which can then provide information about the pulsar geometry. The RVM describes ψ as a function of the pulse phase, ϕ, depending on the magnetic inclination angle, α, and the viewing angle, ζ, which is the angle between the line-of-sight vector and the pulsar spin and can be written as

$$\begin{array}{c}\hfill \psi ={\psi }_0+\arctan \left(\frac{\sin \alpha \sin ({\varphi }_0-\varphi )}{\sin \zeta \cos \alpha -\cos \zeta \sin \alpha \cos ({\varphi }_0-\varphi )}\right),\end{array}$$(1)

where we have modified the equation to follow the observer’s convention.

Many studies have shown that deviations from the RVM model are typical especially for MSPs (e.g. Yan et al. 2011; Dai et al. 2015), and we therefore do not expect good agreement with the RVM model for PSR J1933−6211. However, in particular cases, such as for MSP PSR J1811−2405, the PA values follow an RVM model, which has proven effective in breaking the x_p sin i degeneracy in Shapiro delay measurements to obtain an accurate orbital inclination (Kramer et al. 2021).

Our obtained PA values for PSR J1933−6211, shown in the top and middle panels of Fig. 1, clearly exhibit more complex variations than the simple RVM S-shaped swing described above. The sharp change in slope of the PA points, especially towards the rear end of the profile, suggests that the orthogonally polarised modes are mixed, which makes some PA points unreliable.

Within our plotted PA values, S-shaped curves are discernible, and we therefore attempt to fit Eq. (1) to select PA values after removing points that deviated from an RVM-like swing. A blind fit of the remaining points after accounting for a PA jump from orthogonally polarised modes (at phase 0.43) and with a flat prior on α and ζ, provides a surprisingly precise value for ζ of 34 (1) deg. The posterior distribution of α is bimodal at both 36 (1) deg and 34 (1) deg. These values are consistent with the similar analysis of Kramer et al. (2021), according to the DT92 convention.

For systems in which the spin of the pulsar is expected to be aligned with the orbital angular momentum (e.g. PSR J1933−6211), ζ ≡ i ≡ 180 − i_DT92. However, our timing measurement of the inclination angle (see Sect. 4) is inconsistent with ζ (or 180 − ζ). This confirms that the PA swing indeed follows more complex variations than can be explained by RVM.

As an additional check, we set ζ ≡ i and ζ ≡ 180 − i (where i was obtained from timing) as the prior and performed constrained fits. These provided an α value of 41.66 (4) and 121.99 (4) deg, respectively. The corresponding RVM curves are shown in the middle and top panels of Fig. 1. However, a Bayes factor test between the blind and constrained fit shows that the blind fit is strongly favoured (BFs > 400). We conclude that even for the curated PA points that seem to follow an S-type curve, their variations do not follow the RVM.

4. Results: Timing analysis

The complete set of spin and astrometric timing parameters using the T2 model is provided in Table 2, while the measured binary parameters are contained in Table 3. We also present a set of derived quantities in both tables, which include the values for timing parameters used in ELL1-type orbital models (T_asc, ϵ₁ and ϵ₂). These are useful to record fold-mode data for the pulsar.

We used the CHAINCONSUMER library (Hinton 2016) to visualise the TEMPONEST T2 posterior distributions, with corner plots showing the 1D and 2D posterior distributions of the parameters. Figure 3 shows the resulting output for a subset of timing parameters of interest. Here, we obtained the red-noise model in TEMPONEST with 5000 live points to produce well-sampled distributions. Parameter error bars are 1σ uncertainties following the default smoothing as applied through CHAINCONSUMER.

Fig. 3. Output posterior distributions for PSR J1933−6211 parameters. These were obtained from TEMPONEST sampling using the T2 orbital model and a red-noise model, and were generated using the CHAINCONSUMER package. We present the orbital and post-Keplerian parameters. Full details are provided in Table 2. The obtained pulsar mass (Mp) distribution was computed using the mass function and the posterior distributions on Mc, i, xp, and Pb.

The solution presented in Tables 2 and 3 provides a good description of the timing data. In the top panel of Fig. 4, we show the timing residuals having implemented this best-fit model as a function of the observing date and observing system (see the figure caption for a description of the colouring). The timing residuals show the difference between the observed barycentric ToA value (obtained using the techniques described in Sect. 2.4) and the predicted barycentric arrival time for that particular pulsar rotation based on the single best-fit timing model as obtained above. The validity of the timing model is evident from the low weighted rms (w_rms = 1.372 μs) of the residuals, the obtained reduced χ² value of 0.99. There appears to be no unmodelled trends in the residuals.

Fig. 4. Timing residuals across epochs (top) and orbital phase (bottom). Data include historic Parkes ToAs (green and purple; PKS Multibeam CPSR2 1341 MHz and 1405 MHz) recent and historic Parkes CASPSR data (red; PKS Multibeam CASPSR), as well as recent Parkes UWL/Medusa ToAs (blue; PKS UWL Medusa) and MeerKAT L-band ToAs (gold; MKT L/1K PTUSE). Top: after applying the best-fit values as obtained using the timing and noise model described in Sect. 2.5, we obtain a weighted rms of 1.259 μs. Bottom: residuals as a function of orbital phase, measured from superior conjunction (the ascending node is at a phase of −0.25). A Shapiro delay signal is discernible when setting the companion mass to zero (Mc = 0) while keeping all other parameters fixed, including xp sin i. We overplot the expected theoretical signal based on the best-fit inclination and companion mass values in black. The line width indicates combined 1σ deviations in companion mass and inclination angle.

In subsequent sections, we highlight a few of the physically interesting parameter results obtained from the timing analyses and resulting posterior distribution, in particular, some of the astrometric parameters (including parallax and distance estimates) and the PK parameters that allow for estimates of the component masses and orbital orientation of the system.

4.1. Proper motion

Our updated position and proper motion values for PSR J1933−6211 provide an improvement in precision by a factor of 8 compared to the values published in Graikou et al. (2017). From the measured proper motion values in right ascension and declination (μ_α, μ_δ), we obtain a total proper motion magnitude value of 12.42 (3) mas yr⁻¹. The corresponding position angle of the proper motion, Θ_μ = tan⁻¹ (μ_α/μ_δ) is 333.1 (1) deg in the observer’s convention (see Sect. 3 and Fig. 2).

4.2. Kinematic effects on the pulsar timing parameters

The moderate distance of PSR J1933−6211 and the combination of relatively large proper motion, large projected semi-major axis of its orbit (x_p), relatively low orbital inclination of i ∼ 55 deg and high timing precision provide a rare combination of criteria that enable the detection of subtle kinematic effects that help constrain the 3D geometry of the system. These effects, first described in detail by Kopeikin (1995, 1996), must be modelled precisely; otherwise, the unaccounted-for delays will pollute our measurement of the weak Shapiro delay in this system. This is a consequence of the small i.

We now describe these effects in more detail. They depend on the absolute orientation of the system, which is given by the position angle of the line of nodes (the intersection of the orbital plane with the plane of the sky), Ω, and the orbital inclination, i.

4.2.1. Proper motion contributions to $\dot{\omega }$ and ẋ

The high proper motion of the PSR J1933−6211 binary leads to a constant change in the viewing angle of the pulsar, which manifestd as a constantly changing longitude of periastron (ω) and orbital inclination i; the latter might measurably change x, which is given by x_p sin i/c, even if the semi-major axis of the pulsar orbit (x_p) does not actually change. In the observer’s convention, these kinematic contributions to ẋ and $\dot{\omega }$ are given by

$$\begin{array}{cc}& {\dot{x}}_{\mu }=1.54\times {10}^{-16}xcoti(-{\mu }_{\alpha }cos\mathrm{\Omega }+{\mu }_{\delta }sin\mathrm{\Omega }){\mathrm{s}\mathrm{s}}^{-1}\hfill \end{array}$$(2)

$$\begin{array}{cc}& \dot{{\omega }_{\mu }}=2.78\times {10}^{-7}csci({\mu }_{\alpha }cos\mathrm{\Omega }+{\mu }_{\delta }sin\mathrm{\Omega })\mathrm{deg}{\mathrm{yr}}^{-1}.\hfill \end{array}$$(3)

The expression in Eq. (2) is identically to Eq. (1) in Guo et al. (2021), for example, which provides a convention-independent alternative. We note that for Eq. (2) to become valid in the DT92 convention, the angles need to be transformed accordingly, with Ω = 90° −Ω_DT92 and i = 180° −i_DT92.

Given the low orbital eccentricity, we do not measure a significant $\dot{\omega }$. However, we measure a highly significant ẋ, 1.60 (5)×10⁻¹⁴ s s⁻¹ assuming the DD model. A detailed analysis of all possible contributions to ẋ (e.g. Lorimer & Kramer 2012) shows that this must be almost exclusively caused by the proper motion according to Eq. (2). For this reason, the measured ẋ leads to constraints on the orbital orientation of the system, that is, Ω and i (see Sect. 4.3 for details).

4.2.2. Annual orbital parallax

The variation in the Earth’s position as it orbits the Sun causes small annual changes to x_p (from the apparent change in the orbital inclination caused by the Earth’s motion) and ω. This effect, termed the AOP, is generally very small. However, it is the key for determining the absolute orbital orientation of the system. The reason is that if we measure sin i from Shapiro delay and ẋ for a particular binary, Eq. (2) yields four possible solutions for Ω and i. Measuring the impact of the annular orbital parallax through ẋ allows us to ultimately break this degeneracy.

This cyclic effect of AOP, which has variations at both orbital and annual timescales, imprints itself on x and ω and can be expressed as in Kopeikin (1995),

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\pi }=-\frac{1}{\mathit{cd}}(\mathit{r}·{\mathit{r}}_{\mathrm{p}}-({\mathit{K}}_0\mathit{r})({\mathit{K}}_0{\mathit{r}}_{\mathrm{p}})),\end{array}$$(4)

with c the speed of light, and d the distance between the binary and the SSB. The vectors r = (X, Y, Z) and r_p describe the Earth’s position with respect to the SSB and the pulsar position with respect to the SSB, respectively. The unit normal vector K₀, points from the SSB to the barycentre of the binary. The values of r and r_p will depend on the Solar System ephemeris model that is employed, and they vary with time.

Following the expressions in Kopeikin (1995), we simplify Eq. (4) to obtain an estimate on the expected peak-to-peak amplitude of the AOP. In doing so, we make the simplifying assumption that both the pulsar’s binary orbit and Earth’s orbit are circular (e = 0) and find

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\pi }=& \frac{x_{\mathrm{p}}}{d}[({\mathrm{\Delta }}_{{\mathit{I}}_0}sin\mathrm{\Omega }-{\mathrm{\Delta }}_{{\mathit{J}}_0}cos\mathrm{\Omega })sin({\omega }_{\mathrm{PB}}t)coti\hfill \\ \hfill & -({\mathrm{\Delta }}_{{\mathit{I}}_0}cos\mathrm{\Omega }+{\mathrm{\Delta }}_{{\mathit{J}}_0}sin\mathrm{\Omega })cos({\omega }_{\mathrm{PB}}t)csci],\hfill \end{array}$$(5)

where x_p, Ω, and i are as before (and given in Table 3), and ω_PB = 2π/P_b is the binary orbital frequency.

The unit vectors (I₀, J₀, K₀), describe the coordinate system of the pulsar reference frame, with its origin at the binary system barycentre. Following Kopeikin (1995),

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{{\mathit{I}}_0}& =(\mathit{r}·{\mathit{I}}_0)=-Xsin\alpha +Ycos\alpha \hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{{\mathit{J}}_0}& =(\mathit{r}·{\mathit{J}}_0)=-Xsin\delta cos\alpha -Ysin\delta sin\alpha ,\hfill \end{array}$$(7)

with (α, δ) the right ascension and declination of PSR J1933−6211, and r = (X, Y, Z) as before. Using the same current JPL solar ephemeris as in our timing results (DE440), which is contained within the jplephem package and implemented in astropy, we obtain the Earth’s (X, Y, Z) coordinates as a function of our observing MJD range.

We next use Eq. (5) to compute the resulting Δ_π oscillatory trend as a function of MJD and find a peak-to-peak orbital parallax of PSR J1933−6211 of ∼100 ns.

Table 2 shows that the precision of x_p, following our timing timing analysis, is of the order of 50 ns, so that an AOP contribution per ToA ranging from approximately −50 to 50 ns will have a measurable and time-dependent impact on x. The importance of using the T2 model to account for this AOP and its contribution to ẋ is re-emphasised by this comparison.

4.2.3. Distance estimates from Ṗ_b and ϖ

We measure a decay of the orbital period of 7 (1)×10⁻¹³ s s⁻¹, as presented in Table 3. This Ṗ_b measurement can arise from a number of contributing effects,

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b}}={\dot{P}}_{\mathrm{b}}^{\mathrm{GR}}+{\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}+{\dot{P}}_{\mathrm{b}}^{\dot{m}}+{\dot{P}}_{\mathrm{b}}^{\mathrm{T}},\end{array}$$(8)

where the terms indicate contributions due to gravitational wave decay (GR), kinematic contributions due to changing Doppler shift (kin), mass loss in the system (ṁ), and tidal dissipation of the orbit. We find that the only non-negligible contribution for PSR J1933−6211 arises from the kinematic contributions, which consist of two secular acceleration effects,

$$\begin{array}{c}\hfill {\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{kin}}={\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{Shk}}+{\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{Gal}}·\end{array}$$(9)

Here, ${\dot{P}}_{\mathrm{b}}^{\mathrm{shk}}$ is the acceleration due to transverse motion, also known as the Shklovskii effect, and ${\dot{P}}_{\mathrm{b}}^{\mathrm{Gal}}$ is the acceleration of the binary in the gravitational field of the Milky Way due to differential rotation. ${\dot{P}}_{\mathrm{b}}^{\mathrm{shk}}$ depends on the transverse proper motion of the pulsar (μ_T) and the distance to the pulsar (d_PSR) and is related by

$$\begin{array}{c}\hfill {\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{Shk}}=2.43\times {10}^{-21}{\left(\frac{{\mu }_{\mathrm{T}}}{\mathrm{mas}{\mathrm{yr}}^{-1}}\right)}^2\left(\frac{d_{\mathrm{PSR}}}{\mathrm{kpc}}\right)·\end{array}$$(10)

${\dot{P}}_{\mathrm{b}}^{\mathrm{Gal}}$ also depends on d_PSR, along with a rotation model for the Galaxy that provides the position of the Solar System and the pulsar with respect to the Galactic barycentre, and their relative accelerations. To compute the planar and azimuthal Galactic contribution to Ṗ_b, we use

$$\begin{array}{cc}& {\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{Gal}\mathrm{pl}}=-\frac{1}{c}\frac{{\mathrm{\Omega }}_{\odot }^2}{R_{\odot }}(cosl+\frac{\beta }{{\beta }^2+{sin}^{2}l})cosb\hfill \end{array}$$(11)

$$\begin{array}{cc}& {\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{Gal}\mathrm{az}}=-\frac{K_z\mid sinb\mid }{c},\hfill \end{array}$$(12)

as in Lazaridis et al. (2009), and implemented in the GALDYNPSR library (Pathak & Bagchi 2018). Here, (l, b) are the Galactic coordinates of the pulsar, and β ≡ (d/R_⊙) cos b − cos l, R_⊙ = 8.275 kpc, and Ω_⊙ = 240.5 km s⁻¹ are the Galactic distance of Earth and the orbital velocity. Current estimates of these parameters can be obtained from McMillan (2017), where K_z/c is the vertical component of Galactic acceleration,

$$\begin{array}{c}\hfill \frac{K_z}{c}=-1.08\times {10}^{-19}[0.58+\frac{1.25}{{(z_{\mathrm{kpc}}^2+0.0324)}^{1/2}}]z_{\mathrm{kpc}},\end{array}$$(13)

with z_kpc the vertical height of the pulsar in kiloparsec (Pathak & Bagchi 2018).

Since both the Shklovskii and the Galactic acceleration effects depend linearly on d_PSR, we can use the Ṗ_b measurement to provide a constraint on the pulsar distance independent of the more standard distance constraint obtained from the parallax measurements (Bell & Bailes 1996).

From our timing analysis, we have a direct measurement of the pulsar parallax of ϖ = 1.0 (3). A simple inversion of this measurement provides a distance estimate of 1.0 (3) kpc. However, given the low ∼5σ significance of the measurement, this simple inversion is prone to the Lutz-Kelker bias (Lutz & Kelker 1973) of exponentially increasing stellar density with distance. We corrected for this bias using a scaled probability density function following Antoniadis (2021; see also Verbiest et al. 2012; Bailer-Jones et al. 2018; Jennings et al. 2018),

$$\begin{array}{c}\hfill P(d_{\mathrm{PSR}}|\varpi )=\frac{1}{2L^3}{e}^{-d_{\mathrm{PSR}}/L}{d}^2{e}^{{(-1/d_{\mathrm{PSR}}-\varpi )}^2/2{\sigma }_{\varpi }^2}.\end{array}$$(14)

Here, we adopted an exponentially decreasing space density (EDSD) prior to avoid the divergence issues implicit in the original Lutz-Kelker correction (see Bailer-Jones et al. 2018, for details). L can be thought of as a characteristic length scale, which we set equal to 1.35 kpc, following Antoniadis (2021). The estimate of the distance corrected for the L-K bias from the timing parallax is $1.2_{-0.3}^{+0.5}$ kpc.

As described above, we can also obtain an additional distance estimate from the kinematically dominated Ṗ_b value. We used model C of the GALDYNPSR library (Pathak & Bagchi 2018) to evaluate Eqs. (10)–(13) together with the current values of R_⊙ and Ω_⊙ given above, to compute all kinematic contributions to Ṗ_b. We obtain a distance estimate of $1.7_{-0.3}^{+0.2}$ kpc.

This is consistent with the L-K corrected distance estimate from parallax. We also combined the probability densities of the distance estimates from ϖ and Ṗ_b to obtain a more constraining distance of 1.4 (2) kpc and $1.6_{-0.3}^{+0.2}$ kpc without and with correcting for the L-K bias, respectively. Figure 5 provides the PDF of the distance constraints for all these considerations.

Fig. 5. Constraints on the distance to PSR J1933−6211. The dash-dotted grey and blue lines indicate the constraints from the measurement of timing parallax without and with corrections for Lutz-Kelker bias using an EDSD prior, respectively. The dashed brown line indicates constraints from the kinematic contribution to the measured orbital period decay (${\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}$). The solid grey line shows the combined constraints from a simple inversion of the timing parallax and ${\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}$, with grey hatched lines marking the 68% C.L. This corresponds to a distance of 1.4 (2) kpc. The solid blue line shows the combined constraints from the timing parallax with the EDSD prior and ${\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}$, with the solid blue shaded region denoting its 68% C.L. This corresponds to a distance of $1.6_{-0.3}^{+0.2}$ kpc.

Comparing these distance estimates to the DM-based distance estimates of the NE2001 and YMW16 electron density models, which predict 510 pc and 650 pc, respectively, we find that both electron density models significantly underestimate the distance along this line of sight. We note that discrepancies between DM-estimated distances and parallax-inferred distances are common (e.g. Stovall et al. 2019), especially for high Galactic latitudes, and that independent distance measurements serve to improve electron density models for particular lines of sight.

4.3. Shapiro delay, masses, and orbital orientation

As Fig. 4 shows, the Shapiro delay signal in this pulsar has a maximum of only 7.16 μs, a consequence of the far-from-edge-on configuration, and the reason why this delay was not detected until now. From the DDH model, we can estimate ς and h₃ from this signal. Combining this with the measurement of ẋ, we obtain constraints on the orbital inclination i and Ω. These are depicted graphically in Fig. 6, where the constraints from ς are presented by the dotted black lines and the constraints from ẋ are presented by the brown lines. According to these DDH obtained values, two possible solutions exist, one solution with i ∼ 55 deg, and a second solution with i ∼ 125 deg.

Fig. 6. Allowable orbital orientation for PSR J1933−6211 binary. In the main panel, we display the orbital orientation space in the observer’s convention (the Ω according to the DT92 convention is shown by the scale on the right, the cos i in that convention is simply −cos i). The lines represent constraints from two parameters measured with a DDH solution: in black we show the constraint from ς, and in brown we show the constraint from ẋ (1σ solid, 3σ dashed) assuming it is solely caused by the proper motion of the system; the position angle of the latter is shown by the dot-dashed purple line. The contours, enclosing progressively darker shades of red, include 98%, 86% and 39% C.L. of the 2D probability density function (pdf) and represent the T2 timing model solutions calculated by TEMPONEST and as presented in Table 3. This model takes into account all kinematic effects, including the AOP. The marginalised constraints on cos i and Ω are shown as 1D histograms (in red) in the top and side panels, respectively. For this model, with T0 ∼ 53 004, the AOP is clearly detected, as we can see from the fact that the degeneracy between the two possible cos i − Ω solutions has been lifted: no probability remains in the solution at Ω ∼ 60 deg, cos i ∼ −0.6. We also include the cos i − Ω probability density contours as obtained with T0 set to 55 734 (dashed grey) and 58 836 (dot-dashed blue), respectively. As described in Sect. 4.3.1, we unexpectedly observe i and Ω estimates to depend on T0.

Similarly, the T2 binary timing model can be used to obtain Shapiro delay estimates (see Table 3). However, this model simultaneously takes into account the effect of the AOP, the importance of which will become clear below.

We performed a full TEMPONEST analysis of the parameters in the T2 model. The resulting parameter uncertainties and their correlations are shown in Fig. 3.

The associated values for Ω and cos i are shown as blue shaded contours in Fig. 6. Comparing these with the constraints derived from the DDH model, we see that the degeneracy between the two possible i − Ω solutions is lifted: The solution at i ∼ 125 deg is excluded by the measurement of the AOP within the T2 model, which constrains Ω to ${255}_{-14}^{+8}$ deg. Within this narrow window, the inclination range is better constrained by ẋ than by ς; a consequence of this is that the uncertainty on cos i in the top panel is significantly narrower than the 1σ uncertainty of ς, resulting in an unusually precise measurement of i, 55 (1) deg. These estimates are for T₀ ≃ 53004.16. The reason we specify this will become evident below.

4.3.1. Dependence of the orbital orientation estimates on T₀

We observe an unexpected dependence of the constraints derived for the 3D orientation of the pulsar (i.e. i and Ω) on the fiducial orbit that we chose to measure T₀ (or equivalently, T_asc). While we describe the changes only with respect to T₀ in the following, we observe a similar dependence using T_asc, using the ELL1 formulation within the T2 model.

What appears to be a significant detection of AOP at T₀ ≃ 53 004.13, strong enough to entirely rule out the other Ω − cos i island in Fig. 6, becomes less significant for T₀ values set to the later epochs of the data set. Of the three distinct observing campaigns on the pulsar (see Fig. 4), the T₀ value in Table 3 is roughly in the middle of the first campaign with the CPSR2 backend, conducted soon after discovery. We repeated all the analyses with T₀ ≃ 55 734.66 at the centre of the CASPSR data taken around 2011 and T₀ ≃ 58 836.96, which is the middle of our latest, largest, and most sensitive dataset from the MeerKAT L-band and the Parkes UWL receivers. The corresponding posteriors of i and Ω are also shown in Fig. 6, where the reduction in our sensitivity to AOP is evident. We rigorously tested whether these dependences were due to our software implementations by performing simulations that we detail in Appendix A. We also repeated our analysis of the data with twice the number of TEMPONEST live points (i.e. 10 000) for T₀ ≃ 53 004.13 and 55 734.66 to understand whether we sufficiently sampled the global minima. While doing this, we extended the initial prior range for the parameters without physically motivated custom priors from ±40σ to ±100σ (the parameters with custom priors already had liberal prior distributions; see Sect. 2.6). This ensured that we sampled a larger parameter space and that our solutions were indeed the global minima. We find results consistent with Fig. 6, and for T₀ ≃ 53 004.13 consistent with Table 3. Based on these results and the simulations, we conclude that we do not find strong evidence that the T₀ dependence is caused by the timing software or the analysis method.

This leaves the tantalising possibility that this is indeed physical, which we do not fully understand. The fact that regardless of T₀, we obtain probability islands in the same quadrants as in Fig. 6 validates the robustness of our measurement of ẋ. All other parameters are seen to be almost identical across all the T₀ runs. The nominal proper motion of the system combined with a long orbital period negates the need for any additional or higher-order corrections to the astrometric and relativistic parameters other than what is already modelled by the T2 model. Hence, the physical origin of the dependence of AOP on T₀ is currently unclear. However, because our simulations suggest that we might be able to consistently obtain the 3D position for a similar data set with the same cadence and noise properties (see Appendix A for more details), we chose the T₀ value that provided the most constraints on the 3D geometry for Table 3.

4.3.2. Self-consistent mass measurements

Regardless of the sense of i, its precise measurement means that the weak Shapiro delay signal is used solely to determine the companion mass, M_c = 0.43 (5) M_⊙, that is, no precision is lost because of the correlation between i and M_c. From the mass function and the precise i and M_c, we find a pulsar mass of $M_{\mathrm{p}}=1.4_{-0.2}^{+0.3}{M}_{\odot }$. This mass measurement is consistent for all T₀ values.

The orbital models used in these analyses are independent of theory; however, we know from many other experiments (Bertotti et al. 2003; Freire et al. 2011; Guo et al. 2021) that for weakly gravitating objects such as the Sun or WD stars, the Shapiro delay constraints from the Shapiro delay parameters can be translated directly into the constraints on M_c and i. In addition, the constraints from ẋ and the AOP are purely geometric, such that our TEMPONEST analysis with the T2 model yields M_c, i, and Ω directly, without the need for further assumptions on the theory of gravity used, as would be required if additional PK parameters had been measured.

4.4. Testing the time variation of the gravitational constant with PSR J1933−6211

The fact that our measurement of Ṗ_b is consistent with almost entirely resulting from kinematic contributions together with an independent pulsar distance measured by the timing parallax allowed us to perform a test of the rate of change in the (local) gravitational constant ($\dot{G}$) over the time span of our observations. This change in G is predicted by several classes of alternative theories of gravity, including scalar-tensor gravity. This would produce an additional contribution to Ṗ_b that we can assume to be the residual measurement,

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b}}^{\mathrm{res}}={\dot{P}}_{\mathrm{b}}^{\mathrm{obs}}-{\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}.\end{array}$$(15)

Using the nominal 1σ uncertainty of the distance from the L-K corrected estimate of ϖ, we obtain ${\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}=4_{-1}^{+2}\times {10}^{-13}$ s s⁻¹ and hence ${\dot{P}}_{\mathrm{b}}^{\mathrm{res}}=-2_{-3}^{+2}\times {10}^{-13}$ s s⁻¹.

This residual Ṗ_b can be compared (to leading order and assuming zero contribution from the companion because it is a WD) with the expected Ṗ_b from $\dot{G}$,

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b}}^{\dot{G}}=-2\frac{\dot{G}}{G}[1-\frac{2M_{\mathrm{p}}+3M_{\mathrm{c}}}{2(M_{\mathrm{p}}+M_{\mathrm{c}})}{s}_{\mathrm{p}}]P_{\mathrm{b}},\end{array}$$(16)

where s_p is the sensitivity of the NS, which is defined as

$$\begin{array}{c}\hfill s_{\mathrm{p}}\equiv -\frac{\partial M_{\mathrm{p}}}{\partial G}{\mid }_N,\end{array}$$(17)

where N is the fixed number of baryons in the NS (Lazaridis et al. 2009). This sensitivity of an NS depends on the mass, the equation of state (EoS), and the theory of gravity considered. Rewriting Eq. (16) as

$$\begin{array}{c}\hfill \frac{\dot{G}}{G}=\frac{\mathbb{A}}{(1-\mathbb{B}{s}_{\mathrm{p}})},\end{array}$$(18)

we obtain $\mathbb{A}=-1.1_{-0.7}^{+1.1}\times {10}^{-19}$ s⁻¹ and 𝔹 = 1.115 ± 0.015. Similar to Zhu et al. (2019), we considered Jordan–Fierz–Brans–Dicke (JFBD) theory and AP4 EoS as an example and find s_p ≃ 0.17. This provides a limit on $\dot{G}/G=-4_{-3}^{+4}\times {10}^{-12}$ yr⁻¹, consistent with the prediction of $\dot{G}/G=0$ by General Relativity. Similar tests have been conducted using PSRs J0437−4715, J1713+0747 and J1738+0333 (Verbiest et al. 2008; Zhu et al. 2019; Freire et al. 2012), for instance, the most constraining of which is J1713+0747, which is about four times more sensitive than our results here. Future timing measurements that increase the significance of the timing parallax will aid in performing more stringent tests of $\dot{G}$.

5. Discussion and conclusions

We have presented the results of our timing of PSR J1933−6211, which combined recent Parkes and MeerKAT timing measurements with earlier Parkes measurements, for a total timing baseline of about 19 years. Because of the high timing precision provided by MeerKAT, the results include precise astrometry, in particular, the first measurement of the parallax of this system, the measurement of several kinematic effects on the binary orbit (including AOP), and a first measurement of its Shapiro delay. The measurement of the AOP is noteworthy, as this effect has only been detected in four pulsar binaries, namely PSRs J0437−4715 (the closest and brightest MSP in the sky, van Straten et al. 2001; J2234+0611 Stovall et al. 2019; J1713+0747 Zhu et al. 2019; and J2222−0137 Guo et al. 2021).

A detailed analysis of the above effects allowed us for the first time to measure the component masses: $M_{\mathrm{p}}=1.4_{-0.2}^{+0.3}{M}_{\odot }$ and M_c = 0.43 (5) M_⊙ and the full orbital orientation of the system (i = 55 (1) deg, $\mathrm{\Omega }={255}_{-14}^{+8}$ deg), although the robustness of the latter measurements is seen to depend on the fiducial T₀, as seen in Sect. 4.3.1. The root cause of this dependence is currently unclear. An independent measurement of i and Ω will allow a better understanding of this problem. This independent measurement is possible using scintillation velocity measurements, as has been demonstrated by Reardon et al. (2019), although the current data set does not have the necessary frequency resolution needed for the analysis.

Nevertheless, the mass measurements are robust; the companion mass is significantly more massive than the Tauris & Savonije (1999) prediction for He WDs, indicating that the companion is most likely a CO WD.

We note that the estimated characteristic age of ∼24 Gyr of the pulsar exceeds the Hubble time. This emphasises that for the life cycles of recycled millisecond pulsars, the characteristic age tends to lose its meaning as the underlying assumptions are no longer valid. This implies, for instance, that after recycling, the spin period of this pulsar was close to its current spin period. Nevertheless, we expect this recycled MSP to have a real age of several billion years, such that the WD companion is likely old and cool. Hence, optical observations of PSR J1933−6211, combined with the mass and distance estimates derived herein, can be used to test WD cooling models (Bhalerao & Kulkarni 2011; Kaplan et al. 2014; Bassa et al. 2016; Bergeron et al. 2022). Similarly, optical and infrared photometry can constrain the atmospheric composition of the WD, and using the cooling models that survive the tests above, can provide an estimate of its cooling age. These will allow constraining the spin of the pulsar at birth, and placing additional constraints on the accretion history and origin of the system (Bhalerao & Kulkarni 2011; Tauris et al. 2011).

By analogy with PSR J1614−2230, it is possible that these fast-spinning pulsars with CO WD companions evolved via case A RLO. The very long accretion episode associated with case A RLO is consistent with the very old characteristic age and low B-field of PSR J1933−6211 (see Table 2). Despite this, the mass of PSR J1933−6211 implies that it has not gained more than ∼0.2 M_⊙. This suggests that accretion is generally extremely inefficient. These conclusions agree with the conclusions of Tauris et al. (2011), who pointed out that PSR J1614−2230 is massive mainly because it was born this way, with mass transfer accounting for at most 0.2 − 0.3 M_⊙. It also agrees with the wider range of MSP masses, where no obvious correlation with spin or orbital parameters has been observed; even the eccentric MSPs, which have a rather uniform set of orbital parameters that suggest a uniform evolutionary mechanism, seem to have a wide range of masses (e.g. Serylak et al. 2021 and references therein). This provides additional evidence that NS masses are in general acquired at birth, and are not much affected by their subsequent evolution, instead being a product of supernova physics.

Finally, the measurements presented in this work highlight the capabilities of MeerKAT for precise timing and detailed investigations of pulsar binaries. Without the great sensitivity of MeerKAT, most of these results would not have been obtainable. For example, continuing a monthly campaign on PSR J1933−6211 for the next five years should lead to an increase in the detection significance of Ṗ_b by a factor of 3, and consequently, in equal fashion, improve our distance measurements and constraints on $\dot{G}$. Within the next few years, many other southern binaries will not only have their masses measured accurately, but several of them will also yield new tests of gravity theories from the measurement of multiple PK parameters as part of the MeerTime/RelBin project.

Acknowledgments

We thank the referee for valuable comments on the manuscript. We thank Norbert Wex, Kuo Liu and Matthew Miles for valuable discussions and Robert Main for comments on the manuscript. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time space reference systems department department of the Paris Observatory. MeerTime data is housed on the OzSTAR supercomputer at Swinburne University of Technology maintained by the Gravitational Wave Data Centre and ADACS via NCRIS support. The Parkes radio telescope (Murriyang) is part of the Australia Telescope National Facility (https://ror.org/05qajvd42) which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Wiradjuri people as the traditional owners of the Observatory site. This research has made extensive use of NASA’s Astrophysics Data System (https://ui.adsabs.harvard.edu/) and includes archived data obtained through the CSIRO Data Access Portal (http://data.csiro.au). Parts of this research were conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. V.V.K., P.C.C.F., M.K., J.A., M.Ci.B., D.J.C., and A.P. acknowledge continuing valuable support from the Max-Planck Society. J.A. acknowledges support from the European Commission (Grant Agreement number: 101094354), the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of “Science and Society – Action Always strive for excellence – “Theodoros Papazoglou” (Project Number: 01431). A.Po. and M.Bu. acknowledge the support from the Ministero degli Affari Esteri e della Cooperazione Internazionale – Direzione Generale per la Promozione del Sistema Paese – Progetto di Grande Rilevanza ZA18GR02. M.Bu. and A.Po. acknowledge support through the research grant “iPeska” (PI: Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016. R.M.S. acknowledges support through Australian Research Council Future Fellowship FT190100155. J.P.W.V. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039). This publication made use of open source python libraries including Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration 2018) and Chain Consumer (Hinton 2016), galpy (Bovy 2015), GalDynPsr (Pathak & Bagchi 2018) along with pulsar analysis packages: PSRCHIVE (Hotan et al. 2004), TEMPO2 (Hobbs et al. 2006), TEMPONEST (Lentati et al. 2014).

References

Appendix A: Simulations

We conducted data simulations using the TOASIM software suite, a plugin to TEMPO2. This allowed us to test the validity of the obtained orbital inclination (i) and longitude of the ascending node (Ω) values as presented in Table 3. We observed that the probabilities associated with (i, Ω) and their symmetric solutions (180° −i, 360° −Ω) depend on T₀, as shown in Fig. 6 and described in Sect. 4.3.1.

We used the formIdeal function in TOASIM to obtain idealised ToAs that have the exact same cadence (including the gaps in observing), backends, and jumps as the original data. To do this, we added Gaussian noise (using the addGaussian function) that is statistically equivalent to the noise found to be present in the actual data. To our simulated ToAs we also added stochastic achromatic red noise (using the addRedNoise function) with the same amplitude and spectral index as obtained from the actual data (computed from TEMPONEST and presented in Table 3). We executed these simulating steps for two input ephemerides, one ephemeris with its fiducial T₀ value at MJD 53004.13 (simulation 1), and the other with T₀ at MJD 58836.96 (simulation 2).

After creating these realistic datasets, we ran a TEMPONEST analysis on each, identical to what was done for real data in Sect. 2.5, including the same number of live points and other MULTINEST configurations. We performed the TEMPONEST analysis on both simulations for a range of associated input ephemerides with varying T₀ values: MJDs {53004.13, 55734.66 and 58831.96}. For simulation 2, we also ran it with two additional input T₀ values: MJD 54260.46 and MJD 57503.78. This was done to investigate whether we observe any trends in obtained (i, Ω) as a function of where we placed T₀ and to confirm whether there is any consequence when we place the T₀ values at the gaps in the PSR J1933−6211 timing baseline. The TEMPONEST input ephemerides for both simulated data sets had i and Ω set to 55.3 and 254 degrees, respectively. As for the real data, the sampling priors on i and Ω covered all possible values, that is, 0° ≤i ≤ 180° and 0° ≤Ω ≤ 360°.

The results of the simulations are shown in Fig. A.1. In none of these simulations do we find significant probabilities for the alternative solution (i > 90°, Ω < 180°), in contrast to what we find in the real data. This likely rules out our software as the main reason for this behaviour. With the simulations, we can confirm that for a data set that has the same cadence, noise properties, and residual rms timing precision, we are able to break the degeneracy for the angles i and Ω.

Fig. A.1. Orbital orientation for the PSR J1933−6211 binary as presented in Fig. 6, but using simulated data. In the top panel, the input ephemerides used to generate the fake data set had T0 set to MJD 53004.16 (simulation 1) and in the bottom panel to MJD 58836.96 (simulation 2). The TEMPONEST analyses were run with varying input ephemerides T0 values, as shown in the plot legends.

All Tables

All Figures

Fig. 1. Averaged flux-calibrated and RM-corrected PSR J1933−6211 pulse profile (bottom) and associated polarisation position angles (middle and top). Bottom: this averaged flux-calibrated and RM-corrected profile was obtained from adding all 24.5 h of MeerKAT L-band data, providing an S/N ≈ 3000, a peak flux density of 8.7 mJy, and a mean flux density of 1.1 mJy. The profile shows significant linear (red, L) and circular (blue, V) polarisation, the sharped-edged features of which make it an ideal candidate for MTM timing techniques when detecting it brightly. Middle and top: the linear polarisation PA are plotted for all phase bins with significant (> 5σ) linear polarisation. To visualise PA swings, we present the PA values, together with PA+180 deg and PA−180 deg (grey data points) in both panels. The solid (brown) line shows the best-fit RVM model when using a prior of ζ = 125 deg and obtaining α = 41.66 (4) deg (top), or a prior of ζ = 55 deg and obtaining α = 121.99 (4) deg (middle) as described in Sect. 3, with their ±180 deg equivalent as dashed (brown) lines. This fit was obtained using only the data points in black and removing the fainter data points of the centrally plotted PA values. The dotted lines (blue) in both panels are PA ± 90 deg to show the potential PA values of orthogonally polarised modes.

In the text

Fig. 2. Schematic of the 3D geometry for the PSR J1933−6211 binary in the observer’s convention. The orbital plane is defined by the unit vectors i and j, with the orbital angular momentum perpendicular to the plane, along the direction of k. The spin of the pulsar S is assumed to be aligned with the orbital angular momentum. The plane of the sky is defined by the east and north unit vectors, with n being the unit vector of our line of sight that is defined from the pulsar to the observer. The angle between n and k by definition is the orbital inclination angle, i. In the observer’s convention, this is equal to the polar angle of the pulsar ζ. The projection of S on the sky plane forms. The angle between north and i is the longitude of the ascending node. ρ denotes the opening angle of the emission cone, with β the impact angle of n on the cone, as measured from the magnetic axis μ. αp is the inclination of μ with respect to S.

In the text

Fig. 3. Output posterior distributions for PSR J1933−6211 parameters. These were obtained from TEMPONEST sampling using the T2 orbital model and a red-noise model, and were generated using the CHAINCONSUMER package. We present the orbital and post-Keplerian parameters. Full details are provided in Table 2. The obtained pulsar mass (Mp) distribution was computed using the mass function and the posterior distributions on Mc, i, xp, and Pb.

In the text

Fig. 4. Timing residuals across epochs (top) and orbital phase (bottom). Data include historic Parkes ToAs (green and purple; PKS Multibeam CPSR2 1341 MHz and 1405 MHz) recent and historic Parkes CASPSR data (red; PKS Multibeam CASPSR), as well as recent Parkes UWL/Medusa ToAs (blue; PKS UWL Medusa) and MeerKAT L-band ToAs (gold; MKT L/1K PTUSE). Top: after applying the best-fit values as obtained using the timing and noise model described in Sect. 2.5, we obtain a weighted rms of 1.259 μs. Bottom: residuals as a function of orbital phase, measured from superior conjunction (the ascending node is at a phase of −0.25). A Shapiro delay signal is discernible when setting the companion mass to zero (Mc = 0) while keeping all other parameters fixed, including xp sin i. We overplot the expected theoretical signal based on the best-fit inclination and companion mass values in black. The line width indicates combined 1σ deviations in companion mass and inclination angle.

In the text

Fig. 5. Constraints on the distance to PSR J1933−6211. The dash-dotted grey and blue lines indicate the constraints from the measurement of timing parallax without and with corrections for Lutz-Kelker bias using an EDSD prior, respectively. The dashed brown line indicates constraints from the kinematic contribution to the measured orbital period decay (${\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}$). The solid grey line shows the combined constraints from a simple inversion of the timing parallax and ${\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}$, with grey hatched lines marking the 68% C.L. This corresponds to a distance of 1.4 (2) kpc. The solid blue line shows the combined constraints from the timing parallax with the EDSD prior and ${\dot{P}}_{\mathrm{b}}^{\mathrm{kin}}$, with the solid blue shaded region denoting its 68% C.L. This corresponds to a distance of $1.6_{-0.3}^{+0.2}$ kpc.

In the text

Fig. 6. Allowable orbital orientation for PSR J1933−6211 binary. In the main panel, we display the orbital orientation space in the observer’s convention (the Ω according to the DT92 convention is shown by the scale on the right, the cos i in that convention is simply −cos i). The lines represent constraints from two parameters measured with a DDH solution: in black we show the constraint from ς, and in brown we show the constraint from ẋ (1σ solid, 3σ dashed) assuming it is solely caused by the proper motion of the system; the position angle of the latter is shown by the dot-dashed purple line. The contours, enclosing progressively darker shades of red, include 98%, 86% and 39% C.L. of the 2D probability density function (pdf) and represent the T2 timing model solutions calculated by TEMPONEST and as presented in Table 3. This model takes into account all kinematic effects, including the AOP. The marginalised constraints on cos i and Ω are shown as 1D histograms (in red) in the top and side panels, respectively. For this model, with T0 ∼ 53 004, the AOP is clearly detected, as we can see from the fact that the degeneracy between the two possible cos i − Ω solutions has been lifted: no probability remains in the solution at Ω ∼ 60 deg, cos i ∼ −0.6. We also include the cos i − Ω probability density contours as obtained with T0 set to 55 734 (dashed grey) and 58 836 (dot-dashed blue), respectively. As described in Sect. 4.3.1, we unexpectedly observe i and Ω estimates to depend on T0.

In the text

	Fig. A.1. Orbital orientation for the PSR J1933−6211 binary as presented in Fig. 6, but using simulated data. In the top panel, the input ephemerides used to generate the fake data set had T0 set to MJD 53004.16 (simulation 1) and in the bottom panel to MJD 58836.96 (simulation 2). The TEMPONEST analyses were run with varying input ephemerides T0 values, as shown in the plot legends.
In the text