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 I. Cognard L. Guillemot M. E. Lower A. Possenti G. Theureau
 


Mass measurements and orbital dynamics of PSR J1618-3921


Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany
Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans, CNRS, F-45071 Orléans, France 
Nançay Radio Astronomy Observatory, Observatoire de Paris, Université PSL, CNRS, Université d’Orléans, 18330 Nançay, France, 
Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris Cité, CNRS, F-92190 Meudon, France
INAF-Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius, Italy
Australia Telescope National Facility, CSIRO, Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia
South African Radio Astronomy Observatory, 2 Fir Street, Black River Park, Observatory 7925, South Africa
Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, VIC 3122, Australia
ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Swinburne University of Technology, PO Box 218, VIC 3122, Australia


   


 
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 observation 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. 
   We aim to measure the masses of the pulsar and its companion, as well as constraining the orbital configuration of PSR J1618-3921. This facilitates understanding similarities among eMSPs and could offer hints on their formation mechanism. 
   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 perform a full analysis on this joint dataset with a timing baseline of 23 years. We also use 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.
   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 ω = 0.00145(10).
   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^+2.07_-1.47 and ς =  0.68^+0.13_-0.09. Under the assumption of the validity of General Relativity (GR), the self-consistent combination of these three parameters lead to mass estimates of the total and individual masses in the binary of M_tot= 1.42^+0.20_-0.19, M_c = 0.20^+0.11_-0.03, and M_p = 1.20^+0.19_-0.20. We detect an unexpected change in the orbital period of Ṗ_ b =-2.26^+0.35_-0.33 ×  10^-12, which 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 Ṗ_ b. We also detect a significant second derivative of the spin frequency, f̈. The RVM fit revealed a viewing angle of ζ = 111(1). Furthermore, we report an unexpected, abrupt change of the mean pulse profile in June 2021 with unknown origin.
   We propose that the anomalous P_b and f̈ we measure for J1618-3921 indicate an additional varying acceleration due to a nearby mass, i.e., 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 P_b-M_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 Ṗ_ b and f̈. This will help us further understand the nature of this system, and perhaps improve our understanding of eMSPs in general.
   
   Triple trouble with PSR J1618-3921: Mass measurements and orbital dynamics of an eccentric millisecond pulsar
    K. Grunthal 1kgrunthal@mpifr-bonn.mpg.de
          V. Venkatraman Krishnan 1
          P. C. C. Freire 1
          M. Kramer 1
          M. Bailes 8,9
          S. Buchner 7
          M. Burgay 5
          A. D. Cameron 8,9
          C.-H.R. Chen 1
          I. Cognard 2,3
          L. Guillemot 2,3
          M. E. Lower 6
          A. Possenti 5
          G. Theureau 2,3,4
          
    September 9, 2024
=======================================================================================================================================================================================================================================================================================================================================================================================================








§ INTRODUCTION



As so-called lighthouses in the sky, pulsars are a peerless species of astronomical objects. These highly magnetised neutron stars emit a beam of an electromagnetic radiation along their magnetic poles, which is visible as a steady train of pulses at a radio telescope as the beam periodically sweeps across the observer's line-of-sight. Due to the high accuracy of atomic reference clocks and low-noise receivers in modern radio telescopes, the times-of-arrival (ToAs) of the pulses at the telescope's location are precisely recorded. The motion of the pulsar, the radio emission propagation through the interstellar medium (ISM) as well as 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 which 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. Sec. <ref>), allows uniquely precise measurements of the spin, astrometric and orbital parameters because these pulsars exhibit a uniquely stable rotational behaviour <cit.>. 

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) of the L-band receiver, its wide spectral coverage from (856 to 1712, thus a bandwidth of 856 MHz) and the high aperture efficiency of its 64×13.5 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 <cit.>. Furthermore, the ultrawide-low band receiver (UWL) of the Murriyang radio telescope have also significantly increased its spectral coverage and sensitivity.

This work was conducted as part of the “RelBin” project <cit.>, which is one of the core sub-projects of the MeerTime project, a five-year Large Survey Project <cit.>, aiming to use the precision of the MeerKAT telescope to explore fundamental physics via pulsar timing. As outlined in <cit.>, 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 <cit.>, but also improves pulsar population studies by yielding a continuously growing catalogue of precise NS mass measurements and constraining binary evolution theories <cit.>. 


