# Radio timing constraints on the mass of the binary pulsar PSR J1528-3146

A Berthereau, L Guillemot, P.C.C Freire, M Kramer, V Venkatraman

Krishnan, I Cognard, G Theureau, M Bailes, M.C Bernadich, M.E Lower

► To cite this version:

A Berthereau, L Guillemot, P.C.C Freire, M Kramer, V Venkatraman Krishnan, et al., Radio timing constraints on the mass of the binary pulsar PSR J1528-3146. Astronomy and Astrophysics - A&A, 2023, 674, pp.A71. 10.1051/0004-6361/202346228. hal-04085015

HAL Id: hal-04085015

https://hal.science/hal-04085015

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radio timing constraints on the mass of the binary pulsar PSR J1528-3146

A. Berthereau<fnref n="1" />,<fnref n="2" />, L. Guillemot<fnref n="1" />,<fnref n="2" />, P. C. C. Freire<fnref n="3" />, M. Kramer<fnref n="3" />, V. Venkatraman Krishnan<fnref n="3" />, I. Cognard<fnref n="1" />,<fnref n="2" />, G. Theureau<fnref n="1" />,<fnref n="2" />,<fnref n="4" />, M. Bailes<fnref n="5" />,<fnref n="6" />, M. C. i Bernadich<fnref n="3" />, and M. E. Lower<fnref n="7" />

Received 23 February 2023 / Accepted 12 April 2023

## ABSTRACT

*Context.* PSR J1528-3146 is a 60.8 ms pulsar orbiting a heavy white dwarf (WD) companion, with an orbital period of 3.18 d. The pulsar was discovered in the early 2000s in a survey at 1.4 GHz of intermediate Galactic latitudes conducted with the Parkes radio telescope. The initial timing analysis of PSR J1528-3146, using data recorded from 2001 and 2004, did not reveal any relativistic perturbations to the orbit of the pulsar or to the propagation of its pulses. However, with an orbital eccentricity of ~0.0002 and a large companion mass on the order of $1 M_\odot$, this system has been deemed likely to exhibit measurable perturbations.

*Aims.* This work is aimed at characterizing the pulsar's astrometric, spin, and orbital parameters by analyzing timing measurements conducted at the Parkes, MeerKAT, and Nançay radio telescopes over nearly two decades. The measurement of post-Keplerian perturbations to the pulsar's orbit can be used to constrain the masses of the two component stars of the binary and, in turn, to offer insights into the history of the system.

*Methods.* We analyzed timing data from the Parkes, MeerKAT, and Nançay radio telescopes collected over about 16 yr, obtaining a precise rotation ephemeris for PSR J1528-3146. A Bayesian analysis of the timing data was carried out to constrain the masses of the two components and the orientation of the orbit. We further analyzed the polarization properties of the pulsar to constrain the orientation of the magnetic axis and of the line of sight with respect to the spin axis.

*Results.* We measured a significant rate of advance of periastron, for the first time, and we set constraints on the Shapiro delay in the system and on the rate of change of the projected semi-major axis of the pulsar's orbit. The Bayesian analysis yielded measurements for the pulsar and companion masses of $M_p = 1.61^{+0.14}_{-0.13} M_\odot$ and $M_c = 1.33^{+0.08}_{-0.07} M_\odot$ (68% C.L.), respectively, confirming that the companion is indeed massive. This companion mass as well as other characteristics of PSR J1528-3146 indicate that this pulsar is very similar to PSR J2222-0137, a 32.8 ms pulsar orbiting a WD whose heavy mass (~1.32 $M_\odot$) has been considered unique among pulsar-WD systems until now. Our measurements suggest common evolutionary scenarios for PSRs J1528-3146 and J2222-0137.

**Key words.** ephemerides – pulsars: individual: J1528–3146

## 1. Introduction

Pulsars are highly magnetized, rapidly rotating neutron stars that emit periodic trains of radio pulses. The pulsar population can be divided into two main groups: “normal” pulsars and “millisecond” pulsars (MSPs). The latter are defined as pulsars with small rotational periods ($P \le 30$ ms) and very small period increases ($\dot{P} \le 10^{-17}$ s s⁻<fnref n="1" />). In most cases, MSPs are found to be members of binary systems, and are thought to have been “recycled” by the accretion of matter and transfer of angular momentum from their binary companion, spinning up their rotation to millisecond periods (Bisnovatyi-Kogan & Komberg 1974; Alpar et al. 1982). In some cases, this accretion process can stop before the pulsar gets fully recycled, leading to so-called “mildly recycled” pulsars with rotational periods between 20 and 100 ms. This process happens mainly if the companion stars are more massive:

such stars evolve more rapidly and therefore any accretion episodes will generally be much shorter.

PSR J1528-3146 is a member of the latter sub-class of mildly-recycled pulsars. This $P \sim 60.8$ ms pulsar was discovered in a survey conducted at 1.4 GHz with the Parkes radio telescope (also known as ‘Murriyang’) in Australia, as part of an extension of the Swinburne Intermediate Latitude survey, which discovered a total of 26 pulsars (Jacoby et al. 2007). Timing measurements conducted from 2001 to 2004 revealed that the pulsar is in a binary system with an orbital period of $P_b \sim 3.18$ days, a projected semi-major axis $x \equiv a \sin i/c \sim 11.4$ s (where $a$ is the semi-major axis of the pulsar’s orbit, $i$ is the orbital inclination, and $c$ is the speed of light), and an orbital eccentricity $e \sim 0.0002$ which is large for systems with comparable orbital periods (e.g., Berezina et al. 2017). The orbital period and semi-major axis imply a relatively massive companion: assuming a

<footnote n="1">¹ Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans/CNRS, 45071 Orléans Cedex 02, France
e-mail: anais.berthereau@cnrs-orleans.fr</footnote>

<footnote n="2">² Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France</footnote>

<footnote n="3">³ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="4">⁴ LUTH, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92195 Meudon, France</footnote>

<footnote n="5">⁵ Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, Vic 3122, Australia</footnote>

<footnote n="6">⁶ ARC Centre of Excellence for Gravitational Wave Discovery, OzGrav, Mail H29, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="7">⁷ Australia Telescope National Facility, CSIRO, Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia</footnote>
pulsar mass of $1.4 M_\odot$, the mass function of $0.159 M_\odot$ indicates that the minimum companion mass (calculated for an edge-on orbit, i.e., with an inclination of $90^\circ$) is ~0.94 $M_\odot$. With such properties the system can be classified as an intermediate-mass binary pulsar (IMBP, see e.g., Camilo 1996; Tauris et al. 2012); the large mass of the companion suggests that it must be a CO or ONeMg white dwarf (WD). Optical observations of the field around PSR J1528-3146 with the 6.5 m Baade telescope at the Magellan Observatory led to the detection of a potential counterpart, with an R-band magnitude of about 24.2 and a B-band magnitude of about 24.5 (Jacoby et al. 2006). If the counterpart is indeed associated with the J1528-3146 system, then the optical emission very likely originates from the WD.

The initial timing study by Jacoby et al. (2007) over approximately three years did not reveal any relativistic perturbations to the pulsar’s orbit, such as the rate of periastron advance. Nevertheless, the non-negligible orbital eccentricity and large companion mass make it a promising system for measuring such relativistic perturbations, that can be described in pulsar timing analyses by so-called “post-Keplerian” (PK) parameters (e.g., Damour & Taylor 1992).

Two decades after the discovery of PSR J1528-3146, the accumulation of high-quality pulsar timing data by the Nançay radio telescope (NRT) in France and the advent of highly sensitive new-generation instruments such as the MeerKAT 64-dish array in South Africa (Jonas 2009) and the ultra-wideband low receiver at the Parkes telescope (UWL, see Hobbs et al. 2020) motivated a new study of this pulsar and of its timing properties. In particular, a new timing study of PSR J1528-3146 should enable more precise measurements of its parameters and could yield first detections of PK parameters in this system and possibly measurements of the masses of the pulsar and of its companion.

In this article we report on the results from the analysis of a ~16 yr timing data set on PSR J1528-3146 with the Parkes, NRT, and MeerKAT radio telescopes. The article is organized as follows: in Sect. 2, we describe the details of the observations of PSR J1528-3146 and the extraction of pulsar timing data from these observations. In Sect. 3, we present an analysis of the linear polarization properties of the pulsar and the implications of the results in terms of emission geometry. In Sect. 4, we present the main timing results from our study and the constraints on the mass of the pulsar and of its companion. We finally summarize our results and discuss on the main implications from this study in Sect. 5.

