© The Authors 2025

1. Introduction

Globular clusters (GCs) are old, gravitationally bound stellar systems that are roughly spherical and contain some of the oldest stars in our galaxy. The large number density of stars in their cores (∼10³–10⁶ pc⁻³, Baumgardt & Hilker 2018) result in a high probability of close stellar encounters. Some of these encounters result in exchange interactions, where a main sequence (MS) star might exchange an MS companion with an old neutron star (NS) in the cluster. As the MS star evolves, it starts transferring matter to the NS, resulting in the formation of a low-mass X-ray binary (LMXB). For this reason, GCs are known to have ∼10³ times more X-ray binaries per unit of stellar mass than the Galactic disk (Clark 1975). This abundance is a clear sign that these X-ray binaries form dynamically.

In an LMXB, the NS is ‘recycled’ and ‘spun up’ to high spin frequencies by accretion of matter from the evolving companion, leading to the formation of a millisecond pulsar (MSP) – a class of fully recycled pulsars with spin periods of a few milliseconds and very small spin period derivatives (Tauris & van den Heuvel 2023). In this process, the orbit is circularised through tidal interactions. As the companion evolves to a white dwarf (WD) and the NS to an MSP, the low eccentricity of the system is normally retained, as observed in the vast majority of MSPs in the Galactic disk.

However, there are exceptions to this rule. In multiple star systems, things can go very differently: One exceptional system observed in the Galactic disk is PSR J1903+0327: this MSP (P = 2.15 ms) has a large orbital eccentricity (e = 0.43) and a massive companion (Champion et al. 2008), which was later found to be a 1.03 M_⊙ main sequence star (Freire et al. 2011). Detailed observational (Freire et al. 2011) and theoretical (Portegies Zwart et al. 2011; Pijloo et al. 2012) studies of this system indicate that it started its life as a triple system, which later became unstable, leading to the ejection of the donor to the pulsar.

Globular clusters are also multiple star systems. There are 344 pulsars currently known in 45 GCs¹; its MSP population represents about 40% of the total known MSP population. Although many of these MSPs have circular orbits and WD companions, some GCs have a majority of isolated pulsars and slow pulsars with higher magnetic fields (e.g. Terzan 1, Singleton et al. 2024, NGC 6517, Yin et al. 2024, NGC 6522, Abbate et al. 2023, NGC 6624, Freire et al. 2011; Abbate et al. 2022, NGC 6752, Corongiu et al. 2024, M15, Wu et al. 2024; Zhou et al. 2024); additionally, some very dense GCs have MSPs with massive companions in eccentric orbits (e.g. Terzan 5, Padmanabh et al. 2024, NGC 6544, Lynch et al. 2012, NGC 6624, Ridolfi et al. 2021, NGC 6652, DeCesar et al. 2015 and M30, Balakrishnan et al. 2023).

That these unusual pulsars are almost exclusively found in the denser clusters was explained by Verbunt & Freire (2014) using the concept of the encounter rate per binary (γ_i). A GC can have a large total stellar encounter rate and form many LMXBs and binary MSPs, but if the number of encounters per star over the age of the GC is small, then an LMXB will likely evolve undisturbed into a ‘normal’ circular MSP – low-mass WD system, as generally observed in the Galactic disk and the lower-density GCs. However, if the stellar density is very high, then the probability of subsequent encounters for any particular system increases. Only then are we likely to see a recycled pulsar – the end product of an LMXB formed in an exchange interaction – go into subsequent (‘secondary’) exchange interactions. These can form exotic systems like the aforementioned MSPs with eccentric orbits and massive companions.

The GC NGC 1851 is located in the southern constellation of Columba at a distance of 11.66 ± 0.25 kpc (Libralato et al. 2022). It has a very high central density of about 3 × 10⁶ M_⊙ pc⁻³ (Baumgardt & Hilker 2018) and has a moderately large encounter rate per binary (γ_i = 12.4 γ_M4, Verbunt & Freire 2014). The effect on the pulsar population is clear: of the 14 known pulsars, eight are in binaries (Ridolfi et al. 2022); and of these, three have eccentric orbits and unusually massive companions (Ridolfi et al. 2019; Barr et al. 2024), the characteristics of secondary exchange encounter products. These systems are especially interesting targets for further analysis.

One of these eccentric binaries, PSR J0514−4002A, was found in 2003 using the Giant Metrewave Radio Telescope (GMRT). It was the first pulsar discovered in NGC 1851 (Freire et al. 2004); for this reason it is also known (and will be designated in this work) as NGC 1851A. The timing solution was obtained from data taken with the Green Bank Telescope (GBT) (Freire et al. 2007) and later extended with data from the GMRT (Ridolfi et al. 2019). This binary pulsar is unlike anything found in the Galactic disk: its fast spin (P ∼ 5 ms), high orbital eccentricity (e ∼ 0.89) and large (∼ 1 M_⊙) companion mass made it the second clear case of a secondary exchange encounter product (the first case was PSR B2127+11C in the globular cluster M15, see Prince et al. 1991).

In this work, we aim to address the unresolved questions about NGC 1851A raised by Ridolfi et al. (2019) (hereafter R19). First, their estimate of the proper motion suggested that the pulsar was not bound to the cluster. Second, the observed large positive change in the orbital period of the system could not be explained by the acceleration caused by the gravitational field of the GC, as this would also cause a large positive spin period derivative, which is not observed. Third, their mass measurements left open the question about the nature of the companion. Finally, and most importantly, they measured a large second spin frequency derivative, which is indicative of a very large change of acceleration of the system (a “jerk”), but the implications of this were not discussed in detail.

These puzzles motivated an extension of the timing dataset with sensitive observations from the MeerKAT radio telescope over a long time period, with a few dedicated observations of the system around the periastron and superior conjunction. This extended timing analysis allowed for a robust measurement of the proper motion, the orbital period derivative, and improvements of the constraints on the other relativistic parameters of the system, which in turn allowed more precise measurements of the component masses and a better measurement of the jerk and its change with time.

In this paper, we present the up-to-date timing solution of NGC 1851A, obtained with observations made using the upgraded GMRT (henceforth uGMRT) and MeerKAT, in addition to the earlier GBT and GMRT dataset from R19, over a total time span of ∼18 years. The observations, the subsequent data reduction, and how the resulting pulse times of arrival (TOAs) were analysed are discussed in Section 2. The system parameters and other derived results benefit largely from the much longer timing baseline and the inclusion of new datasets with good timing precision; these are discussed in Section 3. Apart from the timing parameters of NGC 1851A, this section addresses the unresolved issues from R19. Next comes the analysis of the large observed jerk, which we argue is caused by the presence of a third object near the NGC 1851A system. In Section 4, we discuss the implications of this jerk for understanding the variation of the orbital period and the spin-down of the pulsar. In Section 5, we explore the characteristics of this third object. In Section 6, we update the mass measurements for the NGC 1851A system. In Section 7 we briefly summarise the results and arrive at a few general conclusions about the nature of the system.

2. Observations and data analysis

The new data analysed in this work include the observations taken with the uGMRT between 2018 and 2020 and the ones later carried out using MeerKAT, spanning over a total time of more than 4 years. A full description of the dataset, including a breakdown of the individual datasets is given in Table 1.

2.1. GMRT observations

The observing campaigns for NGC 1851A with the GMRT have been running since 2017 April (MJD 57872), during which the data were taken in the Phased Array (PA) mode with a sampling time of 81.92 μs. Since September 2017, besides the PA mode, all observations are also taken simultaneously using the coherent de dispersion (CDP) mode of the GMRT Wideband Backend (GWB, De & Gupta 2016; Reddy et al. 2017). For all observations made in the CDP mode, the data were coherently de dispersed at the dispersion measure (DM) of NGC 1851A of ∼52.14 pc cm⁻³; for a detailed description of these data, see R19. This work looks into all CDP data taken for NGC 1851A between 2019 April (MJD 58601) to 2020 March (MJD 58928) with the 250–500 MHz Band-3 receiver of the uGMRT, a substantially longer dataset than had been analysed by R19 (Table 2).

2.2. MeerKAT observations

NGC 1851 has been observed by MeerKAT, as part of both the MeerTime (Bailes et al. 2020) and TRAPUM (Stappers & Kramer 2016) large science programs. Most of MeerKAT data used in this paper were taken as part of the MeerTime GC pulsar timing programme. We observed the pulsar with the MeerKAT on 27 different epochs, from January 2021 to April 2023 (Table 1), and the data were recorded in coherently dedispersed, full-Stokes mode using the Pulsar Timing User Supplied Equipment (PTUSE) backend (Bailes et al. 2020). The first seven timing observations with MeerTime, which targeted precisely NGC 1851A, were taken in parallel to the observations performed for TRAPUM. This includes 3 observations done with the L-band receiver at a central frequency of 1284 MHz and a bandwidth of 856 MHz, and the rest using the Ultra High Frequency (UHF) receivers at a central frequency of 816 MHz and a bandwidth of 544 MHz. These data were sampled every 7.53 μs and possessed 4096 channels.

Further UHF observations, carried out during the MeerKAT GC survey, were recorded with a sampling time of 9.41 μs across 1024 channels, each with a bandwidth of ∼0.53 MHz. Later, the sampling time was increased to 15.06 μs to reduce the acquired data volume. The different channelisations of the UHF data resulted in a 0.2 MHz difference in the central frequencies. All observations taken with MeerKAT ranged between 2 hours to 4 hours in duration, amounting to a total of 64 hours, and the details for all observations are presented in Table 1. A detailed description of the setups used for these observations can be found in Ridolfi et al. (2022) and Barr et al. (2024).

2.3. Data analysis

For all GMRT observations taken in the PA and the CDP mode with 200 MHz bandwidth, offsets of 1.34217728 s and 2.01326592 s² were included in the analysis (Appendix A) for the respective set of times of arrival (TOAs). Given the 4.99 ms period of this pulsar, these corrections accounted for hundreds of rotations, which would be missed if these offsets were not taken into account.

The initial step of analysis for all search-mode observations taken with the uGMRT was done using gptool³, a tool used for data reduction and bandpass normalisation of the beamformer data of the uGMRT (Susobhanan et al. 2021). The ephemeris obtained from the analysis of R19 was used to fold the pulsar at its topocentric period for every observation, after an initial radio frequency interference (rfi) removal, using the prepfold and rfifind routines of PRESTO⁴ (Ransom 2001) respectively. These folded archives were then converted into PSRFITS format, using the psrconv routine of PSRCHIVE, and the number of frequency channels in the resulting folded archives were downsampled to 256. For the two frequency bands, a high signal-to-noise ratio analytic profile (a ‘template’) was created by adding up all the archives using the psradd routine of PSRCHIVE⁵ and fitting two Gaussian components to it.

The pulse profiles for each of the detections of the pulsar were then cross-correlated with the corresponding templates, and TOAs were extracted using the pat routine of PSRCHIVE. Depending on the scintillation and the brightness of the pulsar for individual observations, TOAs were extracted every ∼10–60 minutes, and with 2 and 4 sub-bands for the GMRT and the MeerKAT datasets respectively. In addition to the offsets for the GMRT datasets described above, a time difference of 1.88 μs was also included for the TOAs created from the latter UHF archives to account for the 0.2 MHz difference in the central frequency⁶. This corresponds to the difference in the sampling time between the two sets of the UHF observations. A careful phase alignment of the standard profiles for the CDP observations of the GMRT and the L-band and UHF observations of MeerKAT was done to a phase precision of 0.01 for the removal of any arbitrary timing offset (Guo et al. 2021).

2.4. Timing procedure

In this work, the observations taken with MeerKAT and uGMRT were added to the earlier data used by R19. For the CDP data from the GMRT discussed in R19, all the observations were re-folded after performing the standard procedure of bitshift conversion (downsampling the data from a 16-bit to an 8-bit format, Gautam et al. 2022), bandpass normalisation and rfi mitigation, as is used for the rest of the uGMRT dataset as described above. This procedure enabled the use of a single template for all GMRT CDP data and ensured a better timing precision. The final addition of the latest observations taken in April 2023 with MeerKAT has extended the timing baseline of the pulsar to almost two decades.