The known pulsar population can be split into two large sub-groups based on their rotational behaviour and spin evolution. The so-called millisecond pulsars exhibit a rotational period of less than 30, as well as a relatively low inferred magnetic field strength (∼ 10^8-9). Additionally, about 80% of MSPs are found in binary systems, with main sequence (MS) stars, other neutron stars (NS) or white dwarfs (WD) as their companions, among the latter the Helium white dwarfs (HeWDs) are the most numerous.
 
In the current binary evolution models, these systems originate from a stellar binary, in which the more massive star already evolved into a 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 𝒪(), 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 a MSP <cit.>. At the end of this stage, the binary consists of a MSP and a stripped stellar core, which depending on its mass evolves either into a NS or a WD <cit.>. 

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, <cit.> derived a relation between the binary orbital period P_b and the mass of a HeWD companion M_WD (which we will 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. <cit.> <cit.>) and usually holds for these binaries. 

The tidal interactions accompanying the RLO lead to a circularization of the binary orbit <cit.>, as well as to an alignment of both the pulsar's spin axis with the angular momentum axis of the orbit <cit.>. 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 a 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 <cit.>.

In globular clusters, interactions with passing-by external stars can disturb the circular orbits of MSP - HeWDs, which is confirmed by the large number of eccentric binary MSPs in globular clusters[For a list of pulsars in globular clusters, see <https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html>.].
Apart from these cases, the majority of MSP - WD systems in the Galactic disk exhibit the expected small residual eccentricities <cit.>: 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, ), with 0.027 < e < 0.13 and 22 < P_b < 32d have been discovered <cit.>. These systems clearly do not follow the e-P_b relation predicted by <cit.> and became known as eccentric millisecond pulsar binaries (eMSPs). These systems are puzzling; their formation mechanism has not yet been fully understood <cit.>.



A possibility could be the formation in a triple system which became unstable, ejecting one of the components, as proposed for PSR J1903+0327 <cit.> by <cit.> and <cit.>.
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), which is seen by <cit.> and <cit.> 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 <cit.>. Lately, the timing analysis of J0955-6150 <cit.> 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 towards new formation models.


The discovery of PSR J1618-3921 (henceforth J1618-3921, similarly all other J2000 object names refer to pulsars if not indicated otherwise) was reported by <cit.> as part of a 1.4-GHz survey of the intermediate Galactic latitudes with the Parkes radio telescope. It is a recycled Galactic-disk pulsar in a binary orbit with a period of 22.7 and a low-mass companion, presumably accompanied by a low-mass HeWD. With a rotational period of 11.98, 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 <cit.>. It is now thought to belong to the eMSP class <cit.>; 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, <cit.> 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 Tab. 3 of <cit.>; for completeness also shown in the second column Tab. <ref>. This shows that the pulsar is a 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 program, 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 will start with a brief summary of the observations of J1618-3921 in Section <ref>. Section <ref> will cover the profile analysis; Section <ref> contains the timing analyses, where we report our new timing solution, including constraints of additional parameters compared to these reported by <cit.>, which include the constraints on the mass of the system. This will be followed in Section <ref> by a thorough discussion of the current state of knowledge on eMSPs in Section <ref>, 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 <ref>.








§ OBSERVATIONS AND DATA PROCESSING







 §.§ Parkes

The first observations of J1618-3921 at the Parkes Radio telescope date back to the 1999 project P309 <cit.>, 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 "multi-beam" receiver <cit.>, with a central frequency
of 1374 and a bandwidth of 288. 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 bands, with central frequencies at 1341 and 1405 respectively. Further technical details of these observations are described in <cit.>.

Making use of the ultra-wide band receiver together with the Medusa-backend <cit.>, 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 centred around a frequency of 2368. When used in pulsar folding mode, the data have a typical sub-integration length of 30 with a resolution of 128 channels per each of the 26 sub-bands, i.e. each channel has a bandwidth of 1, 1024 phase bins and full polarisation information <cit.>.



 §.§ Nançay