## 2. Data

### 2.1. Observations

The NRT dataset consists of observations carried out with the Berkeley-Orléans-Nançay (BON) backend until August 2011 and the NUPPI (a version of the Green Bank Ultimate Pulsar Processing Instrument, GUPPI, made for the NRT) backend afterwards.

The BON observations of PSR J1528-3146 started in November 2006, and were carried out at a central frequency of 1398 MHz. These observations were coherently dedispersed and carried out over a frequency bandwidth of 64 MHz until July 2008, and over a bandwidth of 128 MHz after that date.

The NUPPI instrument, which became the primary pulsar timing instrument at Nançay after August 2011, increased the bandwidth of NRT pulsar observations to 512 MHz that are also coherently dedispersed. NUPPI observations were conducted at

a central frequency of 1484 MHz. In this study we considered NUPPI data recorded until May 2022. Although pulsar observations with the BON backend were maintained for another few years after 2011 for testing purposes, we discarded BON observations made simultaneously with NUPPI observations, as they represented duplicated astrophysical information.

The Parkes radio telescope dataset included the timing data presented in Jacoby et al. (2007), recorded between September 2001 and September 2004, and more recent data taken with the UWL recorded from October 2019 to January 2021. The former observations were conducted with the 21-cm multibeam receiver (Staveley-Smith et al. 1996) at a central frequency of 1390 MHz, using the $2 \times 0.5 \times 512$ MHz analogue filterbank backend (AFB) and were thus incoherently dedispersed. Observations with the UWL receiver were carried out using the MEDUSA backend, that performed coherent dedispersion over a much larger bandwidth of 3328 MHz, centered at 2368 MHz.

The MeerKAT dataset we used for the study is comprised of observations made between March 2019 and April 2021, as part of the RelBin program of the MeerTime large science project. The RelBin program aims at improving the measurement and detecting new relativistic effects in known pulsars, increasing the number of pulsars with mass measurements, and testing gravity theories, with the excellent timing capabilities of MeerKAT (Kramer et al. 2021a). Pulsar observing set-ups, verifications and early science results with MeerKAT are presented in Bailes et al. (2020). Observations were calibrated in polarization as described in Serylak et al. (2021). MeerKAT observations were conducted at a central frequency of 1284 MHz, and covered a frequency bandwidth of about 800 MHz. While the bulk of the observations had a typical duration of about 30 min, an intensive orbital campaign was conducted around 2019 December 12 (MJD 58829). On that day, the system was observed for about four hours to cover a wide range of orbital phases centered on superior conjunction in order to significantly increase the sensitivity of our subsequent timing analyses to a faint Shapiro delay (Shapiro 1964) signal that had already been detected in Nançay data. Additionally, three other 30 min observations of the pulsar were conducted within three days centered on MJD 58829, in order to sample orbital phase ranges away from superior conjunction.

In Fig. 1 we show the MeerKAT pulse profile for PSR J1528-3146 from the latter observation, which represented the highest S/N observation in our total dataset. The main radio component in the pulse profile is at rotational phase ~0.54 in this graph. As can be seen from Fig. 1, the unprecedented sensitivity of MeerKAT observations has enabled us to detect a previously unseen weak component separated from the main pulse by approximately half a rotation. The presence of this interpulse component provides useful information about the pulsar's emission geometry, as discussed further in Sect. 3.

### 2.2. Data analysis

All the data reduction steps for the datasets described above were done using the PSRCHIVE analysis library (Hotan et al. 2004). Observations were cleaned of radio frequency interference (RFI) and also polarization-calibrated.

The typically ~1 h long Nançay observations with the BON and NUPPI backends were split into sub-integrations of 10 min each. For the BON dataset, these sub-integrations were split into two frequency sub-bands of 32 MHz each for the data recorded before July 2008, and two sub-bands of 64 MHz for the data taken after the bandwidth increase. Sub-integrations for
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** MeerKAT pulse profiles at 1284 MHz for PSR J1528−3146, from a ~4-h observation at MJD 58829, corresponding to the highest S/N observation in our dataset. Left: total intensity as a function of observing frequency and rotational phase. Right: integrated pulse profile, with the main radio pulse centered at rotational phase ~0.54 and a previously-unknown weak radio component at phase ~0.03, indicated with an arrow. The inset at the top left shows a zoom-in of the pulse profile around the weak radio component.</figcaption></figure>

<table><caption>**Table 1.** Summary of observation parameters and TOA properties for each of the datasets considered in the timing study.</caption><thead><tr><th>Dataset</th><th>Parkes AFB</th><th>Parkes UWL</th><th>MeerKAT</th><th>Nançay BON</th><th>Nançay NUPPI</th></tr></thead><tbody><tr><td>Start of observations (MJD)</td><td>52 179</td><td>58 759</td><td>58 557</td><td>54 065</td><td>55 848</td></tr><tr><td>End of observations (MJD)</td><td>53 262</td><td>59 216</td><td>59 319</td><td>55 819</td><td>59 725</td></tr><tr><td>Observing frequency (MHz)</td><td>1398</td><td>2368</td><td>1284</td><td>1398</td><td>1484</td></tr><tr><td>Frequency bandwidth (MHz)</td><td>256</td><td>3328</td><td>776/856 <sup>(c)</sup></td><td>64/128 <sup>(d)</sup></td><td>512</td></tr><tr><td>Number of TOAs</td><td>169</td><td>53</td><td>936</td><td>866</td><td>11 454</td></tr><tr><td>Bandwidth per TOA (MHz)</td><td>128</td><td>3328 <sup>(b)</sup></td><td>97/107 <sup>(c)</sup></td><td>32/64 <sup>(d)</sup></td><td>64</td></tr><tr><td>Time per TOA (s)</td><td>— <sup>(a)</sup></td><td>600</td><td>600</td><td>600</td><td>600</td></tr><tr><td>Median TOA uncertainty (μs)</td><td>19.3</td><td>3.2</td><td>3.5</td><td>30.4</td><td>29.0</td></tr><tr><td>Weighted residual rms (μs)</td><td>12.8</td><td>4.3</td><td>3.0</td><td>16.3</td><td>16.3</td></tr><tr><td>EFAC</td><td>1.13</td><td>1.69</td><td>1.22</td><td>0.98</td><td>0.98</td></tr></tbody></table>

**Notes.** (a) These observations were crunched in time to single sub-integrations. See Jacoby et al. (2007) for details on the corresponding observations. (b) For this dataset the TOAs were generated using the wide-band template matching technique implemented in *PulsePortraiture* (see Pennucci et al. 2014). (c) The bulk of the MeerKAT observations covered a total of 856 MHz of bandwidth. For these observations, we grouped the channels into eight sub-bands of 107 MHz each. Other observations covered a slightly shorter bandwidth of 776 MHz – for these, we also grouped the channels in eight sub-bands of 97 MHz. (d) The original frequency bandwidth for the BON backend of 64 MHz increased to 128 MHz in July 2008. For observations taken prior to July 2008, the bandwidth of 64 MHz was split into 32 MHz sub-bands; for the other observations, the bandwidth of 128 MHz was split into 64 MHz sub-bands.

the NUPPI data were divided into eight frequency sub-bands of 64 MHz each. For the BON and NUPPI data, we combined 10 observations with the highest signal-to-noise ratios (S/N) that we then smoothed to produce noise-free standard profiles at 1.4 GHz with the same frequency resolution. Times of arrival (TOAs) to be used in the timing analysis (described in Sect. 4) were derived from the Nançay data using the Matrix Template Matching (MTM) method of the pat routine of PSRCHIVE (van Straten 2006).

The TOAs were extracted from the “historic” Parkes dataset by comparing observations with a high S/N template profile, as described in Jacoby et al. (2007). Each observation was split into two frequency sub-bands of 128 MHz, so that two TOAs were extracted from each observation. On the other hand, the very large frequency bandwidth covered by Parkes observations with the UWL prompted us to consider a different approach for extracting TOAs. Instead of splitting the very large bandwidth into several sub-bands (which often results in sub-bands

containing no clear detections at the highest frequencies), we used the wide-band template matching technique implemented in the *PulsePortraiture* library<fnref n="1" /> (Pennucci et al. 2014) and extracted one TOA for the entire bandwidth for every 10 min of observation.