We used the TEMPO⁷ timing package for the analysis of the TOAs and the Jet Propulsion Laboratory’s DE 421 Solar System ephemeris (Folkner et al. 2009) to take into account the motion of the radio telescopes relative to the Solar System Barycentre (SSB). The parameters in the timing solution are specified in dynamical barycentric time (TDB).

For the initial description of the orbit we used the DDH and the DDGR orbital models, which are derived from the ‘DD’ orbital model (Damour & Deruelle 1986). Like the DD model, the DDH model is a theory-independent model, where each relativistic effect is quantified by a ‘post-Keplerian’ (PK) parameter. The main difference relative to the DD model is that it uses the orthometric parametrisation of the Shapiro delay (Freire & Wex 2010); this was chosen because the two PK parameters that describe the Shapiro delay (h₃ and ς) are much less correlated than the ‘range’ and ‘shape’ parameters that describe the Shapiro delay in the DD model (especially for orbital inclinations that are far from edge-on, as is the case for NGC 1851A). These parameters can be related to the Keplerian parameters and the masses of the components of the system using a gravity theory. The DDGR model is a modification of the “DD” model that fits directly for the total mass and the companion mass of the system under the assumption that all relativistic effects are as predicted by general relativity (GR).

The significant change in the line-of-sight acceleration (the jerk) inferred from the spin period derivatives required the inclusion of higher order orbital frequency derivatives for the system (see further details in Section 4). In principle, this can be done easily using the BTX model (Shaifullah et al. 2016), which is a reimplementation of the BT model (Blandford & Teukolsky 1976) that includes higher order derivatives for the orbital frequency. However, the BTX model is not suitable for eccentric systems like NGC 1851A since it estimates that the longitude of periastron increases linearly with time. The DD-based timing models, on the other hand, correctly account for the evolution of the longitude of periastron, which increases linearly with the true anomaly. This feature significantly improves the quality of the fit for NGC 1851A. However, as currently implemented in TEMPO, these models are unable to account for higher orbital frequency derivatives. For this reason, we implemented modified versions of these models that can take into account the jerk of the system in a self-consistent way (see details in Appendix E). All updated parameters and results reported in this work were derived using these modified DD models.

3. Results

The extended timing baseline of ∼18 years, inclusion of the 2-s offset in the GMRT dataset (see Appendix A) and modified analysis procedure has been the key in obtaining precise measurements of a handful of astrometric and spin parameters for NGC 1851A. These values and their interpretations are discussed in Sect. 3. The updated timing parameters for NGC 1851A obtained after this analysis are presented in Table 2. All the parameters and their uncertainties are derived using the TEMPO timing software and the modified DDH model mentioned above.

The residuals between the predicted and measured times of arrival of the pulsar, as a function of time and orbital phase, are given in Fig. 1. There are no apparent systematic trends in the residuals, and the reduced χ² = 1.0044 of the overall fit shows that the timing model provides a good estimation of the TOAs. This was obtained after individual adjustments for each of the GBT, GMRT, and MeerKAT datasets separately, and these adjustment factors (EFACs) are noted in Table 1.

Fig. 1. Top to bottom: Post-fit timing residuals for NGC 1851A obtained with the DDH timing model, as a function of (a) epoch and (b) orbital phase. The orbital phase 0 denotes periastron and the residual 1-σ uncertainties are indicated by vertical error bars. The various colours denote the different telescopes and backends used for the timing. Orange indicates the early GBT data, and purple and red indicate the uGMRT data taken in the PA and CDP modes respectively. Grey and green represent the data taken with the L-band and the UHF receivers of the MeerKAT.

3.1. Astrometric parameters

We found that the position of the pulsar at the reference epoch is $\alpha ={05}^{\mathrm{h}}{14}^{\mathrm{m}}06\stackrel{\text{s}}{.}697099(37)$ and $\delta =-40°{02}^{\prime }48\stackrel{″}{.}90709(23)$, which is 1$\stackrel{″}{.}$58 from the centre of the GC at ${\alpha }^{\prime }={05}^{\mathrm{h}}{14}^{\mathrm{m}}06\stackrel{\text{s}}{.}755$ and ${\delta }^{\prime }=-40°{02}^{\prime }47\stackrel{″}{.}47$ (Miocchi et al. 2013). The latter study also includes an estimate of the core radius, 5$\stackrel{″}{.}$4, this means that the pulsar is well within the core of the cluster.

The quantity that has changed the most notably from the introduction of the additional 2 s offset in the GMRT data and the continued timing analysis, is the proper motion of the system. We obtained μ_α = 2.61 ± 0.13 mas yr⁻¹ in right ascension and μ_δ = −0.90 ± 0.11 mas yr⁻¹ in declination. The total proper motion of NGC 1851A is 2.76 (12) mas yr⁻¹ and the corresponding position angle is 109 (2) deg. Using the images obtained from the Hubble Space Telescope (HST) Legacy Survey of Galactic Globular Clusters and Gaia astrometric data, Libralato et al. (2022) have provided an updated measurement of the proper motion of the host cluster in right ascension and declination as μ_α = 2.128 ± 0.031 mas yr⁻¹, μ_δ = − 0.646 ± 0.032 mas yr⁻¹. This means that the motion of the pulsar relative to the cluster is Δμ_α = 4.48 ± 0.13 mas yr⁻¹, Δμ_δ = −0.25 ± 0.11 mas yr⁻¹. The corresponding relative transverse velocity of NGC 1851A to the centre of the cluster, using the distance estimate of d = 11.66 ± 0.25 kpc (Libralato et al. 2022), is v_T = 30 ± 7 km s⁻¹. While the three-dimensional velocity for any binary pulsar in a GC is not known, the transverse velocity of NGC 1851A is smaller than the cluster’s central escape velocity (v_esc) of 42.9 km s⁻¹ (Baumgardt & Hilker 2018), being therefore consistent with the system’s association to the cluster. This was already indicated by the close sky proximity to the centre of the cluster and the similarity of the DMs of all the other pulsars in the cluster (Ridolfi et al. 2022).

3.2. Keplerian orbital parameters

With an orbital eccentricity of e  ≃  0.89, NGC 1851A is the third most eccentric known binary system in a GC, after NGC 6652A (e  ≃  0.95, DeCesar et al. 2015) and Terzan 5 ap (e  ≃  0.905, Padmanabh et al. 2024). Using the values of the orbital period of the system (P_b) and the projected semi-major axis of the pulsar orbit in time units (x ≡ a_p sin i/c) from Table 2, we obtained for the mass function

$$\begin{array}{c}\hfill f({\overline{M}}_{\mathrm{p}},{\overline{M}}_{\mathrm{c}})=\frac{{({\overline{M}}_{\mathrm{c}}sini)}^3}{{({\overline{M}}_{\mathrm{p}}+{\overline{M}}_{\mathrm{c}})}^2}=\frac{4{\pi }^2}{T_{\odot }}\frac{x^3}{{P_{\mathrm{b}}}^2}=0.1454280(16),\end{array}$$(1)

where the exact quantity $T_{\odot }\equiv {(\mathcal{G}\mathcal{M})}_{\odot }^{\mathrm{N}}/c^3=4.925490947641266978\cdots \mathrm{\mu }\mathrm{s}$ is the nominal solar mass parameter ${(\mathcal{G}\mathcal{M})}_{\odot }^{\mathrm{N}}$ (Prša et al. 2016) in time units, and M_p and M_c are the masses of the pulsar and the companion respectively. In equations that contain T_⊙ explicitly, one has to use the adimensional mass parameters ${\overline{M}}_i\equiv G{M}_i/{(\mathcal{GM})}_{\odot }^{\mathrm{N}}$, (i = p, c, …), which are the numerical values of M_i when expressed in units of solar mass (M_⊙).

This high mass function value is indicative of a massive companion: assuming that the pulsar mass is larger than the smallest NS mass known (M_p  =  1.17 M_⊙, Martinez et al. 2015) and an edge-on orbit (i = 90deg), we derived M_c > 0.84 M_⊙. As discussed previously by Freire et al. (2004, 2007) and R19, these combined factors reinforce the conclusion that NGC 1851A is the product of a secondary exchange encounter: the previous low-mass companion to the pulsar was replaced by its current massive counterpart. The chaotic nature of such exchange encounters results in high orbital eccentricities. If the new companion to the pulsar is degenerate, then there is not much tidal circularisation.

3.3. Rate of advance of periastron

The updated timing model for the pulsar includes a very precise measurement of the rate of advance of periastron: ${\dot{\omega }}_{\mathrm{obs}}=0.012961(1)\mathrm{deg}{\mathrm{yr}}^{-1}$. This is fully consistent with the value presented by R19 but twice as precise and is 40 times more precise than the first value quoted by Freire et al. (2007); this improvement is due to the extension of the dataset and the usage of pulse profile templates with consistent phase definition (as described earlier). This parameter is a quantification of the average change in the longitude of periastron of the system due to various relativistic and classical effects.

Under the assumption that the ${\dot{\omega }}_{\mathrm{obs}}$ is purely relativistic, we derived the total mass of the system from this measurement:

$$\begin{array}{c}\hfill {\overline{M}}_{\mathrm{tot}}=\frac{1}{T_{\odot }}{\left[\frac{{\dot{\omega }}_{\mathrm{Rel}}}{3}(1-e^2)\right]}^{3/2}{\left(\frac{P_{\mathrm{b}}}{2\pi }\right)}^{5/2}=2.47346(28).\end{array}$$(2)

This mass value is slightly larger and more than twice as precise than the last published work (R19); the stability of the observed $\dot{\omega }$ values (see Appendix B) suggest that it is robust.

However, there are additional effects that contribute to ${\dot{\omega }}_{\mathrm{obs}}$. The potentially largest is due to spin-orbit coupling (${\dot{\omega }}_{\mathrm{SO}}$), if the companion WD is fast rotating. This includes the contributions due to Lense-Thirring precession and the one caused by the spin-induced quadrupole moment of the rotating WD companion (Venkatraman Krishnan et al. 2020). For the 1.09 M_⊙ companion of NGC 1851A (see next section), the maximum contribution, when the WD is rotating close to break-up, is ${\dot{\omega }}_{\mathrm{SO}}\approx 0.03{\dot{\omega }}_{\mathrm{obs}}$. As it becomes clear below, even such an extreme rotation of the companion would have no influence on the conclusions of the present work.

The second largest is due to the interaction with a previously undetected third body, and is calculated in detail in Sect. 5. Additional contributions include ${\dot{\omega }}_{\mathrm{PM}}$, arising due to the proper motion of the system (Kopeikin 1996). For the measured proper motion of NGC 1851A as given in Table 2, its maximum value is ${\dot{\omega }}_{\mathrm{PM}}=8.72\times {10}^{-7}\mathrm{deg}{\mathrm{yr}}^{-1}$. This is comparable to the current measurement uncertainty in ${\dot{\omega }}_{\mathrm{obs}}$, and thus can be further neglected for all analysis.

The nominal measurement of ${\dot{\omega }}_{\mathrm{obs}}$ and its constraints are displayed in red in Figure 2. From the point where it intersects the constraint from the mass function (which excludes all points in gray), we derived a minimum companion mass of M_c = 0.96 M_⊙ and a maximum pulsar mass of M_p = 1.51 M_⊙.

Fig. 2. Orbital inclination and mass constraints for NGC 1851A. The solid and dotted lines represent the nominal values and 1-σ constraints for the post-Keplerian parameters (red: $\dot{\omega }$, cyan: γ, purple: h3) obtained with the DDH binary model. Bottom left plot: cos i versus the companion mass (Mc); Bottom right plot: Pulsar mass (Mp) versus Mc. The grey portions on the two plots indicate the regions excluded for Mp < 0 and sin i > 1 respectively. The normalised 1-D pdfs for cos i, Mp and Mc are shown in the top left, top right, and right panels.

3.4. Einstein delay

Owing to the long timing baseline of NGC 1851A, we precisely measured another relativistic parameter, the amplitude of the Einstein delay (γ) for the system. In GR, it is given by:

