© The Authors 2024

1. Introduction

As so-called lighthouses in the sky, pulsars are peerless astronomical objects. These highly magnetised neutron stars (NSs) emit a beam of an electromagnetic radiation along their magnetic poles, which is visible as a steady train of pulses at a radio telescope as the beam periodically sweeps across the observer’s line of sight. Due to the high accuracy of atomic reference clocks and low-noise receivers in modern radio telescopes, the times of arrival (ToAs) of the pulses at the telescope’s location are precisely recorded. The motion of the pulsar, the radio emission propagation through the interstellar medium (ISM), and the motion of the radio telescope through the Solar System causes the ToAs to deviate from a purely periodic behaviour. Measuring the ToAs and fitting a model to them that accounts for all these possible effects is known as pulsar timing. In particular, the timing of millisecond pulsars (MSPs), a certain sub-population of pulsars (cf. Sect. 1), allows uniquely precise measurements of the spin, astrometric, and orbital parameters because these pulsars exhibit a uniquely stable rotational behaviour (Lorimer & Kramer 2005).

Timing of pulsars in the southern hemisphere experienced a step change in precision with the arrival of the MeerKAT telescope: the low system temperature (∼18 K) of the L-band receiver, its wide spectral coverage (from 856 to 1712 MHz, and thus a bandwidth of 856 MHz) and the high aperture efficiency of its 64 × 13.5 m offset Gregorian dishes (which improve upon the Parkes Murriyang radio telescope gain by a factor of four) make MeerKAT a powerful addition to other existing radio observatories, significantly increasing the radio sensitivity in the southern hemisphere (Jonas & MeerKAT Team 2016). Furthermore, the ultrawide low band receiver (UWL) of the Murriyang radio telescope has also significantly increased its spectral coverage and sensitivity.

This work was conducted as part of the ‘RelBin’ project (Kramer et al. 2021), which is one of the core sub-projects of the MeerTime project, a five-year Large Survey Project (Bailes et al. 2020), aiming to use the precision of the MeerKAT telescope to explore fundamental physics via pulsar timing. As is outlined in Kramer et al. (2021), the main aim of RelBin is detecting or improving on the measurement of timing parameters related to relativistic effects in the orbital motion of binary systems. Due to the high precision of observations with MeerKAT, this project offers not only a wide range of tests of gravity theories (e.g. Hu et al. 2022), but also improves pulsar population studies by yielding a continuously growing catalogue of precise NS mass measurements and constraining binary evolution theories (e.g. Serylak et al. 2022).

The known pulsar population can be split into two large sub-groups based on their rotational behaviour and spin evolution. The so-called millisecond pulsars exhibit a rotational period of less than 30 ms, as well as a relatively low inferred magnetic field strength (∼10^8 − 9 G). Additionally, about 80% of MSPs are found in binary systems, with main sequence (MS) stars, other NS, or white dwarfs (WDs, of which the helium white dwarfs (HeWDs) are the most numerous) as their companions.

In the current binary evolution models, these systems originate from a stellar binary, in which the more massive star already evolved into an NS. As the companion star leaves the MS and becomes a red giant, it fills its Roche lobe and overflows it. Some of the matter in this so-called Roche-lobe-overflow (RLO) accretes onto the NS. During this period of 𝒪(Gyr), the system is detectable as a low-mass X-ray binary (LMXB). The mass transfer from the red giant to the NS also transfers orbital angular momentum to the NS, leading to a significant spin-up of the NS, such that it becomes an MSP (Radhakrishnan & Srinivasan 1982; Alpar et al. 1982). At the end of this stage, the binary consists of an MSP and a stripped stellar core, which depending on its mass evolves either into an NS or a WD (Tauris & van den Heuvel 2023).

In most observed cases, the companion is a WD; the pulsars in these systems have significantly shorter spin periods, owing to the slower evolution and longer accretion episodes associated with lighter companions. By means of detailed numerical simulations, Tauris & Savonije (1999) derived a relation between the binary orbital period, P_b, and the mass of an HeWD companion, M_WD (which we shall refer to as the TS99 relation). Using catalogues of known MSP-HeWD system masses and comparing them to the latest stage of simulation results, this relation has been reviewed intensely over past decades (see e.g. Smedley et al. 2013, Hui et al. 2018) and usually holds for these binaries.

The tidal interactions accompanying the RLO lead to a circularisation of the binary orbit (Phinney 1992), as well as to an alignment of both the pulsar’s spin axis with the angular momentum axis of the orbit (Phinney & Kulkarni 1994). Since the companion then evolves slowly into a WD, this low-eccentricity orbit and the spin alignment should be retained at later stages. In systems where the companion becomes an NS, the mass loss and the kick associated with the supernova event that forms the second NS will cause a significant increase in the eccentricity of the orbit (e), if not outright disruption, and in many cases a misalignment of the spin of the recycled pulsar with the angular momentum of the post-SN orbit (for a review, see Tauris et al. 2017).

In globular clusters, interactions with external stars passing by can disturb the circular orbits of MSP-HeWDs, which is confirmed by the large number of eccentric binary MSPs in globular clusters¹. Apart from these cases, the majority of MSP-WD systems in the Galactic disc exhibit the expected small residual eccentricities Phinney (1992): there are no nearby stars to perturb them, and the evolution of the companion to a WD does not increase e. Nevertheless over the last decade, six systems with low-mass companions (which in one case are confirmed as HeWDs, Antoniadis et al. 2016), with 0.027 < e < 0.13 and 22 < P_b < 32d have been discovered (see Table 1 and Table 1 in Serylak et al. 2022). These systems clearly do not follow the e-P_b relation predicted by Phinney (1992) and became known as eccentric millisecond pulsar binaries (eMSPs). These systems are puzzling; their formation mechanism has not yet been fully understood (Serylak et al. 2022).

One possibility could be the formation in a triple system that became unstable, ejecting one of the components, as proposed for PSR J1903+0327 (Champion et al. 2008) by Freire et al. (2011) and Portegies Zwart et al. (2011). Intuitively, such a chaotic process should lead to a diversity of orbital configurations and companion types. However, eMSPs do not only have similar orbits, but also similar companion masses (all consistent with being HeWDs). This is seen by Freire & Tauris (2014) and Knispel et al. (2015) as a strong indicator in favour of a deterministic process with a fixed outcome.

For this reason, five competing theories were put forward in order to explain the formation of Galactic eMSPs. They commonly rely on the TS99 relation, but describe various perturbative mechanisms capable of introducing an eccentricity of the binary orbit. A broader introduction to these can be found in Serylak et al. (2022). Lately, the timing analysis of J0955−6150 (Serylak et al. 2022) revealed that this system violates the TS99 relation, which is not compatible with all five theories.

The following analysis of the eMSP PSR J1618−3921 aims to broaden the knowledge about these systems, to find any similarities that could pave the way toward new formation models.

The discovery of PSR J1618−3921 (henceforth J1618−3921, similarly all other J2000 object names refer to pulsars if not indicated otherwise) was reported by Edwards & Bailes (2001) as part of a 1.4-GHz survey of the intermediate Galactic latitudes with the Parkes radio telescope. It is a recycled Galactic-disc pulsar in a binary orbit with a period of 22.7 days and a low-mass companion, presumably accompanied by a low-mass HeWD. With a rotational period of 11.98 ms, but unmeasured period derivative, it was suspected of being an MSP. As a result of the first observations, J1618−3921 stood out from the pulsar population in the Galactic Plane due to its anomalously large orbital eccentricity of 0.027 (Bailes 2007). It is now thought to belong to the eMSP class (Bailes 2007; Serylak et al. 2022); it is however the pulsar with by far the lowest eccentricity and longest spin period within that sub-population.

After a decade of sporadic observations with Parkes, Octau et al. (2018) aimed to precisely measure the pulsar’s spin, astrometric and orbital parameters via a set of dense observations of the pulsar with the Nançay radio telescope (NRT): 51 h of regular observations spread over three observing campaigns. This resulted in the first ever timing solution for this system, its parameters are given in Table 3 of Octau et al. (2018); for completeness also shown in the second column Table 3. This shows that the pulsar is an MSP (from the small period derivative) and confirm the unusual orbital eccentricity. Due to limited precision (this means, a comparably large mean uncertainty in the Nançay ToAs) and timing baseline, the observations were not sufficient to reveal additional timing parameters such as the pulsar’s proper motion, the rate of advance of periastron or the Shapiro delay.

After the addition of J1618−3921 to the RelBin programme, it has been regularly observed with the MeerKAT radio telescope. In addition to that we have also started observing it regularly with the Parkes Radio Telescope and continued observations at NRT. Using all extent data on this pulsar – adding up to a total baseline of more than 23 years – we derived an updated timing solution that improves on both numerical precision and the number of measured relativistic effects of the binary orbit, including the first estimates of the component masses.

In the course of the paper, we start with a brief summary of the observations of J1618−3921 in Section 2. Section 3 covers the profile analysis; Section 4 contains the timing analyses, where we report our new timing solution, including constraints of additional parameters compared to these reported by Octau et al. (2018), which include the constraints on the mass of the system. This is followed in Section 5 by a thorough discussion of the current state of knowledge on eMSPs in Section 3, with special focus on the combined results from the timing of other eMSPs and our J1618−3921 timing parameters. Finally, we conclude by summarising our results in Section 6.