Finally, the MeerKAT observations were divided into eight frequency sub-bands 97 MHz (resp. 107 MHz) for observations with a 759 MHz (resp. 856 MHz) bandwidth and the sub-integrations were grouped into segments of 10 min. For each resulting element of time and frequency, we extracted a TOA using the standard Fourier Domain Markov chain Monte Carlo (FDM) method of pat, comparing the profiles with a smoothed version of the highest S/N MeerKAT observation with the same number of frequency sub-bands. A summary of the TOA dataset produced for the timing analysis is presented in Table 1.

<footnote n="1">¹ https://github.com/pennucci/PulsePortraiture</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Polar representation of the MeerKAT pulse profile at 1284 MHz displayed in Fig. 1, with radio flux densities shown in logarithmic scale. The black dashed line shows the position of the interpulse component and the red dashed line shows the same position shifted by 180°.</figcaption></figure>

### 3. Pulse polarization analysis

As mentioned in Sect. 2, the unprecedented sensitivity of MeerKAT radio observations of PSR J1528–3146 enabled us to detect a weak profile component that is well separated from the main emission pulse, for the first time. This newly detected weak component can be seen in the polar representation of the pulse profile from the MeerKAT observation at MJD 58829 shown in Fig. 2. The dashed lines show the position of this profile component and that of the same position shifted by half a rotation. The phase location was determined from a fit of the pulse profile displayed in Fig. 1 with a Lorentzian line for the highest peak and two Gaussian lines for the other two components, on top of a constant baseline. Denoting $\phi_{M,1}$ and $\phi_{M,2}$ as the phases of the two components of the main pulse, and $\phi_{IP}$ that of the weak component, we found $\phi_{M,1} = 0.469(2)$, $\phi_{M,2} = 0.54410(2)$, and $\phi_{IP} = 0.033(8)$ from the fit of the pulse profile. The new component is therefore separated from the main emission cone by almost exactly half a rotation and it can thus be considered as interpulse emission, likely originating from the magnetic pole opposite to the one responsible for the main emission pulse. The presence of this interpulse component indicates that the angle between the pulsar’s spin axis and its magnetic axis, $\alpha$, must be close to 90°.

Tighter constraints on the value of $\alpha$ and that of $\zeta$, namely, the angle between the pulsar’s spin axis and the observer’s line of sight, may be placed by fitting the variation of the polarization position angle (PPA) $\Psi$ of the linearly polarized component of the pulsar as a function of the rotational phase, with a rotating vector model (RVM, Radhakrishnan & Cooke 1969). The position angles $\Psi$ are related to the Q and U Stokes parameters of the pulsar’s emission as $\Psi = \frac{1}{2} \arctan(\frac{U}{Q})$. In the following, we use the modified form introduced in Johnston & Kramer (2019):

$$ \Psi = \Psi_0 + \arctan \frac{\sin \alpha \sin (\phi - \phi_0 - \Delta)}{\sin \zeta \cos \alpha - \cos \zeta \sin \alpha \cos (\phi - \phi_0 - \Delta)}. \quad (1) $$

In the above expression, $\phi_0$ denotes the pulse phase for which $\Psi = \Psi_0$ (assuming $\Delta = 0$), and the term $\Delta$ can take account of a potential difference in emission heights between the two magnetic poles of the pulsar. We used PPA values as a function of rotational phase extracted from the 1284 MHz MeerKAT observation at MJD 58829. We first formed a time-averaged and dispersed version of the observation, keeping the initial frequency resolution of 928 channels over the total bandwidth of 776 MHz. We used the “Bayesian Lomb-Scargle Periodogram” (BGLSP) method of RMcalc<fnref n="2" /> to fit for the pulsar’s rotation measure, RM, which quantifies the amount of Faraday rotation occurring along the line of sight caused by the presence of magnetic field in the interstellar medium. Denoting $\lambda$ as the radio wavelength, we have:

$$ \Delta\Psi_{\text{PPA}} = \text{RM} \times \lambda^2. \quad (2) $$

Our BGLSP analysis yielded an RM value of –19 (1) rad m⁻¹, consistent with the value published in Kramer et al. (2021a). We finally corrected the values of $\Psi$ across the frequency bandwidth of the MeerKAT observation for the above-described Faraday rotation term and built a frequency-averaged pulse profile for PSR J1528–3146 with the four Stokes parameters. Using this profile, we computed the linearly and circularly polarized intensity as shown as red and blue lines in Fig. 3a (top panel), together with the total intensity as a black line. We then performed several RVM fits to the derived PPA data, shown in the bottom panel of Fig. 3a. The fits were obtained using the Ultra-Nest sampler by Buchner (2021), using different priors for the angles in different runs. In each case, we used a standard expression, $\ln(\mathcal{L}) = -\chi^2/2$, as a likelihood function. In the first fit, we assumed our random $\alpha$ and $\zeta$ angles to form spin and magnetic field vectors that are distributed uniformly over a unit-sphere, while $\Delta = 0$. The prior for $\Psi_0$ was a uniform distribution over [–90°, 90°] and, similarly, $\phi_0$ had a uniform prior over [0, 1]. A corner plot of the obtained posterior distributions of the joint RVM model parameters, with the off-diagonal elements representing the correlations between parameters and the diagonal elements denoting the marginalized histograms is shown in Fig. 3b. The PPA values of the discovered interpulse component significantly help to constrain the derived parameter range. In addition to a solution at $(\alpha, \zeta) = (\sim 103°, \sim 92°)$, a second solution (less favoured by the log-likelihood statistic) is seen for $(\alpha, \zeta) = (\sim 80°, \sim 72°)$. Either case suggests an orthogonal rotor, as expected when an interpulse is present. The solid black line in the lower panel of Fig. 3a is the derived RVM corresponding to the preferred solution (with the gray band indicating 68% uncertainties). The dashed black lines show the RVM shifted by ±90°. If this fit is accurate, it indeed suggests that the main pulse and interpulse are emitting in different orthogonal emission states, as is found quite frequently in interpulse-pulsars (see e.g., Johnston & Kramer 2019).

As we discuss further in Sect. 4.3, the Bayesian analysis of the timing data suggests an orbital inclination angle of either $i \sim 57°$ or $i \sim 123°$. For a fully recycled pulsar, we can expect the spin axis of the pulsar to be aligned with the orbital angular momentum vector, so that $i \approx \zeta$, and also $\Psi_0 \approx \Omega$, where $\Omega$ is the longitude of the ascending node (see Sect. 4.3). The usage of polarization data and RVM fits to determine the orbital configuration was discussed recently, for instance, by Kramer et al. (2021a) or Guo et al. (2021). Clearly, the above $\zeta$-values are not consistent with the results from timing. In a second and third fit, we therefore applied uniform priors for $\zeta$ restricted to a

<footnote n="2">² https://gitlab.mpifr-bonn.mpg.de/nporayko/RMcalc</footnote>
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Rotating vector model analysis of PSR J1528−3146. Left: top panel shows the polarization profile of PSR J1528−3146 as observed at 1284 MHz. The black line represents total power, the red shows linearly and the blue line circularly polarized intensity, respectively. Bottom panel displays the derived PPA values and a number of RVM fits. The black lines represents the resulting RVM for an unconstrained fit covering the whole parameter space. The dashed lines indicate vertical shifts of the RVM by ±90° to represent possible orthogonal polarized modes. The red and blue curves are RVM fits for a restricted range of ζ, forcing it to be consistent with values of the orbital inclination angle as derived from timing in this paper. The magenta line displays a possible shift in the emission height between the main pulse and the interpulse. See text for details. Right: corner plot of the posterior distributions of the joint RVM model parameters, with the off-diagonal elements representing the correlations between parameters and the diagonal elements denoting the marginalized histograms. The posteriors shown here are derived for the unconstrained fit leading to the black line on the left.

6° interval centred on the respective inclination angle values. The red curves display the results for ζ ~ 57° (α = 69(1)°, ζ = 59(1)°, Ψ₀ = -23(1)°), while the blue lines are for ζ ~ 123° (α = 128(1)°, ζ = 120(1)°, and Ψ₀ = -10(1)°), respectively. The latter case, results in an “outer” line of sight, whereas the former is an “inner” one (see e.g., Lorimer & Kramer 2012). It is clear that the red curve, ζ ~ 57°, is a better fit to the data, but it is still far from ideal.