$$\begin{array}{c}\hfill \gamma =T_{\odot }^{2/3}e{\left(\frac{P_{\mathrm{b}}}{2\pi }\right)}^{1/3}\frac{{\overline{M}}_{\mathrm{c}}({\overline{M}}_{\mathrm{c}}+{\overline{M}}_{\mathrm{tot}})}{{{\overline{M}}_{\mathrm{tot}}}^{4/3}}·\end{array}$$(3)

This delay quantified by γ results from the combination of the varying special relativistic time dilation due to the changing velocity of the pulsar in its orbit and the varying gravitational redshift due to the companion. The Einstein delay was earlier measured by R19, and this allowed for the determination of the individual masses for the system. However, in this work, we report an updated value, γ = 0.01925 (48) s. The constraints on γ are indicated by the cyan lines in Fig. 2. This value is 1.8 times more precise than the value determined by R19, but it is also 2.6-σ smaller. This difference is discussed in more detail in Appendix B.

With this measurement and the knowledge of M_tot as derived above, we derived the masses using Equations (21) and (22) from R19. For NGC 1851A, under the assumption of GR, we obtained M_c = 1.09 (2) M_⊙, M_p = 1.38 (2) M_⊙ and i = 61.5deg or i = 118.5deg. These differ from the mass estimates of R19 because of the difference in the measured γ.

As discussed by R19, for wide binaries the measurement of γ is covariant with secular variations of the projected semi-major axis (ẋ). In addition, they assumed that the ẋ of NGC 1851A originates from the proper motion of the system; this is an effect that must be taken into account for any wide binary system. However, the realization that the NGC 1851A system has a nearby companion means that perturbations from that third object can also contribute to ẋ (ẋ_tidal), and therefore affect γ and the mass measurements. The magnitude of this effect is estimated in Sect. 5, and the individual masses are estimated in Sect. 6.

3.5. Shapiro delay

With an aim to quantify the relativistic light-propagation delay in the system, we took extensive amounts of data at periastron and superior conjunction. In our DDH orbital model, we fixed an estimate of the orthometric ratio ς derived from i, $\varsigma =\frac{sini}{1+|cosi|}=0.59$; from this we obtained a measurement of the orthometric amplitude h₃ = 0.52 (53) μs. This indicates that we do not detect any Shapiro delay for the system. This is in agreement with the expected non-significant measurement due to the face-on inclination of the system derived from the other PK parameters.

4. Period derivatives of a pulsar in a globular cluster

Many of the new parameters we measure in this work are the derivatives of the spin and orbital frequencies. Given the common origin of these terms, and some of the subtleties in their analysis, we discuss them at length in this dedicated section. This includes new material that is relevant for the analysis of, in principle, all pulsars in GCs where statistically significant changes in the acceleration are measurable, i.e., the vast majority of GC pulsars with long timing baselines. In the following section, we discuss the physical causes of the acceleration of this system and its variation over time.

4.1. A theoretical interpretation

As we cannot measure radial motion, we observe modified ‘Doppler-shifted’ values of the intrinsic spin period of the pulsar (P_int) and orbital period of the binary (P_b, int). The observed parameters are related to these values as:

$$\begin{array}{cc}\hfill P_{\mathrm{obs}}& =D^{-1}{P}_{\mathrm{int}}\simeq P_{\mathrm{int}}[1+V_{\mathrm{l}}/c],\hfill \end{array}$$(4a)

$$\begin{array}{cc}\hfill P_{\mathrm{b},\mathrm{obs}}& =D^{-1}{P}_{\mathrm{b},\mathrm{int}}\simeq P_{\mathrm{b},\mathrm{int}}[1+V_{\mathrm{l}}/c],\hfill \end{array}$$(4b)

where D is the Doppler shift factor (Damour & Taylor 1992), the last term is a leading order term and V_l is given by

$$\begin{array}{c}\hfill V_{\mathrm{l}}={\mathit{V}}_{\mathrm{CM}}·{\widehat{\mathit{K}}}_0,\end{array}$$(5)

where V_CM is the velocity of the centre of mass of the binary (CM) relative to the Solar System barycentre (SSB) and ${\widehat{\mathit{K}}}_0$ is the unit vector pointing from the SSB to CM. Although this radial motion also affects other distance and mass parameters of the binary, these contributions are negligible compared to the observed values and hence ignored in further analysis.

A change of V_l with time, given by a_l = dV_l/dt, causes a secular change in the Doppler shift. Differentiating Equations (4a) and (4b), we obtain (Phinney 1993):

$$\begin{array}{cc}\hfill {\left(\frac{\dot{P}}{P}\right)}_{\mathrm{obs}}& ={\left(\frac{\dot{P}}{P}\right)}_{\mathrm{int}}+\frac{a_{\mathrm{l}}}{c},\hfill \end{array}$$(6a)

$$\begin{array}{cc}\hfill {\left(\frac{\dot{P_{\mathrm{b}}}}{P_{\mathrm{b}}}\right)}_{\mathrm{obs}}& ={\left(\frac{\dot{P_{\mathrm{b}}}}{P_{\mathrm{b}}}\right)}_{\mathrm{int}}+\frac{a_{\mathrm{l}}}{c}·\hfill \end{array}$$(6b)

For a pulsar in a GC, a_l has contributions from multiple factors, and can be explicitly written as

$$\begin{array}{c}\hfill a_{\mathrm{l}}=a_{\mathrm{Shk}}+a_{\mathrm{Gal}}+a_{\mathrm{GC}}+a_{\mathrm{NS}},\end{array}$$(7)

where a_Shk indicates the apparent acceleration due to Shklovskii effect (Shklovskii 1970), a_Gal is the difference between the accelerations of the Solar System and NGC 1851 in the field of the Galaxy, a_GC is the line-of-sight acceleration of the pulsar in the gravitational field of the GC and a_NS is the acceleration in the field of nearby stars. A detailed investigation of the various accelerations is presented in Sect. 5.

To derive the second order spin period and orbital period derivatives of the system, we differentiated Equations (6a) and (6b). For a pulsar in a globular cluster, the change in the line-of-sight acceleration is much larger than the change in the first term on the right of these equations and thus dominates the contribution to the second order derivatives (Phinney 1992). Under this consideration, and using Equation (7) from Freire et al. (2017), we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{\ddot{P}}{P}\right)}_{\mathrm{obs}}& \simeq \frac{{\dot{a}}_{\mathrm{l}}}{c}\simeq -\left(\frac{\ddot{f}}{f}\right),\hfill \\ \hfill {\left(\frac{\ddot{P_{\mathrm{b}}}}{P_{\mathrm{b}}}\right)}_{\mathrm{obs}}& =\frac{{\dot{a}}_{\mathrm{l}}}{c}\simeq -\left(\frac{\ddot{f_{\mathrm{b}}}}{f_{\mathrm{b}}}\right),\hfill \end{array}\end{array}$$(8)

where f and f_b denotes the spin frequency of the pulsar and the orbital frequency of the binary, respectively. Note that in these equations, any ‘intrinsic’ terms are very likely to be negligible, a point of importance for what follows.

Equation (8) can be differentiated to get quite accurate expressions for further higher order derivatives, for k > 2, as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{1}{P}\frac{d^{k}P}{d{t}^k}& =\frac{1}{c}\frac{d^{k-1}{a}_{\mathrm{l}}}{d{t}^{k-1}}\simeq \frac{f^{(n)}}{f},\hfill \\ \hfill \frac{1}{P_{\mathrm{b}}}\frac{d^k{P}_{\mathrm{b}}}{d{t}^k}& =\frac{1}{c}\frac{d^{k-1}{a}_{\mathrm{l}}}{d{t}^{k-1}}\simeq \frac{f_{\mathrm{b}}^{(n)}}{f_{\mathrm{b}}}·\hfill \\ \hfill \end{array}\end{array}$$(9)

The values of all the spin and orbital period derivatives for NGC 1851A are given in Table 3, and a detailed discussion for each of them is given below.

4.2. Measured spin period (frequency) derivatives

The period derivatives of pulsars in GCs, unlike the ones in the Galactic disk, are highly influenced by the gravitational field of the globular cluster and nearby stars. In this section, we discuss the spin period (equivalently, spin frequency) derivatives for NGC 1851A. These are measured for two different reference epochs: Epoch 1 (MJD = 53623, 10 September 2005) was used by Freire et al. (2007) and R19, Epoch 2 (MJD = 59728, 29 May 2022) is used in Table 2. This latter epoch is within the recent phase of intense MeerKAT timing; for this reason the timing parameters in Table 2 have smaller correlations than if they were measured for Epoch 1.

4.2.1. First spin frequency derivative

We report an updated and much more precise measurement of the first spin frequency derivative for NGC 1851A (see Table 3), which translates into a spin period derivative of Ṗ_obs = 2.3538 (20) × 10⁻²⁰ ss⁻¹. This is ∼34 times larger than the value reported by R19, but measured for Epoch 2. This change is predominantly due to a large variation in the acceleration of the system: if we instead used Epoch 1, we obtained a value of ḟ that is consistent with that presented in those works, but still ∼27 times smaller than the value measured for Epoch 2. We now discuss this large change of Ṗ.

4.2.2. Higher spin frequency derivatives

In the extremely dense core of a GC, a pulsar experiences variations of its acceleration due to its movement in the gravitational potential of the cluster and its interaction with the nearby stars (Phinney 1992). Equation (8) shows that this change in the line-of-sight acceleration of the system, also known as the line-of-sight ‘jerk’ (ȧ_l), is manifested by the second period (frequency) derivative of the pulsar (Table 3). From its value at Epoch 2, we estimated the jerk for NGC 1851A as

$$\begin{array}{c}\hfill {\dot{a}}_l=3.20\times {10}^{-18}\mathrm{m}{\mathrm{s}}^{-3}.\end{array}$$(10)

Although a similar value for $\ddot{f}$ was measured by R19, they did not comment on its unusual size, which is orders of magnitude larger than the observed values for other MSPs in globular clusters (Freire et al. 2017; Prager et al. 2017), and only a factor of 5 smaller than observed in PSR B1620−26, the triple system in the globular cluster M4 (Thorsett et al. 1999). This is so large that it causes the large variation of Ṗ_obs between Epochs 1 and 2. Inevitably, this jerk affects the observed orbital period and its derivatives, as discussed in Sect. 4.3. As mentioned above, we discuss the implications of this unusually large jerk in Sect. 5.

The long timing baseline and the precise timing solution allowed for a 19.4-σ measurement of the third spin frequency derivative for NGC 1851A as given in Table 3. Its negative value implies that the large jerk (with a negative sign in $\ddot{f}$) is becoming even larger in magnitude. For a pulsar in the core of such a dense GC, such higher order derivatives of the spin period are fairly common: unlike the lower spin frequency derivatives, which can be caused by the mean field of the GC, the higher order derivatives are caused by the gravitational fields of nearby stars (Phinney 1992).

4.2.3. Evolution of spin frequency derivatives

We now look in detail at the evolution of the spin frequency (f) and its time derivatives (ḟ, $\ddot{f}$, and $\stackrel{⃛}{f}$), using them as the coefficients in a Taylor expansion (Lorimer & Kramer 2004),

$$\begin{array}{c}\hfill f(t)=f_0+{\dot{f}}_0(t-t_0)+\frac{1}{2}{\ddot{f}}_0{(t-t_0)}^2+\frac{1}{6}{\stackrel{⃛}{f}}_0{(t-t_0)}^3+\dots ,\end{array}$$(11)

around Epoch 1 and 2. Figure 3 shows the evolution of the spin frequency and its first and second derivatives over the total timing baseline of this pulsar. As discussed previously, the high jerk for the system is clearly evident from the evolution of the first spin frequency derivative with time. In the last plot, we can also appreciate how fast the jerk is increasing with time, having increased by about 65% between Epochs 1 and 2.

Fig. 3. Top to bottom: Evolution of the observed spin frequency (f), its first derivative (ḟ), and its second derivative ($\ddot{f}$) as a function of time (MJD). The vertical yellow and green lines indicate Epochs 1 and 2 respectively. The values from R19 are shown in red, and the results from this work are shown in blue. The large change in the ḟ value between the two epochs can be clearly seen.