2. Observations and data processing

2.1. Parkes

The first observations of J1618−3921 at the Parkes Radio telescope date back to the 1999 project P309 (Edwards & Bailes 2001), followed by observations in 2001 during P360. In total, the pulsar was observed on six days in August 1999 and on three days in 2001, spanning the orbital phase from 0 to 0.3 and 0.5 to 0.7, respectively. Both runs use the central beam of the 13-beam 21 cm ‘multi-beam’ receiver (Staveley-Smith et al. 1996), with a central frequency of 1374 MHz and a bandwidth of 288 MHz. After a change to the CPSR-2 (Caltech-Parkes-Swinburne-Recorder) backend, J1618−3921 was monitored again in the first half of 2003 with a monthly cadence (covering the orbital phase between 0.2 and 0.7) and twice in 2005 with a gap of five days. These observations were now to made simultaneously using two different 64 MHz bands, with central frequencies at 1341 MHz and 1405 MHz, respectively. Further technical details of these observations are described in Edwards & Bailes (2001), Manchester et al. (2001).

Making use of the ultra-wide band receiver together with the Medusa-backend (Hobbs et al. 2020), observations of J1618−3921 with Parkes resumed in 2019, and continue at the time of writing on a regular basis. The UWL receiver has a bandwidth of 3328 MHz centred around a frequency of 2368 MHz. When used in pulsar folding mode, the data have a typical sub-integration length of 30 s with a resolution of 128 channels for each of the 26 sub-bands; that is, each channel has a bandwidth of 1 MHz, 1024 phase bins and full polarisation information (Hobbs et al. 2020).

2.2. Nançay

As pointed out in Octau et al. (2018), J1618−3921 was first observed at Nançay in May 2009 with the Berkeley-Orléans-Nançay (BON) instrument. Due to a lack of detailed information on spin, orbital parameters and the dispersion measure (DM), these first observations were conducted using the ‘survey’ mode. The incoherent de-dispersion and coarse time resolution associated with this mode lead to very large ToA uncertainties. After the change of the Nançay instrumentation to the NUPPI, a clone of the Green Bank Ultimate Pulsar Processing Instrument (GUPPI) in August 2011, J1618−3921 was still observed in survey mode, with a total bandwidth of 512 MHz divided in 1024 channels with a 64 μs sampling. When a coherent timing solution for J1618−3921 was found, observations were continued using the ‘timing’ mode of NUPPI from December 2014 on. In this mode, NUPPI is able to coherently de-disperse the data and also samples with higher time resolution. This leads to a significant improvement in the quality of the observations, which is visible in the decrease in the mean ToA uncertainty. The observation lengths vary between 1500 and 3400 s, with sub-integrations that vary between 15 and 30 s. In the other axes, all data files have the same resolution of 128 frequency bins, 2048 phase bins and full polarisation information.

2.3. MeerKAT

As part of the RelBin programme (Kramer et al. 2021) at the MeerKAT telescope, J1618−3921 has been observed since March 2019, yielding a total observation time of 28.85 hours. All observations use the L-band receiver (central frequency of 1284 MHz and an effective bandwidth of 776 MHz) together with the PTUSE backend. All technical set-up details can be found in Bailes et al. (2020), Serylak et al. (2021) give a thorough description of the polarisation and flux calibration. The typical sub-integration length is 8 s, and each observation contains usually 2048 sub-integrations at a frequency resolution of 1024 channels over the full bandwidth, 1024 phase bins and the full polarisation information.

Comparing the details of the MeerKAT observations with the Parkes UWL observations, clearly the former have exceptionally low noise, resulting in outstanding quality of profile measurements. This is evident from the mean ToA uncertainty, that is almost a factor of six lower for the MeerKAT observations than for the Parkes (a full discussion of the timing procedure and ToA derivation will be given in Sect. 4). However, the Parkes observations do reveal the structure of the pulse profile at higher frequencies. A summary of all observations is presented in Table 2.

2.4. Data processing

Following standard data reduction procedures in pulsar timing, we used the PSRCHIVE (Hotan et al. 2004) software package. If not explicitly indicated otherwise, all programs or commands referred to in this section are part of this package.

The early Parkes data sets were manually cleaned from radio frequency interference (RFI) using PAZI and PSRZAP. We used the PSRPYPE pipeline² for the data reduction of the UWL observations, that have observing lengths between 2048 and 14 402 seconds. PSRPYPE uses the CLFD software package³ (Morello et al. 2018) RFI cleaning and flux calibration measurements of the Hydra A radio galaxy, returning cleaned and flux calibrated pulsar archives. In order to polarisation calibrate the observations, METM (Measurement Equation Template Matching) (van Straten 2013) was performed on the observations, using off-target calibration observations with injected pulses from a noise diode. The calibrated and cleaned UWL-data was folded into 13 frequency sub-bands.

By default, all pulsar fold-mode observations conducted with MeerKAT as part of the RelBin programme are put through the MEERPIPE pipeline, that performs the RFI excision and polarisation calibration. MEERPIPE is a modified version of COASTGUARD (Lazarus et al. 2016). For the polarisation calibration, a calibration observation is performed before each pulsar observation session, from which the Jones matrices used to calibrate the pulsar observations are obtained. For more details, see Kramer et al. (2021). The cleaned and calibrated files are then decimated in time, frequency and polarisation to the desired resolution; in the case of this work means a scrunching factor of 116 in frequency, 128 in time and a full scrunch in polarisation. This leaves observations containing 8 frequency channels across the 775 MHz.

The NRT data archives went through the full data reduction scheme described in Octau et al. (2018). For the final analysis, we re-installed our latest ephemeris to the data and folded each observation completely in time and polarisation. These archives had a sufficient S/N to keep a resolution of four frequency channels across all observations. We used frequency-resolved templates to account for the strong profile evolution across frequency. These were generated by iteratively running PAAS on the four frequency channels. Then we obtain frequency resolved ToAs via the PAT command.

3. Radio emission properties

3.1. Change in profile with frequency

If not otherwise indicated, for all analyses of the pulse profile, the integrated profile was obtained by summing up all observations of J1618−3921 on a backend-wise basis and summing them along the time, frequency and polarisation axes. The left part of Fig. 1 shows the profile as seen by MeerKAT’s L-band receiver after ∼26 hours of integration, the middle part shows the equivalent for Parkes with the UWL receiver after ∼29 observing hours and the right part corresponds to the ∼50 hours of observations with the Nançay radio telescope. The pulse profile shows a main pulse with a duty cycle of roughly 20%. It consists of two sharper peaks, where the first one exhibits a small sub-peak on its right side. For the MeerKAT observations, the first sub-pulse peaks at ∼6/7 of the peak intensity of the second pulse. The main pulse is preceded by a low-intensity pulse with 1/7 of the main pulse amplitude, that is located at ∼110° beforehand. The shape of that secondary pulse is somewhat different than that of the main pulse, with a plateau-like feature on its left-hand side and a wider peak. Although it has a duty cycle of only around 15%, due to its low amplitude and shape plateau it appears more smeared out than the main pulse.

Fig. 1. Frequency resolved intensity profiles from observations with the MeerKAT L-band receiver (left, Tobs ∼ 28.85 h), the Parkes UWL receiver (middle, Tobs ∼ 28.95 h) and the NRT L-band receiver (right, Tobs ∼ 48.9 h). The top panel of each plot shows the total intensity profile across one period, in case of the MeerKAT and Parkes observations it is flux calibrated. The NRT observations are not flux-calibrated. The bottom panels show the frequency resolved dynamic spectra across adjacent frequency bands. The intensity scale of the MeerKAT and NRT dynamic spectra were adjusted to fit the range of the Parkes spectrum for ease of comparison. The MeerKAT observations were frequency scrunched to 8 channels, those from Parkes down to 13 channels and the NRT data was decimated to 4 channels. The number of channels was chosen such that the frequency resolution is kept as large as possible while providing a S/N in each band that allows for a sufficient ToA precision.

In all plots in Fig. 1, the heat map in the lower sub-figure resolves the pulse into the different frequency bands, a brighter colour indicating a larger intensity. Clearly, the intensity of the pulse decreases with increasing frequency, meaning that the pulsar has a steep spectrum. Spiewak et al. (2022) found a spectral index of −2.28 ± 0.04. At the same time the profile is broader at lower frequencies. For the main pulse this means that the two sharp peaks almost merge into one single broad peak at the lowest frequencies. In light of the template matching used in pulsar timing to create the ToAs, this might be a significant impairment of the ToA precision in the lower frequency bands.

3.2. Polarisation properties

Fig. 2 shows the polarisation profile of J1618−3921 as recorded with the MeerKAT L-band receiver and corrected for the Rotation Measure given in Spiewak et al. (2022), as well as the evolution of the position angle (PA) across the pulsar’s phase. The PAs are measured in the so-called ‘observer’s convention’. The PA displays sudden jumps at the edges of the main pulse that are coincident with the sharp drops in the total linear polarisation. These features are consistent with arising from orthogonal polarisation modes (OPMs; Manchester et al. 1975; Manchester 1975), a phenomena that is either intrinsic to the emission of the pulsar (e.g. Gangadhara 1997), or result from propagation effects in the pulsar magnetosphere (e.g. Blandford & Scharlemann 1976; Melrose & Stoneham 1977).