We can improve the fit by allowing for Δ ≠ 0. This was done in our last fit, the result of which is shown as the magenta line in the PPA panel, resulting in Δ = 0.07 (2) and a significant difference in emission height: (α = 67 (2)°, ζ = 59 (2)°, Ψ₀ = -25 (1)°). Regardless of a non-zero Δ or otherwise, overall, the solution for ζ ~ 57° describes the data much better than the solution with i = 123°. It is interesting to also study the values for Ω derived from timing (see Sect. 4 and Table 2) and listed in Table 3. Our value derived for an RVM with ζ close to 57°, Ψ₀ = -23° is consistent with a value for Ω in the range derived for Solution 4 in Table 3. Nevertheless, the overall quality of the fit when forcing ζ to lie close to the derived orbital inclination angle is clearly limited. An unconstrained fit prefers a value much closer to 90°. As shown in Sect. 4.3, an orbital inclination value close to 90° is excluded by the timing analysis of PSR J1528−3146. Whether the applicability of the RVM is at fault here or whether the assumption of aligned spin-vectors is wrong, these are interesting questions that merit further investigation and discussion.

## 4. Timing analysis

The timing analysis was carried out by analyzing the topocentric TOA data presented in Sect. 2, using the TEMPO2 pulsar timing package (Hobbs et al. 2006). The measured TOAs

were first converted to Barycentric Coordinate Time (TCB), taking the known clock corrections for the different telescopes and using the JPL DE438 solar system ephemeris<fnref n="3" />. The orbital motion of the pulsar was described using the DDH orbital model, which uses the orthometric parameterization of the Shapiro delay (Freire & Wex 2010). This parameterization is well-suited for systems with low orbital inclinations (which is the case of PSR J1528−3146, as shown later in this section). The reason is that it describes the Shapiro delay using two parameters: the orthometric amplitude ($h_3$) and the orthometric ratio ($\zeta$), which are, in such systems, much less correlated than the range ($r$) and shape ($s$) parameters of the DD orbital model (Damour & Deruelle 1985, 1986).

We used TEMPO2 to fit for the pulsar’s sky position and proper motion, its spin parameters (spin frequency and first and second time derivatives), orbital parameters, and PK parameters describing perturbations to the classical Keplerian orbit. The TOAs from the different datasets were aligned with each other using time offsets, which were fitted for as part of the timing analysis. The dispersion measure (DM) and its variations were modeled as a piecewise-constant function, using the “DMX” model. The TOA dataset was divided into segments of 100 days, and constant DM values were fitted in each of the intervals containing timing data. In a first iteration of the analysis, we obtained a preliminary timing model for PSR J1528−3146 which we then used to determine “EFAC” factors for each of the individual datasets, namely, scaling factors applied to the TOA uncertainties such that the individual reduced $\chi^2$ values become equal to unity. The EFAC factors we determined through this

<footnote n="3">³ https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de438s.bsp.lbl</footnote>
<table><caption>**Table 2.** Parameters for PSR J1528−3146 from the timing analysis, and derived quantities.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td colspan="2">Observation and data reduction parameters</td></tr><tr><td>Span of timing data (MJD)</td><td>52179-59725</td></tr><tr><td>Number of TOAs</td><td>13 478</td></tr><tr><td>Reference epoch (MJD)</td><td>56 000</td></tr><tr><td>Time scale</td><td>TCB</td></tr><tr><td>Solar system ephemeris</td><td>DE438</td></tr><tr><td>Solar wind electron number density, n<sub>0</sub> (cm<sup>-3</sup>)</td><td>4</td></tr><tr><td>Weighted rms residual (µs)</td><td>9.245</td></tr><tr><td>Reduced χ<sup>2</sup></td><td>0.993</td></tr><tr><td colspan="2">Astrometric, spin and dispersion parameters</td></tr><tr><td>Right ascension, α (J2000)</td><td>15:28:34.95312 (6)</td></tr><tr><td>Declination, δ (J2000)</td><td>-31:46:06.866 (2)</td></tr><tr><td>Proper motion in α, µ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-1.39 (9)</td></tr><tr><td>Proper motion in δ, µ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-3.7 (2)</td></tr><tr><td>Spin frequency, ν (Hz)</td><td>16.4413566796379 (11)</td></tr><tr><td>Spin frequency first derivative, ˙ν (10<sup>-17</sup> Hz s<sup>-1</sup>)</td><td>-6.7175 (9)</td></tr><tr><td>Spin frequency second derivative, ˙˙ν (10<sup>-28</sup> Hz s<sup>-2</sup>)</td><td>2.8 (6)</td></tr><tr><td>Dispersion measure, DM (cm<sup>-3</sup> pc)</td><td>18.159 <sup>(a)</sup></td></tr><tr><td>Rotation measure, RM (rad m<sup>-2</sup>)</td><td>-19 (1) <sup>(b)</sup></td></tr><tr><td colspan="2">Binary parameters</td></tr><tr><td>Orbital model</td><td>DDH</td></tr><tr><td>Orbital period, P<sub>b</sub> (days)</td><td>3.1803502 (2)</td></tr><tr><td>Projected semi-major axis of the pulsar orbit, x (lt-s)</td><td>11.452301 (11)</td></tr><tr><td>Epoch of periastron, T<sub>0</sub> (MJD)</td><td>56000.7860 (3)</td></tr><tr><td>Orbital eccentricity, e</td><td>0.0002137 (2)</td></tr><tr><td>Longitude of periastron, ω (°)</td><td>297.30 (4)</td></tr><tr><td>Rate of advance of periastron, ω̇ (° yr<sup>-1</sup>)</td><td>0.057 (3)</td></tr><tr><td>Rate of change of orbital semi-major axis, ˙x (10<sup>-15</sup> lt-s s<sup>-1</sup>)</td><td>-2.6 (1.9)</td></tr><tr><td>Orthometric amplitude of the Shapiro delay, h<sub>3</sub> (µs)</td><td>1.2 (1)</td></tr><tr><td>Orthometric ratio of the Shapiro delay, ζ</td><td>0.44 (8)</td></tr><tr><td colspan="2">Derived parameters</td></tr><tr><td>Galactic longitude, l (°)</td><td>337.9386445</td></tr><tr><td>Galactic latitude, b (°)</td><td>20.2167420</td></tr><tr><td>Total transverse proper motion, µ<sub>T</sub> (mas yr<sup>-1</sup>)</td><td>3.9 (2)</td></tr><tr><td>Position angle of proper motion (J2000), Θ<sub>μ</sub> (°)</td><td>200.8 (1.8)</td></tr><tr><td>DM-derived distance (NE2001), d (kpc)</td><td>0.8(2) <sup>(c)</sup></td></tr><tr><td>DM-derived distance (YWM16), d (kpc)</td><td>0.8(2) <sup>(c)</sup></td></tr><tr><td>Mass function, f (M<sub>☉</sub>)</td><td>0.1594451 (5)</td></tr><tr><td>Spin period, P (ms)</td><td>60.822231369658 (4)</td></tr><tr><td>Apparent spin period derivative, ˙P (10<sup>-19</sup> s s<sup>-1</sup>)</td><td>2.4850 (3)</td></tr><tr><td>Total kinematic contribution to spin period derivative, ΔP (10<sup>-21</sup> s s<sup>-1</sup>)</td><td>1.2 (3) <sup>(d)</sup></td></tr><tr><td>Intrinsic spin period derivative (Ṗ<sub>int</sub>) (10<sup>-19</sup> s s<sup>-1</sup>)</td><td>2.473 (3)</td></tr><tr><td>Spin-down power, Ė (10<sup>31</sup> erg s<sup>-1</sup>)</td><td>4.340 (5) <sup>(e)</sup></td></tr><tr><td>Surface magnetic field strength, B<sub>surf</sub> (10<sup>9</sup> G)</td><td>3.928 (2) <sup>(e)</sup></td></tr><tr><td>Characteristic age, τ<sub>c</sub> (10<sup>9</sup> yr)</td><td>3.896 (5) <sup>(e)</sup></td></tr><tr><td>Pulsar mass, M<sub>p</sub> (M<sub>☉</sub>)</td><td>1.61<sup>+0.14</sup> <sup>(f)</sup><sub>-0.13</sub></td></tr><tr><td>Companion mass, M<sub>c</sub> (M<sub>☉</sub>)</td><td>1.33<sup>+0.08</sup> <sup>(f)</sup><sub>-0.07</sub></td></tr></tbody></table>