Despite these changes between Epochs 1 and 2, and between our Epoch 1 and the values measured by R19, we can see that the spin evolution we observe is mostly consistent with that described by them, at least for the period of time covered by their earlier timing solution.

4.3. Orbital period derivatives

The extension of the timing baseline was hugely beneficial for a refined measurement of the change in the orbital period for NGC 1851A. The inclusion of orbital frequency derivatives higher than first order is, as explained below, necessary for a correct description of the orbital evolution of the system, and has far-reaching implications.

4.3.1. First orbital period derivative

Another quantity that has changed considerably relative to R19 is the measurement of the orbital period derivative of the system. Equation (4b) shows that the observed change in the orbital period, like the spin period, is also affected significantly by a secular change in the Doppler shift; thus the acceleration should affect the Ṗ and Ṗ_b of a pulsar in a globular cluster in a similar way. The unexpectedly large positive value of Ṗ_b observed by R19 could not be explained by a large acceleration contribution from the cluster’s field: such an acceleration would result in a value of Ṗ that is much larger than that observed. The expected intrinsic variation in the orbital period (Ṗ_b,int) from orbital decay due to gravitational wave emission is orders of magnitude smaller than the observed value, and also of negative sign, and thus could not account for the large positive value of Ṗ_b.

In this work, we revisited this problem with a much longer dataset and planned observations to check the robustness of the measurement. The updated value of the Ṗ_b,obs for this system, derived using the DDH model, is given in Table 3. This measurement is 13 times more precise than the estimate by R19, and 2.6-σ smaller. The reasons for the difference are in part related to the discovery of the 2-second clock offset between the earlier GMRT dataset, but also to the decrease in correlation between different parameters. Furthermore, as we see below, the jerk should also cause this value to change with time.

Knowing the masses (from $\dot{\omega }$ and γ), we used GR (using the equation of Peters 1964) to estimate the orbital decay due to the emission of gravitational waves (GW), which we assumed is the intrinsic variation of the orbital period. This gave Ṗ_b,int = −0.153 × 10⁻¹² ss⁻¹, which is one order of magnitude smaller than Ṗ_b,obs; this means that the latter is dominated by the contribution from a_l.

4.3.2. Higher orbital frequency derivatives

Higher order orbital frequency derivatives are generally not measurable for MSPs in globular clusters, unless they are in ‘spider’ systems (Freire et al. 2017); in these cases the variations are not related to variations in the accelerations of those systems and arise due to the changing gravitational quadrupole moment of their companion (e.g. Ridolfi et al. 2016). The large line-of-sight jerk for NGC 1851A, which changes the Ṗ_obs by a factor of 27 between the two reference epochs, should cause a similar correlated change in the Ṗ_b,obs for the pulsar. Even though the latter effect is not clearly measurable, it should be taken into consideration; as we’ll see below this makes a real difference.

A precise description of the variation of the orbital phase requires the inclusion of higher order derivatives in the orbital frequency (f_b = 1/P_b), described as orbital frequency derivatives (f_b⁽ⁿ⁾) for the system. Since the higher derivatives of f and f_b are affected solely by the derivatives of the acceleration – as mentioned above, there should be no intrinsic contributions to them – the f_b⁽ⁿ⁾ (n  ≥  2) values can be estimated directly from the corresponding values of f⁽ⁿ⁾ by rewriting Eq. (9):

$$\begin{array}{c}\hfill f_{\mathrm{b}}^{(n)}=f^{(n)}\frac{f_{\mathrm{b}}}{f},\end{array}$$(12)

which assumes that these derivatives are measured at the same time. Such an estimate is intrinsically much more precise than any direct measurement of f_b⁽ⁿ⁾ given the much higher precision in the measurement of f⁽ⁿ⁾, i.e., in the measurement of the spin phase relative to the orbital phase. We modified our timing models to take into account the jerk and its variation in a self-consistent way using these equations (see Appendix E). The values of f_b⁽ⁿ⁾ derived from the values of f⁽ⁿ⁾ are quoted in the second part of Table 2.

As we see next, a direct argument for the necessity of taking the f_b⁽ⁿ⁾ into account comes from the inconsistent estimates of intrinsic spin period derivative, Ṗ_int (and therefore of quantities derived from it assuming rotating dipole model, like characteristic age, surface magnetic field and spin-down power) for Epochs 1 and 2, that we obtain when the change in the orbital period is estimated by a single derivative. The agreement of the values of Ṗ_int after the addition of ${\ddot{f}}_{\mathrm{b}}$ and ${\stackrel{⃛}{f}}_{\mathrm{b}}$ can be seen in the second part of Table 3.

4.4. Intrinsic spin period derivative and the characteristic age

The similarity of Equations (4a) and (4b) means that the difference of these two can be used to estimate the intrinsic change in the spin period of the pulsar as

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

In Table 3, we calculate this parameter and the characteristic age for Epochs 1 and 2, using the equation in Lorimer & Kramer (2004), from the orbital and spin period derivatives measured at each epoch. In the first part of the table, we see how neglecting any higher orbital period derivatives produced completely different ‘local estimates’ of Ṗ_int. In the second part of the table, we see how taking the higher order estimates of ${\ddot{f}}_{\mathrm{b}}$ into account leads to consistent estimates for Ṗ_int. This quantity is intrinsic to the pulsar and (for MSPs) should not change significantly with the time at which it is estimated. This highlights the importance of taking these higher order orbital frequency derivatives into account.

We can obtain the characteristic age (τ_c) and surface magnetic field (B) from the intrinsic spin period derivative (Ṗ_int) of the pulsar. For NGC 1851A, we estimate τ_c = 4.8 Gyr, and B = 2.8  × 10⁸ G, making this pulsar a rather typical MSP.

5. Investigating the external influence on the binary

5.1. On the need for a third object

It is quite evident from Fig. 3 that the pulsar’s apparent spin-up changed to an apparent spin-down just a few years before Epoch 1: this, and the continued large change in the apparent spin-down between Epochs 1 and 2 hint towards the fact that the system is being substantially perturbed in the globular cluster environment.

As pulsars slow down due to the loss of rotational energy, the intrinsic spin period derivative (first term on the right of Equation (6a)) has a small, positive and constant value. The measured value of Ṗ_obs thus implies an upper limit on the line-of-sight acceleration of the pulsar, and the latter estimated at Epoch 2 using Equation (6a), is $a_{\mathrm{l},\max }=c{(\dot{P}/P)}_{\mathrm{obs}}=1.41\times {10}^{-9}\mathrm{m}{\mathrm{s}}^{-2}$.

The various contributions to the acceleration are given explicitly in Equation (7), and the values are discussed here. In this work, we have precisely measured the total proper motion of the system, μ_T = 2.76 (12) mas yr⁻¹. This, and the distance to the cluster (d ≈ 11.7 kpc, Libralato et al. 2022) yield a_Shk = μ²d = 6.5 (6)  × 10⁻¹¹ m s⁻². Using the Milky Way potential from McMillan (2017) and the Galactic coordinates of the globular cluster, we obtained a_Gal = −1.174  × 10⁻¹¹ m s⁻². The sum of a_Shk and a_Gal is about 27 times smaller than a_l, max at Epoch 2 given above. For most pulsars, the explanation for this is that a_GC is much larger, and the effect of the passing stars (a_NS) is negligible in comparison (Phinney 1993). However, the large change in the line of sight acceleration between Epochs 1 and 2 suggests otherwise.

The a_GC term can also change with time. NGC 1851 has a high velocity dispersion of 10.2 km s⁻¹ and a small core radius of 5$\stackrel{″}{.}$4 (Miocchi et al. 2013); using these values and Equation (10) from Freire et al. (2017), we estimated the maximum line-of-sight jerk that could be induced by the gravitational potential of the cluster. We used the escape velocity of the cluster, 42.9 km s⁻¹ (Baumgardt & Hilker 2018), as the maximum velocity of the pulsar relative to the cluster, and the projected separation of NGC 1851A from the core of the cluster as given in Table 2. We derived

$$\begin{array}{c}\hfill |{\dot{a}}_{\mathrm{l},\mathrm{GC},\max }|=1.40\times {10}^{-19}\mathrm{m}{\mathrm{s}}^{-3},\end{array}$$(14)

which is 23 times smaller than the observed jerk for the system. This suggests that ȧ_l is not due to the mean field of the cluster, but from the varying accelerations due to close encounters with nearby stars, which can be significantly larger (Phinney 1993; Prager et al. 2017).

Not only does this make it evident that the most likely cause for the large observed jerk is an external mass, but it also opens up two scenarios. If the third mass is bound to the inner binary, we have a system similar to PSR B1620−26 (Thorsett et al. 1993, 1999); such triple systems could be formed through dynamical encounters in these dense environments. In the alternative scenario, the third body is not bound to the inner binary and is merely flying by. Such a phenomenon has not been previously observed for any pulsar in a globular cluster. The probability of this for the M4 triple system was negligible (Thorsett et al. 1999); however, given the extremely dense and compact core of NGC 1851, this possibility cannot be ruled out for NGC 1851A.

We determined the closest distance of approach (b_star) of a nearby mass to induce the observed jerk for NGC 1851A (Equation (4.4) from Phinney 1993), and the probability of a passing star being so close is given by (Phinney 1993)

$$\begin{array}{c}\hfill p_{\mathrm{star}}(<b_{\mathrm{star}})=1-\exp [-(4\pi /3)n{b}_{\mathrm{star}}^3],\end{array}$$(15)

where n is the central density of the globular cluster. For a star of mass 1 M_⊙, the estimated closest approach is 498 astronomical units (AU), and the corresponding probability is 0.196. The variation of the b_star and p_star values as a function of masses are shown in Fig. 4 below. These values suggest that for NGC 1851A, the probability of a flyby is a fair possibility, owing to the high stellar density of NGC 1851. In the next section, we discuss further implications for these values and investigate the nature of the external influence in bound and unbound orbits.

Fig. 4. Closest distance of approach and probability, for a nearby star of mass M in the globular cluster, to induce the observed jerk in the NGC 1851A system.

5.2. Developing a three-body model

The precise measurements of ḟ, $\ddot{f}$, and $\stackrel{⃛}{f}$ for NGC 1851A allowed for a characterisation of the external dynamical influence on the system. An attempt to fit for f⁽⁴⁾ and f⁽⁵⁾ in our timing solution revealed that they are not yet precisely measured and have large uncertainties, and this has prevented a complete model for this three-body problem. However, for our analysis, we used the preliminary measurements to obtain further constraints on the solutions. An added uncertainty comes from the observed line-of-sight acceleration of the system, which is well within the maximum possible acceleration due to the globular cluster (as discussed above), and thus the exact contribution of the nearby star to ḟ is not determinable.

We developed a model to use the first three time derivatives of the pulse frequency for NGC 1851A and the limits on the fourth and fifth derivatives to derive the mass and orbital parameters of the external mass. With the caveats mentioned above, the dynamically induced frequency derivatives solely due to the third body can be given by

$$\begin{array}{c}\hfill f^{(n)}=-f\frac{{\mathit{a}}^{(n-1)}·{\widehat{\mathit{K}}}_0}{c},\end{array}$$(16)

and a can be explicitly expressed using the orbital orientation of the binary and the eccentricity of the external orbit. We followed the method outlined in Joshi & Rasio (1997) (hereafter JR97) and Perera et al. (2017), and the equations used in our analysis are presented in Appendix C. The inferred orbital period of the external mass is much larger than that of the inner binary, and this allows us to make a good approximation of the system as the combination of two non-interacting Keplerian orbits (Thorsett et al. 1999). The orbital elements for this three-body model are illustrated in Fig. 5. The eccentric binary, currently the system NGC 1851A, is treated as a single object of mass M_tot = M_p + M_c, in an orbit with period P_b, M, eccentricity e_b, M, and semi-major axis x_M. The true anomaly, longitude of periastron, and orbital phase of the centre of mass of this binary with respect to the entire triple configuration are λ_M, ω_M, and ψ_M. The outer companion has a mass m₃, and the orbital elements are P_b, 3 (for the case of an elliptical orbit), e₃, and x₃. The corresponding angles for the external mass are λ₃ and ω₃, where λ₃ = λ_M and ω₃ = ω_M + 180deg. This external mass is moving at a velocity v₃ at a radial distance r₃ with respect to the centre of mass of the inner binary.