Fig. 2. Polarisation profile of J1618−3921 obtained from integrating 29 hours of observations with the MeerKAT L-band receiver. The upper part shows the total intensity (light blue), as well as the linear (red) and circular polarisation (dark blue) fraction. The lower part shows the evolution of the position angle (PA) across the pulsar’s phase. The PA exhibits the characteristic swing as well as some phase jumps. The solid red line corresponds to the RVM fit to the PA, while the narrow grey band indicates the uncertainties of the fit result. The dashed line marks the RVM solution separated by 90 deg from the main one in order to include the jumped PA values. Details on the fit and the PA behaviour are discussed in Sect. 3.2.

At the right edge of the pulse we find a jump of clearly less than 90°, with an offset of only 60°–70° from the nominal PA swing. This indicates that these jumps do not originate purely from linear modes, but most likely from magnetospheric propagation effects creating circular modes as well (Petrova 2001, 2006; Melrose et al. 2006; Dyks 2020).

We can draw information about the geometry of J1618−3921 from the highly resolved swing of the polarisation angle across the main pulse. This can be explained by means of the RVM (Radhakrishnan & Cooke 1969): The emitted electromagnetic waves are polarised along the magnetic field lines, that point radially outward along the pulsar’s cone. As the beam moves across the line of sight, the observer sees these field lines under an ever changing angle (Lorimer & Kramer 2005). Exploiting basic geometric considerations, the RVM yields, for the position angle, ψ:

$$\begin{array}{c}\hfill tan(\psi -{\psi }_0)=\frac{sin\alpha sin(\varphi -{\varphi }_0)}{sin(\alpha +\beta )cos\alpha -cos(\alpha +\beta )sin\alpha cos(\varphi -{\varphi }_0)}\end{array}$$(1)

where α is inclination angle of the magnetic axis relative to the spin axis and ζ is the angle between the line of sight and the spin axis of the pulsar. This is connected to β (the minimum distance between the magnetic axis and the line of sight) via ζ = α + β (Lorimer & Kramer 2005). This minimum distance happens at spin phase ϕ₀; this is where ψ has the steepest slope, the corresponding PA of the linear polarisation is ψ = ψ₀. The angles in Eq. (1) are defined as in Radhakrishnan & Cooke (1969), that is, the ‘RVM/DT92’ convention (Damour & Taylor 1992). With the polarisation angle measurements from the MeerKAT observations (all data points in Fig. 2), we determine the RVM parameter posteriors in their joint parameter space following the method outlined in Johnston & Kramer (2019). The model also accounts for the possibility of OPM jumps and includes the corrected values in the fit. Keeping in mind the caveats associated with the RVM model (see e.g. Johnston & Kramer 2019), the results from the best fit model are shown in terms of corner plots in Fig. 3. Following Everett & Weisberg (2001), Johnston & Kramer (2019), Kramer et al. (2021), the results are presented using the RVM/DT92 convention. We obtain $\alpha =62.{27}_{-0.25}^{+0.26°}$ and $\zeta =110.{63}_{-0.93}^{+1.02°}$, quoting the 68% confidence levels on the posteriors.

Fig. 3. Corner plot showing the posterior distributions from fitting the RVM to the position angle variation observed by the MeerKAT telescope.

3.3. Change in profile with time

While inspecting the timing residuals we encountered an intriguing feature in the MeerKAT observations: starting with the observation from 2021-07-06, all residuals are offset by about 1 μs with respect to all residuals before that in the data set, while the MeerKAT residuals from observations before July 2021 align with the residuals from the other telescopes after fitting for a jump between them.

We found a change in the mean pulse profile to be the reason for the jump in the residuals. In Fig. 4, we show the summed profiles from all MeerKAT observations before the jump occurred, with a total of 26 hours, and from the 7 hours of observations since July 2022 that lead to the jumped ToAs, respectively. In the following we refer to the first one as the ‘pre-change profile’ and to the latter one as the ‘post-change profile’. Both profiles are generated by integrating the respective archives in time, frequency and polarisation. The first panel in Fig. 4 contains their difference (‘residual profile’), calculated by matching the pre- and post-change profile with the ProfileShiftFit subroutine from the PYTHON interface of PSRCHIVE (Hotan et al. 2004) and subtracting the re-scaled version of the latter one from the first one⁴. The underlying method of alignment is a χ²-fit of the Fourier-transformed profiles to each other to determine the respective phase shift and scale offset. The re-scaling process consists of applying the phase-rotation and overall intensity scaling of the fitting process to the latter profile. It is clearly visible that the profiles significantly differ from each other.

Fig. 4. Phase evolution of the total profile pre- (middle panel) and post (lower panel) change. The upper panel shows the difference between both.

To reassure ourselves that the change we see was actually occurring in June 2022, we performed a set of control analyses. To this end, we split the frequency and polarisation scrunched data from the pre- and post-change archives into two observations each. Then we repeated the subtraction procedure with these observations for all possible combinations. As was expected, the fitting amongst each own data set (pre with pre and post with post) yielded flat residual profiles in both cases. When cross-correlated (pre with post and vice versa), the shape of the deviation was reproduced when correlating the profiles between the two data sets. These results indicate that we are dealing with a genuine change in the mean pulse profile from July 2021.

A few of these profile changes have been reported in the literature over the past years. One prominent example of a DM-related profile change is found in the observations of J1713+0747 (Lin et al. 2021), that was originally associated with a DM-change. A characteristic for a DM-related profile change is a f⁻² frequency dependence; that means, this effect should dominate in the lower frequency bands. In contrast to that, the frequency dependence of the profile change of J1643−1224 (Shannon et al. 2016) excluded a DM-origin. Here, Shannon et al. (2016) point to changes in the emission region of the pulsar as being accountable for a change in the emission profile. These changes in the pulsar itself are responsible for profile changes. As a DM or magnetospheric origin of the change are difficult to distinguish, we investigated the MeerKAT observations further.

We performed a qualitative analysis of the frequency dependence by repeating the fitting and subtraction procedure on a per-sub-band basis. In doing so, we are unfortunately limited by the S/N of the observations. As we split all MeerKAT observations into eight sub-bands, we chose to display the frequency dependence at the same resolution as in Fig. 5. Evidently the deviation dominates in the lower frequency bands, but the nature of the change and the available S/N prevent us to confirm or refute a f⁻² dependence. The maximum frequency resolution feasible was sixteen sub-bands, where the deviations were most strongly visible in bands 0 to 2, weaker in bands 3 and 4, and absent from band 5 onward.

Fig. 5. Difference between the “pre-change” and “post-change” profiles on a per-channel basis. The profiles were created by integrating the observations of the respective time span in time and splitting them into eight frequency sub-bands. Before subtracting, both profiles were aligned using the ProfileShiftFit subroutine from the PYTHON interface of PSRCHIVE.

If the profile change were purely DM-related, we should be able to reproduce it by suitably altering the DM on the total pre-jump archive with the highest frequency resolution (928 channels) accordingly. After we scrunch this archive in frequency, it should give a similar residual profile as seen in Fig. 5 when compared to the pre-jump profile with the original DM. By fitting for DM and spin frequency on a per-observation basis, we retrieve the effective change from variations in the profile. By visually inspecting the resulting DM evolution, the profile change caused an alteration of around −0.01 pc cm⁻³ in the dispersion measure. We interpret this change as not physical, but caused by the impact on the fit of the profile change. Surprisingly, a reduction of the DM in the archive header by 0.01 pc cm⁻³ in the reverse engineering scheme laid out above, did not reproduce the profile change we show in Fig. 5. This is a strong indicator that the profile change is caused by magnetospheric changes, rather than by the ISM.

A change in the magnetosphere or the viewing geometry might also alter the polarisation properties of the radio beam. Thus, we assessed the difference of the PA across the total profile prior to and after the jump. We did not find any indications of a change.

Putting everything together, the frequency-resolved analysis of the jump points toward a non-ISM-related profile change, as we were not able to reproduce the profile change by introducing an artificial DM change for the pre-change observations (before July 2022). We point out that we could not investigate if the change could be caused by a strong scattering event, as our spectral analysis is limited by the steep spectral index and the subsequently low S/N in the upper bands.

Since July 2021 observations of J1618−3921 were not only conducted at MeerKAT, but also with the Parkes and Nançay radio telescopes; thus we inspected the other data sets for further traces of the timing jump. With only one observation from NRT in that time span we cannot make a meaningful statement concerning any impact of the profile change. In contrast, we have several observations at the Parkes radio telescope before and after the profile change. The summed profile resulting from the Parkes observations after July 2021 does not show any significant differences to the summed profile of the observations before that date. However, the mean ToA uncertainty of the Parkes observations is much larger than the size of the respective jump needed for the MeerKAT data set. Thus, we treat these ToAs jointly.

4. Timing analysis

