# Triple trouble with PSR J1618-3921: Mass measurements and orbital dynamics of an eccentric millisecond pulsar

K Grunthal, V Venkatraman Krishnan, P.C.C Freire, M Kramer, M Bailes, S Buchner, M Burgay, A.D Cameron, C.-H.R Chen, Ismaël Cognard, et al.

► To cite this version:

K Grunthal, V Venkatraman Krishnan, P.C.C Freire, M Kramer, M Bailes, et al. Triple trouble with PSR J1618-3921: Mass measurements and orbital dynamics of an eccentric millisecond pulsar. Astronomy & Astrophysics - A&A, 2024, 691, pp.A22. 10.1051/0004-6361/202347482. hal-04710212

HAL Id: hal-04710212

https://hal.science/hal-04710212v1

Submitted on 13 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Triple trouble with PSR J1618-3921: Mass measurements and orbital dynamics of an eccentric millisecond pulsar

K. Grunthal<fnref n="1" />,*, V. Venkatraman Krishnan<fnref n="1" />, P. C. C. Freire<fnref n="1" />, M. Kramer<fnref n="1" />, M. Bailes<fnref n="8" />,<fnref n="9" />, S. Buchner<fnref n="7" />, M. Burgay<fnref n="5" />, A. D. Cameron<fnref n="8" />,<fnref n="9" />, C.-H.R. Chen<fnref n="1" />, I. Cognard<fnref n="2" />,<fnref n="3" />, L. Guillemot<fnref n="2" />,<fnref n="3" />, M. E. Lower<fnref n="6" />, A. Possenti<fnref n="5" />, and G. Theureau<fnref n="2" />,<fnref n="3" />,<fnref n="4" />

Received 17 July 2023 / Accepted 21 August 2024

## ABSTRACT

*Context.* PSR J1618-3921 is one of five known millisecond pulsars (MSPs) in eccentric orbits (eMPSs) located in the Galactic plane, whose formation is poorly understood. Earlier studies of these objects revealed significant discrepancies between observations and predictions from standard binary evolution scenarios of pulsar-helium white dwarf (HeWD) binaries, especially in the case of PSR J0955-6150, for which mass measurements ruled out most eMSP formation models.

*Aims.* We aim to measure the masses of the pulsar and its companion, and constrain the orbital configuration of PSR J1618-3921. This facilitates understanding similarities among eMSPs and could offer hints as to their formation mechanism.

*Methods.* We conducted observations with the L-band receiver of the MeerKAT radio telescope and the UWL receiver of the Parkes Murriyang radio telescope between 2019 and 2021. These data were added to archival Parkes and Nançay observations. We performed a full analysis on this joint data set with a timing baseline of 23 years. We also used the data from recent observations to give a brief account of the emission properties of J1618-3921, including a rotating vector model (RVM) fit of the linear polarisation position angle of the pulsar.

*Results.* From the timing analysis, we measure a small but significant proper motion of the pulsar. The long timing baseline allowed for a highly significant measurement of the rate of advance of periastron of $\dot{\omega} = (0.00145 \pm 0.00010)^\circ \text{yr}^{-1}$. Despite the tenfold improvement in timing precision from MeerKAT observations, we can only report a low-significance detection of the orthometric Shapiro delay parameters, $h_3 = 2.70_{-1.47}^{+2.07} \mu\text{s}$ and $\zeta = 0.68_{-0.09}^{+0.13}$. Under the assumption of the validity of general relativity (GR), the self-consistent combination of these three parameters leads to mass estimates of the total and individual masses in the binary of $M_{\text{tot}} = 1.42_{-0.19}^{+0.20} \text{M}_\odot$, $M_c = 0.20_{-0.03}^{+0.11} \text{M}_\odot$, and $M_p = 1.20_{-0.20}^{+0.19} \text{M}_\odot$. We detect an unexpected change in the orbital period of $\dot{P}_b = -2.26_{-0.33}^{+0.35} \times 10^{-12}$, that is an order of magnitude larger and carries an opposite sign to what is expected from the Galactic acceleration and the Shklovskii effect, which are a priori the only non-negligible contributions expected for $\dot{P}_b$. We also detect a significant second derivative of the spin frequency, $\dot{f}$. The RVM fit reveals a viewing angle of $\zeta = (111 \pm 1)^\circ$. Furthermore, we report an unexpected, abrupt change in the mean pulse profile in June 2021 of unknown origin.

*Conclusions.* We propose that the anomalous $\dot{P}_b$ and $\dot{f}$ we measure for J1618-3921 indicate an additional varying acceleration due to a nearby mass. The J1618-3921 binary system is likely part of a hierarchical triple, but with the third component much farther away than the outer component of the MSP in a triple star system, PSR J0337+1715. This finding suggests that at least some eMSPs might have formed in triple star systems. Although the uncertainties are large, the binary companion mass is consistent with the $\dot{P}_b - M_{\text{WD}}$ relation, which has been verified for circular HeWD binaries and also for the two HeWDs in the PSR J0337+1715 system. Future regular observations with the MeerKAT telescope will, due to the further extension of the timing baseline, improve the measurement of $\dot{P}_b$ and $\dot{f}$. This will help us further understand the nature of this system, and perhaps improve our understanding of eMSPs in general.

**Key words.** binaries: general – stars: neutron – pulsars: individual: PSR J1618-3921

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

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