**Notes.** Numbers in parentheses are the nominal 1σ statistical uncertainties on the last quoted digits. <sup>(a)</sup>DM variations were modeled using piecewise-constant deviations around this value, as enabled by the “DMX” model of Tempo2. <sup>(b)</sup>Obtained using the RMcalc program, see Sect. 3 for details. <sup>(c)</sup>Distances estimated using the NE2001 (Cordes & Lazio 2002) and YMW16 (Yao et al. 2017) models, assuming a 20% uncertainty. <sup>(d)</sup>Contributions from the Shklovskii effect (Shklovskii 1970) and from the differential Galactic acceleration, as calculated using the GalPot model (McMillan 2017; Dehnen & Binney 1998). The uncertainty quoted here corresponds to the maximum excursion from the central value, when varying the pulsar distance by ±20%. <sup>(e)</sup>See Lorimer & Kramer (2012) for definitions. <sup>(f)</sup>Determined from a Bayesian analysis of the timing data, see Sect. 4.3 for details.

procedure are listed in Table 1. As can be seen from the table, the EFACs we obtained were all close to 1. A final iteration of the timing analysis with scaled TOA uncertainties yielded the parameters for PSR J1528−3146 listed in Table 2. The differences between measured TOAs and those predicted by the

best-fit timing model (the so-called timing residuals) are presented in Figs. 4 and 5. The residuals are normally-distributed around 0, demonstrating that the parameters in Table 2 describe the timing data appropriately. Figure 6 shows the best-fit piecewise-constant DM values obtained from the timing analysis.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Timing residuals as a function of time (top) and orbital phase (bottom), with the datasets considered in our analysis shifted by different amounts for ease of comparison. In the bottom plot, the orbital phase corresponding to superior conjunction (i.e., the moment of the orbit when the pulsar is on the opposite side of the companion) was shifted to occur at 0.25.</figcaption></figure>

The DM displays long-term variations around the value of 18.163 (6) measured by Jacoby et al. (2007).

The analysis of ~16 yr of timing data, which includes data of unprecedented quality from the MeerKAT telescope, enabled us to measure the pulsar's proper motion and several PK parameters for the first time. In the following subsections, we present these measurements in detail and discuss their implications with respect to PSR J1528-3146.

## 4.1. Proper motion, distance and spin parameters

The timing analysis yielded a precise proper motion measurement. We find the total transverse proper motion to be 3.9 (2) mas yr⁻¹, with a position angle of Θµ = 200.8 (1.8)°. Given that 0° indicates north in the “observer’s” convention we used, the position angle of the proper motion of PSR J1528-3146 indicates a southward motion. Our analysis did not enable us to measure the timing parallax, which otherwise provides a direct constraint on the distance. However, the NE2001 (Cordes & Lazio 2002) and YMW16 (Yao et al. 2017)

Galactic free-electron density models both place the pulsar at a distance of d = 0.8 (2) kpc, assuming conservative uncertainty levels of 20%; we use this distance in the estimates given below.

With these proper motion and distance values, we can estimate the intrinsic spin-down rate of the pulsar, Pint, as:

$$ \frac{\dot{P}_{\text{int}}}{P} = \frac{\dot{P}}{P} - \frac{\mu_{\text{T}}^2 d}{c} - \frac{a}{c}, \qquad (3) $$

where P is the rotational period of the pulsar and P is its apparent spin-down rate listed in Table 2, µT is the total transverse proper motion of the pulsar, d its distance, c is the speed of light, and a represents the difference of the accelerations of the pulsar and the Solar System in the gravitational field of the Galaxy, projected along the line of sight from the Solar System to the pulsar.

The term (µT²d)/c represents the contribution from the Shlovskii effect (Shklovskii 1970) and it is relatively small for PSR J1528-3146 due to its relatively small distance and transverse proper motion. The last term, a/c, can be calculated using the GalPot model of gravitational potential of the Galaxy (McMillan 2017; Dehnen & Binney 1998); it is also small given
<figure><img src="image_6.png" /><figcaption>Fig. 5. Timing residuals as a function of time (top) and orbital phase (bottom), with no time offsets applied to the timing residuals for the different TOA datasets. No trends are visible in the timing residuals plotted as a function of orbital phase or as a function of time, indicating that the timing model describes the data appropriately.</figcaption></figure>

the small distance of the pulsar. The latter two terms combined represent a total kinematic contribution to the spin derivative $\Delta\dot{P} = 1.2(3) \times 10^{-21} \text{ s s}^{-1}$ that is small compared to the apparent spin-down rate $\dot{P} = 2.4850(3) \times 10^{-19} \text{ s s}^{-1}$, so that our estimate for the pulsar's intrinsic spin period derivative $\dot{P}_{\text{int}} = \dot{P} - \Delta\dot{P}$ is very close to $\dot{P}$ (see Table 2).

The derived intrinsic spin-down rate was then used to determine the pulsar's spin-down power, $\dot{E}$, magnetic field at the stellar surface, $B_{\text{surf}}$, and characteristic age, $\tau_{c}$. Although PSR J1528-3146 is located at a small distance to the Earth, its spin-down power, $\dot{E} \sim 4.3 \times 10^{31} \text{ erg s}^{-1}$, is very low compared to that of known $\gamma$-ray pulsars (see e.g., Abdo et al. 2013), and no detection of $\gamma$-ray pulsations has been reported to date for this pulsar. The nearest object in the Fermi Large Area Telescope (LAT) 12-Year Point Source Catalog (4FGL-DR3, Abdollahi et al. 2022) is located 45' away from the pulsar, meaning that PSR J1528-3146 does not, in fact, have any counterpart in the Fermi LAT catalog of $\gamma$-ray sources.

## 4.2. Keplerian and post-keplerian parameters

### 4.2.1. Rate of advance of periastron

From our timing analysis, we measured $\dot{\omega}_{\text{obs}} = 0.057(3)^{\circ} \text{yr}^{-1}$. In the absence of a third object, the observed advance of periastron is given by Lorimer & Kramer (2012):

$$ \dot{\omega}_{\text{obs}} = \dot{\omega}_{\text{GR}} + \dot{\omega}_{\text{kin}} + \dot{\omega}_{\text{SO}}, \quad (4) $$

where the first term is caused by the leading order relativistic effect, which, assuming GR, only depends on the total mass of the system, $M_{\text{tot}}$ (Robertson 1938; Taylor & Weisberg 1982):

$$ M_{\text{tot}} = \frac{1}{T_{\odot}} \left( \frac{\dot{\omega}_{\text{GR}}}{3} (1 - e^2) \right)^{3/2} \left( \frac{P_{\text{b}}}{2\pi} \right)^{5/2} M_{\odot}, \quad (5) $$

where $T_{\odot} = GM_{\odot}^N/c^3 = 4.925490947... \text{ ms}$ is an exact quantity, the solar mass parameter ($GM_{\odot}^N$, Prša et al. 2016) in time units. Assuming $\dot{\omega}_{\text{GR}} \sim \dot{\omega}_{\text{obs}}$, we find a total mass of $M_{\text{tot}} = 2.8(2) M_{\odot}$.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Best-fit piecewise-constant DM values from the timing analysis (see Sect. 4). We used segments of 100 days and excluded time intervals that contained no timing data.</figcaption></figure>

This total mass constraint is displayed graphically in Fig. 7 (red lines) and forms the basis of the comparison with the estimated mass from the Bayesian analysis in Sect. 4.3.

We now go on to evaluate the above assumption that $\dot{\omega}_{GR} \sim \dot{\omega}_{obs}$. The second term in Eq. (4) is caused by kinematic effects and is given by Kopeikin (1996):

$$ \dot{\omega}_{\text{kin}} = \frac{\mu_T}{\sin i} \cos (\Theta_{\mu} - \Omega), \qquad (6) $$