4.1. Generating times of arrival

We produced the ToAs for all data sets using the standard template matching technique employed in pulsar timing: The ToAs are calculated by correlating a standard profile against a profile the actual each observation archive over polarisation and a suitable amount of time and frequency channels. The time and frequency resolution for each telescope is chosen in a trade-off against the resulting ToA precision, resulting in the number of frequency channels specified in Table 2. For frequency-resolved ToAs we created a frequency resolved standard profile by iteratively running PAAS on the integrated profile in each frequency channels. The ToAs were obtained via the PAT command. The significant decrease in intensity in the higher frequency channels for the MeerKAT and Parkes observations result in large ToA uncertainties in these bands. For the timing analysis, we carefully discarded these ToAs in order to reduce to computational load of the analysis without altering the fit results. At most MJDs we are still left with a frequency resolution of up to 9 (7) channels for the Parkes (MeerKAT) data, which is a large improvement to the previous work (Octau et al. 2018). Due to the low S/N of the earlier data from Parkes and the NRT, those ToAs were generated using the fully integrated observations, in other words, one frequency channel per observation.

4.2. Fitting timing models

To analyse the final data set containing 1535 ToAs we use the timing software package TEMPO2 (Hobbs et al. 2006), that performs a least-square minimisation of the residuals based on the χ² statistic as well as a Bayesian noise analysis using the TEMPONEST plugin. In contrast to the standard TEMPO2 usage, the TEMPONEST plugin relies on Bayesian parameter estimation, that (among other features) enables the fit for stochastic noise processes such as white noise, red timing noise and changes in dispersion measure using power law based models (Lentati et al. 2014).

The different data sets were combined by introducing a jump between each of them, with the MeerKAT data set before July 2022 as the reference data set. These jumps were treated as free fitting parameters in the TEMPO2 fit, while usually being marginalised over in the TEMPONEST analysis. Additionally, parts of the MeerKAT data set were corrected with known jumps.

By default, all ToA timestamps were recorded with an on-site reference clock. To be able to combine measurements from different telescopes, these are then converted to Coordinated Universal Time (UTC). Furthermore, UTC is converted to the main realisation of the terrestrial time (TT), the high-precision coordinate time standard called ‘International Atomic Time’ (TAI, temps atomique international). It is defined via the theoretically elapsed proper time on the Earth’s geoid and thus not prone to Earth’s rotational variations as UTC is. Finally the ToAs are transformed to the Solar System Barycentre (SSB), by accounting for the relative motion between each telescope and the SSB with JPL’s Solar System Ephemeris DE436.

For the binary orbit, TEMPO2 provides several models based on the calculations by Damour & Deruelle (1986) which provided a standard orbital model (henceforth “DD”). In this model, the orbital motion is parameterised by five Keplerian parameters (binary orbital period, P_b, longitude of periastron, ω, time of periastron passage, T₀, orbital eccentricity, e, and orbital semi-major axis projected along the line-of-sight, x) and a few additional “post-Keplerian” parameters that quantify, in a theory-independent way, the relativistic deviations from the Keplerian orbital motion. Relevant here are: the rate of change of the orbital period, Ṗ_b, the rate of periastron advance, $\dot{\omega }$, the Einstein delay, γ, (that quantifies the effects of the the varying gravitational redshift and special-relativistic time dilation) as well as the Shapiro delay, that affects the propagation time of the radio waves to Earth. In the DD model, the latter effect is parameterised using the ‘range’, r and ‘shape’, s, parameters. In GR, these are related to the companion mass, M_c, and the sine of the orbital inclination angle, ι, respectively (Damour & Taylor 1992).

Upon deriving a timing solution for J1618−3921 we analysed the ToAs with the theory-independent DDH model developed by Freire & Wex (2010), that differs from the DD model only in the parameterisation of the Shapiro delay: the new PK parameters (h₃ and ς) are less correlated than r and s, especially for systems with small orbital inclinations like J1618−3921. In addition, in a later stage of the analysis, we used the ‘DDGR’ model. Unlike the ‘DD’ and ‘DDH’ models, it is not theory-independent, but assumes that general relativity is the correct gravity theory, where no PK parameters are fit, only the total system mass and the companion mass. Due to the geometry of the system, the DDH model allowed for a more stable fit than the DDGR model.

After obtaining a first timing solution, that phase-connected the ToAs across the complete timing baseline, we updated the ephemeris in all available observations. With the new ephemeris installed, we repeated the entire process to obtain better profiles and standard templates. With these updated standards we then re-calculated the ToAs.

4.3. Bayesian timing and noise models

After deriving a final stable fit in TEMPO2 with the DDH model, we performed a Bayesian non-linear fit of the timing model by means of the TEMPONEST software package. Using the parameters from the TEMPO2 output ephemeris as the input for TEMPONEST, we derived a timing solution that additionally accounted for the commonly known noise parameters: Unrecognised systematics in the ToA uncertainties are modelled by the white noise parameters EFAC F and EQUAD Q on a per-backend basis. Therefore the uncertainty σ_ToA, old of each ToA is re-scaled as ${\sigma }_{\mathrm{ToA},\mathrm{new}}=\sqrt{Q^2+F^2{\sigma }_{\mathrm{ToA},\mathrm{old}}}$ (Lentati et al. 2014). For the chromatic models, we obtained an amplitude, A, and a spectral index, γ (Lentati et al. 2014).

In order to find the best-fitting chromatic noise model, we proceeded in a two-fold way: On the one hand we compared the evidence returned by the sampler MULTINEST (Feroz et al. 2019) for different combinations of noise models (Red noise (RN) only, DM noise only, Red and DM noise). On the other hand we also varied the number of noise model coefficients between 45, 60 and 100, and compared the resulting time-domain realisation between the different models. The realisations were produced using the methods of the LA FORGE github repository⁵ (Hazboun 2020) adapted for the relevant models at hand. The most favoured models were the 60 and 100 coefficient DM-only models, with a difference in the log-evidences of 29. From comparing 100 averaged realisations of both noise models to the ToAs, we decided to chose the 100 coefficient model, as it visibly reflected the ToA changes more precise than the 60 coefficient model. The respective time-domain noise realisations are shown as the blue lines in the lower plot of Fig. 6. TEMPONEST accounts for the DM noise in terms of a power law model (Lentati et al. 2014), where for the chosen model we have an amplitude of A_DM = −10.37 and a slope of γ_DM = 0.94. This slope is exceptionally shallow for a noise process whose slope is usually expected to be of the order of 2. From Fig. 6 we can deduce that the residuals exhibit some significant small-timescale variations that might give rise to the shallow slope. Nevertheless, the time-domain noise realisation in the lower plot of Fig. 6 shows that the noise model seems to match the visible trends in the data; hence, we regard the noise model as satisfactory.

Fig. 6. Post-fit timing residuals of PSR J1618−3921. Upper two plots: Post-fit timing residuals as a function of time for the best-fit timing solution (TEMPONEST fit without removing the 100 coefficient DM noise model) of PSR J1618−3921 given in Table 3. Lower plot: Post-fit timing residuals after subtracting the DM noise model. The colours denote the different backends and systems as listed in Table 2. In all plots, the uncertainties are re-scaled with the White Noise parameters EFAC and EQUAD (cf. Table 2). The middle plot also contains the time domain realisation of the 100 parameter DM noise model: the blue lines show the 100 randomly created model realisations, and the black dots indicate the median across all these at each ToA.

As the data set exhibits large gaps in the beginning of the observations, we also investigated the covariance between the jumps and the timing parameters setting up a TEMPONEST analysis, where the jumps are also treated as free parameters. We did not find a significant change in any of the timing parameters.

The timing parameters of the best-fit solution from TEMPONEST using the DDH model are presented in the third column in Table 3. Each parameter is quoted as the maximum of the marginalised posterior together with the respective left and right 39% confidence limits. The timing residuals achieved from this solution are shown in Fig 6. Table 3 also shows the corresponding parameters reported by Octau et al. (2018), with blank entries when the parameter was fit for for the first time in the scope of this work. In Fig. 7 we show both the 2D-correlation contours and the 1D posterior distributions resulting from the TEMPONEST analysis for a chosen subset of fitted parameters we attribute a higher relevance in this work.

Fig. 7. Corner plot for the relevant subset of timing parameters from the DDH model for J1618−3921 derived by applying the DDH binary models to the ToA data set applying non-linear Bayesian timing techniques using the software TEMPONEST to the ToAs set shown in Fig. 6. The diagonal elements show the 1D marginalised posterior distributions for each parameter, the shaded region indicates the 1σ credibility interval. The 2D contours populating the off-diagonal elements show the correlation between pairs of parameters, where the lines mark the 39%, 86% and 98% credibility regions, going from dark to light shaded.

In the following, we present the individual timing parameters in greater detail and and discuss their implications for the binary system based on the numeric values derived from the best-fit TEMPONEST solution.

4.4. Position and proper motion