<footnote n="2">² Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans, CNRS, F-45071 Orléans, France</footnote>

<footnote n="3">³ Nançay Radio Astronomy Observatory, Observatoire de Paris, Université PSL, CNRS, Université d'Orléans, 18330 Nançay, France</footnote>

<footnote n="4">⁴ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris Cité, CNRS, F-92190 Meudon, France</footnote>

<footnote n="5">⁵ INAF-Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius, Italy</footnote>

<footnote n="6">⁶ Australia Telescope National Facility, CSIRO, Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia</footnote>

<footnote n="7">⁷ South African Radio Astronomy Observatory, 2 Fir Street, Black River Park, Observatory 7925, South Africa</footnote>

<footnote n="8">⁸ Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, VIC 3122, Australia</footnote>

<footnote n="9">⁹ ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Swinburne University of Technology, PO Box 218, VIC 3122, Australia</footnote>

<footnote>* Corresponding author; kgrunthal@mpifr-bonn.mpg.de</footnote>
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 subgroups 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 $O(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<fnref n="1" />. 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

<footnote n="1">¹ For a list of pulsars in globular clusters, see https://www3.mpi-fr-bonn.mpg.de/staff/pfreire/GCpsr.html</footnote>
<table><caption>**Table 1.** Binary parameters for all currently known Galactic-disc eccentric MSPs.</caption><thead><tr><th>Pulsar</th><th>P (ms)</th><th>P<sub>b</sub> (d)</th><th>e</th><th>M<sub>PSR</sub> (M<sub>☉</sub>)</th><th>M<sub>c</sub> (M<sub>☉</sub>)</th><th>M<sub>c,theo</sub> (M<sub>☉</sub>)</th><th>Presumable nature of system</th><th>Refs.</th></tr></thead><tbody><tr><td>J1903+0327</td><td>2.1499</td><td>95.1741</td><td>0.4367</td><td>1.667(21)</td><td>1.029(8)</td><td>–</td><td>Disrupted triple</td><td>i, j</td></tr><tr><td>J1618–3921</td><td>11.9873</td><td>22.7456</td><td>0.0274</td><td>1.20<sup>+0.19</sup><sub>−0.20</sub></td><td>0.20<sup>+0.11</sup><sub>−0.03</sub></td><td>0.269–0.297</td><td>Triple</td><td>e, f, here</td></tr><tr><td>J1950+2414</td><td>4.3048</td><td>22.1914</td><td>0.0798</td><td>1.496(23)</td><td>0.2795<sup>+0.0046</sup><sub>−0.0038</sub></td><td>0.268–0.296</td><td>Binary</td><td>a, b</td></tr><tr><td>J2234+0611</td><td>3.5766</td><td>32.0014</td><td>0.1293</td><td>1.353<sup>+0.014</sup><sub>−0.017</sub></td><td>0.298<sup>+0.015</sup><sub>−0.012</sub></td><td>0.281–0.310</td><td>Binary</td><td>c, d</td></tr><tr><td>J1946+3417</td><td>3.1701</td><td>27.0199</td><td>0.1345</td><td>1.827(13)</td><td>0.2654(13)</td><td>–</td><td>Binary</td><td>g, h</td></tr><tr><td>J1146–6610</td><td>3.7223</td><td>62.7712</td><td>0.0074</td><td>–</td><td>–</td><td>0.307–0.339</td><td>Binary</td><td>k</td></tr><tr><td>J0955–6150</td><td>1.9993</td><td>24.5784</td><td>0.1175</td><td>1.71(2)</td><td>0.254(2)</td><td>0.271–0.300</td><td>Binary</td><td>l</td></tr></tbody></table>

**Notes.** The first columns show the pulsar period, *P*, in milliseconds, their orbital period, *P*<sub>b</sub>, in days and their orbital eccentricity, *e*. In case mass measurements are available, the pulsar and companion mass, *M*<sub>PSR</sub> and *M*<sub>c</sub>, respectively, are given as well. For comparison, we also calculated the companion mass, *M*<sub>c,theo</sub>, theoretically expected from the *P*<sub>b</sub>*M*<sub>WD</sub> relation by Tauris & Savonije (1999). Refs: (a) Knispel et al. (2015), (b) Zhu et al. (2019), (c) Deneva et al. (2013), (d) Stovall et al. (2019), (e) Edwards & Bailes (2001), (f) Octau et al. (2018), (g) Barr et al. (2013), (h) Barr et al. (2017), (i) Champion et al. (2008), (j) Freire et al. (2011), (k) Lorimer et al. (2021), (l) Serylak et al. (2022).

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<fnref n="2" /> for the data reduction of the UWL observations, that have observing lengths between 2048 and 14 402 seconds. PSRPYPE uses the CLFD software package<fnref n="3" /> (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.

<table><caption>**Table 2.** Summary of all observations of J1618–3921 used in this work.</caption><thead><tr><th rowspan="3">Telescope<br>Receiver<br>Backend<br>start<br>finish</th><th colspan="4">Parkes</th></tr><tr><th colspan="2">Multi-beam</th><th colspan="2">UWL</th></tr><tr><th>FB 1<br>Aug. 1999<br>Sep. 2001</th><th>BIT<br>Aug. 1999<br>Oct. 2001</th><th>CPSR-2<br>Jan. 2003<br>Jul. 2005</th><th>Medusa<br>Oct. 2019<br>Apr. 2022</th></tr></thead><tbody><tr><td>T<sub>obs,tot</sub></td><td>51 min</td><td>1 h14 min</td><td>1 h49 min</td><td>28 h57 min</td></tr><tr><td>f<sub>0</sub> (MHz)</td><td>1374</td><td>1374</td><td>1374</td><td>2368</td></tr><tr><td>BW (MHz)</td><td>288</td><td>288</td><td>288</td><td>3328</td></tr><tr><td>N<sub>chn</sub></td><td>2</td><td>2</td><td>2</td><td>13</td></tr><tr><td>σ&#773;<sub>ToA</sub></td><td>93 µs</td><td>23 µs</td><td>37 µs</td><td>32 µs</td></tr><tr><td>EFAC</td><td>0.98</td><td></td><td>0.28</td><td>0.98</td></tr><tr><td>log<sub>10</sub>(EQUAD)</td><td>-9.77</td><td></td><td>-9.17<br>-4.79</td><td>-4.29</td></tr><tr><th colspan="5">Telescope</th></tr><tr><th rowspan="3">Receiver<br>Backend<br>start<br>finish</th><th colspan="3">Nançay</th><th>MeerKAT</th></tr><tr><th colspan="2">L-band</th><th colspan="2">L-band</th></tr><tr><th>BON<br>May 2009<br>Mar. 2011</th><th>NUPPI<br>Oct. 2013<br>Oct. 2014</th><th>PTUSE<br>Dec. 2014<br>Mar. 2022</th><th>Mar. 2019<br>Jun. 2022</th></tr><tr><td>T<sub>obs,tot</sub></td><td>6 h56 min</td><td>3 h26 min</td><td>48 h55 min</td><td>28 h51 min</td></tr><tr><td>f<sub>0</sub> (MHz)</td><td>1398</td><td>1484</td><td>1484</td><td>1284</td></tr><tr><td>BW (MHz)</td><td>128</td><td>512</td><td>512</td><td>776</td></tr><tr><td>N<sub>chn</sub></td><td>1</td><td>1</td><td>4</td><td>8</td></tr><tr><td>σ&#773;<sub>ToA</sub></td><td>227 µs</td><td>193 µs</td><td>30 µs</td><td>6 µs</td></tr><tr><td>EFAC</td><td>0.82</td><td>0.50</td><td>0.86</td><td>1.09</td></tr><tr><td>log<sub>10</sub>(EQUAD)</td><td>-13.82</td><td>-5.25</td><td>-4.57</td><td>-8.43</td></tr></tbody></table>

**Notes.** For each campaign, the receiver and backend are listed, as well as its data span. The data set is characterised by the total observation time $T_{\text{obs,tot}}$ and the centre frequency and bandwidth, $f_0$ and BW of the receiver. We also list the number of frequency channels, $N_{\text{chn}}$, each observed bandwidth was subdivided into, the mean ToA uncertainty, $\bar{\sigma}_{\text{ToA}}$, and the derived white noise parameters EFAC and EQUAD, that were determined with TEMPONEST. The CPSR-2 recorder independently records two frequency bands (1341 MHz and 1405 MHz) and is thus fit with two sets of EFAC and EQUAD values, one for each band.

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.

<footnote n="2">² Publicly available under https://github.com/vivekvenkris/psrpype</footnote>

<footnote n="3">³ Publicly available under https://github.com/v-morello/clfd</footnote>
# 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.

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

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, $\psi$:

$$ \tan(\psi - \psi_0) = \frac{\sin \alpha \sin(\phi - \phi_0)}{\sin(\alpha + \beta) \cos \alpha - \cos(\alpha + \beta) \sin \alpha \cos(\phi - \phi_0)} \quad (1) $$

where $\alpha$ is inclination angle of the magnetic axis relative to the spin axis and $\zeta$ is the angle between the line of sight and the spin axis of the pulsar. This is connected to $\beta$ (the minimum distance between the magnetic axis and the line of sight) via $\zeta = \alpha + \beta$ (Lorimer & Kramer 2005). This minimum distance happens at spin phase $\phi_0$; this is where $\psi$ has the steepest slope, the corresponding PA of the linear polarisation is $\psi = \psi_0$. 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.26}_{-0.25}$ and $\zeta = 110.63^{+1.02}_{-0.93}$, quoting the 68% confidence levels on the posteriors.

## 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<fnref n="4" />. The underlying method of alignment is a $\chi^2$-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.

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

<footnote n="4">⁴ This numerical output and graphical display is similar to running the pat -t -s <standard> <archive> command.</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Frequency resolved intensity profiles from observations with the MeerKAT L-band receiver (left, $T_{\text{obs}}$ ~ 28.85 h), the Parkes UWL receiver (middle, $T_{\text{obs}}$ ~ 28.95 h) and the NRT L-band receiver (right, $T_{\text{obs}}$ ~ 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.

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

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

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

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

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

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
<figure><img src="image_5.png" /></figure>

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

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

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

change is a $f^{-2}$ 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^{-2}$ 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.

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 \text{ pc cm}^{-3}$ 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 \text{ pc cm}^{-3}$ 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 itera-
tively 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 care-
fully 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) chan-
nels for the Parkes (MeerKAT) data, which is a large improve-
ment 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 per-
forms a least-square minimisation of the residuals based on
the $\chi^2$ 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 dif-
ferent telescopes, these are then converted to Coordinated Uni-
versal Time (UTC). Furthermore, UTC is converted to the main
realisation of the terrestrial time (TT), the high-precision coor-
dinate 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 account-
ing 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 pro-
vided a standard orbital model (henceforth “DD”). In this model,
the orbital motion is parameterised by five Keplerian parame-
ters (binary orbital period, $P_b$, longitude of periastron, $\omega$, time of
periastron passage, $T_0$, orbital eccentricity, $e$, and orbital semi-
major axis projected along the line-of-sight, $x$) and a few addi-
tional “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, $\dot{P}_b$, the rate of periastron advance, $\dot{\omega}$, the Einstein delay,
$\gamma$, (that quantifies the effects of the varying gravitational red-
shift 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, $i$, respectively (Damour & Taylor 1992).

Upon deriving a timing solution for J1618–3921 we anal-
ysed the ToAs with the theory-independent DDH model devel-
oped by Freire & Wex (2010), that differs from the DD model
only in the parameterisation of the Shapiro delay: the new PK
parameters ($h_3$ and $\zeta$) 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: Unrecog-
nised 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 $\sigma_{ToA,old}$ of each ToA is
re-scaled as $\sigma_{ToA,new} = \sqrt{Q^2 + F^2\sigma_{ToA,old}$ (Lentati et al. 2014).
For the chromatic models, we obtained an amplitude, *A*, and a
spectral index, $\gamma$ (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 reali-
sation between the different models. The realisations were pro-
duced using the methods of the LA FORGE github repository<fnref n="5" />
(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 com-
paring 100 averaged realisations of both noise models to the
ToAs, we decided to chose the 100 coefficient model, as it vis-
ibly reflected the ToA changes more precise than the 60 coeffi-
cient model. The respective time-domain noise realisations are
shown as the blue lines in the lower plot of Fig. 6. TEMPON-
EST 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 $\gamma_{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

<footnote n="5">⁵ Freely accessible via https://github.com/nanograv/la_forge</footnote>
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

trends in the data; hence, we regard the noise model as satisfac-
tory.

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 anal-
ysis, 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 corre-
sponding 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.

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<sup>6</sup> (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 com-
panion. 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).

We are able to measure both the proper motion in Right
Ascension $\mu_{\alpha} = 1.24^{+0.14}_{-0.13}$ mas yr$^{-1}$ and Declination $\mu_{\delta} = (-2.37 \pm 0.35)$ mas yr$^{-1}$. This leads to a total proper motion of $(-2.5 \pm 0.3)$ mas yr$^{-1}$. 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_{\text{T}} = (33 \pm 4)$ km s$^{-1}$.

<footnote n="6">⁶ https://aladin.cds.unistra.fr/</footnote>
### 4.4.1. Spin-down and higher frequency derivatives

An important quantity describing a pulsar's properties is the intrinsic spin down, $\dot{P}_{\text{int}}$. For a pulsar at a distance, $d$, moving with a relative proper motion, $\mu$, 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

$$ \frac{\dot{P}_{\text{int}}}{P} = \frac{\dot{P}_{\text{obs}}}{P} + \frac{\dot{D}}{D}, \qquad (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:

$$ \frac{\dot{D}}{D} = -\frac{1}{c} \left[ \mathbf{K}_0 \cdot (\mathbf{a}_{\text{PSR}} - \mathbf{a}_{\text{SSB}}) + \frac{V_T^2}{d} \right] = -\frac{a}{c} - \frac{\mu^2 d}{c}, \qquad (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, $\mathbf{a}_{\text{PSR}}$, and the solar system barycenter (SSB), $\mathbf{a}_{\text{SSB}}$, onto the unit vector, $\mathbf{K}_0$, 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_{\odot} = 8.275 \pm 0.034$ kpc and $V_{\odot} = 240.5 \pm 4.1$ km s$^{-1}$ (GRAVITY Collaboration 2021).

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

Furthermore, with $\ddot{f} = -1.0(2)10^{-27}$ 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 ($O(10^{-33})$, 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.

## 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.00008}_{-0.00010}$ yr$^{-1}$. If this effect is purely relativistic, it yields a direct measurement of the total mass of the system, $M_{\text{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)

$$ \dot{\omega}_{\mu} = \frac{\mu}{\sin \iota} \cos(\Theta_{\mu} - \Omega), \qquad (4) $$

where $\Theta_{\mu}$ is the proper motion position angle and $\Omega$ the position angle of the line of nodes. Assuming an optimal alignment ($\cos(\Theta_{\mu} - \Omega) = 1$), it contributes at the order of $\dot{\omega}_{\mu} \sim 8 \times 10^{-7} \text{° 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:

$$ \dot{\omega}_{\text{triple}} = \left(\frac{\dot{x}}{x}\right)_{\text{triple}} \frac{2[\sin^2 \theta_3(5 \cos^2 \Phi_3 - 1) - 1]}{\cot \iota \sin 2\theta_3 \cos(\omega + \Phi_3)}. \qquad (5) $$

Including $\dot{x}$ 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 $O(1)$, the fit value of $\dot{x}$ gives an upper limit to the contribution of the periastron advance from the putative third body of $\dot{\omega}_{\text{triple}} < 3 \times 10^{-7} \text{° 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

$$ M_{\text{tot}} = \frac{1}{T_{\odot}} \left[ \frac{\dot{\omega}}{3} (1 - e^2) \right]^{3/2} \left( \frac{P_b}{2\pi} \right)^{5/2}, \qquad (6) $$

where $T_{\odot} \equiv \mathcal{G}M_{\odot}^N/c^3 = 4.92549094764126...$ $\mu$s is an exact quantity that follows from the exact definitions of the speed of light, $c$, and the solar mass parameter $\mathcal{G}M_{\odot}^N$ (Prša et al. 2016). From the best-fit parameters, we derive a total mass estimate of $1.42^{+0.20}_{-0.19} M_{\odot}$. Comparing this result with the mass measurements for similar NS-WD binaries$^7$, 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^{+2.07}_{-1.47}$ $\mu$s and $\zeta = 0.68^{+0.13}_{-0.09}$. 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 $\chi^2$-grid analysis of the $M_{\text{PSR}} - \cos \iota$ 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.

<footnote>$^7$ For instance those listed under https://www3.mpifr-bonn.mpg.de/staff/pfreire/NS_masses.html</footnote>
<table><caption>**Table 3.** Timing parameters from Octau et al. (2018) and the TEMPONEST fit performed in this work.</caption><thead><tr><th>Parameter</th><th>Octau et al. (2018)</th><th>This work</th></tr></thead><tbody><tr><td>Right ascension, α (J2000)</td><td>16:18:18.8248(3)</td><td>16:18:18.82500(3)</td></tr><tr><td>Declination, δ (J2000)</td><td>-39:21:01.815(10)</td><td>-39:21:01.832(1)</td></tr><tr><td>Reference epoch (MJD)</td><td>56000</td><td>59000</td></tr><tr><td>Frequency, f (s<sup>-1</sup>)</td><td></td><td>83.421562665386(3)</td></tr><tr><td>Frequency derivative, ḟ (10<sup>-16</sup> s<sup>-2</sup>)</td><td></td><td>-3.7437(6)</td></tr><tr><td>Second frequency derivative, ḟ̇ (10<sup>-27</sup> s<sup>-3</sup>)</td><td></td><td>-1.0(2)</td></tr><tr><td>Dispersion measure, DM (cm<sup>-3</sup> pc)</td><td>117.965(11)</td><td>117.950<sup>+0.003</sup><sub>-0.002</sub></td></tr><tr><td>Dispersion measure derivative, DM1 (cm<sup>-3</sup> pc s<sup>-1</sup>)</td><td></td><td>-0.0062(5)</td></tr><tr><td>Second Dispersion measure derivative, DM2 (cm<sup>-3</sup> pc s<sup>-2</sup>)</td><td></td><td>-0.0008<sup>+0.0002</sup><sub>-0.0001</sub></td></tr><tr><td>Right ascension proper motion, μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td></td><td>1.24<sup>+0.14</sup><sub>-0.13</sub></td></tr><tr><td>Declination proper motion, μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td></td><td>-2.5(3)</td></tr><tr><td>Orbital period, P<sub>b</sub> (d)</td><td>22.74559403(19)</td><td>22.7455991<sup>+0.0000003</sup><sub>-0.0000004</sub></td></tr><tr><td>Orbital period derivative, Ṗ<sub>b</sub> (10<sup>-11</sup>)</td><td></td><td>-2.26<sup>+0.35</sup><sub>-0.33</sub></td></tr><tr><td>Projected semi-major axis of orbit, x (lt-s)</td><td>10.278300(5)</td><td>10.278285<sup>+0.000001</sup><sub>-0.000002</sub></td></tr><tr><td>Epoch of periastron, T<sub>0</sub> (MJD)</td><td>56012.21639(15)</td><td>59014.635117<sup>+0.000021</sup><sub>-0.000015</sub></td></tr><tr><td>Longitude of periastron, ω (°)</td><td>-6.717(3)</td><td>353.2919<sup>+0.0002</sup><sub>-0.0003</sub></td></tr><tr><td>Longitude of periastron derivative, ω̇ (°yr<sup>-1</sup>)</td><td></td><td>0.00142<sup>+0.00008</sup><sub>-0.00010</sub></td></tr><tr><td>Orbital eccentricity, e</td><td>0.0274133(10)</td><td>0.02741231(1)</td></tr><tr><td>Shapiro delay amplitude, h<sub>3</sub> (10<sup>-7</sup> s)</td><td></td><td>2.70<sup>+2.07</sup><sub>-1.47</sub></td></tr><tr><td>Orthometric ratio, ζ</td><td></td><td>0.68<sup>+0.13</sup><sub>-0.09</sub></td></tr><tr><td>Span of timing data (MJD)</td><td>54963.0-57869.1</td><td>51395.2-55553.4</td></tr><tr><td>Number of ToAs</td><td>70</td><td>1535</td></tr><tr><td>Weighted residual rms (μs)</td><td>25.3</td><td>8.11</td></tr><tr><td>Reduced χ<sup>2</sup> value</td><td>1.2</td><td>0.91</td></tr><tr><td colspan="3"><strong>Derived parameters</strong></td></tr><tr><td>Galactic longitude, l (°)</td><td>340.72</td><td>340.724887</td></tr><tr><td>Galactic latitude, b (°)</td><td>7.89</td><td>7.888043</td></tr><tr><td>DM-derived distance (NE2001), d (kpc)</td><td>2.7</td><td>2.7</td></tr><tr><td>DM-derived distance (YMW16), d (kpc)</td><td>5.5</td><td>5.5</td></tr><tr><td>Rotational period, P (ms)</td><td>11.987308585310(22)</td><td>11.98730841341(1)</td></tr><tr><td>Period derivative, Ṗ (10<sup>-20</sup>)</td><td>5.408(18)</td><td>5.3796(9)</td></tr><tr><td>Total proper motion, μ (mas yr<sup>-1</sup>)</td><td>&lt;6.0</td><td>2.8(3)</td></tr><tr><td>Heliocentric transverse velocity, v<sub>T</sub> (km s<sup>-1</sup>)</td><td></td><td>36(4)</td></tr><tr><td>Total mass, M<sub>tot</sub> (M<sub>☉</sub>)</td><td></td><td>1.42<sup>+0.20</sup><sub>-0.19</sub></td></tr><tr><td>Pulsar mass, M<sub>p</sub> (M<sub>☉</sub>)</td><td></td><td>1.20<sup>+0.19</sup><sub>-0.20</sub></td></tr><tr><td>Companion mass, M<sub>c</sub> (M<sub>☉</sub>)</td><td></td><td>0.20<sup>+0.11</sup><sub>-0.03</sub></td></tr><tr><td>Intrinsic spin period derivative, Ṗ<sub>int</sub> (10<sup>-20</sup>)</td><td></td><td>18(3)</td></tr><tr><td>Surface magnetic field, B (10<sup>9</sup> G)</td><td>0.814</td><td>1.5</td></tr><tr><td>Characteristic age, τ<sub>c</sub> (Gyr)</td><td>3.5</td><td>1.1</td></tr><tr><td>Spin-down luminosity, Ė (10<sup>33</sup> erg s<sup>-1</sup>)</td><td>1.24</td><td>4.1</td></tr></tbody></table>

**Notes.** All uncertainties are quoted to the left and right 39% confidence limits. We used the DDH model to fit for the Shapiro delay. The second column quotes the fitted and derived timing parameters from Octau et al. (2018) to the precision as given in Table 3 of their work. The numbers missing in the second column of the table have not been fit for by Octau et al. (2018). The reference epoch used for position and for period differs between the previous work and this work. In the second half of the table we present quantities derived from the fit values. Opposite to Octau et al. (2018) we measure the rotational frequency and its derivatives; hence, we quote to their period derivative in the second section. For the mass estimates see Sect. 4.5.3, for the equations to derive $B, τ_c$ and $Ė$ see Lorimer & Kramer (2005). The last three values are derived from $Ṗ_{int}$, meaning they are corrected for the kinematic effects.

### 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 $\chi^2$ value of the TEMPO2 fit, where a lower $\chi^2$ value describes a better fit. The result is a map of $\chi^2$ 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
<figure><img src="image_8.png" /><figcaption>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\sigma$ 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.</figcaption></figure>

method we constrain the companion mass to $0.20^{+0.11}_{-0.03} M_{\odot}$, the pulsar mass to $1.20^{+0.19}_{-0.20} M_{\odot}$ and the total mass to $1.42^{+0.20}_{-0.19} 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,

$$ \left(\frac{\dot{P}_b}{P_b}\right)^{\text{obs}} = \left(\frac{\dot{P}_b}{P_b}\right)^{\text{GW}} + \left(\frac{\dot{P}_b}{P_b}\right)^{\text{int}} - \frac{\dot{D}}{D}, \qquad (7) $$

where apart from the kinematic contributions ($\dot{D}/D$), also emission of gravitational waves (GW) and mass loss from the system ($\dot{m}$) might significantly contribute to the measured value. Evaluating the expressions given in Lorimer & Kramer (2005) for the latter two effects, we find $(\dot{P}_b/P_b)^{GW} \sim -1 \times 10^{-23} s^{-1}$ and $(\dot{P}_b/P_b)^{int} \sim 4 \times 10^{-28} 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}_b = -\dot{D}/DP_b \sim +0.05 \times 10^{-12}$. Surprisingly, the best-fit timing model reveals a measured orbital period change of $-2.2^{+0.35}_{-0.33} \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
<figure><img src="image_9.png" /><figcaption>Fig. 8. Excerpt from the DECAPSS 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.</figcaption></figure>

of magnitude too small to provide an explanation for the large observed value of $\dot{P}_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, $\gamma$, or any variation in the projected semi-major axis, $\dot{x}$, 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 $\dot{x}$ and $\gamma$ 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 $\iota$ coincides with the viewing angle $\zeta$. 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 \iota$. Hence, we cannot distinguish if the corresponding inclination angle is $\iota$ or $180^{\circ} - \iota$. In case of a reliable RVM fit, this ambiguity can be solved by comparing $\iota$ to $\zeta$. This can also not be done directly, since the above RVM equation assumes that $\psi$ increases clockwise on the sky, opposite to the astronomical convention, where $\psi$ 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 $\zeta$ with $180^{\circ} - \iota$ or $\iota$ with $180^{\circ} - \zeta$, 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^{\circ} - \zeta = 69.37^{+1.02^{\circ}}_{-0.93^{\circ}}$, that $\sin \iota$ translates into $\iota = (66 \pm 14)^{\circ}$. This is also confirmed by performing two further RVM fits in which we restricted the variation in $\zeta$ to one of the ranges allowed by the timing results (cf. Sect. 4) on $\sin \iota$, 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 $\dot{P}_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 $\dot{P}_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 $\dot{P}_b$, and still obtain a positive value. Subtracting Eq. (7) from Eq. (2), and neglecting the GW emission terms, we obtain:

$$ \dot{P}_{\text{int}} = \dot{P}_{\text{obs}} - P \left( \frac{\dot{P}_{\text{b}}}{\dot{P}_{\text{b}}} \right)_{\text{obs}} \quad (8) $$

since $\dot{P}_{\text{b,obs}}$ is negative, this has the effect of increasing our estimate of $\dot{P}_{\text{int}}$ to ~1.8(3) × 10⁻¹⁹, that is ~3.4 times larger than the observed $\dot{P}$. 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 $\dot{P}_{\text{int}}$ between this work and the work from Octau et al. (2018) lead to significant differences in the values of $\tau_c$ and $\dot{E}$.

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^{\circ}$, 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;
<figure><img src="image_10.png" /><figcaption>**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 $\dot{\Omega}$. 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 $P_b$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**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 $\chi^2$ map. The left of the two inner plots shows the $M_c-\cos i$ diagram, that was sampled evenly. The grey area is excluded because $M_{PSR}$ must be positive. The plot to its right shows the projection of the $M_c-\cos i$ pdf on the $M_{PSR}-M_c$ space using the mass function. The grey area is excluded because $\sin i \le 1$. The outer three plots display the projected distributions for $\cos i$, $M_{PSR}$, and $M_c$, the hatched area corresponds to the $1\sigma$ intervals. The unconstrained inclination angle shows that we have a non-detection of the Shapiro delay; thus, we do not provide the confidence intervals.</figcaption></figure>
in this case, we would not expect large variations in the line-of-sight acceleration and thus on the *f*. Instead, we measure a large *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*<sub>b</sub> > 200 yr). This was confirmed by the detection of a nearby star with the same proper motion. Comparing the measured value *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 *P*<sub>b</sub> 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 synchronization 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)<fnref n="8" />. With the knowledge of fewer derivatives, we can only put a few constraints on the orbit (Bassa et al. 2016): *P*<sub>b</sub> relates to the corresponding acceleration from the third body, *a*, as *a/c* ~ *P*<sub>b</sub>*P*<sub>b</sub>. Similarly *ḟ* relates to the change in the acceleration as *ā/c* ~ *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_{b,3} \sim a/\dot{a} \sim 300$ 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 \sim 0.2$ mag<fnref n="9" /> and adopting a distance of 5.5 kpc, a 0.56 M<sub>☉</sub> 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<sub>☉</sub> or a compact object.

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

We consulted the Gaia DR3 (Gaia Collaboration 2023)<fnref n="10" /> 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 faint 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

<footnote n="8">⁸ The first derivative of *f* generally cannot be used as intended by these authors, because of the a priori unknowable pulsar spin-down, but also because, in the system studied in these works (PSR B1620-26), the acceleration caused by the host globular cluster (M4) is also hard to estimate, given the lack of a precise 3-D position of the pulsar relative to M4. However, as was mentioned before, in the case of PSR J1618-3921, we have direct access to the acceleration of the system via *P*<sub>b</sub>, that means that the equations of Joshi & Rasio (1997) can indeed be used.</footnote>

<footnote n="9">⁹ Estimated via Galactic Dust Reddening and Extinction https://irsa.ipac.caltech.edu/applications/DUST/</footnote>

<footnote n="10">¹⁰ https://gea.esac.esa.int/archive/</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** LOS acceleration $a_{LOS}$ as a function of the projection angle $\alpha$ 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⊙.</figcaption></figure>

those for a 0.6 M⊙ (or K9V) star with $A_V \sim 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_{LOS} = \frac{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 $\alpha$ 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 \alpha$. The resulting acceleration curves calculated under the previously outlined assumptions for the three objects marked in Fig. 8 are shown in Fig. 11.

All these objects cause accelerations that are roughly two orders of magnitude smaller than the acceleration $a_{LoS, P_\odot} = -3.45 \times 10^{-9} \text{ m s}^{-2}$ obtained from $\dot{P}_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 \pm 1)^\circ$. 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_3$ and $\zeta$. 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 $\dot{x}$ (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.

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

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

Antoniadis, J., Kaplan, D. L., Stovall, K., et al. 2016, ApJ, 830, 36

Archibald, A. M., Gusinskaia, N. V., Hessels, J. W. T., et al. 2018, Nature, 559, 73

Bailes, M. 2007, ArXiv e-prints [arXiv:astro-ph/0702698]

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

Barr, E. D., Champion, D. J., Kramer, M., et al. 2013, MNRAS, 435, 2234

Barr, E. D., Freire, P. C. C., Kramer, M., et al. 2017, MNRAS, 465, 1711

Bassa, C. G., Janssen, G. H., Stappers, B. W., et al. 2016, MNRAS, 460, 2207

Baumann, M., Boch, T., Pineau, F.-X., et al. 2022, ASP Conf. Ser., 532, 7

Blandford, R. D., & Scharlemann, E. T. 1976, MNRAS, 174, 59

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Champion, D. J., Ransom, S. M., Lazarus, P., et al. 2008, Science, 320, 1309

Cordes, J., Lazio, T., Chatterjee, S., Arzoumanian, Z., & Chernoff, D. 2002, in 34th COSPAR Scientific Assembly, 34, 2305

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

Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840

Deneva, J. S., Stovall, K., McLaughlin, M. A., et al. 2013, ApJ, 775, 51

Dyks, J. 2020, MNRAS, 495, L118

Edwards, R. T., & Bailes, M. 2001, ApJ, 553, 801

Everett, J. E., & Weisberg, J. M. 2001, ApJ, 553, 341

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, Open J. Astrophys., 2, 10

Freire, P. C. C., & Tauris, T. M. 2014, MNRAS, 438, L86

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

Freire, P. C. C., Bassa, C. G., Wex, N., et al. 2011, MNRAS, 412, 2763

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Gangadhara, R. T. 1997, ArXiv e-prints [arXiv:astro-ph/9707168]

GRAVITY Collaboration (Abuter, R., et al.) 2021, A&A, 647, A59

Hazboun, J. S. 2020, https://doi.org/10.5281/zenodo.4152550

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

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

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

Hu, H., Kramer, M., Champion, D. J., et al. 2022, A&A, 667, A149

Hui, C. Y., Wu, K., Han, Q., Kong, A. K. H., & Tam, P. H. T. 2018, ApJ, 864, 30

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

Jonas, J., & MeerKAT Team 2016, in MeerKAT Science: On the Pathway to the SKA, 1

Joshi, K. J., & Rasio, F. A. 1997, ApJ, 479, 948

Knispel, B., Lyne, A. G., Stappers, B. W., et al. 2015, ApJ, 806, 140

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

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

Lazaridis, K., Wex, N., Jessner, A., et al. 2009, MNRAS, 400, 805

Lazarus, P., Karuppusamy, R., Graikou, E., et al. 2016, MNRAS, 458, 868

Lentati, L., Alexander, P., Hobson, M. P., et al. 2014, MNRAS, 437, 3004

Lin, F. X., Lin, H.-H., Luo, J., et al. 2021, MNRAS, 508, 1115

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

Lorimer, D. R., Kawash, A. M., Freire, P. C. C., et al. 2021, MNRAS, 507, 5303

Manchester, R. N. 1975, PASA, 2, 334

Manchester, R. N., Taylor, J. H., & Huguenin, G. R. 1975, ApJ, 196, 83

Manchester, R., Lyne, A., Camilo, F., et al. 2001, MNRAS, 328, 17

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993

Melrose, D. B., & Stoneham, R. J. 1977, PASA, 3, 120

Melrose, D., Miller, A., Karastergiou, A., & Luo, Q. 2006, MNRAS, 365, 638

Morello, V., Barr, E. D., Cooper, S., et al. 2018, MNRAS, 483, 3673

Nguyen, C. T., Costa, G., Girardi, L., et al. 2022, A&A, 665, A126

Octau, F., Cognard, I., Guillemot, L., et al. 2018, A&A, 612, A78

Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9

Petrova, S. A. 2001, A&A, 378, 883

Petrova, S. A. 2006, MNRAS, 366, 1539

Phinney, E. S. 1992, Phil. Trans. R. Soc. London Ser. A, 341, 39

Phinney, E. S., & Kulkarni, S. R. 1994, ARA&A, 32, 591

Portegies Zwart, S., van den Heuvel, E. P. J., van Leeuwen, J., & Nelemans, G. 2011, ApJ, 734, 55

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

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

Radhakrishnan, V., & Srinivasan, G. 1982, Curr. Sci., 51, 1096

Ransom, S. M., Stairs, I. H., Archibald, A. M., et al. 2014, Nature, 505, 520

Rasio, F. A. 1994, ApJ, 427, L107

Ridolfi, A., Freire, P. C. C., Gupta, Y., & Ransom, S. M. 2019, MNRAS, 490, 3860

Schlafly, E. F., Green, G. M., Lang, D., et al. 2018, ApJS, 234, 39

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

Serylak, M., Venkatraman Krishnan, V., Freire, P. C. C., et al. 2022, A&A, 665, A53

Shannon, R. M., Lentati, L. T., Kerr, M., et al. 2016, ApJ, 828, L1

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

Smedley, S. L., Tout, C. A., Ferrario, L., & Wickramasinghe, D. T. 2013, MNRAS, 437, 2217

Spiewak, R., Bailes, M., Miles, M. T., et al. 2022, PASA, 39, e027

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

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

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

Tauris, T. M., & van den Heuvel, E. P. J. 2023, Physics of Binary Star Evolution: From Stars to X-ray Binaries and Gravitational Wave Sources (Princeton University Press)

Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170

Toonen, S., Hamers, A., & Portegies Zwart, S. 2016, Comput. Astrophys. Cosmol., 3, 6

van Straten, W. 2013, ApJS, 204, 13

van Straten, W., Manchester, R. N., Johnston, S., & Reynolds, J. E. 2010, PASA, 27, 104

Voisin, G., Cognard, I., Freire, P. C. C., et al. 2020, A&A, 638, A24

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

Zhu, W. W., Freire, P. C. C., Knispel, B., et al. 2019, ApJ, 881, 165