which is at most $1.2 \times 10^{-6} \text{ yr}^{-1}$, given the constraints on $i$ and $\Omega$ from the Bayesian analysis presented in Sect. 4.3. This is four orders of magnitude smaller than the observed advance of periastron. The last term, $\dot{\omega}_{SO}$, comes from the spin-orbit coupling. This term is currently only measurable in the case of the double pulsar, PSR J0737−3039A/B (Kramer et al. 2021b). Given PSR J1528−3146’s much wider orbit compared to those of PSRs J0737−3039A and B, this term about two orders of magnitude smaller. Therefore, the initial $\dot{\omega}_{GR} \sim \dot{\omega}_{obs}$ assumption is warranted.

### 4.2.2. Shapiro delay

The DDH solution uses the orthometric parameterization of Freire & Wex (2010) to describe the Shapiro delay. The orthometric amplitude ($h_3$) and the orthometric ratio ($\zeta$) parameters can be written as a function of the Shapiro range ($r$) and shape ($s$) of the DD model, as follows:

$$ h_3 = r\zeta^3, \qquad \zeta = \frac{s}{1+\bar{c}} = \left(\frac{1-\bar{c}}{1+\bar{c}}\right)^{1/2}, \qquad (7) $$

where $\bar{c} \equiv \sqrt{1-s^2} = |\cos i|$. From our timing, we obtain significant measurements of both parameters: $h_3 = 1.2(1)\mu$s and $\zeta = 0.44(8)$. These values can thus be used to constrain the individual masses of the components of the system. In Fig. 7, these constraints are displayed as blue solid (respectively, dashed) lines for $h_3$ (resp., $\zeta$).

### 4.2.3. Variation of the orbital period

Our timing analysis of PSR J1528−3146 using the DDH model did not yield a significant measurement of the variation of the system’s orbital period, $\dot{P}_b$. The timing model presented in Table 2 therefore does not include a measurement of $\dot{P}_b$. Fitting for this parameter leads to $\dot{P}_{b,obs} = -0.7(2.2) \times 10^{-14} \text{ s s}^{-1}$, which is consistent with zero given the present uncertainties.

The observable orbital period change can be written as follows (Lorimer & Kramer 2012):

$$ \left(\frac{\dot{P}_b}{P_b}\right) = \left(\frac{\dot{P}_b}{P_b}\right)^{\text{GW}} + \left(\frac{\dot{P}_b}{P_b}\right)^{\text{m}} + \left(\frac{\dot{P}_b}{P_b}\right)^{\text{T}} + \frac{\mu_T^2 d + a}{c}, \qquad (8) $$

where the first term corresponds to the intrinsic orbital period variation caused by gravitational wave emission. Assuming the masses for the pulsar and its companion derived from the Bayesian analysis presented in Sect. 4.3, we estimate this term to be about $-4.8 \times 10^{-15} \text{ s s}^{-1}$ (Peters 1964), namely, an order of magnitude smaller than the uncertainties.

The second term is the orbital period change caused by mass loss of the system:

$$ \left(\frac{\dot{P}_b}{P_b}\right)^{\text{m}} = \frac{8\pi^2 G I}{T_{\odot}c^5 M_{\text{tot}}} \frac{\dot{P}}{P^3}, \qquad (9) $$

where $G$ denotes Newton’s gravitational constant, and $I \sim 10^{45} \text{ g cm}^2$ is the moment of inertia of the pulsar. Assuming $M_{\text{tot}} = 2.97 M_{\odot}$ (see Sect. 4.3), we find that the expected contribution from mass loss to the orbital period change is four orders of magnitude smaller than $\dot{P}_{b,obs}$, and is thus negligible. Similarly, the third term in Eq. (8), which is due to tidal torques, can be neglected since the companion star is not deformed and is not interacting with the pulsar.

Finally, the last term corresponds to the kinematic contribution, that includes the Shklosvskii effect and the Galactic differential acceleration. We estimate it to be $\dot{P}_{b,kin} = 5.4(1.4) \times 10^{-15} \text{ s s}^{-1}$. This is of the same order as $\dot{P}_{b,GW}$, so that the two terms dominate in the calculation of $\dot{P}_{b,obs}$. However, they mostly cancel each other: we find a prediction for $\dot{P}_b$ of $\sim 0.7(1.4) \times 10^{-15} \text{ s s}^{-1}$. Our measurement of the apparent orbital period variation is compatible with the predicted $\dot{P}_b$, but its uncertainty is still several times larger than that prediction.

### 4.2.4. Variation of the projected semi-major axis

The projected semi-major axis $x$ and its secular change $\dot{x}$ have been measured using the DDH model (see Table 2) to be $\dot{x} = 2.6(1.9) \times 10^{-15} \text{ lt-s s}^{-1}$. Assuming the absence of a third body in the system, the change in $x$ can be written as the sum of various contributions (Lorimer & Kramer 2012):

$$ \left(\frac{\dot{x}}{x}\right)^{\text{obs}} = \left(\frac{\dot{x}}{x}\right)^{\text{PM}} + \left(\frac{\dot{x}}{x}\right)^{\text{GW}} + \frac{d\epsilon_{\text{A}}}{dt} + \frac{d\mu_{\text{T}}^2 + a}{c} + \left(\frac{\dot{x}}{x}\right)^{\text{m}} + \left(\frac{\dot{x}}{x}\right)^{\text{SO}}. \qquad (10) $$

A detailed evaluation of all these terms was done by Guo et al. (2021) for PSR J2222−0137, whose orbital properties are similar to those of PSR J1528−3146. Guo et al. (2021) concluded that all terms except for the first one are negligible. Similar calculations for PSR J1528−3146 enable us to conclude that the same applies to this system. The leading term in Eq. (10) describes the secular change in $\dot{x}$ caused by the proper motion of the system, as follows (modified from Kopeikin 1996):

$$ \left(\frac{\dot{x}}{x}\right)^{\text{PM}} = \mu_T \cot i \sin (\Theta_{\mu} - \Omega). \qquad (11) $$

We find that $\dot{x}^{\text{PM}}$ is at most $3.5 \times 10^{-15} \text{ lt-s s}^{-1}$, comparable to $\dot{x}_{\text{obs}}$. This means that the latter measurement already constrains the orbital orientation of the system, in particular, the value of $\sin (\Theta_{\mu} - \Omega)$. These constraints are depicted graphically in Fig. 8.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Constraints on the mass of the companion of PSR J1528−3146 as a function of the cosine of the orbital inclination (left), and as a function of the mass of the pulsar (right). In the left-hand plot, the gray region is excluded by the requirement that the pulsar mass must be greater than zero, and the narrow orange region near cos i = 0 is excluded by the marginally significant constraints on x̄. On the right-hand side of the figure, mass values in the gray region are excluded by the mass function measurement. The solid red lines indicate the contours of the 1σ uncertainty regions for ω, the 1σ uncertainty regions for h₃ (resp., ζ) are displayed as solid blue (resp., dashed blue) lines; these parameters are measured in the DDH orbital model. The solid contours include 1σ and 2σ equivalent percentiles for the joint posterior pdfs for cos i and M_c (left) and M_c−M_c (right). The marginalized posterior pdfs for cos i, M_c and M_p are displayed for each plot axis. In the right-hand panel, the probability density increases from right to left.</figcaption></figure>

4.3. Bayesian analysis

In order to estimate the physical parameters of the system (pulsar mass, M<sub>p</sub>, companion mass, M<sub>c</sub>, orbital inclination angle, i, and longitude of the ascending node, Ω) while self-consistently taking all available constraints into account, we performed a Bayesian analysis of the timing data. Following the method described in Stovall et al. (2019), we constructed a χ<sup>2</sup> map of the three-dimensional (3D) parameter space of cos i−Ω−M<sub>tot</sub>. For randomly aligned orbital angular momenta and systems with unknown total mass, these parameters have a priori uniform probability densities. These uniform distributions are our Bayesian priors.

The orbital model we used to take the kinematic effects into account was the T2 model implemented in TEMPO2 (Edwards et al. 2006), which not only takes into account the secular effects on x̄ and ω (Kopeikin 1996), but also the effects caused by the varying line of sight to the system caused by the orbital motion of the Earth, such as the annual orbital parallax (Kopeikin 1995). The latter terms are generally very small and are not included in any other orbital models. In some cases the annual orbital parallax allows for an elimination of the degeneracy between i and 180°−i (see e.g., Stovall et al. 2019). The T2 model uses Ω and i as parameters, calculating all kinematic effects internally from these two numbers, which are specified in the “observer’s convention”, where Ω is