As usual, the timing solution provides the pulsar’s position with very high accuracy. With a location at RA (J2000) 16h 18′18.824940(38)″ and DEC (J2000) −39° 21′01.815(10)″, we searched the second data release of the DECam Plane Survey (DECaPS2) (Schlafly et al. 2018), a five-band optical and near-infrared survey of the southern Galactic plane, using the Aladin Lite web interface⁶ (Baumann et al. 2022). The corresponding excerpt from the survey image with a field of view of about 17 as around the pulsar’s position is shown in Fig. 8. At the position of the pulsar (indicated by the purple hair-cross on the image), we cannot identify any counterpart for either the pulsar or its companion. This implies that the electromagnetic emission of both bodies is below the detection thresholds of this survey, that are quoted to 23.7, 22.7, 22.2, 21.7, and 20.9 mag in the grizY bands (Schlafly et al. 2018).

Fig. 8. Excerpt from the DECaPS2 survey (Schlafly et al. 2018) visualisation taken from Aladin Lite (Baumann et al. 2022) in a field of view of ∼17 as around the position of J1618−3921. The pulsar’s position is indicated by the purple hair-cross. The nearest sources are enumerated with the numbers 1–3.

We are able to measure both the proper motion in Right Ascension ${\mu }_{\alpha }=1.{24}_{-0.13}^{+0.14}\mathrm{mas}{\mathrm{yr}}^{-1}$ and Declination μ_δ = ( − 2.37 ± 0.35) mas yr⁻¹. This leads to a total proper motion of ( − 2.5 ± 0.3) mas yr⁻¹. Furthermore, combining the timing model value of the dispersion measure DM with models of the electron distribution of the Galaxy, we infer a distance to the pulsar of 2.7 to 5.5 kpc. For the lower boundary to the distance window we apply the NE2001 model (Cordes et al. 2002), the upper boundary is based on the YMW16 model (Yao et al. 2017). Using the distance from the NE2001 model, we translate the measured proper motions into the heliocentric velocity of the binary system of v_T = (33 ± 4) km s⁻¹.

4.4.1. Spin-down and higher frequency derivatives

An important quantity describing a pulsar’s properties is the intrinsic spin down, Ṗ_int. For a pulsar at a distance, d, moving with a relative proper motion, μ, any time-related measurement is influenced by the change in the Doppler shift. Thus we correct the precisely measured period derivative $\dot{P}=-(5.37620\pm 0.00068)\times {10}^{-20}$ to

$$\begin{array}{c}\hfill \frac{{\dot{P}}_{\mathrm{int}}}{P}=\frac{{\dot{P}}_{\mathrm{obs}}}{P}+\frac{\dot{D}}{D},\end{array}$$(2)

where D is the Doppler factor caused by the unknown radial velocity of the pulsar and $\dot{D}$ its derivative. Although neither D or $\dot{D}$ are known, their ratio can be estimated as:

$$\begin{array}{c}\hfill \frac{\dot{D}}{D}=-\frac{1}{c}[{\mathit{K}}_0·({\mathit{a}}_{\mathrm{PSR}}-{\mathit{a}}_{\mathrm{SSB}})+\frac{V_{\mathrm{T}}^2}{d}]=-\frac{a}{c}-\frac{{\mu }^{2}d}{c},\end{array}$$(3)

where the first term holds the contribution of the line-of-sight acceleration, a, by projecting the difference between the Galactic acceleration at the position of the pulsar, a_PSR, and the solar system barycenter (SSB), a_SSB, onto the unit vector, K₀, pointing from the Earth to the pulsar. The second term, that depends on the transverse velocity, V_T, and the distance to the pulsar, d, is the Shklovskii term (Shklovskii 1970).

We obtain suitable values of the Galactic acceleration at the SSB and the position of the pulsar using the Milky Way mass model presented by Lazaridis et al. (2009). For the position and velocity of the solar system barycenter we assumed R_⊙ = 8.275 ± 0.034 kpc and V_⊙ = 240.5 ± 4.1 km s⁻¹ (GRAVITY Collaboration 2021).

Using the NE2001 distance estimate, the Shklovskii effect contributes Pμ²d/c = 6.2 × 10⁻²². The Galactic acceleration field partly compensates for this effect with an excess period change of Pa/c = −1.4 × 10⁻²². We therefore arrive at an intrinsic spin-down of ${\dot{P}}_{\mathrm{int}}=5.33326\times {10}^{-20}$, that is only slightly smaller than Ṗ_obs (Ṗ_int = 0.991 Ṗ_obs).

Furthermore, with $\ddot{f}=-1.0(2){10}^{-27}{\mathrm{s}}^{-3}$ we find a non-zero value of the second derivative of the spin frequency. This value is multiple orders of magnitude larger than what is expected from a pure spin-down (𝒪(10⁻³³), assuming a characteristic age of 10 Gyr and a braking index of 3) and among the very few values of $\ddot{f}$ measured for the 333 pulsars with P < 30 ms. Outside of globular clusters, only 9 measurements of $\ddot{f}$ have been made Manchester et al. (2005), mostly for highly energetic gamma-ray MSPs, where timing noise could be happening, additionally some of these systems are in “black widow” binaries with strong outgassing. In one case (J1024−0719), the pulsar is known to have a distant companion, a K dwarf Bassa et al. (2016), in another case, J1903+0327, the system is thought to have formed in a triple system that later became unstable Freire et al. (2011); perhaps the third object was not fully ejected and is still somewhere in the vicinity of the system. A comparison of the timing residuals for the timing models with and without this parameter is shown in the upper plot of Fig. 9. Higher derivatives of f are likely to originate from a varying acceleration along the line of sight of the binary system. The implications of the measurement of $\ddot{f}$ on the nature of the system and other timing parameters will be discussed in more detail in Sect. 5.3.

Fig. 9. Demonstration of the best-fit timing solution for J1618-3921. Top: Timing residuals as a function of time for the best-fit timing model with (black) and without (red) considering the second derivative of the rotational frequency $\ddot{f}$. Bottom: Timing residuals as a function of orbital phase for the best-fit timing model with (black) and without (red) considering the derivative of the orbital period Ṗb.

4.5. Post-Keplerian parameters

4.5.1. Rate of advance of periastron

The orbital eccentricity of the system and the long timing baseline allow a highly significant measurement of the rate of advance of periastron, despite the wide orbit: $\dot{\omega }=0.{00142}_{-0.00010}^{+0.00008}{\mathrm{yr}}^{-1}$. If this effect is purely relativistic, it yields a direct measurement of the total mass of the system, M_tot.

In order to gauge the reliability and meaning of the measurement of $\dot{\omega }$, we have to consider the possibility of additional non-relativistic effects. The most important of these is a proper motion contribution, ${\dot{\omega }}_{\mu }$. This contribution is given by (Kopeikin 1996)

$$\begin{array}{c}\hfill {\dot{\omega }}_{\mu }=\frac{\mu }{sin\iota }cos({\mathrm{\Theta }}_{\mu }-\mathrm{\Omega }),\end{array}$$(4)

where Θ_μ is the proper motion position angle and Ω the position angle of the line of nodes. Assuming an optimal alignment (cos(Θ_μ − Ω) = 1), it contributes at the order of ${\dot{\omega }}_{\mu }\sim 8\times {10}^{-7}°{\mathrm{yr}}^{-1}$.

As is discussed in Rasio (1994), Joshi & Rasio (1997), a third body in the system can add a contribution to the observed periastron advance:

$$\begin{array}{c}\hfill {\dot{\omega }}_{\mathrm{triple}}={\left(\frac{\dot{x}}{x}\right)}_{\mathrm{triple}}\frac{2[{sin}^2{\theta }_3(5{cos}^2{\mathrm{\Phi }}_3-1)-1]}{cot\iota sin2{\theta }_{3}cos(\omega +{\mathrm{\Phi }}_3)}.\end{array}$$(5)

Including ẋ in the timing model fit yields $\dot{x}=(2\pm 8)\times {10}^{-15}$, which is consistent with a non-detection. Considering that the geometric terms in Eq. (5) contribute at 𝒪(1), the fit value of ẋ gives an upper limit to the contribution of the periastron advance from the putative third body of ${\dot{\omega }}_{\mathrm{triple}}<3\times {10}^{-7}°{\mathrm{yr}}^{-1}$.

Compared to the measured rate of advance of periastron, both contributions are negligibly small, so we conclude that the measured value of $\dot{\omega }$ is within measurement precision, relativistic. The relativistic $\dot{\omega }$ relates to the total mass of the system as

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

where $T_{\odot }\equiv {\mathcal{GM}}_{\odot }^{\mathrm{N}}/c^3=4.92549094764126\cdots \mu \mathrm{s}$ is an exact quantity that follows from the exact definitions of the speed of light, c, and the solar mass parameter ${\mathcal{GM}}_{\odot }^{\mathrm{N}}$(Prša et al. 2016). From the best-fit parameters, we derive a total mass estimate of $1.{42}_{-0.19}^{+0.20}{\mathrm{M}}_{\odot }$. Comparing this result with the mass measurements for similar NS-WD binaries⁷, we find that our measurement lies well within the expected mass range.

4.5.2. Shapiro delay

