© The Authors 2023

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 ≲ 30 ms) and very small period increases (Ṗ ≲ 10⁻¹⁷ s s⁻¹). 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 ∼ 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 ∼ 3.18 days, a projected semi-major axis x ≡ a sin i/c ∼ 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 ∼ 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 pulsar mass of 1.4 M_⊙, the mass function of 0.159 M_⊙ indicates that the minimum companion mass (calculated for an edge-on orbit, i.e., with an inclination of 90°) is ∼0.94 M_⊙. 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 × 0.5 × 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.

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.

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

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 ϕ_M, 1 and ϕ_M, 2 as the phases of the two components of the main pulse, and ϕ_IP that of the weak component, we found ϕ_M, 1 = 0.469 (2), ϕ_M, 2 = 0.54410 (2), and ϕ_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, α, must be close to 90°.

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

Tighter constraints on the value of α and that of ζ, 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) Ψ 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 Ψ are related to the Q and U Stokes parameters of the pulsar’s emission as $\mathrm{\Psi }=\frac{1}{2}arctan\left(\frac{U}{Q}\right)$. In the following, we use the modified form introduced in Johnston & Kramer (2019):

$$\begin{array}{c}\hfill \mathrm{\Psi }={\mathrm{\Psi }}_0+arctan\frac{sin\alpha sin(\varphi -{\varphi }_0-\mathrm{\Delta })}{sin\zeta cos\alpha -cos\zeta sin\alpha cos(\varphi -{\varphi }_0-\mathrm{\Delta })}·\end{array}$$(1)

In the above expression, ϕ₀ denotes the pulse phase for which Ψ = Ψ₀ (assuming Δ = 0), and the term Δ 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 dedispersed 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² 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 λ as the radio wavelength, we have:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{PPA}}=\mathrm{RM}\times {\lambda }^2.\end{array}$$(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 Ψ 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 UltraNest sampler by Buchner (2021), using different priors for the angles in different runs. In each case, we used a standard expression, ln(ℒ) = − χ²/2, as a likelihood function. In the first fit, we assumed our random α and ζ angles to form spin and magnetic field vectors that are distributed uniformly over a unit-sphere, while Δ = 0. The prior for Ψ₀ was a uniform distribution over [ − 90° ,90° ] and, similarly, ϕ₀ 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 (α, ζ) = (∼103°, ∼92°), a second solution (less favoured by the log-likelihood statistic) is seen for (α, ζ) = (∼80°, ∼72°). Either case suggests an orthogonal rotator, 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).

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.

As we discuss further in Sect. 4.3, the Bayesian analysis of the timing data suggests an orbital inclination angle of either i ∼ 57° or i ∼ 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 ≈ ζ, and also Ψ₀ ≈ Ω, where Ω 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 ζ-values are not consistent with the results from timing. In a second and third fit, we therefore applied uniform priors for ζ restricted to a 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³. 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₃) and the orthometric ratio (ς), 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 χ² values become equal to unity. The EFAC factors we determined through this 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-consant DM values obtained from the timing analysis. The DM displays long-term variations around the value of 18.163 (6) measured by Jacoby et al. (2007).

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.

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.

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.

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, Ṗ_int, as:

$$\begin{array}{c}\hfill \frac{{\dot{P}}_{\mathrm{int}}}{P}=\frac{\dot{P}}{P}-\frac{{\mu }_{\mathrm{T}}^{2}d}{c}-\frac{a}{c},\end{array}$$(3)

where P is the rotational period of the pulsar and Ṗ 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 $({\mu }_{\text{T}}^{2}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 the small distance of the pulsar. The latter two terms combined represent a total kinematic contribution to the spin derivative ΔṖ = 1.2(3) × 10⁻²¹ s s⁻¹ that is small compared to the apparent spin-down rate Ṗ = 2.4850 (3) × 10⁻¹⁹ s s⁻¹, so that our estimate for the pulsar’s intrinsic spin period derivative Ṗ_int = Ṗ − ΔṖ is very close to Ṗ (see Table 2).

The derived intrinsic spin-down rate was then used to determine the pulsar’s spin-down power, Ė, magnetic field at the stellar surface, B_surf, and characteristic age, τ_c. Although PSR J1528−3146 is located at a small distance to the Earth, its spin-down power, Ė ∼ 4.3 × 10³¹ erg s⁻¹, is very low compared to that of known γ-ray pulsars (see e.g., Abdo et al. 2013), and no detection of γ-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 γ-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 }}_{\mathrm{obs}}=0.057(3)°$ yr⁻¹. In the absence of a third object, the observed advance of periastron is given by Lorimer & Kramer (2012):