Fig. 5. Definition of the orbital elements for a three-body model with two Keplerian orbits. The ‘inner binary’ is represented by Mtot and the external mass by m3, and both the scenarios of an elliptical orbit and a hyperbolic fly-by encounter are shown.

We left the eccentricity of the external orbit (e₃) as a free parameter, and selected a value for λ_M and ψ_M for each trial (under the assumption of uniform prior probability distributions). Using the measurements of $\ddot{f}$ and $\stackrel{⃛}{f}$ and Equations (C.4) and (C.5), we did an initial filtering of solutions to identify those where the contribution to the ḟ from the external mass could be up to 5 times the observed value (to account for the acceleration of NGC 1851). We then calculated the estimated values for f⁽⁴⁾ and f⁽⁵⁾, and rejected further trials by comparing them with the measured limits. For our calculations, we used 3-σ upper limits on the f⁽⁴⁾ and f⁽⁵⁾ values obtained from our timing analysis. This has allowed us to exclude the solutions for bound orbits with very short orbital periods and unbound orbits with very close distances of approach. For all the allowed solutions, we calculated m₃ and a₃ and derived the relevant parameters for the external mass.

5.3. Case I: a second mass in a hierarchical triple system

We first considered the scenario where the system is in a hierarchical triple. This refers to a binary with the pulsar and its massive WD companion very close to each other, and the external third body as a second companion in a much longer orbit compared to the separation of the inner binary.

Initially, we assumed that the external orbit is circular (e₃ = 0), and attempted to derive a complete solution for the second companion. The measured first and third frequency derivatives of the pulsar can be used in Equations (24) and (25) from JR97 to directly derive λ_M and ${\dot{\lambda }}_M$. However, the same signs of the derivatives produce an imaginary ${\dot{\lambda }}_M$, and straightaway eliminates the possibility of a circular orbit for the second companion.

We next applied our method described in the previous section to elliptical orbits for the second companion, and looked at a discrete range of eccentricities between 0.1 and 0.99, with the allowed space 0  ≤ λ_M,  ω_M,  ψ_M < 360 deg. An initial analysis performed with only the values of ḟ, $\ddot{f}$, and $\stackrel{⃛}{f}$ suggested that the second companion could have a mass as low as 2.08 × 10⁻¹⁰ M_⊙ and an orbital period as small as 98 days. To ensure consistency with the timing baseline of the system and its long-term stability, it is important that the second companion has a significantly larger orbit and does not come too close to the inner binary. This motivated us to further use the limits on f⁽⁴⁾ and f⁽⁵⁾ derived from the timing analysis, and we obtained a three-parameter family of solutions as illustrated in Fig. 6.

Fig. 6. Family of solutions for the mass m3, orbital period Pb, 3, and relative distance R for a second companion of NGC 1851A using the limits on all the frequency derivatives mentioned in the text. The different eccentricities are shown in different colours. The constant lines of ${\dot{\omega }}_{\mathrm{tidal}}/{\dot{\omega }}_{\mathrm{obs}}$ are shown with dotted black lines in the bottom plot, while the dash-dotted crimson lines indicate constant ẋtidal.

For an elliptical outer orbit, we concluded that a second companion to NGC 1851A could have a mass m₃ between 6.6  ×  10⁻⁵ to 2.3 M_⊙, with an orbital period larger than 50 years and going up to several Myr. Given the poor constraints on the current measurements of the higher frequency derivatives, the exact mass of the companion and the parameters of the outer orbit are undetermined. The nature of the objects that satisfy the allowed range of masses is pretty variable, ranging from sub-stellar objects like brown dwarfs and planets to stellar mass objects, neutron stars and black holes.

5.4. Case II: a hyperbolic encounter with a passing star

We also considered the scenario where a hyperbolic encounter with a close-passing star is the cause for the remarkably large ȧ_l observed for NGC 1851A. A hyperbolic, or fly-by encounter, is a scenario when a nearby star passes by the target object in an unbound, hyperbolic trajectory. Such a close fly-by encounter involving an external star and a binary pulsar, capable of inducing such a large jerk in the latter, has not been previously observed. However, it has been suggested that clusters with high stellar density would be favorable hosts (Mukherjee et al. 2021). This proposition for NGC 1851A is motivated by the exceptionally high central density of NGC 1851. The cluster also exhibits an intermediate rate of stellar interactions per binary, as evident by NGC 1851A itself being the product of a secondary exchange encounter. This parameter exponentially increases the likelihood of stellar encounters and close interactions. This makes the dynamical influence from a nearby star a likely possibility, given that the pulsar is located at a separation of 1$\stackrel{″}{.}$58 from the core.

As detailed earlier, we varied both the λ_M and ψ_M values and analysed orbits with eccentricities (e₃) 1.1, 2.0, 3.0, and 5.0. Unlike a bound orbit, the true anomaly for a hyperbolic orbit varies between −π + cos⁻¹(1/e₃) and π − cos⁻¹(1/e₃). For all eccentricities explored in this work, we find that the λ₃ is always less than 0, suggesting that the external mass is approaching the binary in the core of the globular cluster. We implemented an additional criterion on the asymptotic velocity of the external mass (v_∞), constraining it to be less than the central escape velocity of the cluster (v_esc). With the increase in the eccentricity of the hyperbolic orbit, only the solutions with high velocity and larger distance are recovered. The allowed family of solutions are shown in Fig. 7.

5.5. Secular perturbations to the binary orbit

It is very important to discuss the dynamical interactions for such a three-body system. We looked at the secular perturbations induced on the inner orbit due to the external mass, under the approximation that the latter has a fixed position in space with respect to the inner binary. These perturbative effects are mostly evident in the orbital eccentricity (e), the longitude of periastron (ω), and the orbital period (P_b).

The changes in P_b, due to this external mass, have been discussed in detail in the previous sections, and the changes in e are not yet measurable. Here we make predictions for the change in the longitude of periastron of the inner binary due to the external mass (${\dot{\omega }}_{\mathrm{tidal}}$) by calculating the orbital perturbation rates. For all possible orbital orientations, we estimated the maximum tidal contribution using the mass and radial separation of the external mass and the orbital parameters of the inner binary (Equation (D.14)). All detailed equations and calculations are presented in Appendix D.

In the bottom panels of Figs. 6 and 7, we show the tidal contribution to the observed $\dot{\omega }$ value for bound and unbound orbits. For all eccentricities ranging from 0.1 to 5.0, ${10}^{-6}\le {\dot{\omega }}_{\mathrm{tidal}}/{\dot{\omega }}_{\mathrm{obs}}\le 1.3\times {10}^{-3}$, suggesting that the effect of the external body on the ${\dot{\omega }}_{\mathrm{obs}}$ could be an order of magnitude larger than the measurement uncertainty of $\dot{\omega }$, but this changes the total mass measurement by about 10⁻³. We can thus conclude that the ${\dot{\omega }}_{\mathrm{obs}}$ is mostly relativistic, and can be used for inferring the total mass of the system.

Fig. 7. Allowed solutions showing the mass m3, asymptotic velocity v∞, and relative distance R for a nearby star in a fly-by encounter, in the core of the globular cluster. Constant ${\dot{\omega }}_{\mathrm{tidal}}/{\dot{\omega }}_{\mathrm{obs}}$ and constant ẋtidal are shown with black dotted and crimson dash-dotted lines respectively in the bottom plot.

In these panels, we also see that the maximum ẋ_tidal is of the order of 1.4 × 10⁻¹⁴ s s⁻¹, and normally much less than this value. Since, as mentioned above, this value is correlated with the value of γ, it can also affect the mass measurements for the individual components of the binary, as discussed in the following section.

5.6. Optical constraints on the nature of the companions of NGC 1851A

Given the old age of NGC 1851 (11.0 Gyr, VandenBerg et al. 2013), any stars still in the MS must have a mass of ∼ 0.8 M_⊙ or smaller: more massive stars will have long left the MS. The possible exception are blue straggler stars (BSS), which might have masses just up to twice that, ∼ 1.6 M_⊙. From our calculations in the previous section, we concluded that any MS star must have a radial distance (R) from NGC 1851A smaller than about 800 au (bottom panel in Fig. 7). For a BSS with a higher mass, this distance does not change much. Given the parallax of this GC, 0.0857 mas, the region within R will have an angular radius of $\sim 0\stackrel{″}{.}069$.

We inspected archival HST images of NGC 1851; these are listed and described in detail by Barr et al. (2024). To locate the pulsar, we used the same astrometric solution that was used by the latter authors when searching the optical counterpart of NGC 1851E. The vicinity of the pulsar can be seen in Fig. 8, which includes part of an image taken with HST’s Wide Field Camera 3, the blue circle representing the region of interest, within 0$\stackrel{″}{.}$069 of NGC 1851A. Given the pulsar’s proximity to the cluster centre, where the crowdedness is severe, the faintest MS stars that can be detected in these images have masses of ∼0.7 M_⊙.

Fig. 8. Location of NGC 1851A, with the blue circles corresponding to radius of 0$\stackrel{″}{.}$069 around the position of the pulsar as estimated for the time the left frame was taken (MJD = 55558). Given the small proper motion, correcting for it does not induce any noticeable changes in the position of the pulsar. Images from the Wide Field Camera 3, taken with the F275W filter (left) and F814W filter (right).

As we can see, no sources are detected within the region of interest, either in the near-UV image on the left (taken with the F275W filter) or the near-IR image on the right (taken with the F814W survey). This excludes the possibility of the third star being a BSS, a giant or a MS star with a mass above ∼0.7 M_⊙.

Apart from obtaining constraints on the third star, these images also allowed us to obtain definite conclusions regarding the massive binary companion of NGC 1851A. If it were an extended star, we would not be able to rely on the PK parameters to measure its mass. However, from the mass function we derived a lower limit for the companion mass of 0.84 M_⊙ (Sect. 3.2). Any BSS, MS or giant star with such a mass or larger should be clearly detectable at the position of the pulsar. This non-detection confirms that the companion of NGC 1851A is a compact object, either a WD or a NS. Such a conclusion was already expected from the lack of eclipses and DM variations in the timing data (R19).

6. Self-consistent analysis of masses and orbital orientation

The compact nature of the companion of NGC 1851A implies that the measured PK parameters quantify relativistic effects in the orbit, with small perturbations caused by the tides induced on the binary by the third object. This means that, within the uncertainty limits imposed by those perturbations, the GR equations of these PK parameters can be used to infer the masses of the components of the NGC 1851A system. We now proceed to do this.

6.1. Bayesian map

As discussed in detail by R19, for wide binary systems, the effect of the Einstein delay (γ) is indistinguishable from a secular change of projected semi-major axis (ẋ). For this reason, we cannot separate both effects, but only measure their sum. Therefore, when estimating the masses from $\dot{\omega }$ and γ, we must take into account possible ẋ contributions to our measurement of γ, which we can only estimate theoretically.

R19 did these calculations and came to the conclusion that by far the largest contribution to ẋ is caused by the change of viewing angle that originates from the proper motion of the system (Arzoumanian et al. 1996; Kopeikin 1996). This depends on the proper motion (with total magnitude μ and position angle Θ_μ), the unknown position angle of the line of nodes of the binary (Ω) and the orbital inclination. Rewriting the formula given by Kopeikin (1996), we obtain:

$$\begin{array}{c}\hfill {\dot{x}}_{\mu }=x\mu cotisin({\mathrm{\Theta }}_{\mu }-\mathrm{\Omega }).\end{array}$$(17)

Using our new estimate of the proper motion (2.76 mas yr⁻¹), and with i ∼ 60deg, we obtained |ẋ_μ| ≤ 0.9 × 10⁻¹⁴ s s⁻¹. This is smaller than estimated by R19, both because of the smaller proper motion we measured, but also because of the more edge-on inclination. Nevertheless, the uncertainty in the value of γ is also smaller, so that the effect of the proper motion on the measurement of γ and in the mass estimates is still relevant.