With J1618−3921 being a pulsar in the RelBin programme, one of the main aims of this work is achieving a significant Shapiro delay measurement by means of the high timing precision that comes along with MeerKAT observations. The rather low flux density, combined with a low inclination angle made a precise measurement of the Shapiro delay difficult. We were able to stabilise the DDH model based TEMPO2 fit with ToAs gained from a dedicated superior conjunction observation campaign toward a low-significance detection of the Shapiro delay. From the TEMPONEST analysis we found $h_3=2.{70}_{-1.47}^{+2.07}\mathrm{\mu }\mathrm{s}$ and $\varsigma =0.{68}_{-0.09}^{+0.13}$. In order to convert these measurements and the measurement of $\dot{\omega }$ into constraints on the mass and the inclination angle of the system, we perform a χ²-grid analysis of the M_PSR − cosι space (cf. Sect. 4.5.3). The unconstrained inclination angle in the right plot of Fig. 10 resulting from the analysis demonstrates that we did not arrive at a significant measurement of the Shapiro delay.

Fig. 10. Constraints on the orbital inclination and masses of the J1618−3921 binary. The dashed red lines correspond to constraints derived from the TEMPONEST fit for the rate of advance of periastron $\dot{\omega }$ presented in Table 3 assuming the validity of GR and regarding the effect as purely relativistic. The black lines include 68.3 and 95.4% of all probability of the 2-dimensional probability distribution functions (pdfs) derived from the χ2 map. The left of the two inner plots shows the Mc-cosι diagram, that was sampled evenly. The grey area is excluded because MPSR must be positive. The plot to its right shows the projection of the Mc-cosι pdf on the MPSR-Mc space using the mass function. The grey area is excluded because sinι ≤ 1. The outer three plots display the projected distributions for cosι, MPSR, and Mc. the hatched area corresponds to the 1σ intervals. The unconstrained inclination angle shows that we have a non-detection of the Shapiro delay; thus, we do not provide the confidence intervals.

4.5.3. Mass measurement

We now estimate the masses with the highly significant detection of $\dot{\omega }$ and the weak Shapiro delay constraints using the analysis technique outlined in Barr et al. (2017). At each grid point corresponding to a (M_PSR, M_c)-pair we fix the respective values of M_tot and M_c in a DDGR ephemeris adapted from the actual TEMPONEST results, that is then used in a TEMPO2 fit. With the two mass values, the DDGR model self-consistently accounts for all observed relativistic parameters except for the orbital decay, where we know there are large contributions from other causes. The goodness of the fit is quantified by the χ² value of the TEMPO2 fit, where a lower χ² value describes a better fit. The result is a map of χ² values across the M_PSR-M_c-grid, that can be translated into credibility contours by subtracting the global minimum value across the map from all map points. The result is displayed in the mass-mass diagram in Fig. 10, together with the credibility band from the rate of advance of periastron. With this method we constrain the companion mass to $0.{20}_{-0.03}^{+0.11}{\mathrm{M}}_{\odot }$, the pulsar mass to $1.{20}_{-0.20}^{+0.19}{\mathrm{M}}_{\odot }$ and the total mass to $1.{42}_{-0.19}^{+0.20}{\mathrm{M}}_{\odot }$ (68.3% confidence limits).

4.5.4. Change in orbital period

The impact of the change in the orbital period on the timing residuals is shown in the lower plot of Fig. 9. Similar to the measurement of the spin period derivative, the observed rate of change in the orbital period is the sum of various effects,

$$\begin{array}{c}\hfill {\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{obs}}={\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\mathrm{GW}}+{\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}^{\dot{m}}-\frac{\dot{D}}{D},\end{array}$$(7)

where apart from the kinematic contributions ($\dot{D}/D$), also emission of gravitational waves (GW) and mass loss from the system (ṁ) might significantly contribute to the measured value. Evaluating the expressions given in Lorimer & Kramer (2005) for the latter two effects, we find ${({\dot{P}}_{\mathrm{b}}/P_{\mathrm{b}})}^{\mathrm{GW}}\sim -1\times {10}^{-23}{\mathrm{s}}^{-1}$ and ${({\dot{P}}_{\mathrm{b}}/P_{\mathrm{b}})}^{\dot{m}}\sim 4\times {10}^{-28}{\mathrm{s}}^{-1}$. Compared to our measured value, these contributions are negligible.

Thus, the only significant term comes from $-\dot{D}/D$. Using the value calculated in section 4.4.1, we obtain ${\dot{P}}_{\mathrm{b}}=-\dot{D}/D{P}_{\mathrm{b}}\sim +0.05\times {10}^{-12}$. Surprisingly, the best-fit timing model reveals a measured orbital period change of $-2.2_{-0.33}^{+0.35}\times {10}^{-11}$. This is not only two orders of magnitude larger than expected, but also carries an opposite sign. All considered effects are multiple orders of magnitude too small to provide an explanation for the large observed value of Ṗ_b. A possible solution to this tension is the presence of an additional acceleration caused by a third body in the vicinity of the binary, as is discussed in Sect. 5.3.

4.5.5. Other parameters

As for similar systems (cf. Serylak et al. 2022), we are not able to obtain a significant measurement of the Einstein delay amplitude, γ, or any variation in the projected semi-major axis, ẋ, since their contributions to the residuals are beyond the current precision of our ToAs; furthermore, given the orbital periods of these pulsars, the timing effect of ẋ and γ are strongly correlated (Ridolfi et al. 2019). Moreover, we do not detect derivatives of the spin frequency higher than $\ddot{f}$.

5. Discussion

5.1. Comparison to Octau et al.

In comparison to the work by Octau et al. (2018), we use data not only from NRT, but also Parkes and MeerKAT, including observations that span back to 1999. These early observations were available previously, but only the high quality of the MeerKAT observations, together with the observation density achieved by combining three radio telescopes guaranteed a timing solution that was robust enough to extend the timing model back to 1999, through a very sparse set of observations. This large timing baseline, plus the precise recent timing, allows for the measurement of timing parameters that were not previously available: proper motion, of higher order spin and DM derivatives and post-Keplerian parameters. We are also able to significantly improve on the measurement and variation of the DM. In comparison to the four frequency channels obtained from the third NRT observation run, the large-bandwidth observations with the Parkes UWL receiver have a S/N that allows us to separate them into 13 frequency channels with often a reasonable ToA precision. Although we have to discard the ToAs from the high frequency channels, we still achieved a significant refinement in the frequency resolution compared to the previous work.

5.2. Orbital geometry

If the spin of the pulsar in a binary is aligned with the orbital angular momentum, the inclination angle ι coincides with the viewing angle ζ. But upon comparing the timing result for the Shapiro delay parameter to the RVM fit results, there are two major caveats: First, in fitting for the Shapiro delay, we determine sinι. Hence, we cannot distinguish if the corresponding inclination angle is ι or 180° −ι. In case of a reliable RVM fit, this ambiguity can be solved by comparing ι to ζ. This can also not be done directly, since the above RVM equation assumes that ψ increases clockwise on the sky, opposite to the astronomical convention, where ψ increases counter-clockwise from the above equation (Damour & Taylor 1992; Everett & Weisberg 2001; van Straten et al. 2010). Hence we have to identify the RVM fit value for ζ with 180° −ι or ι with 180° −ζ, respectively (Kramer et al. 2021).

Taking both these aspects into account, we first of all find with the reference angle from the RVM fit of $180°-\zeta =69.{37}_{-0.93}^{+1.02°}$, that sinι translates into ι = (66 ± 14)°. This is also confirmed by performing two further RVM fits in which we restricted the variation in ζ to one of the ranges allowed by the timing results (cf. Sect. 4) on sinι, respectively.

Although the viewing angle from the RVM fit is consistent with the inclination angle from the timing solution (Table 3), we cannot make any reasonable statement about an alignment or misalignment of both axes due to the highly unconstrained Shapiro delay.

5.3. The cause of the anomalous Ṗ_b and $\ddot{f}$

The significant deviation between measurement and prediction shows that there is another contribution to the pulsar’s acceleration. This additional acceleration completely dominates the expected Galactic gravitational acceleration. Such a strong gravitational field could be produced by a massive nearby object. We can test this hypothesis in a simple way. If the observed Ṗ_b is caused by an unexpected acceleration (and therefore implying a larger than assumed $-\dot{D}/D$ term), then we should be able to re-compute the spin-down of the pulsar using this term, as measured by Ṗ_b, and still obtain a positive value. Subtracting Eq. (7) from Eq. (2), and neglecting the GW emission terms, we obtain:

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{int}}={\dot{P}}_{\mathrm{obs}}-P{\left(\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}\right)}_{\mathrm{obs}},\end{array}$$(8)

since ${\dot{P}}_{\mathrm{b},\mathrm{obs}}$ is negative, this has the effect of increasing our estimate of Ṗ_int to ∼1.8(3)×10⁻¹⁹, that is ∼3.4 times larger than the observed Ṗ. From this value, we estimate the characteristic age, the spin-down luminosity as well as the surface magnetic field of the pulsar (Lorimer & Kramer 2005). These values can be seen in Table 3, there we see how the change in the value of Ṗ_int between this work and the work from Octau et al. (2018) lead to significant differences in the values of τ_c and Ė.