$$\begin{array}{c}\hfill {\dot{\omega }}_{\mathrm{obs}}={\dot{\omega }}_{\mathrm{GR}}+{\dot{\omega }}_{\mathrm{kin}}+{\dot{\omega }}_{\mathrm{SO}},\end{array}$$(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_tot (Robertson 1938; Taylor & Weisberg 1982):

$$\begin{array}{c}\hfill M_{\mathrm{tot}}=\frac{1}{T_{\odot }}{\left(\frac{{\dot{\omega }}_{\mathrm{GR}}}{3}(1-e^2)\right)}^{3/2}{\left(\frac{P_{\mathrm{b}}}{2\pi }\right)}^{5/2}{M}_{\odot },\end{array}$$(5)

where $T_{\odot }=\mathcal{G}{\mathcal{M}}_{\odot }^{\mathrm{N}}/c^3=4.925490947\dots \mathrm{\mu }$s is an exact quantity, the solar mass parameter ($\mathcal{G}{\mathcal{M}}_{\odot }^{\mathrm{N}}$, Prša et al. 2016) in time units. Assuming ${\dot{\omega }}_{\mathrm{GR}}\sim {\dot{\omega }}_{\mathrm{obs}}$, we find a total mass of M_tot = 2.8 (2) M_⊙. 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.

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 ẋ. 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 $\dot{\omega }$, the 1σ uncertainty regions for h3 (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 Mc (left) and Mc − Mc (right). The marginalized posterior pdfs for cos i, Mc and Mp are displayed for each plot axis. In the right-hand panel, the probability density increases from right to left.

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

$$\begin{array}{c}\hfill {\dot{\omega }}_{\mathrm{kin}}=\frac{{\mu }_{\mathrm{T}}}{sini}cos({\mathrm{\Theta }}_{\mu }-\mathrm{\Omega }),\end{array}$$(6)

which is at most 1.2 × 10⁻⁶ ° yr⁻¹, given the constraints on i and Ω 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 }}_{\mathrm{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 }}_{\mathrm{GR}}\sim {\dot{\omega }}_{\mathrm{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₃) and the orthometric ratio (ς) parameters can be written as a function of the Shapiro range (r) and shape (s) of the DD model, as follows:

$$\begin{array}{c}\hfill h_3=r{\varsigma }^3,\varsigma =\frac{s}{1+\overline{c}}={\left(\frac{1-\overline{c}}{1+\overline{c}}\right)}^{1/2},\end{array}$$(7)

where $\overline{c}\equiv \sqrt{1-s^2}=|cosi|$. From our timing, we obtain significant measurements of both parameters: h₃ = 1.2 (1) μs and ς = 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₃ (resp., ς).

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, Ṗ_b. The timing model presented in Table 2 therefore does not include a measurement of Ṗ_b. Fitting for this parameter leads to Ṗ_b,obs = −0.7 (2.2) × 10⁻¹⁴ s s⁻¹, which is consistent with zero given the present uncertainties.

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

$$\begin{array}{c}\hfill \left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)={\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{GW}}+{\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\dot{m}}+{\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{T}}+\frac{{\mu }_{\mathrm{T}}^{2}d+a}{c},\end{array}$$(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 × 10⁻¹⁵ s s⁻¹ (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:

$$\begin{array}{c}\hfill {\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\dot{m}}=\frac{8{\pi }^{2}G}{T_{\odot }{c}^5}\frac{I}{M_{\mathrm{tot}}}\frac{\dot{P}}{P^3},\end{array}$$(9)

where G denotes Newton’s gravitational constant, and I ∼ 10⁴⁵ g cm² is the moment of inertia of the pulsar. Assuming M_tot = 2.97 M_⊙ (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 Ṗ_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 Shklosvkii effect and the Galactic differential acceleration. We estimate it to be Ṗ_b,kin = 5.4 (1.4) × 10⁻¹⁵ s s⁻¹. This is of the same order as Ṗ_b,GW, so that the two terms dominate in the calculation of Ṗ_b,obs. However, they mostly cancel each other: we find a prediction for Ṗ_b of ∼0.7 (1.4)×10⁻¹⁵ s s⁻¹. Our measurement of the apparent orbital period variation is compatible with the predicted Ṗ_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 ẋ have been measured using the DDH model (see Table 2) to be ẋ = 2.6 (1.9) × 10⁻¹⁵ lt-s s⁻¹. 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):

$$\begin{array}{c}\hfill {\left(\frac{\dot{x}}{x}\right)}^{\mathrm{obs}}={\left(\frac{\dot{x}}{x}\right)}^{\mathrm{PM}}+{\left(\frac{\dot{x}}{x}\right)}^{\mathrm{GW}}+\frac{\mathrm{d}{ϵ}_{\mathrm{A}}}{\mathrm{d}t}+\frac{\mathrm{d}{\mu }_{\mathrm{T}}^2+a}{c}+{\left(\frac{\dot{x}}{x}\right)}^{\dot{m}}+{\left(\frac{\dot{x}}{x}\right)}^{\mathrm{SO}}.\end{array}$$(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 ẋ caused by the proper motion of the system, as follows (modified from Kopeikin 1996):

$$\begin{array}{c}\hfill {\left(\frac{\dot{x}}{x}\right)}^{\mathrm{PM}}={\mu }_{\mathrm{T}}cotisin({\mathrm{\Theta }}_{\mu }-\mathrm{\Omega }).\end{array}$$(11)

We find that ẋ_PM is at most 3.5 × 10⁻¹⁵ lt-s s⁻¹, comparable to ẋ_obs. This means that the latter measurement already constrains the orbital orientation of the system, in particular, the value of sin(Θ_μ − Ω). These constraints are depicted graphically in Fig. 8.

Fig. 8. Constraints on the longitude of the ascending node, Ω, 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, Θμ. The constraints given by ẋobs and its 1-σ uncertainty are displayed as solid orange lines, the constraints given by the measurement of ς 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 Ω, with darker shades corresponding to higher probabilities. The marginal plots display the marginalized posterior pdfs for cos i and Ω. In the right-hand panel, the probability density increases from right to left.

4.3. Bayesian analysis

In order to estimate the physical parameters of the system (pulsar mass, M_p, companion mass, M_c, 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 χ² map of the three-dimensional (3D) parameter space of cos i − Ω − M_tot. 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 ẋ and $\dot{\omega }$ (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 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_tot grid, we fixed Ω and i and then derived the PK parameters M_c, γ and $\dot{\omega }$ 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 χ² value obtained for that point of the cos i − Ω − M_tot 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_tot values between 2 and 4 M_⊙. The ranges and step sizes for the three parameters are given in Table 3. The resulting three-dimensional χ² map was then used to derive the 3D probability density map using the Bayesian likelihood function (Splaver et al. 2002):

$$\begin{array}{c}\hfill p(\mathrm{\Omega },cosi,M)\propto exp[({\chi }_{\min }^2-{\chi }^2)/2],\end{array}$$(12)

where ${\chi }_{\text{min}}^2$ is the minimum χ² value over the entire map. From this three-dimensional joint posterior probability density function (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 − Ω plane (see the main panel of Fig. 8), and marginalized posterior pdfs for cos i, Ω, M_tot, M_p and M_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 χ² and therefore the same associated probabilities. The preferred value for cos i of −0.55 corresponds to an inclination of $i=123.3_{-4.1}^{+4.8}°$, the other solution corresponds to $i=56.7_{-4.8}^{+4.1}°$. The fact that the two maxima have such similar probabilities indicates that the annual orbital parallax 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 − Ω 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 ẋ and ς from the timing analysis, the broad distribution of probability can thus be understood to result from the relatively low precision of ẋ_obs.

For the masses we obtained, from the marginalized posterior pdfs, $M_{\mathrm{tot}}=2.{94}_{-0.20}^{+0.21}{M}_{\odot }$, $M_{\mathrm{p}}=1.{61}_{-0.13}^{+0.14}{M}_{\odot }$, and $M_{\mathrm{c}}=1.{33}_{-0.07}^{+0.08}{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_{\mathrm{tot}}=2.{95}_{-0.41}^{+0.44}{M}_{\odot }$, $M_{\mathrm{p}}=1.{61}_{-0.26}^{+0.28}{M}_{\odot }$, and $M_{\mathrm{c}}=1.{33}_{-0.15}^{+0.17}{M}_{\odot }$. The value of M_tot matches that estimated based on the $\dot{\omega }$ measurement, of 2.8 (2) M_⊙.

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 × 10⁻⁴, we can measure the slow periastron advance of 0.057° yr⁻¹; furthermore, with a low orbital inclination of about 57°, 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 μs represents only 0.00015 of the spin period. The system is thus reminiscent of PSR J2222−0137 (Boyles 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 × 10⁻⁴). 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)°), 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_⊙ and M_c = 1.319 (4) M_⊙; the high total mass M_tot = 3.150 (14) M_⊙ 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_⊙. 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_⊙, and possibly narrower. Finally, PSR J2222−0137 makes a third clump with companion masses around 1.3 M_⊙. 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 ẋ 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 timing 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 ẋ will greatly improve the constraints on the orbital orientation of the system. The measurement of Ṗ_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 Ṗ_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_⊙ (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-σ upper limit (1.89 M_⊙) 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.

Acknowledgments

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

All Tables

All Figures

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.

In the text

	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°.
In the text

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.

In the text

	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.
In the text

	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.
In the text

	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.
In the text

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 ẋ. 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 $\dot{\omega }$, the 1σ uncertainty regions for h3 (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 Mc (left) and Mc − Mc (right). The marginalized posterior pdfs for cos i, Mc and Mp are displayed for each plot axis. In the right-hand panel, the probability density increases from right to left.

In the text

Fig. 8. Constraints on the longitude of the ascending node, Ω, 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, Θμ. The constraints given by ẋobs and its 1-σ uncertainty are displayed as solid orange lines, the constraints given by the measurement of ς 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 Ω, with darker shades corresponding to higher probabilities. The marginal plots display the marginalized posterior pdfs for cos i and Ω. In the right-hand panel, the probability density increases from right to left.

In the text