However, the recognition that NGC 1851A has a nearby companion, and the current limits on its distance and mass imply ẋ_tidal = 1.4 × 10⁻¹⁴ s s⁻¹, which is about 1.5 times larger than ẋ_μ. However, this is an absolute maximum value: for most possible configurations (including the one that yields this maximum value), ẋ_tidal values are significantly smaller. So the assertion that ẋ_μ is the largest contribution to ẋ is very likely correct.

These contributions must be compared to the apparent ẋ that arises from γ and its uncertainty, which according to R19 is given by:

$$\begin{array}{c}\hfill {\dot{x}}_{\gamma }=-\frac{\gamma \dot{\omega }}{\sqrt{1-e^2}}sin\omega ,\end{array}$$(18)

which, for the values in Table 2, results in ẋ_γ = −30.0 ± 0.7 × 10⁻¹⁴ s s⁻¹, which is more than 20 times larger than the maximum possible ẋ_tidal. We thus conclude that our measurement of γ has at least that significance. Furthermore, the uncertainty on the measurement of γ is now smaller than the maximum value of ẋ_μ; this means that the former no longer dominates the mass uncertainty and therefore the individual masses cannot be estimated much more precisely solely by improving γ.

R19 made a fully self-consistent estimate of the component masses and the orbital inclination of the system that takes the ẋ_μ fully into account. For this, they used the DDK orbital model, assumed the validity of GR, and made a map of the quality of fit (χ²) for every point in the orbital orientation space of a binary system with a known total mass M_tot.

Since ẋ_tidal cannot yet be estimated properly, but is likely to be smaller than ẋ_μ, the methodology proposed by R19 should still produce reliable uncertainty estimates. We followed that methodology by creating 2-D grids of the χ² values under the consideration that randomly oriented orbits have a constant probability for cos i. The model also intrinsically estimates all contributions from kinematic effects.

The whole space ranges from −1 to 1 in cos i and 0 to 360 deg in Ω: for each point in this grid, we hold the total mass (M_tot) constant and introduce the corresponding values of i and Ω in our modified DDK model using the KOM and KIN parameters. For each point in the given space, the masses of the pulsar and companion are calculated from the total mass and i using Equation (1); the Einstein delay is calculated using Equation (3) respectively. These values are then added to the model using the M2 and GAMMA parameters respectively (the SINI parameter is defined automatically from KIN). The KOM, KIN, M2, and GAMMA parameters are kept fixed in the fit. Due to the contributions from various kinematic effects, we kept both Ṗ_b and $\dot{\omega }$ as free parameters for generating the χ² grids.

We then ran TEMPO using the DDK model and allowed for a fitting of all parameters not mentioned above. The post-fit χ² value for each (cos i, Ω) combination is assigned to that specific grid point, and this value is then used to calculate the Bayesian 2D probability density function (pdf) at each point using (Splaver et al. 2002):

$$\begin{array}{c}\hfill p(\mathrm{\Omega },cosi)\propto exp\left(\frac{{\chi }_{\min }^2-{\chi }^2}{2}\right)\end{array}$$(19)

where χ_min² is the lowest χ² of the whole grid. Following this, the 2D pdf is projected on the cos i and Ω axes and the corresponding 1D pdfs are calculated, which are shown in the upper left-hand and bottom right-hand panel of Fig. 9. The 2D probability is shown in the central panel of the same figure.

Fig. 9. Central panel: 2-D pdf for the full Ω − cos i space of binary pulsars. The blue and green contours show the probability density for inclination i > 90 deg and i < 90 deg, where increasing probability density is shown with darker shades. The grey regions are excluded based on the requirement that the pulsar mass must be greater than 0, and the dotted red line indicates the position angle of the proper motion of the system. Top panel: Probability density function for cosi, normalised to the maximum. Right-hand panel: Probability density function for Ω.

6.2. Results

It is evident from Fig. 9 that the orbital inclination and Ω are correlated with each other. These variations of the best-fit values of cos i as a function of Ω are a clear indication that the effect of the proper motion-induced ẋ matters for this calculation. The identical probabilities for the two peaks of orbital inclination (i = 62.6 ± 3.0deg, or 117.4 ± 3.0deg) shown in the top plot, is a result of a non-significant detection of the annual orbital parallax, the only effect that can be used to remove this degeneracy; this non-detection is expected given the very large distance to the pulsar. Because of the non-detection of the Shapiro delay, Ω is still not constrained. The derived masses for the pulsar and the companion, to the 68.3 percent confidence limit, are M_p = 1.39(3) M_⊙ and M_c = 1.08(3) M_⊙. This makes the pulsar mass 2.8-σ larger and the companion mass 2.3-σ smaller than what R19 measured, $M_{\mathrm{p}}=1.{25}_{-0.06}^{+0.05}{\mathrm{M}}_{\odot }$ and $M_{\mathrm{c}}=1.{22}_{-0.05}^{+0.06}{\mathrm{M}}_{\odot }$. The resulting lower companion mass now strongly indicates that the companion is a massive white dwarf.

7. Summary and conclusions

In this paper, we have presented the results of our updated timing analysis for NGC 1851A, using data from the GMRT and MeerKAT, over a period of ∼18 years. We accounted for a 2-s offset in the GMRT dataset and extended the timing baseline of the system by more than 4 years since R19. This allowed for a very precise measurement of the position and proper motion of the system. The transverse velocity of the pulsar with respect to the cluster is 30 (7) km s⁻¹, and given the escape velocity of the cluster of 42.9 km s⁻¹, we can conclude that this is consistent with the pulsar’s association to the cluster.

We obtained precise measurements of three spin frequency derivatives for the pulsar. Although the first derivative for pulsars in GCs is generally dominated by the acceleration from the cluster’s mean field, the second derivative implies a line-of-sight jerk for NGC 1851A that is too high to be accounted by the latter. Additionally, the spin and orbital frequency derivatives show very large variations. In this paper we show that to estimate the orbital evolution correctly, we must take into account, higher order (> 1) orbital frequency derivatives, which are in this case estimated directly from the higher order spin frequency derivatives.

The inclusion of these higher derivatives and the aforementioned 2-s offset, in addition to the extension of the timing baseline, offered a precise measurement for Ṗ_b,obs = 2.2 (7) × 10⁻¹² ss⁻¹ at Epoch 2. This measurement, the estimate of the orbital period decay due to GW and the Ṗ_obs measured at Epoch 2 allowed for a robust measurement of the pulsar’s Ṗ_int, from which we derive a characteristic age of 4.8 Gyr and a surface magnetic field of 2.8  ×  10⁸ G. The inclusion of the higher orbital frequency derivatives at Epoch 1 derived from the higher order spin period derivatives measured there leads to the derivation of different values of Ṗ_obs and Ṗ_b,obs at Epoch 1, but from them consistent values of Ṗ_int can be obtained.

The large second spin frequency derivative, we argue, is due to the presence of a third body near the binary, and both elliptical and hyperbolic orbits are possible. Although we cannot yet conclude the exact nature of this object, we derived a family of solutions for the orbital parameters for different orbital eccentricities. This family of solutions indicates that this third object, if a main sequence star, must have a radial distance smaller than 800 au from NGC 1851A. This, in combination with optical images of the cluster, exclude the possibility of the third object being a blue straggler star, or a main sequence star more massive than 0.7 M_⊙. If this system is in a hierarchical triple, it would allow an insight into the formation of such triples in globular clusters. On the other hand, a hyperbolic encounter will make this system the first of its kind and provide a test bed for studying the dynamical interactions in this cluster. Further long-term timing measurements will offer tighter constraints on the higher frequency derivatives, and potentially distinguish between ellipsoidal and hyperbolic orbits.

The absence of any bright star at the position of NGC 1851A in HST images of NGC 1851 confirms that the binary companion to NGC 1851A is a compact object, either a NS or a massive WD. The possibility of an extended companion to NGC 1851A had been advanced earlier because of the long-timescale scintillation seen during different orbital phases of the pulsar, which could have been caused by gas emanating from the companion (Freire et al. 2007). With the discovery of a handful of other pulsars in this cluster by Ridolfi et al. (2022) and the observation of similar scintillation for all of them (including the isolated ones), we can conclude that the observed scintillation of NGC 1851A is likely not due to its companion.

In this work, we obtained improved and robust measurements of the $\dot{\omega }$ and γ parameters for the system. The contribution to the $\dot{\omega }$ from the third body is at least 3 orders of magnitude smaller than the observed value, and thus we concluded that the observed $\dot{\omega }$ is almost purely relativistic. The maximum possible value of ẋ_tidal could change γ by a factor of up to 1/20; however this number is likely much smaller, so we concluded that the observed γ is also mostly relativistic. The masses derived from a self-consistent analysis of the timing, which take into account the contribution of the proper motion to ẋ and γ, are M_tot = 2.47346 M_⊙ (the latter has a relative uncertainty due to ${\dot{\omega }}_{\mathrm{tidal}}$ that could be up to 10⁻³), M_p = 1.39(3) M_⊙ and M_c = 1.08(3) M_⊙. The latter are 2.8 and 2.3-σ different from the masses obtained by R19 respectively, $M_{\mathrm{p}}=1.{25}_{-0.06}^{+0.05}{\mathrm{M}}_{\odot }$ and $M_{\mathrm{c}}=1.{22}_{-0.05}^{+0.06}{\mathrm{M}}_{\odot }$, the reason for this is the decrease in the value of γ that has happened with the inclusion of more data. Given the recent stabilization in the value of γ, we deem these mass values to be more reliable than those published by R19. With the new mass estimate for the companion derived in this work, and considering the currently lightest known NS mass of 1.17 M_⊙ (Martinez et al. 2015), it becomes much more likely that the companion of NGC 1851A is a massive white dwarf.

Acknowledgments

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time– space reference systems department of the Paris Observatory. MeerTime data is housed on the OzSTAR supercomputer at Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government. MeerKAT observations used the PTUSE and APSUSE computing clusters for data acquisition, storage, and analysis. These clusters were funded and installed by the Max- Planck-Institut für Radioastronomie (MPIfR) and the Max- Planck-Gesellschaft. PTUSE was developed with support from the Australian SKA Office and Swinburne University of Technology. The authors also acknowledge Max-Planck-Institut für Radioastronomie funding to contribute to MeerTime infrastructure. We thank the staff of the GMRT that made the observations possible. GMRT is run by the National Centre for Astrophysics of the Tata Institute of Fundamental Research. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. AD, PCCF, TG, NW, DJC, MK, VVK, VB, and PVP acknowledge continuing valuable support from the Max-Planck Society. FA acknowledges that part of the research activities described in this paper were carried out with the contribution of the NextGenerationEU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures, Project IR0000034 – ‘STILES -Strengthening the Italian Leadership in ELT and SKA’. MBa acknowledges support through ARC grant CE170100004. YG acknowledges support from the Department of Atomic Energy, Government of India, under project numbers 12-R&D-TFR-5.02-070 and RTI4002. AP acknowledges that this work was supported in part by the “Italian Ministry of Foreign Affairs and International Cooperation”, grant number ZA23GR03, under the project “RADIOMAP- Science and technology pathways to MeerKAT+: the Italian and South African synergy”. SMR is a CIFAR Fellow and is supported by the NSF Physics Frontiers Center award 1430284. LZ is supported by ACAMAR Postdoctoral Fellowship and the National Natural Science Foundation of China (Grant No. 12103069). Finally, we appreciate the insightful comments and suggestions provided by the anonymous referee on the manuscript.

References

Appendix A: GMRT clock offsets

Two years after its release with 16 antennas, the wideband correlator (GWB) of the GMRT was upgraded to a 32−antenna version in 2017. Following this upgrade, the actual sampling time of the recorded data had a time difference with the accurate reference timestamps of observations ⁸. For the beam data acquired with the GMRT, this difference was constant. The exact values for the PA and CDP beams, for different bandwidths (BWs) and between corresponding dates, are noted in Tables 1a and 1b below respectively. All differences were corrected for beam data acquired after 14 May 2020. We include these values here for easy reference.

Appendix B: Stability of measured post-Keplerian parameters