For pulsars at a low Galactic latitudes, this additional acceleration might be caused by massive molecular clouds in their vicinity. With J1618−3921 located at b = 7.9°, this is unlikely, but not impossible. Another option is that a third body is in a wide orbit around the PSR-WD binary.

The measurement of the second derivative of the spin frequency helps to distinguish between these two scenarios. A molecular cloud would be located at a large distance to the binary and thus its acceleration would appear to be constant; in this case, we would not expect large variations in the line-of-sight acceleration and thus on the $\dot{f}$. Instead, we measure a large $\ddot{f}$ of ( − 1.0 ± 0.2)×10⁻²⁷ s⁻³, that is very likely caused by a variation of the external acceleration. This is a strong indicator that the source of the acceleration is in the vicinity of the binary. Thus we propose that the system is a hierarchical triple system.

This line of arguments is strongly motivated by a similar discussion of the J1024−0719 system (Bassa et al. 2016). Upon its discovery, it was regarded as an isolated pulsar, but the measurement of higher-order spin frequency derivatives led Bassa et al. (2016) to propose a companion in an extremely wide orbit (P_b > 200 yr). This was confirmed by the detection of a nearby star with the same proper motion. Comparing the measured value $\ddot{f}$ for both pulsars, we find that the value for J1618−3921 is a factor of two smaller than for J1024−0719 and thus of a very similar order of magnitude. With the measurement of Ṗ_b we even have the advantage of estimating the acceleration of the inner binary system - this is not possible for J1024−0719, because that pulsar is not already in a binary system.

Keeping in mind that most stars are part of multiple systems, it is no surprise that on rare occasions, binaries with a pulsar are actually part of a higher-order stellar system. Due to stability arguments (Toonen et al. 2016), most of these systems are hierarchical triple systems, that means, they consist of an inner binary in a wider orbit around a third object.

An example is the well-known triple system consisting of the MSP J0337+1715 (Ransom et al. 2014). Detailed timing of this system (Ransom et al. 2014; Archibald et al. 2018; Voisin et al. 2020) revealed that both orbits of the system are co-planar and circular and the WD masses are as predicted by TS99 relation, as was expected from adopting the previously discussed WD-MSP formation scenario (Tauris & van den Heuvel 2023). On the other hand, Toonen et al. (2016) show in a broad study on triple systems that the unique dynamic in these systems also allows for a stable eccentric inner binary. They also point out that mechanisms such as Lidov-Kozai cycles prevent a synchronisation and circularisation of the binary, leading to MSP systems that stand in complete contrast to the formation scenario described by Tauris & van den Heuvel (2023).

If PSR J1618−3921 really has a stellar companion, all derivatives of f are expected to eventually converge on a Keplerian orbit for the outer component (Rasio 1994). Here, J1024−0719 again serves as a precedent; we should consider these MSP companions to be also settled in exceptionally wide orbits. Any associated parameter derivative is therefore expected to show up only in data sets with a combination of a long timing baseline and significant timing precision. Determining the orbital configuration of the outer companion would require the knowledge of at least the first five derivatives of f (Rasio 1994; Joshi & Rasio 1997)⁸. With the knowledge of fewer derivatives, we can only put a few constraints on the orbit (Bassa et al. 2016): Ṗ_b relates to the corresponding acceleration from the third body, a, as $a/c\sim {\dot{P}}_{\mathrm{b}}/P_{\mathrm{b}}$. Similarly $\ddot{f}$ relates to the change in the acceleration as $\dot{a}/c\sim \ddot{f}/f$. From the acceleration and its change, we can place an order-of-magnitude estimate on the orbital period of the third body as $P_{\mathrm{b},3}\sim a/\dot{a}\sim 300\mathrm{yr}$, given the values form our best-fitting timing solution. This is not unexpected, and also highly in line with the findings from Bassa et al. (2016) in the case of J1024−0719.

5.3.1. Optical counterpart

We consulted the DECaPS2 (Schlafly et al. 2018) catalogue to search for a spatially resolved object that could be associated with the PSR J1618−3921 system, and thus be identified as the binary companion or the putative third body. As was mentioned in Section 4.4, no counterparts are identified near the position of PSR J1618−3921. The upper mass limit of any companions (either the binary companion to the pulsar, or the more distant object) can thus be estimated with the depth of the catalogue through comparisons with the expected colours and magnitudes from stellar evolutionary models. We have used the PAdova TRieste Stellar Evolutionary Code (PARSEC v2.0 Bressan et al. 2012; Nguyen et al. 2022) to obtain the grizY magnitudes in the ABmag system to facilitate comparisons with the DECaPS2 catalogue. Applying an extinction A_V ∼ 0.2 mag⁹ and adopting a distance of 5.5 kpc, a 0.56 M_⊙ dwarf star (∼M0V,Pecaut & Mamajek 2013) would have grizY = 23.9, 22.5, 21.7, 21.3, 21.2 mag, respectively. Such a star would be near the detectability limit of the riz bands in the DECaPS2 survey, given its limiting magnitudes 23.7, 22.7, 22.2, 21.7, and 20.9 mag in grizY bands, respectively, and would have been detected in all 5 bands if a smaller distance of 2.7 kpc is adopted. To summarise, any companion at the location of PSR J1618−3921 would have a limiting magnitude detection threshold of 23.5 mag in the G-band, that at the distance to the pulsar of 5.5 kpc corresponds to an absolute magnitude > 9.79 mag. This could be a M-dwarf of mass < 0.56 M_⊙ or a compact object.

5.3.2. Nearby stars, their motions, and their gravitational accelerations

We consulted the Gaia DR3 (Gaia Collaboration 2023)¹⁰ to search for objects that might have a proper motion similar to that of PSR J1618−3921; that is, within the ±3σ error ellipse. This was, incidentally, how the distant binary companion of PSR J1024−0719 (a K7V star) was identified Bassa et al. (2016). No objects with such a proper motion are detected within a radius of 1.4′ around PSR J1618−3921. Using the NE2001 distance for a lower limit, this corresponds to a minimum distance of 0.8 pc.

In the deeper DECaPS2 catalogue (Schlafly et al. 2018) we find three nearby stars; shown in Fig. 8, at a distance of 2″, 4″ and 2″ (following the labels 1 to 3) from PSR J1618−3921. Given the depth of this catalogue, these faints stars are not in the Gaia DR3 catalogue, so an association with PSR J1618−3921 cannot be excluded based on proper motion measurements. Under the assumption that the three objects are stellar type objects and that they are at the same distance as the pulsar, we have extracted their grizY magnitudes from the DECaPS2 catalogue to estimate their masses. We use Star 1 as an example as it has measurements in all 5 bands: 23.3, 22.0, 21.3, 20.7, 20.5 mag, respectively. These magnitudes are in agreement with those for a 0.6 M_⊙ (or K9V) star with A_V ∼ 0.2 and a distance of 5.5 kpc: 23.4, 22.0, 21.3, 21.0, and 20.9 mag, respectively.

For each of the three stars, we make an order-of-magnitude estimate of the line-of-sight (LOS) acceleration, a_LOS, they exert on the pulsar, respectively. Assuming a typical mass of 0.6 M_⊙ for these stars, we derive the estimate via Newton’s law $a_{\mathrm{LOS}}=\frac{\mathit{GM}}{d^2}sin\alpha$, where M denotes the mass of the star, G is the gravitational constant, d is the distance between the pulsar and the star and α is the angle between the vector pointing from the star to the nearest point on the LoS and the vector pointing from the star to the pulsar. Turning the angular distance taken from ALADIN into the physical separation, we use the NE2001 distance, as it gives us an upper limit on the acceleration. This inferred separation is the projected distance r between the pulsar and the star, so d = r/cosα. The resulting acceleration curves calculated under the previously outlined assumptions for the three objects marked in Fig. 8 are shown in Fig. 11.

Fig. 11. LOS acceleration aLOS as a function of the projection angle α onto the LOS for the three objects in the DECaPS2 (Schlafly et al. 2018) catalogue found closest to the position of J1618−3921. The acceleration was calculated from Newton’s first law assuming they are K- or M-stars with a mass of 0.6 M⊙.

All these objects cause accelerations that are roughly two orders of magnitude smaller than the acceleration a_LoS,Pb = −3.45 × 10⁻⁹ m s⁻² obtained from Ṗ_b. Hence the putative wide-orbit companion of J1618−3921 must be closer than these objects, and must have a luminosity below the DECaPS2 (Schlafly et al. 2018) limit: as was mentioned above, it could be an M-dwarf or a compact object.

6. Summary

This paper presents a comprehensive overview of the latest knowledge about the eMSP J1618−3921 using the combined data set from 23 years of observations with Parkes, NRT and MeerKAT radio telescopes and their respective different back-ends.

We present a detailed study on the pulsar’s emission properties with two notable results: First we recorded a profile change that happened around June 2021 with the MeerKAT observations. Our analyses favours an intrinsic profile change over an ISM-related influence, but due to the limited S/N in the upper MeerKAT frequency bands, we cannot finally determine the origin of this change. Furthermore we analysed the behaviour of the position angle of the linear polarisation. Assuming a purely dipolar radio emission, with the PA perfectly following the RVM, we constrained the position of the spin axis of the pulsar to (111 ± 1)°. The uncertainty in the orbital inclination precludes any conclusions on the alignment of the spin axis of the pulsar with the orbital angular momentum.