measured counterclockwise from north through east and i is the
angle between the orbital angular momentum and the line of
sight from the pulsar to the Earth.

For each value of the cos i−Ω−M<sub>tot</sub> grid, we fixed Ω and i and then derived the PK parameters M<sub>c</sub>, γ and ω from the total mass, the inclination and the mass function, using the relevant GR equations (from e.g., Lorimer & Kramer 2012). We then fit a solution with these fixed parameters to the timing data, while allowed all the other timing parameters to vary. When TEMPO2 is done minimizing the weighted sum of the squares of all timing residuals, we record the χ<sup>2</sup> value obtained for that point of the cos i−Ω−M<sub>tot</sub> grid.

We performed the analysis in two distinct cos i regions instead of a continuous range to reduce the computational cost: from −0.8 to −0.4 for the first region and from 0.4 to 0.8 for the second. We explored total mass M<sub>tot</sub> values between 2 and 4 M<sub>☉</sub>. The ranges and step sizes for the three parameters are given in Table 3. The resulting three-dimensional χ<sup>2</sup> map was then used to derive the 3D probability density map using the Bayesian likelihood function (Splaver et al. 2002):

$$
p(\Omega, \cos i, M) \propto \exp\left[\frac{(\chi_{\min}^2 - \chi^2)/2}{2}\right], \quad (12)
$$

where χ²_min is the minimum χ² value over the entire map. From
this three-dimensional joint posterior probability density func-
tion (PDF), we derived 2D joint posterior PDFs for the M_c – cos i
and M_p – M_c planes (see main panels in Fig. 7), for the cos i – Ω
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Constraints on the longitude of the ascending node, $\Omega$, as a function of the cosine of the orbital inclination. The gray shaded regions are excluded by the mass function and the estimate for the total mass of the binary system. The horizontal dashed orange line indicates the position angle of the proper motion, $\Theta_{\mu}$. The constraints given by $\dot{x}_{\text{obs}}$ and its $1-\sigma$ uncertainty are displayed as solid orange lines, the constraints given by the measurement of $\zeta$ are given by the blue dashed lines. These parameters are as measured with the DDH model. The grayscale maps represent the joint posterior pdf for $\cos i$ and $\Omega$, with darker shades corresponding to higher probabilities. The marginal plots display the marginalized posterior pdfs for $\cos i$ and $\Omega$. In the right-hand panel, the probability density increases from right to left.</figcaption></figure>

<table><caption>**Table 3.** Details of the parameters used in the $\chi^2$ sampling.</caption><thead><tr><th>Region</th><th>cos i</th><th>Ω</th><th>Best cos (i)</th><th>Best Ω</th><th>Best M<sub>tot</sub></th><th>χ<sup>2</sup><sub>min</sub></th></tr></thead><tbody><tr><td>1</td><td>-0.8 to -0.4</td><td>20° to 110°</td><td>-0.55</td><td>70°</td><td>2.98 M<sub>☉</sub></td><td>13318.23</td></tr><tr><td>2</td><td>-0.8 to -0.4</td><td>110° to 200°</td><td>-0.55</td><td>150°</td><td>2.98 M<sub>☉</sub></td><td>13318.24</td></tr><tr><td>3</td><td>0.4 to 0.8</td><td>200° to 290°</td><td>0.55</td><td>250°</td><td>2.98 M<sub>☉</sub></td><td>13318.23</td></tr><tr><td>4</td><td>0.4 to 0.8</td><td>290° to 360°</td><td>0.55</td><td>330°</td><td>2.98 M<sub>☉</sub></td><td>13318.24</td></tr></tbody></table>

**Notes.** The $\cos i$ and $\Omega$ parameters were sampled over the ranges given in the second and third columns, with respective steps of 0.01 and 10°. For the $M_{\text{tot}}$ parameter, we explored values ranging from 2 to 4 $M_{\odot}$, with a step of 0.01 $M_{\odot}$. See Sect. 4.3 for additional details regarding this analysis.

plane (see the main panel of Fig. 8), and marginalized posterior pdfs for $\cos i$, $\Omega$, $M_{\text{tot}}$, $M_{\text{p}}$ and $M_{\text{c}}$; most of these appear in the lateral panels of Figs. 7 and 8.

The marginalized posterior pdf of $\cos i$ has two maxima, at -0.55 and 0.55. As can be seen from Table 3, these solutions have virtually the same $\chi^2$ and therefore the same associated probabilities. The preferred value for $\cos i$ of -0.55 corresponds to an inclination of $i = 123.3^{+4.8}_{-4.1}^{\circ}$, the other solution corresponds to $i = 56.7^{+4.1}_{-4.8}^{\circ}$. The fact that the two maxima have such similar probabilities indicates that the annual orbital parallelax is not detected in the considered dataset. As noted in Sect. 3, the polarimetry suggests that the latter may be the preferred solution.

Figure 8 shows four distinct solutions in the $\cos i-\Omega$ space, which are listed in Table 3. Again, none of the solutions is

strongly preferred over the other, as can be seen from Fig. 8. Furthermore, there is a plateau of high probability between these peaks. The high-probability regions are well described by the $\dot{x}$ and $\zeta$ from the timing analysis, the broad distribution of probability can thus be understood to result from the relatively low precision of $\dot{x}_{\text{obs}}$.

For the masses we obtained, from the marginalized posterior pdfs, $M_{\text{tot}} = 2.94^{+0.21}_{-0.20} M_{\odot}$, $M_{\text{p}} = 1.61^{+0.14}_{-0.13} M_{\odot}$, and $M_{\text{c}} = 1.33^{+0.08}_{-0.07} M_{\odot}$; here, the values are the medians of the pdfs and the uncertainties include 68.3% of the total probability around these medians. If we include 95.44% of the probability around the medians, then we obtain $M_{\text{tot}} = 2.95^{+0.44}_{-0.41} M_{\odot}$, $M_{\text{p}} = 1.61^{+0.28}_{-0.26} M_{\odot}$, and $M_{\text{c}} = 1.33^{+0.17}_{-0.15} M_{\odot}$. The value of $M_{\text{tot}}$ matches that estimated based on the $\dot{\omega}$ measurement, of 2.8 (2) $M_{\odot}$.
**5. Discussion and prospects**

In this work, we describe the polarimetric and timing analysis of our radio observations of PSR J1528−3146. This effort has resulted in substantially improved physical parameters for this system as well as vastly improved knowledge of the radio emission properties for the pulsar.

The mass measurements for the PSR J1528−3146 system are remarkable for several reasons: these values depend mostly on the fact that even with an orbital eccentricity as low as $2.1 \times 10^{-4}$, we can measure the slow periastron advance of $0.057^\circ \text{yr}^{-1}$; furthermore, with a low orbital inclination of about $57^\circ$, we were able to measure the Shapiro delay with enough precision to be useful, at least in conjunction with $\dot{\omega}$. Another feature of the system is that the TOA rms of $9.2 \mu\text{s}$ represents only 0.00015 of the spin period. The system is thus reminiscent of PSR J2222−0137 (Boyleς et al. 2013), but with less pronounced characteristics: the latter pulsar spins faster ($P = 32.8$ ms), its orbital period is slightly smaller (2.44 days), and its orbit is slightly more eccentric ($3.8 \times 10^{-4}$). On the other hand, PSR J1528−3146 has a magnetic field that is five times larger. The similarities between the two systems do suggest a close evolutionary relation.

For PSR J2222−0137, the high orbital inclination ($i = 85.27(4)^\circ$), the better timing precision and the slightly more eccentric, faster orbit allow more precise measurements of the Shapiro delay and $\dot{\omega}$. This results in precise mass measurements: $M_p = 1.831(10) M_\odot$ and $M_c = 1.319(4) M_\odot$; the high total mass $M_{\text{tot}} = 3.150(14) M_\odot$ implies that it is the most massive double degenerate binary known in the Galaxy (Guo et al. 2021). Furthermore, as already highlighted by Cognard et al. (2017), the amount of material accreted by PSR J2222−0137 is very small, which implies that the mass of that pulsar is basically its birth mass. This shows that NSs are born with a wide range of masses.