As pointed out in <cit.>, 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 to 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 divided in 1024 channels with a 64 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 of the mean ToA uncertainty. The observation lengths vary between 1500 and 3400, with  sub-integrations that vary between 15 and 30. In the other axes, all data files have the same resolution of 128 frequency bins, 2048 phase bins and full polarisation information.



 §.§ MeerKAT


As part the RelBin programme <cit.> 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 and an effective bandwidth of 776) together with the PTUSE backend. All technical set-up details can be found in <cit.>, <cit.> give a thorough description of the polarisation and flux calibration. The typical sub-integration length is 8, 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, which 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 Sec. <ref>). However, the Parkes observations do reveal the structure of the pulse profile at higher frequencies. A summary of all observations is presented in Table <ref>.










 §.§ Data processing


Following standard data reduction procedures in pulsar timing, we used the psrchive <cit.> 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[publicly available under <https://github.com/vivekvenkris/psrpype>] for the data reduction of the UWL observations, that have observing lengths between 2048 and 14402 seconds. psrpype uses the clfd software package[publicly available under <https://github.com/v-morello/clfd>] <cit.> 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) <cit.> was performed on the observations, using off-target calibration observations with injected pulses from a noise diode. The calibrated and cleaned UWL-data was folded into 13 frequency sub-bands.

By default, all pulsar fold-mode observations conducted with MeerKAT as part of the RelBin program are put through the meerpipe pipeline, which performs the RFI excision and polarisation calibration. meerpipe is a modified version of coastguard <cit.>. 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 <cit.>. The cleaned and calibrated files are then decimated in time, frequency and polarisation to the desired resolution, which 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.

The NRT data archives went through the full data reduction scheme described in <cit.>. 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.






§ RADIO EMISSION PROPERTIES





 §.§ Change of 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. <ref> 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, which 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. <ref>, 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. <cit.> found a spectral index of -2.28(4). 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. 






 §.§ Polarisation properties





Fig. <ref> shows the polarisation profile of J1618-3921 as recorded with the MeerKAT L-band receiver and corrected for the Rotation Measure given in <cit.>, 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 <cit.>, a phenomena that is either intrinsic to the emission of the pulsar <cit.>, or result from propagation effects in the pulsar magnetosphere <cit.>.
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 <cit.>.

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 Rotating Vector Model (RVM) <cit.>: The emitted electromagnetic waves are polarised along the magnetic field lines, which point radially outwards 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 <cit.>.
Exploiting basic geometric considerations, the RVM yields, for the position angle ψ: 

    tan(ψ-ψ_0) = sinαsin(ϕ-ϕ_0)/sin(α+β)cosα - cos(α+β)sinαcos(ϕ-ϕ_0)

where α is inclination angle of the magnetic axis relative to the spin axis and ζ is the angle between the line of sight and the spin axis of the pulsar. This is connected to β (the minimum distance between the magnetic axis and the line of sight) via ζ = α+β <cit.>. This minimum distance happens at spin phase ϕ_0; this is where ψ has the steepest slope, the corresponding PA of the linear polarisation is ψ=ψ_0. The angles in Eq. <ref> are defined as in <cit.>, i.e. “RVM/DT92” convention <cit.>. With the polarisation angle measurements from the MeerKAT observations (all data points in Fig. <ref>), we determine the RVM parameter posteriors in their joint parameter space following the method outlined in <cit.>. 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. <cit.>, the results from the best fit model are shown in terms of corner plots in Fig. <ref>. Following <cit.>, the results are presented using the RVM/DT92 convention.

We obtain α=62.27^+0.26_-0.25 and ζ = 110.63^+1.02_-0.93, quoting the 68% confidence levels on the posteriors.






 §.§ Change of the 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 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 of the mean pulse profile to be the reason for the jump in the residuals. In Fig. <ref>, 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 will 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. <ref> contains their difference ("residual profile"), calculated by matching the pre- and post-change profile with the  subroutine from the python interface of psrchive <cit.> and subtracting the re-scaled version of the latter one from the first one[This numerical output and graphical display is similar to running the  command.]. The underlying method of alignment is a χ^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 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 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 <cit.>, which was originally associated with a DM-change. A characteristic for a DM-related profile change is a f^-2 frequency dependence, i.e. this effect should dominate in the lower frequency bands. In contrast to that, the frequency dependence of the profile change of J1643-1224 <cit.> excluded a DM-origin. Here, <cit.> 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. <ref>. 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. <ref> 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 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 in the reverse engineering scheme laid out above, did not reproduce the profile change we show in Fig. <ref>. 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 towards 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 will treat these ToAs jointly.







§ TIMING ANALYSIS





 §.§ 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 Tab. <ref>. For frequency-resolved ToAs we created a frequency resolved standard profile by iteratively running paas on the integrated profile in each frequency channels. The ToAs were obtained via the pat command. The significant decrease of intensity in the higher frequency channels for the MeerKAT and Parkes observations result in large ToA uncertainties in these bands. For the timing analysis, we carefully discarded these ToAs in order to reduce to computational load of the analysis without altering the fit results. At most MJDs we are still left with a frequency resolution of up to 9 (7) channels for the Parkes (MeerKAT) data, which is a large improvement to the previous work <cit.>. Due to the low S/N of the earlier data from Parkes and the NRT, those ToAs were generated using the fully integrated observations, i.e. one frequency channel per observation.



 §.§ Fitting timing models


To analyse the final data set containing 1535 ToAs we use the timing software package tempo2 <cit.>, which does a least-square minimisation of the residuals based on the χ^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, which (among other features) enables the fit for stochastic noise processes such as red timing noise and changes in dispersion measure using power law based models <cit.>.

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

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

For the binary orbit, tempo2 provides several models based on the calculations by <cit.> which provided a standard orbital model (henceforth "DD"). In this model, the orbital motion is parameterised by five Keplerian parameters (binary orbital period P_b, longitude of periastron ω, time of periastron passage T_0, orbital eccentricity e and orbital semi-major axis projected along the line-of-sight x) and a few additional "post-Keplerian" parameters that quantify, in a theory-independent way, the relativistic deviations from the Keplerian orbital motion. Relevant here are: the rate of change of the orbital period Ṗ_b, the rate of periastron advance ω̇, the Einstein delay γ (which quantifies the effects of the the varying gravitational redshift and special-relativistic time dilation) as well as the Shapiro delay, which 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" parameters s. In GR, these are related to the companion mass M_c and the sine of the orbital inclination angle ι respectively <cit.>.


Upon deriving a timing solution for J1618-3921 we analysed the ToAs with the theory-independent DDH model developed by <cit.>, which differs from the DD model only in the parameterisation of the Shapiro delay: the new PK parameters (h_3 and ς) are less correlated than r and s, especially for systems with small orbital inclinations like J1618-3921.
In addition, in a later stage of the analysis, we used the "DDGR" model, which unlike the "DD" and "DDH" 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, which 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.




 §.§ 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. This plugin relies on Bayesian parameter estimation, which (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 <cit.>. Using the parameters from the tempo2 output ephemeris as the input for temponest, we derived a timing solution which additionally accounted for the commonly known noise parameters: Unrecognised systematics in the ToA uncertainties are modelled by the white noise parameters EFAC F and EQUAD Q on a per-backend basis. Therefore the uncertainty σ_ToA, old of each ToA is re-scaled as σ_ToA,new = √(Q^2 + F^2σ_ToA, old) <cit.>. For the chromatic models, we obtained an amplitude A and a spectral index γ <cit.>.

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 <cit.> for different combinations of noise models (Red noise (RN) only, DM noise only, Red and DM noise). On the other hand we also varied the number of noise model coefficients between 45, 60 and 100, and compared the resulting time-domain realisation between the different models. The realisations were produced using the methods of the La Forge github repository[Freely accessible via <https://github.com/nanograv/la_forge>] <cit.> adapted for the relevant models at hand. The most favoured models were the 60 and 100 coefficient DM-only models, with a difference in the log-evidences of 29. From comparing 100 averaged realisations of both noise models to the ToAs, we decided to chose the 100 coefficient model, as it visibly reflected the ToA changes more precise than the 60 coefficient model. The respective time-domain noise realisations are shown as the blue lines in the lower plot of Fig. <ref>. temponest accounts for the DM noise in terms of a power law model <cit.>, where for the chosen model we have an amplitude of A_DM=-10.37 and a slope of γ_DM=0.94. This slope is exceptionally shallow for a noise process whose slope is usually expected to be of the order of 2. From Fig. <ref> we can deduce that the residuals exhibit some significant small-timescale variations which might give rise to the shallow slope. Nevertheless, the time-domain noise realisation in the lower plot of Fig. <ref> shows that the noise model seems to match the visible trends in the data, hence we regard the noise model as satisfactory. 

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

The timing parameters of the best-fit solution from temponest using the DDH model are presented in the third column in Tab. <ref>. 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 <ref>. Tab. <ref> also shows the corresponding parameters reported by <cit.>, with blank entries when the parameter was fit for for the first time in the scope of this work. In Fig. <ref> 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 will 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.












 §.§ 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) <cit.>, a five-band optical and near-infrared survey of the southern Galactic plane, using the Aladin Lite web interface[<https://aladin.cds.unistra.fr/>] <cit.>. The corresponding excerpt from the survey image with a field of view of about 17 around the pulsar's position is shown in Fig. <ref>. At the position of the pulsar (indicated by the purple hair-cross on the image), we cannot identify any counterpart for either the pulsar or its companion. This implies that the electromagnetic emission of both bodies is below the detection thresholds of this survey, which are quoted to 23.7, 22.7, 22.2, 21.7, and 20.9 mag in the grizY bands <cit.>.



We are able to measure both the proper motion in Right Ascension μ_α=1.24^+0.14_-0.13 and Declination μ_δ=-2.37(35). This leads to a total proper motion of -2.5(3).
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 . For the lower boundary to the distance window we apply the NE2001 model <cit.>, the upper boundary is based on the YMW16 model <cit.>.  Using the distance from the NE2001 model, we translate the measured proper motions into the heliocentric velocity of the binary system of v_T = 33(4). 



  §.§.§ Spin-down and higher frequency derivatives



An important quantity describing a pulsar's properties is the intrinsic spin down Ṗ_int. For a pulsar at a distance d moving with a relative proper motion μ, any time-related measurement is influenced by the change in the Doppler shift. Thus we correct the precisely measured period derivative Ṗ = -5.37620(68)e-20 to

    Ṗ_int/P = Ṗ_obs/P + Ḋ/D,

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

    Ḋ/D = - 1/c[K⃗_0·(a⃗_PSR - a⃗_SSB) + V_T^2/d] = - a/c - μ^2d/c,

where the first term holds the contribution of the line-of-sight acceleration a by projecting the difference between the Galactic acceleration at the position of the pulsar a⃗_PSR and the solar system barycenter (SSB) a⃗_SSB onto the unit vector K⃗_0 pointing from the Earth to the pulsar. The second term, which depends on the transverse velocity V_T and the distance to the pulsar d, is the Shklovskii term <cit.>.

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 <cit.>. For the position and velocity of the solar system barycenter we assumed R_⊙ = 8.275±0.034 and V_⊙ = 240.5±4.1<cit.>.
 
Using the NE2001 distance estimate, the Shklovskii effect contributes Pμ^2d/c=6.2e-22. The Galactic acceleration field partly compensates this effect with an excess period change of Pa/c=-1.4e-22. We therefore arrive at an intrinsic spin-down of Ṗ_int=5.33326e-20, which is only slightly smaller than Ṗ_obs (Ṗ_int = 0.991 Ṗ_obs). 

Furthermore, with f= -1.0(2)e-27 we find a non-zero value of the second derivative of the spin frequency. This value is multiple orders of magnitude larger than what is expected from a pure spin-down (𝒪(10^-33), assuming a characteristic age of 10 and a braking index of 3) and among the very few values of f measured for the 333 pulsars with P < 30. Outside of globular clusters, only 9 measurements of f̈ have been made <cit.>, 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 <cit.>, in another case, J1903+0327, the system is thought to have formed in a triple system that later became unstable <cit.>; 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. <ref>. 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 f on the nature of the system and other timing parameters will be discussed in more detail in Sec. <ref>.




 §.§ Post-Keplerian parameters




  §.§.§ 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: ω̇ =0.00142^+0.00008_-0.00010.
If this effect is purely relativistic, it yields a direct measurement of the total mass of the system, M_ tot.

In order to gauge the reliability and meaning of the measurement of ω̇, we have to consider the possibility of additional non-relativistic effects. The most important of these is a proper motion contribution ω̇_μ. This contribution is given by <cit.>

    ω̇_μ = μ/sinιcos(Θ_μ - Ω),

where Θ_μ is the proper motion position angle and Ω the position angle of the line of nodes. Assuming an optimal alignment (cos(Θ_μ - Ω)=1), it contributes at the order of ω̇_μ∼8e-7.

As discussed in <cit.>, a third body in the system can add a contribution to the observed periastron advance: 

    ω̇_triple = (ẋ/x)_triple2[sin^2θ_3 (5cos^2Φ_3 - 1) -1]/ιsin2θ_3cos(ω+Φ_3).

Including x in the timing model fit yields x=(2±8)×10^-15, which is consistent with a non-detection. Considering that the geometric terms in Eq. <ref> contribute at 𝒪(1), the fit value of x gives an upper limit to the contribution of the periastron advance from the putative third body of ω̇_triple < 3e-7.

Compared to the measured rate of advance of periastron, both contributions are negligibly small, so we conclude that the measured value of ω̇ is within measurement precision, relativistic. The relativistic ω̇ relates to the total mass of the system as

    M_tot = 1/T_⊙[ω̇/3(1-e^2)]^3/2(P_b/2π)^5/2,

where T_⊙≡ GM_⊙^N/c^3= 4.9254909476412675… μ s is an exact quantity that follows from the exact definitions of the speed of light c and the solar mass parameter GM_⊙^N<cit.>. From the best-fit parameters, we derive a total mass estimate of 1.42^+0.20_-0.19. Comparing this result with the mass measurements for similar NS-WD binaries [e.g. those listed under <https://www3.mpifr-bonn.mpg.de/staff/pfreire/NS_masses.html>], we find that our measurement lies well within the expected mass range.




  §.§.§ 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 towards a low-significance detection of the Shapiro delay. From the temponest analysis we found h_3 = 2.70^+2.07_-1.47 and ς = 0.68^+0.13_-0.09. In order to convert these measurements and the measurement of ω̇ into constraints on the mass and the inclination angle of the system, we perform a  χ^2-grid analysis of the M_PSR-cosι space (cf. Sec. <ref>). The unconstrained inclination angle in the right plot of Fig. <ref> resulting from the analysis demonstrates that we did not arrive at a significant measurement of the Shapiro delay.




  §.§.§ Mass measurement



We now estimate the masses with the highly significant detection of ω̇ and the weak Shapiro delay constraints using the analysis technique outlined in <cit.>. 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, which 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 χ^2 value of the tempo2 fit, where a lower χ^2 value describes a better fit. The result is a map of χ^2 values across the M_PSR-M_c-grid, which 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. <ref>, together with the credibility band from the rate of advance of periastron. With this method we constrain the companion mass to 0.20^+0.11_-0.03, the pulsar mass to 1.20^+0.19_-0.20 and the total mass to 1.42^+0.20_-0.19 (68.3 % confidence limits).





  §.§.§ Change of orbital period


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

    ( Ṗ_b/P_b)^obs =  ( Ṗ_b/P_b)^GW + ( Ṗ_b/P_b)^ṁ
        -Ḋ/D,

where apart from the kinematic contributions (Ḋ/D), also emission of gravitational waves (GW) and mass loss from the system (ṁ) might significantly contribute to the measured value. Evaluating the expressions given in <cit.> for the latter two effects, we find (Ṗ_b/P_b)^GW∼-1e-23 and (Ṗ_b/P_b)^ṁ∼4e-28. Compared to our measured value, these contributions are negligible.

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



  §.§.§ Other parameters



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







§ DISCUSSION





 §.§ Comparison to Octau et al.


In comparison to the work by <cit.>, 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.



 §.§ Orbital geometry


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

Taking both these aspects into account, we first of all find with the reference angle from the RVM fit of 180-ζ = 69.37^+1.02_-0.93, that sinι translates into ι=66(14). This is also confirmed by performing two further RVM fits in which we restricted the variation of ζ to one of the ranges allowed by the timing results (cf. Sec. <ref>) on sinι respectively.

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




 §.§ What causes the anomalous Ṗ_b and f̈?



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

    Ṗ_ int = Ṗ_ obs - P ( Ṗ_ b/P_ b)_ obs,

since Ṗ_ b, obs is negative, this has the effect of increasing our estimate of Ṗ_ int to ∼ 1.8(3) × 10^-19, which is ∼3.4 times larger than the observed Ṗ. From this value, we estimate the characteristic age, the spin-down luminosity as well as the surface magnetic field of the pulsar <cit.>. These values can be seen in Table <ref>, there we see how the change in the value of  Ṗ_ int between this work and the work from <cit.> lead to significant differences in the values of τ_c and 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, 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, thus its acceleration would appear to be constant; in this case, we would not expect large variations in the line-of-sight acceleration and thus on the ḟ. Instead, we measure a large f of -1.0(2)e-27, which 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 <cit.>. Upon its discovery, it was regarded as an isolated pulsar, but the measurement of higher-order spin frequency derivatives led <cit.> to propose a companion in an extremely wide orbit (P_b > 200). 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, thus of a very similar order of magnitude. With the measurement of P_b we even have the advantage of estimating the acceleration of the inner binary system - this is not possible for J1024-0719, because that pulsar is not already in a binary system.

Keeping in mind that most stars are part of multiple systems, it is no surprise that on rare occasions, binaries with a pulsar are actually part of a higher-order stellar system. Due to stability arguments <cit.>, most of these systems are hierarchical triple systems, i.e. they consist of an inner binary, which is in a wider orbit around a third object.

An example is the well-known triple system consisting of the MSP J0337+1715 <cit.>. Detailed timing of this system <cit.> revealed that both orbits of the system are co-planar and circular and the WD masses are as predicted by TS99 relation, as expected from adopting the previously discussed WD-MSP formation scenario <cit.>. On the other hand, <cit.> showed in a broad study on triple systems, that the unique dynamic in these systems also allows for a stable eccentric inner binary. They also point out, that mechanisms such as Lidov-Kozai cycles prevent a synchronisation and circularization of the binary, leading to MSP systems that stand in complete contrast to the formation scenario described by <cit.>.

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 <cit.>. 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 <cit.>[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 mentioned before, in the case of PSR J1618-3921, we have direct access to the acceleration of the system via Ṗ_b, which means that the equations of <cit.> can indeed be used.]. With the knowledge of fewer derivatives, we can only put a few constraints on the orbit <cit.>: P_b relates to the corresponding acceleration from the third body a as a/c∼P_b/P_b. Similarly f relates to the change of the acceleration as a/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∼ a/a∼300, given the values form our best-fitting timing solution. This is not unexpected, and also highly in line with the findings from <cit.> in the case of J1024-0719.



  §.§.§ Optical counterpart


We consulted the DECaPS2 <cit.> 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 mentioned in  Section <ref>, 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 <cit.> to obtain the grizY magnitudes in the ABmag system to facilitate comparisons with the DECaPS2 catalogue.
Applying an extinction A_V∼0.2 mag [estimated via Galactic Dust Reddening and Extinction <https://irsa.ipac.caltech.edu/applications/DUST/>.] and adopting a distance of 5.5 kpc, a 0.56M_⊙ dwarf star (∼M0V,)
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  in the G-band, which at the distance to the pulsar of 5.5 corresponds to an absolute magnitude  >9.79. This could be a M-dwarf of mass <   0.56  M_⊙ or a compact object.



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


We consulted the Gaia DR3 <cit.>[<https://gea.esac.esa.int/archive/>] to search for objects that might have a proper motion similar to that of PSR J1618-3921, i.e., within the ± 3σ error ellipse.
This was, incidentally, how the distant binary companion of PSR J1024-0719 (a K7V star) was identified <cit.>.
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.

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

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

All these objects cause accelerations which are roughly two orders of magnitude smaller than the acceleration a_LoS, P_b=-3.45e-09□ obtained from 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 <cit.> limit: as mentioned above, it could be an M-dwarf or a compact object.









§ SUMMARY



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

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

While in previous publications <cit.>, 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 ς. 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 to 3.5) 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 of the orbital period and the large second derivative of the spin frequency, which 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 x (which will constrain the orbital orientation of the system) or higher derivatives of f, which will provide additional information on the companion mass and its orbit.







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, which 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 <cit.>, tempo2 <cit.> and temponest <cit.>, 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.



aa