While in previous publications (Bailes 2007; Octau et al. 2018), orbital and then phase-coherent timing solutions were already published, here we not only report the old timing with significantly improved precision, but we provide the first solution including a binary model with an increased number of post-Keplerian parameters. The stability of the solution is mainly provided by the dense accumulation of data points from joint MeerKAT, Parkes and NRT observations in the recent past. This allowed us to include all available observations up to the very first observations from 1999. This large timing baseline significantly improved the measurement of rate of advance of periastron.

Although the ToAs obtained from monthly observations with the MeerKAT L-band receiver exhibit an outstanding precision compared to ToAs resulting from concurrent observations at the Parkes and Nançay radio telescopes, the low S/N of the pulsar as well as the shallow inclination angle impeded a high-significance detection of the Shapiro delay. Nevertheless we are able to present a first constraint on the orthometric parameters h₃ and ς. Combining the low-significance Shapiro delay detection with the precise measurement of the rate of advance of periastron we are able to present the first ever mass estimates of this system. Unfortunately, the steep spectral index prevents us from obtaining more precise ToAs using the S-band (1.75–3.5 GHz) receiver at MeerKAT. With a factor of two to three improvement in timing precision, the Shapiro delay should be measurable with useful precision, but this will only be possible with future radio telescopes of even higher sensitivity.

The most remarkable result of the timing analysis is the amount of change in the orbital period and the large second derivative of the spin frequency, that indicate that the pulsar is actually part of a triple system. The possibility of the evolution of J1618-3921 as a triple system opens the door for similar evolution of other eMSPs. However, there is at the moment no clear evidence that other eMSPs are part of hierarchical triple systems.

Our long-term plan for this pulsar consists of regular observations of J1618−3921 with the L-band receiver at MeerKAT and the UWL receiver at the Murriyang Parkes radio telescope. We expect that the increased timing baseline will significantly improve all currently measured parameters, but also enable the detection of additional parameters such as ẋ (that will constrain the orbital orientation of the system) or higher derivatives of f, that will provide additional information on the companion mass and its orbit.

Acknowledgments

The authors thank Aurélien Chalumeau, Michael Keith and Aditya Parthasarathy for the fruitful discussions and support regarding the noise model comparison and time domain realisation. All authors affiliated with the Max-Planck-Gesellschaft (MPG) acknowledge its constant support. VVK acknowledges financial support provided under the European Union’s Horizon Europe 2022 Starting Grant “COMPACT” (Grant agreement number: 101078094, PI: Vivek Venkatraman Krishnan). The MeerKAT telescope is operated by the South African Radio Astronomy Observatory (SARAO), a facility of the National Research Foundation, that is 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), as well as the time space reference systems department of the Paris Observatory. MeerTime data is housed on the OzSTAR supercomputer at Swinburne University of Technology, on which significant parts of the data reduction was performed. Parts of this research were supported by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. The Parkes radio telescope (Murriyang) 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. 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. Parts of the data set used in this work include archived data obtained through the CSIRO Data Access Portal (http://data.csiro.au). It also made broad use of the NASA Astrophysics Data System (https://ui.adsabs.harvard.edu/). The analysis done in this publication made use of the open source pulsar analysis packages PSRCHIVE (Hotan et al. 2004), TEMPO2 (Hobbs et al. 2006) and TEMPONEST (Lentati et al. 2014), as well as open source PYTHON libraries including Numpy, Matplotlib, Astropy and Chainconsumer. APo and MBu acknowledge the support from the research grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016. APo and MBu acknowledge that part of this work has been funded using resources from the INAF Large Grant 2022 “GCjewels” (P.I. Andrea Possenti) approved with the Presidential Decree 30/2022.

References

All Tables

All Figures

Fig. 1. Frequency resolved intensity profiles from observations with the MeerKAT L-band receiver (left, Tobs ∼ 28.85 h), the Parkes UWL receiver (middle, Tobs ∼ 28.95 h) and the NRT L-band receiver (right, Tobs ∼ 48.9 h). The top panel of each plot shows the total intensity profile across one period, in case of the MeerKAT and Parkes observations it is flux calibrated. The NRT observations are not flux-calibrated. The bottom panels show the frequency resolved dynamic spectra across adjacent frequency bands. The intensity scale of the MeerKAT and NRT dynamic spectra were adjusted to fit the range of the Parkes spectrum for ease of comparison. The MeerKAT observations were frequency scrunched to 8 channels, those from Parkes down to 13 channels and the NRT data was decimated to 4 channels. The number of channels was chosen such that the frequency resolution is kept as large as possible while providing a S/N in each band that allows for a sufficient ToA precision.

In the text

Fig. 2. Polarisation profile of J1618−3921 obtained from integrating 29 hours of observations with the MeerKAT L-band receiver. The upper part shows the total intensity (light blue), as well as the linear (red) and circular polarisation (dark blue) fraction. The lower part shows the evolution of the position angle (PA) across the pulsar’s phase. The PA exhibits the characteristic swing as well as some phase jumps. The solid red line corresponds to the RVM fit to the PA, while the narrow grey band indicates the uncertainties of the fit result. The dashed line marks the RVM solution separated by 90 deg from the main one in order to include the jumped PA values. Details on the fit and the PA behaviour are discussed in Sect. 3.2.

In the text

	Fig. 3. Corner plot showing the posterior distributions from fitting the RVM to the position angle variation observed by the MeerKAT telescope.
In the text

	Fig. 4. Phase evolution of the total profile pre- (middle panel) and post (lower panel) change. The upper panel shows the difference between both.
In the text

	Fig. 5. Difference between the “pre-change” and “post-change” profiles on a per-channel basis. The profiles were created by integrating the observations of the respective time span in time and splitting them into eight frequency sub-bands. Before subtracting, both profiles were aligned using the ProfileShiftFit subroutine from the PYTHON interface of PSRCHIVE.
In the text

Fig. 6. Post-fit timing residuals of PSR J1618−3921. Upper two plots: Post-fit timing residuals as a function of time for the best-fit timing solution (TEMPONEST fit without removing the 100 coefficient DM noise model) of PSR J1618−3921 given in Table 3. Lower plot: Post-fit timing residuals after subtracting the DM noise model. The colours denote the different backends and systems as listed in Table 2. In all plots, the uncertainties are re-scaled with the White Noise parameters EFAC and EQUAD (cf. Table 2). The middle plot also contains the time domain realisation of the 100 parameter DM noise model: the blue lines show the 100 randomly created model realisations, and the black dots indicate the median across all these at each ToA.

In the text

Fig. 7. Corner plot for the relevant subset of timing parameters from the DDH model for J1618−3921 derived by applying the DDH binary models to the ToA data set applying non-linear Bayesian timing techniques using the software TEMPONEST to the ToAs set shown in Fig. 6. The diagonal elements show the 1D marginalised posterior distributions for each parameter, the shaded region indicates the 1σ credibility interval. The 2D contours populating the off-diagonal elements show the correlation between pairs of parameters, where the lines mark the 39%, 86% and 98% credibility regions, going from dark to light shaded.

In the text

	Fig. 8. Excerpt from the DECaPS2 survey (Schlafly et al. 2018) visualisation taken from Aladin Lite (Baumann et al. 2022) in a field of view of ∼17 as around the position of J1618−3921. The pulsar’s position is indicated by the purple hair-cross. The nearest sources are enumerated with the numbers 1–3.
In the text

Fig. 9. Demonstration of the best-fit timing solution for J1618-3921. Top: Timing residuals as a function of time for the best-fit timing model with (black) and without (red) considering the second derivative of the rotational frequency $\ddot{f}$. Bottom: Timing residuals as a function of orbital phase for the best-fit timing model with (black) and without (red) considering the derivative of the orbital period Ṗb.

In the text

Fig. 10. Constraints on the orbital inclination and masses of the J1618−3921 binary. The dashed red lines correspond to constraints derived from the TEMPONEST fit for the rate of advance of periastron $\dot{\omega }$ presented in Table 3 assuming the validity of GR and regarding the effect as purely relativistic. The black lines include 68.3 and 95.4% of all probability of the 2-dimensional probability distribution functions (pdfs) derived from the χ2 map. The left of the two inner plots shows the Mc-cosι diagram, that was sampled evenly. The grey area is excluded because MPSR must be positive. The plot to its right shows the projection of the Mc-cosι pdf on the MPSR-Mc space using the mass function. The grey area is excluded because sinι ≤ 1. The outer three plots display the projected distributions for cosι, MPSR, and Mc. the hatched area corresponds to the 1σ intervals. The unconstrained inclination angle shows that we have a non-detection of the Shapiro delay; thus, we do not provide the confidence intervals.

In the text

	Fig. 11. LOS acceleration aLOS as a function of the projection angle α onto the LOS for the three objects in the DECaPS2 (Schlafly et al. 2018) catalogue found closest to the position of J1618−3921. The acceleration was calculated from Newton’s first law assuming they are K- or M-stars with a mass of 0.6 M⊙.
In the text