For the same reasons, the mass of PSR J1528−3146 is also very likely its birth mass – its slower spin and higher magnetic field indicate it has been even less recycled than PSR J2222−0137, and therefore it is likely to have accreted less mass. Thus, the precise measurement of the mass of this pulsar will be important – at the very least, it will add to the distribution of known NS birth masses.

The mass of the companion WD is also interesting. First, it could be the most massive WD around any pulsar measured to date. Second, because, as found by McKee et al. (2020), the masses of WD stars orbiting pulsars do seem to clump around a few preferred values. For He WD companions around MSPs, the measured masses increase with orbital period (Tauris & Savonije 1999), reaching a maximum of about $0.4 M_\odot$. Then, for several systems with slower spin periods (tens of ms), there is a second clump in companion masses, between 0.7 and $0.9 M_\odot$, and possibly narrower. Finally, PSR J2222−0137 makes a third clump with companion masses around $1.3 M_\odot$. In this respect, the current mass measurement of the WD companion to PSR J1528−3146 is already important because it confirms the existence of this third clump, and highlights the general pattern. There are very few WD masses outside these ranges: in two cases, they orbit pulsars that are not recycled, and in very few other cases (e.g., PSR J1614−2230, Crawford et al. 2006; Arzoumanian et al. 2018; Shamohammadi et al. 2023), the system is thought to have evolved through case A Roche lobe overflow (Tauris et al. 2011), which is a likely explanation for their very short spin periods.

In the near future, the precision of $\dot{\omega}$ and $\dot{x}$ values for this system will keep improving relatively fast, with the uncertainties decreasing faster than the nominal $T^{-3/2}$ rate (where T is the tim-

ing baseline; see Damour & Taylor 1992), because of the recent improvement in timing precision achieved by MeerKAT. A better $\dot{\omega}$ will improve the mass measurements, and along with better $\dot{x}$ will greatly improve the constraints on the orbital orientation of the system. The measurement of $\dot{P}_b$ will improve faster, with uncertainties decreasing faster than the nominal rate of $T^{-5/2}$ – again, because of the timing precision of MeerKAT.

The improved $\dot{P}_b$ that can be achieved in the near future has the potential for a sensitive search for dipolar gravitational wave (DGW) emission, similar to that done with PSR J2222−0137 (Guo et al. 2021). The non-detection of DGW emission in the latter system was important for ruling out the phenomenon of spontaneous scalarization (Zhao et al. 2022), which was predicted by a class of scalar-tensor gravity theories of Damour & Esposito-Farèse (1993). This happened because of the specific mass of PSR J2222−0137, which is in the so-called “mass gap” of spontaneous scalarization (Shibata et al. 2014), which extends from $1.6$ to $1.9 M_\odot$ (Shao et al. 2017). Whether a similar test for PSR J1528−3146 is useful or not will depend on the mass of the pulsar, which is not yet known with sufficient precision: the lower half of the uncertainty region is below the mass gap, any value between the median and the $2-\sigma$ upper limit ($1.89 M_\odot$) would place the pulsar in that gap. Thus, continued timing of this pulsar will be essential in order to improve the mass measurements and determine its potential as a gravitational laboratory.

**Acknowledgements.** 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 of the Paris Observatory. MeerTime data is housed on the OzSTAR supercomputer at Swinburne University of Technology. The Parkes radio telescope is funded by the Commonwealth of Australia 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 NASAs 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 and the Laureate fellowship number FL150100148. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from the “Programme National de Cosmologie et Galaxies” (PNCG) and “Programme National Hautes Energies” (PNHE) of CNRS/INSU, France. The authors acknowledge the constant support from the Max-Planck-Society. M.B. acknowledges support through ARC grant CE170100004 (OzGrav).

**References**

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53

Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Nature, 300, 728

Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37

Bailes, M., Jameson, A., Abbate, F., et al. 2020, PASA, 37, e028

Berezina, M., Champion, D. J., Freire, P. C. C., et al. 2017, MNRAS, 470, 4421

Bisnovatyi-Kogan, G. S., & Komberg, B. V. 1974, Soviet Astron., 18, 217

Boyleς, J., Lynch, R. S., Ransom, S. M., et al. 2013, ApJ, 763, 80

Buchner, J. 2021, J. Open Source Softw., 6, 3001

Camilo, F. 1996, in IAU Colloq. 160: Pulsars: Problems and Progress, eds. S. Johnston, M. A. Walker, & M. Bailes, ASP Conf. Ser., 105, 539

Cognard, I., Freire, P. C. C., Guillemot, L., et al. 2017, ApJ, 844, 128

Cordes, J. M., & Lazio, T. J. W. 2002, ArXiv e-prints [arXiv:astro-ph/0207156]

Crawford, F., Roberts, M. S. E., Hessels, J. W. T., et al. 2006, ApJ, 652, 1499

Damour, T., & Deruelle, N. 1985, Ann. Inst. Henri Poincaré Phys. Théor, 43, 107

Damour, T., & Deruelle, N. 1986, Ann. Inst. Henri Poincaré Phys. Théor, 44, 263

Damour, T., & Esposito-Farèse, G. 1993, Phys. Rev. Lett., 70, 2220
Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840

Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429

Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549

Freire, P. C. C., & Wex, N. 2010, MNRAS, 409, 199

Guo, Y. J., Freire, P. C. C., Guillemot, L., et al. 2021, A&A, 654, A16

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

Hobbs, G., Manchester, R. N., Dunning, A., et al. 2020, PASA, 37, e012

Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302

Jacoby, B. A., Chakrabarty, D., van Kerkwijk, M. H., Kulkarni, S. R., & Kaplan, D. L. 2006, ApJ, 640, L183

Jacoby, B. A., Bailes, M., Ord, S. M., Knight, H. S., & Hotan, A. W. 2007, ApJ, 656, 408

Johnston, S., & Kramer, M. 2019, MNRAS, 490, 4565

Jonas, J. L. 2009, IEEE Proc., 97, 1522

Kopeikin, S. M. 1995, ApJ, 439, L5

Kopeikin, S. M. 1996, ApJ, 467, L93

Kramer, M., Stairs, I. H., Venkatraman Krishnan, V., et al. 2021a, MNRAS, 504, 2094

Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2021b, Phys. Rev. X, 11, 041050

Lorimer, D. R., & Kramer, M. 2012, Handbook of Pulsar Astronomy
(Cambridge: Cambridge University Press)

McKee, J. W., Freire, P. C. C., Berezina, M., et al. 2020, MNRAS, 499, 4082

McMillan, P. J. 2017, MNRAS, 465, 76

Pennucci, T. T., Demorest, P. B., & Ransom, S. M. 2014, ApJ, 790, 93

Peters, P. C. 1964, Phys. Rev., 136, 1224

Prša, A., Harmanec, P., Torres, G., et al. 2016, AJ, 152, 41

Radhakrishnan, V., & Cooke, D. J. 1969, Astrophys. Lett., 3, 225

Robertson, H. P. 1938, Ann. Math., 39, 101

Serylak, M., Johnston, S., Kramer, M., et al. 2021, MNRAS, 505, 4483

Shamohammadi, M., Bailes, M., Freire, P. C. C., et al. 2023, MNRAS, 520, 1789

Shao, L., Sennett, N., Buonanno, A., Kramer, M., & Wex, N. 2017, Phys. Rev. X, 7, 041025

Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789

Shibata, M., Taniguchi, K., Okawa, H., & Buonanno, A. 2014, Phys. Rev. D, 89, 084005

Shklovskii, I. S. 1970, Soviet Astron., 13, 562

Splaver, E. M., Nice, D. J., Arzoumanian, Z., et al. 2002, ApJ, 581, 509

Staveley-Smith, L., Wilson, W. E., Bird, T. S., et al. 1996, PASA, 13, 243

Stovall, K., Freire, P. C. C., Antoniadis, J., et al. 2019, ApJ, 870, 74

Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928

Tauris, T. M., Langer, N., & Kramer, M. 2011, MNRAS, 416, 2130

Tauris, T. M., Langer, N., & Kramer, M. 2012, MNRAS, 425, 1601

Taylor, J. H., & Weisberg, J. M. 1982, ApJ, 253, 908

van Straten, W. 2006, ApJ, 642, 1004

Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29

Zhao, J., Freire, P. C. C., Kramer, M., Shao, L., & Wex, N. 2022, Class. Quant. Grav., 39, 11LT01