We also did self-consistent checks for the stability of the observed $\dot{\omega }$, γ, and Ṗ_b,obs with the extension of the dataset; this is shown in Fig. B.1. The measured values of the parameters are consistent throughout the dataset and the fluctuations are all consistent within the 1-σ uncertainties. The stability of the measurements indicates that the constraints on the masses are reliable and there are no bad observations that might influence a biased measurement, suggesting that they will not change significantly in the future.

An observing campaign for NGC 1851 during July-August 2023 (Barr et al. 2024) produced a dense dataset for NGC 1851A, consisting of 18 hours of data across six observations with the UHF receiver of MeerKAT. These high-cadence observations improved the sampling and orbital coverage of the TOAs, resulting in precise and stable measurements of the parameters. As evident from the bottom panel of Fig. B.1, a robust measurement for Ṗ_b,obs was obtained at the end of the observing campaign, which remains stable with the addition of more data.

Fig. B.1. Top to bottom: Cumulative weighted TOA count, observed rate of advance of periastron ${\dot{\omega }}_{\mathrm{obs}}$, Einstein delay γ, and observed orbital period derivative Ṗb,obs as a function of time. The values and their residual uncertainties (given by the vertical errorbars) are as measured using the DDH timing model presented in Table 2.

Appendix C: Three-body orbital model

We used equations (11) - (13) given by JR97, and an extension of the same by Perera et al. 2017 (Appendix A) to derive the parameters of the external mass using the explicit expressions for the frequency derivatives of the pulsar. Unlike them, we used the $\ddot{f}$ as our primary measurement and propagated the dependence of the other derivatives on the same. All the equations have a dependence on the sum of λ and ω, the true anomaly and the longitude of periastron respectively. We defined the parameter ψ as ψ = λ + ω, and a schematic of the geometry used is shown in Fig. 5.

The orbital parameters of the centre of mass of the binary are given with a subscript ‘M’ and those of the external mass with a ‘3’, with respect to the centre of mass of the entire three-body configuration. In addition to the three frequency derivatives precisely measured for NGC 1851A, we also obtained constraints on the fourth and fifth derivatives. These derivatives, as purely induced by the external dynamics, after the removal of kinematic and intrinsic contributions, can be given as

$$\begin{array}{cc}& \dot{f}=-f\frac{\mathit{a}·{\widehat{\mathit{K}}}_0}{c},\hfill \end{array}$$(C.1)

$$\begin{array}{cc}& \ddot{f}=-f\frac{\dot{\mathit{a}}·{\widehat{\mathit{K}}}_0}{c},\hfill \end{array}$$(C.2)

$$\begin{array}{cc}& ⋮\hfill \\ \hfill & f^{(5)}=-f\frac{{\mathit{a}}^{(4)}·{\widehat{\mathit{K}}}_0}{c}.\hfill \end{array}$$(C.3)

The terms on the right can be explicitly expressed using the orbital parameters of the binary and the angular separations. Here we present these expressions, derived using the method outlined by JR97, and transformed to a function of ψ and $\dot{\psi }$. These equations are valid for both elliptical and hyperbolic orbits and can be solved to derive the unknown parameters for the third body.

The altered equations used in our code for obtaining the family of solutions presented earlier are

$$\begin{array}{cc}\hfill f^{(2)}& =\frac{B{\dot{\psi }}_M\dot{f}}{A^{2}sin{\psi }_M},\stackrel{⃛}{f}=\frac{C{{\dot{\psi }}_M}^2\dot{f}}{A^{2}sin{\psi }_M},\hfill \end{array}$$(C.4)

$$\begin{array}{cc}\hfill f^{(4)}& =\frac{D{{\dot{\psi }}_M}^3\dot{f}}{A^{2}sin{\psi }_M},f^{(5)}=\frac{E{{\dot{\psi }}_M}^4\dot{f}}{A^{2}sin{\psi }_M},\hfill \end{array}$$(C.5)

where a prime indicates a derivative with respect to ψ_M, and A, B, C, D, and E (Joshi & Rasio 1997; Perera et al. 2017) are defined as

$$\begin{array}{cc}\hfill A& =1+e_{3}cos{\lambda }_M,\hfill \end{array}$$(C.6)

$$\begin{array}{cc}\hfill B& =2A{A}^{\prime }sin{\psi }_M+A^{2}cos{\psi }_M,\hfill \end{array}$$(C.7)

$$\begin{array}{cc}\hfill C& =B^{\prime }+\frac{2B{A}^{\prime }}{A},D=C^{\prime }+\frac{4C{A}^{\prime }}{A},E=D^{\prime }+\frac{6D{A}^{\prime }}{A}.\hfill \end{array}$$(C.8)

We derived $\dot{\psi }$ using the equations C.4, and use it to derive ḟ, f⁽⁴⁾, and f⁽⁵⁾ from f⁽²⁾, using the equations

$$\begin{array}{cc}\hfill {\dot{\psi }}_M& =\left(\frac{\stackrel{⃛}{f}}{f^{(2)}}\right)\left(\frac{B}{C}\right),\hfill \end{array}$$(C.9)

$$\begin{array}{cc}\hfill \dot{f}& =\frac{A^2}{B{\dot{\psi }}_M}{f}^{(2)}sin{\psi }_M,f^{(4)}=\frac{\mathit{DB}}{{{\dot{\psi }}_M}^2},f^{(5)}=\frac{\mathit{EB}}{{{\dot{\psi }}_M}^3}.\hfill \end{array}$$(C.10)

The calculated values of the frequency derivatives are compared with the observed limits, and the allowed solutions for λ, ψ_M, and ${\dot{\psi }}_M$ are used to determine the parameters of the companion using equations (14)-(18) from JR97.

For an external third body in an elliptical orbit, the mass (m₃) and semi-major axis (a₃) are given by

$$\begin{array}{cc}& m_{3}sin{i}_3\approx -\frac{\dot{f}c}{fsin{\psi }_M}{\left[\frac{M_{\mathrm{tot}}^2{A}^2}{G{\dot{\psi }}_M^4}\right]}^{1/3}\hfill \end{array}$$(C.11)

$$\begin{array}{cc}& a_3=\left(\frac{M_{\mathrm{tot}}}{m_3}\right)a_M=\left(\frac{M_{\mathrm{tot}}}{m_3}\right)\frac{h}{1-e_3^2},\hfill \end{array}$$(C.12)

$$\begin{array}{cc}& h=-\frac{\dot{f}c{\mathrm{A}}^2}{fsin{i}_{3}sin{\psi }_M{{\dot{\psi }}_M}^2},\hfill \end{array}$$(C.13)

where h is the semi-latus rectum of the orbit of the binary’s centre of mass, and the semimajor axis of the relative orbit between M_tot and m₃ is a_R = a_M + a₃. The relative separation between the two (R), the orbital frequency and orbital period of the orbit of the external mass (n_b, M, P_b, M), and its (relative) velocity in this orbit (v) are given by (Roy 2005)

$$\begin{array}{cc}& R=\frac{a_R(1-e_3^2)}{1+e_{3}cos{\lambda }_M},\hfill \end{array}$$(C.14)

$$\begin{array}{cc}& n_{\mathrm{b},M}={\left[\frac{G(M_{\mathrm{tot}}+m_3)}{a_R^3}\right]}^{1/2},P_{\mathrm{b},M}=\frac{2\pi }{n_{\mathrm{b},M}},\hfill \end{array}$$(C.15)

$$\begin{array}{cc}& v={\left[G(M_{\mathrm{tot}}+m_3)\left(\frac{2}{R}-\frac{1}{a_R}\right)\right]}^{1/2}.\hfill \end{array}$$(C.16)

For the external body in a hyperbolic orbit, the mass m₃ can be derived using equation C.11; however, the semi-major axis a₃, velocity v and relative separation R of this body are given as (Roy 2005)

$$\begin{array}{cc}& a_3=\left(\frac{M_{\mathrm{tot}}}{m_3}\right)a_M=\left(\frac{M_{\mathrm{tot}}}{m_3}\right)\frac{h}{e_3^2-1},\hfill \end{array}$$(C.17)

$$\begin{array}{cc}& v={\left[G(M_{\mathrm{tot}}+m_3)\left(\frac{2}{R}+\frac{1}{a_R}\right)\right]}^{1/2}\hfill \end{array}$$(C.18)

$$\begin{array}{cc}& R=\frac{a_R(e_3^2-1)}{1+e_{3}cos{\lambda }_M}.\hfill \end{array}$$(C.19)

The asymptotic velocity of this body at infinity can be further derived as

$$\begin{array}{c}\hfill v_{\infty }={\left[\frac{G(M_{\mathrm{tot}}+m_3)}{a_R}\right]}^{1/2}.\end{array}$$(C.20)

Since the system is present in the core of the GC, we additionally ensure, for each trial, that the asymptotic velocity of the external mass, v_∞, is less than twice the central escape velocity of the cluster. As discussed earlier, we also put limits on ψ_M according to the eccentricity e₃ of the external mass’ trajectory.

Appendix D: Perturbations of the binary orbit by the third body

The presence of a distant third body leads to tidal perturbations of the pulsar binary. It is important to quantify these perturbations and estimate their relevance for the conclusions of the main part of the paper, most importantly the total mass of the system derived under the assumption that the observed $\dot{\omega }$ can be interpreted as the relativistic periastron advance of a binary system. Here we summarise the calculations needed to estimate the secular changes in the orbital parameters due to the perturbations by a distant third mass m₃. The orbital parameters for the orbit of this third body are taken from Sect. C.

Using the same convention as in Sect. 3, we describe the geometry of the binary orbit (relative motion r = r_p − r_c) by the following Keplerian parameters:

• a: semimajor axis of the relative binary motion

• e: orbital eccentricity

• i: orbital inclination with respect to the plane of the sky, i.e. the angle between the orbital angular momentum (direction $\widehat{\mathit{k}}$) and the line of sight from the observer to the pulsar binary (direction ${\widehat{\mathit{K}}}_0$)

• Ω: longitude of the ascending node (direction $\widehat{\mathit{i}}$), measured from a fixed direction in the plane of the sky, e.g. North (clockwise as seen by the observer).

• ω: angular distance of the periastron from the ascending node. In pulsar astronomy this angle is called the longitude of periastron.⁹

In the following, we used a coordinate system XYZ where the X-Y plane is the plane of the sky and Z is in the direction away from the observer (${\widehat{\mathit{K}}}_0$). Without loss of generality, we can assume the X-axis to be aligned with the ascending node of the binary orbit at the epoch t₀ (practically identical to T₀). Consequently, the unit vectors related to the orbital plane as defined by Damour & Taylor (1992) are given by

$$\begin{array}{c}\hfill \widehat{\mathit{i}}=(1,0,0),\widehat{\mathit{j}}=(0,cosi,sini),\widehat{\mathit{k}}=(0,-sini,cosi).\end{array}$$(D.1)

As a first approximation, which is sufficient for this paper, we can assume in the following that the perturbing third body has a fixed direction in space. The direction from the centre of mass of the binary to the third body at a distance R is given by the unit vector $\widehat{\mathit{N}}$. Using standard spherical coordinates (φ, ϑ) with respect to XYZ one has

$$\begin{array}{c}\hfill \widehat{\mathit{N}}=(sin{\vartheta }_{3}cos{\phi }_3,sin{\vartheta }_{3}sin{\phi }_3,cos{\vartheta }_3).\end{array}$$(D.2)

From Sect. C one has

$$\begin{array}{c}\hfill R=\frac{a_R(1-e_3^2)}{1+e_{3}cos{\lambda }_M},\end{array}$$(D.3)

where a_R is the semimajor axis of the relative orbit between the third body and the centre of mass of the binary system. Furthermore, one has (cf. Sect. C)

$$\begin{array}{c}\hfill {\vartheta }_3=arccos(N_Z)=arccos(sin{\psi }_{3}sin{i}_3).\end{array}$$(D.4)

where ψ₃ = ω₃ + λ_M. Note that the azimuth-angle φ₃ of the third body is unconstrained by the frequency derivatives (cf. Sect. C) and therefore needs to be treated as a free parameter.

In the following we need the direction cosines

$$\begin{array}{cc}& {\xi }_{\widehat{\mathit{i}}}:=\widehat{\mathit{N}}·\widehat{\mathit{i}}=sin{\vartheta }_{3}cos{\phi }_3,\hfill \end{array}$$(D.5)

$$\begin{array}{cc}& {\xi }_{\widehat{\mathit{j}}}:=\widehat{\mathit{N}}·\widehat{\mathit{j}}=cos{\vartheta }_{3}sini+sin{\vartheta }_{3}sin{\phi }_{3}cosi,\hfill \end{array}$$(D.6)

$$\begin{array}{cc}& {\xi }_{\widehat{\mathit{k}}}:=\widehat{\mathit{N}}·\widehat{\mathit{k}}=cos{\vartheta }_{3}cosi-sin{\vartheta }_{3}sin{\phi }_{3}sini.\hfill \end{array}$$(D.7)

The perturbation of a binary orbit by a distant third body can be calculated using the formalism of osculating elements (see for instance Sect. 3.4.1 in the book by Poisson & Will (2014)). Integrating over one full binary orbit one obtains the following shifts in the orbital parameters:

$$\begin{array}{cc}& \mathrm{\Delta }a=0,\hfill \end{array}$$(D.8)

$$\begin{array}{cc}& \mathrm{\Delta }e=\frac{15\pi \eta }{2}e\sqrt{1-e^2}\hfill \\ \hfill & \times [({\xi }_{\widehat{\mathit{i}}}^2-{\xi }_{\widehat{\mathit{j}}}^2)sin(2\omega )-2{\xi }_{\widehat{\mathit{i}}}{\xi }_{\widehat{\mathit{j}}}cos(2\omega )],\hfill \end{array}$$(D.9)

$$\begin{array}{cc}& \mathrm{\Delta }i=\frac{3\pi \eta }{2}\frac{{\xi }_{\widehat{\mathit{k}}}}{\sqrt{1-e^2}}\hfill \\ \hfill & \times [(2+3e^2){\xi }_{\widehat{\mathit{i}}}+5e^2{\xi }_{\widehat{\mathit{i}}}cos(2\omega )+5e^2{\xi }_{\widehat{\mathit{j}}}sin(2\omega )],\hfill \end{array}$$(D.10)

$$\begin{array}{cc}& \mathrm{\Delta }\mathrm{\Omega }=\frac{3\pi \eta }{2}\frac{{\xi }_{\widehat{\mathit{k}}}}{\sqrt{1-e^2}sini}\hfill \\ \hfill & \times [(2+3e^2){\xi }_{\widehat{\mathit{j}}}+5e^2{\xi }_{\widehat{\mathit{i}}}sin(2\omega )-5e^2{\xi }_{\widehat{\mathit{j}}}cos(2\omega )],\hfill \end{array}$$(D.11)

$$\begin{array}{cc}& \mathrm{\Delta }\omega =\frac{3\pi \eta }{2}\sqrt{1-e^2}\hfill \\ \hfill & \times [1-3{\xi }_{\widehat{\mathit{k}}}^2+5({\xi }_{\widehat{\mathit{i}}}^2-{\xi }_{\widehat{\mathit{j}}}^2)cos(2\omega )+10{\xi }_{\widehat{\mathit{i}}}{\xi }_{\widehat{\mathit{j}}}sin(2\omega )]\hfill \\ \hfill & -cosi\mathrm{\Delta }\mathrm{\Omega },\hfill \end{array}$$(D.12)

where

$$\begin{array}{c}\hfill \eta =\frac{m_3}{M_{\mathrm{tot}}}{\left(\frac{a}{R}\right)}^3\ll 1.\end{array}$$(D.13)

These equations agree with the ones obtained by Lidov (1962) (cf. his Eqs. 17).¹⁰ The corresponding (secular) rates of change of the orbital parameters are simply obtained by dividing the parameter shifts by the orbital period, for example $\dot{\omega }\equiv \mathrm{\Delta }\omega /P_{\mathrm{b}}$.

As equation (D.11) suggests, the secular change in the longitude of periastron depends on the orientation of the third body and the inclination, eccentricity, and longitude of periastron of the inner binary. For the purpose of this work, it is sufficient to estimate the maximum expected tidal perturbation to the observed rate of advance of periastron, ${\dot{\omega }}_{\mathrm{tidal},\max }$, over all possible directions of the external third body. Using the estimated orbital inclination of ∼63 deg for NGC 1851A from Table 2 and leaving η as a free parameter, we derive

$$\begin{array}{c}\hfill |\dot{\omega }{|}_{\mathrm{tidal},\max }\approx 8.1\times {10}^{-5}\left(\frac{m_3}{1{\mathrm{M}}_{\odot }}\right){\left(\frac{R}{100\mathrm{AU}}\right)}^{-3}\mathrm{deg}{\mathrm{yr}}^{-1}.\end{array}$$(D.14)

Similarly, from equation (D.10), using Δx = x coti Δi, one finds

$$\begin{array}{c}\hfill |\dot{x}{|}_{\mathrm{tidal},\max }\approx 1.8\times {10}^{-13}\left(\frac{m_3}{1{\mathrm{M}}_{\odot }}\right){\left(\frac{R}{100\mathrm{AU}}\right)}^{-3}.\end{array}$$(D.15)

On a final note, in the limit for small eccentricities, when neglecting terms 𝒪(e²), Eqs. (D.9) and (D.10) agree with the corresponding results obtained by Rasio (1994) and used by Joshi & Rasio (1997). This is not the case for Eq. (D.12). It appears that Rasio (1994) has overlooked the vertical (out of plane) tidal component (Eq. 4 in Rasio 1994) when calculating $\dot{\omega }$. Nevertheless, even for orbits with a low eccentricity, this contribution is relevant. It results (generally) in a tilting of the orbital plane, which leads to a change of the line of the nodes within the plane of the sky and, more importantly, within the orbital plane, which in turn contributes to $\dot{\omega }$, since ω is defined with respect to the ascending node. This effect gives rise to the ‘−cos i ΔΩ’ term in Eq. (D.12), which is absent in Eq. (5) of Rasio (1994).

Appendix E: Modifications to orbital models

The large line-of-sight jerk for a pulsar in a GC is a measure of the dynamical influence from an external body, and this affects both the spin period and the orbital period of the pulsar. This is evident in the resulting measurement of the large second period derivative and suggests considerable orbital variability in the system. In this scenario, the variations in the orbital period (P_b) can be modelled with multiple orbital frequency derivatives, and the orbital frequency f_b is the rate of change of the number of orbits. At any given time t, the number of orbits can be expressed in a Taylor expansion as

$$\begin{array}{c}\hfill N_{\mathrm{b}}={\displaystyle \sum _{k=0}^n}\frac{1}{(k+1)!}{f}_b^{(k)}{(t-T_0)}^{k+1},\end{array}$$(E.1)

where n is the number of orbital frequency derivatives, and f_b^(k) is the k^th orbital frequency derivative at the reference epoch of periastron passage T₀.

However, if we’re confident that the ‘intrinsic’ contributions to these higher spin and orbital frequency derivatives are negligible (for instance, when the expected contributions to $\ddot{f}$ and $\stackrel{⃛}{f}$ from timing noise are much smaller than the observations), the orbital frequency derivatives can be obtained from the higher order spin frequency derivatives according to eq. 12. Replacing those in equation (E.1), and using the P_b and Ṗ_b used in those models (PB and PBDOT), we obtain:

$$\begin{array}{cc}& N_{\mathrm{b}}=\frac{1}{P_{\mathrm{b}}}[(t-T_0)-\frac{1}{2}\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}{(t-T_0)}^2\hfill \\ \hfill & +\frac{1}{6}\frac{\ddot{f}}{f}{(t-T_0)}^3+\frac{1}{24}\frac{\stackrel{⃛}{f}}{f}{(t-T_0)}^4+\cdots ],\hfill \end{array}$$(E.2)

where the reference epoch for the measurement of the ḟ^(b) must be close to T₀. We thus modified the DD, DDK, DDH, and DDGR orbital models, where the second part of the equation is used if a flag (JERK) is set to 1 and any higher-order spin frequency derivatives are being fit. These models then produce self-consistent values for f⁽ⁿ⁾ (with n > 1), where the effect of the jerk is taken into account not only for the spin phase but also for the orbital phase. If the JERK flag is set to 0, or not specified at all, then only the first part of the equation is used, as in the current implementation.

All Tables

All Figures

Fig. 1. Top to bottom: Post-fit timing residuals for NGC 1851A obtained with the DDH timing model, as a function of (a) epoch and (b) orbital phase. The orbital phase 0 denotes periastron and the residual 1-σ uncertainties are indicated by vertical error bars. The various colours denote the different telescopes and backends used for the timing. Orange indicates the early GBT data, and purple and red indicate the uGMRT data taken in the PA and CDP modes respectively. Grey and green represent the data taken with the L-band and the UHF receivers of the MeerKAT.

In the text

Fig. 2. Orbital inclination and mass constraints for NGC 1851A. The solid and dotted lines represent the nominal values and 1-σ constraints for the post-Keplerian parameters (red: $\dot{\omega }$, cyan: γ, purple: h3) obtained with the DDH binary model. Bottom left plot: cos i versus the companion mass (Mc); Bottom right plot: Pulsar mass (Mp) versus Mc. The grey portions on the two plots indicate the regions excluded for Mp < 0 and sin i > 1 respectively. The normalised 1-D pdfs for cos i, Mp and Mc are shown in the top left, top right, and right panels.

In the text

Fig. 3. Top to bottom: Evolution of the observed spin frequency (f), its first derivative (ḟ), and its second derivative ($\ddot{f}$) as a function of time (MJD). The vertical yellow and green lines indicate Epochs 1 and 2 respectively. The values from R19 are shown in red, and the results from this work are shown in blue. The large change in the ḟ value between the two epochs can be clearly seen.

In the text

	Fig. 4. Closest distance of approach and probability, for a nearby star of mass M in the globular cluster, to induce the observed jerk in the NGC 1851A system.
In the text

	Fig. 5. Definition of the orbital elements for a three-body model with two Keplerian orbits. The ‘inner binary’ is represented by Mtot and the external mass by m3, and both the scenarios of an elliptical orbit and a hyperbolic fly-by encounter are shown.
In the text

Fig. 6. Family of solutions for the mass m3, orbital period Pb, 3, and relative distance R for a second companion of NGC 1851A using the limits on all the frequency derivatives mentioned in the text. The different eccentricities are shown in different colours. The constant lines of ${\dot{\omega }}_{\mathrm{tidal}}/{\dot{\omega }}_{\mathrm{obs}}$ are shown with dotted black lines in the bottom plot, while the dash-dotted crimson lines indicate constant ẋtidal.

In the text

	Fig. 7. Allowed solutions showing the mass m3, asymptotic velocity v∞, and relative distance R for a nearby star in a fly-by encounter, in the core of the globular cluster. Constant ${\dot{\omega }}_{\mathrm{tidal}}/{\dot{\omega }}_{\mathrm{obs}}$ and constant ẋtidal are shown with black dotted and crimson dash-dotted lines respectively in the bottom plot.
In the text

Fig. 8. Location of NGC 1851A, with the blue circles corresponding to radius of 0$\stackrel{″}{.}$069 around the position of the pulsar as estimated for the time the left frame was taken (MJD = 55558). Given the small proper motion, correcting for it does not induce any noticeable changes in the position of the pulsar. Images from the Wide Field Camera 3, taken with the F275W filter (left) and F814W filter (right).

In the text

Fig. 9. Central panel: 2-D pdf for the full Ω − cos i space of binary pulsars. The blue and green contours show the probability density for inclination i > 90 deg and i < 90 deg, where increasing probability density is shown with darker shades. The grey regions are excluded based on the requirement that the pulsar mass must be greater than 0, and the dotted red line indicates the position angle of the proper motion of the system. Top panel: Probability density function for cosi, normalised to the maximum. Right-hand panel: Probability density function for Ω.

In the text

	Fig. B.1. Top to bottom: Cumulative weighted TOA count, observed rate of advance of periastron ${\dot{\omega }}_{\mathrm{obs}}$, Einstein delay γ, and observed orbital period derivative Ṗb,obs as a function of time. The values and their residual uncertainties (given by the vertical errorbars) are as measured using the DDH timing model presented in Table 2.
In the text