Some compact objects observed in gravitational wave events have masses in the gap between known neutron stars NSs and black holes BHs. The nature of these mass gap objects is unknown as is the formation of their host binary systems. We report pulsar timing observations made with the Karoo Array Telescope MeerKAT of PSR J05144002E an eccentric binary millisecond pulsar in the globular cluster NGC 1851. We found a total binary mass of 3.887  0.004 solar masses MSUBSUB and multiwavelength observations show that the pulsars binary companion is also a compact object. The companions mass 2.09 to 2.71 MSUBSUB 95 confidence interval is in the mass gap indicating either a very massive NS or a lowmass BH. We propose that the companion formed in a merger between two earlier NSs.
['Astrophysics - High Energy Astrophysical Phenomena', 'General Relativity and Quantum Cosmology', 'Nuclear Theory']
A pulsar in a binary with a compact object in the mass gap between neutron stars and black holes
{'No Header': "A pulsar in a binary with a compact object in the mass gap between neutron stars and black holes  \n4  \nEwan D. Barr, 1 ∗† Arunima Dutta, 1 ∗† Paulo C. C. Freire, 1 Mario Cadelano, 3 , 4 Tasha Gautam, 1 Michael Kramer, 1 Cristina Pallanca, 3 , 4 Scott M. Ransom, 9 Alessandro Ridolfi, 1 , 7 Benjamin W. Stappers, 2 Thomas M. Tauris, 1 , 13 Vivek Venkatraman Krishnan, 1 Norbert Wex, 1 Matthew Bailes, 5 , 6 Jan Behrend, 1 Sarah Buchner, 8 Marta Burgay, 7 Weiwei Chen, 1 David J. Champion, 1 C. -H. Rosie Chen, 1 Alessandro Corongiu, 7 Marisa Geyer, 8 , 12 ‡ Y. P. Men, 1 Prajwal V. Padmanabh, 1 , 10 , 11 Andrea Possenti 7  \n- 1 Max-Planck-Institut fur Radioastronomie, 53121 Bonn, Germany\n- 2 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,\n- The University of Manchester, Manchester, M13 9PL, UK\n- Dipartimento di Fisica e Astronomia 'Augusto Righi', Universit'a degli Studi di Bologna, 40129 Bologna, Italy\n- 3  \nIstituto nazionale di astrofisica-Osservatorio di Astrofisica e Scienze dello Spazio di Bologna,  \nI-40129 Bologna, Italy  \n- 5 Centre for Astrophysics and Supercomputing, Swinburne University of Technology,  \nPO Box 218 Hawthorn, Vic, 3122, Australia  \n- 6 Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav)\n- 7 Istituto nazionale di astrofisica-Osservatorio Astronomico di Cagliari, I-09047, Selargius, Italy\n- 8 South African Radio Astronomy Observatory, Observatory 7925, South Africa\n- 9 National Radio Astronomy Observatory, Charlottesville, VA, 22903, USA\n- 10 Max Planck Institute for Gravitational Physics (Albert Einstein Institute),  \nD-30167 Hannover, Germany  \n- 11 Leibniz Universitat Hannover, D-30167 Hannover, Germany\n- 12 Department of Astronomy, University of Cape Town, Rondebosch, Cape Town, 7700,", 'South Africa': '13 Department of Materials and Production, Aalborg University  \nDK-9220 Aalborg Øst, Denmark  \n∗ Corresponding author;  \nE-mail: ebarr@mpifr-bonn.mpg.de, adutta@mpifr-bonn.mpg.de  \n† These authors contributed equally to this work.  \n‡ Present address: High Energy Physics, Cosmology & Astrophysics Theory (HEPCAT) Group,  \nDepartment of Mathematics & Applied Mathematics, University of Cape Town,  \nCape Town 7700, South Africa  \nAmong the compact objects observed in gravitational wave merger events a few have masses in the gap between the most massive neutron stars (NSs) and least massive black holes (BHs) known. Their nature and the formation of their merging binaries are not well understood. We report on pulsar timing observations using the Karoo Array Telescope (MeerKAT) of PSR J0514 -4002E, an eccentric binary millisecond pulsar in the globular cluster NGC 1851 with a total binary mass of 3 . 887 ± 0 . 004 solar masses ( M ⊙ ). The companion to the pulsar is a compact object and its mass (between 2 . 09 and 2 . 71 M ⊙ , 95% confidence interval) is in the mass gap, so it either is a very massive NS or a low-mass BH. We propose the companion was formed by a merger between two earlier NSs.  \nGlobular clusters (GCs) are dense, gravitationally bound, stellar clusters. They have been observed to host a large number of low-mass X-ray binaries (LMXBs) consisting of a compact object accreting material from a donor star. LMXBs are ∼ 10 3 times more abundant per unit of stellar mass in GCs than in the disk of the Milky Way galaxy (Galactic plane) (1). This is due to the high stellar densities at the centre of GCs, which increase the rate of exchange encounters in which neutron stars (NSs) acquire low-mass main sequence (MS) companions. The MS stars evolve until they start transferring mass to the NS, at which point an LMXB is formed.  \nThese X-ray binaries are expected to produce millisecond pulsars (MSPs, radio-emitting neutron stars with spin periods P < 10 ms) in almost circular orbits around low-mass companions(2, 3). There are a total of 305 pulsars known in 40 GCs(4), the vast majority of which are MSPs. Most of the systems in GCs are similar to the MSP population found in the Galactic plane, although their orbital eccentricities are often higher, which is thought to be a result of close encounters with other stars(5).  \nIn GCs with the densest cores, any particular star - or MSP - is likely to experience multiple  \nexchange encounters over its lifetime(6). A possible outcome is the exchange of a low-mass companion of an MSP for either a massive white dwarf (WD) or another NS, resulting in a massive, eccentric MSP binary (7-10). Observing such systems allows their component masses to be measured and can test theories of gravity(11). The same process could also produce a MSP-black hole (BH) system (12).', 'The millisecond pulsar binary PSR J0514 -4002E': "A survey searching for MSPs in GCs (10) has been carried out using the MeerKAT radio telescope array in South Africa (13, 14). The results of the survey (15) included 13 MSPs in NGC 1851, a GC located in the Southern constellation of Columba(16). These include three massive, eccentric MSP binaries: PSR J0514 -4002A (8, 17), PSR J0514 -4002D and PSR J0514 -4002E(16).  \nThe latter has a spin period ( P ) of 5 . 6 ms, an orbital period ( P b ) of 7 . 44 days and an orbital eccentricity ( e ) of 0 . 71 (16). The projected semi-major axis of the pulsar's orbit ( x ≡ a p sin i/c , where a p is the semi-major axis of the pulsar orbit, i the orbital inclination and c is the speed of light in vacuum) is 27.8 seconds (16). The mass function ( f ) is thus:  \nf ( m p , m c ) ≡ ( m c sin i ) 3 ( m p + m c ) 2 = 4 π 2 c 3 G x 3 P b 2 = 0 . 41672 ± 0 . 00022 M ⊙ , (1)  \nwhere m c is the mass of the companion, m p is the mass of the pulsar and G is the gravitational constant. The unit M ⊙ is the mass of the Sun; we use the nominal solar mass adopted by the International Astronomical Union (18). All uncertainties are a confidence interval (C. I.) corresponding to a 68.3% confidence level, unless otherwise stated. Assuming m p ≥ 1 . 17 M ⊙ (corresponding to the lowest NS mass measured (19) and a theoretical lower limit (20)) and an edge-on orbit ( i = 90 · ), the mass function alone indicates m c ≥ 1 . 40 M ⊙ .", 'Radio timing observations': "To determine the spin, astrometric and orbital parameters for all the pulsars in NGC 1851, we conducted 24 observations of this globular cluster using MeerKAT. The observations used either the L-band (856 to 1712 MHz) or Ultra High Frequency (UHF, 544 to 1088 MHz) receivers(10) and were performed between January 2021 and August 2022. Data acquisition and initial reduction were performed using the Pulsar Timing User Supplied Equipment (PTUSE) instrument(21). We analysed the resulting times of arrival (ToAs) of the pulsed signal to determine an initial phase-coherent timing model for PSR J0514 -4002E(22). Using this model, we recovered the previously undetected signals from this pulsar in six archival observations of NGC1851 made with the 800 MHz (795 to 845 MHz) and S-band (1.73 to 2.60 GHz) receivers on the Robert C. Byrd Green Bank Telescope (GBT) (23) between December 2005 and August 2006. We combined the MeerKAT and GBT ToAs (22) and re-fitted the timing parameters to determine a refined timing model. The results are listed in Table 1 and the fitting residuals between this model and the observed ToAs are shown in Figure 1.  \nThe timing model includes a precise measurement of the binary's rate of periastron advance, ˙ ω = 0 . 03468 ± 0 . 00003 · yr -1 . We obtain a consistent value (but slightly higher uncertainty) if we consider the MeerKAT data alone (Figure S2). As discussed in detail below, this effect is relativistic (with additional contributions being small relative to the measurement uncertainty (see Supplementary Text)) and its magnitude implies a high total system mass. We therefore performed an additional dense observing campaign (22) to search for Shapiro delay, a relativistic light propagation delay in the system (24). The longer time span provided by the GBT data also enabled us to search for the Einstein delay, an effect caused by the varying time dilation experienced by the pulsar at different orbital phases. We do not detect either the Shapiro delay or Einstein delay in our timing data, setting 95% confidence upper limits of h 3 < 1 . 48 µs and  \nFigure 1: Timing model fitting residuals for PSR J0514 -4002E . Shown are the residuals between the observed ToAs and the timing model presented in Table 1 as a function of observing epoch (A) , in units of Modified Julian Date (MJD), and orbital phase (B) . An orbital phase of 0 corresponds to periastron, and superior conjunction occurs at an orbital phase of 0.008 (shown with the brown dashed line). The vertical error bars indicate the 1σ uncertainties. The green points indicate data taken with the 800 MHz receiver on the GBT. Orange and blue points indicate data taken with the L-Band and UHF receivers of MeerKAT, respectively.  \n<!-- image -->  \nγ E < 25 ms, where h 3 is the orthometric amplitude of the Shapiro delay (25) and γ E is the Einstein delay parameter. These non-detections further constrain the companion mass, as discussed below.  \nTable 1: Timing model for PSR J0514 -4002E. The model was derived from the MeerKAT and GBT data with the single GBT observation at S-band excluded during fitting (22). Because the proper motion and parallax cannot be measured from the timing data, we adopt the bulk proper motion and parallax for NGC 1851 derived from HST and Gaia observations(26). Values are reported in the Dynamic Barycentric Time timescale and uncertainties are 68.3% C. I.  \n- ∗ Assuming the observed ˙ ω is due to relativistic effects (see Supplementary Text)", 'Near-ultraviolet and optical observations': 'If the ≥ 1 . 40 M ⊙ companion of PSR J0514 -4002E were a main-sequence star, it should be detectable at optical wavelengths. We searched for an optical source using archival Hubble Space Telescope (HST) observations with the Wide Field Camera 3 (WFC3) in the F275W and F336W filters. This filter combination is particularly sensitive to blue stars such as blue straggler stars (BSSs) and WDs (27, 28). Bright BSSs are common in GCs, formed through stellar collisions or mass-transfer in a binary system; they can in principle have masses compatible with those predicted for the companion of PSR J0514 -4002E.  \nNo optical source is detected at the position of PSR J0514 -4002E (Figures 2A and 2B). The closest stellar source to the pulsar position is a star with a color-magnitude position consistent with those of BSSs (Figure 2C) and is offset by 90 milliarcseconds (mas), which is more than 6 times the astrometric precision (14 mas). It excludes a physical association, because the orbital separation between the pulsar and companion calculated from our timing model and located at the cluster distance of 11.66 kpc (26) is < 10 -3 mas. Assuming the optical source is a hydrogenburning star, we estimate its mass as ∼ 1 . 2 M ⊙ (22), which is lower than the lower limit on m c from the mass function. Another star, located 100 mas from the pulsar position, is a red-giant. Red-giant branch stars in old GCs such as NGC 1851 have masses of about 0 . 7 -0 . 8 M ⊙ (29), also lower than the minimum companion mass. We therefore conclude that the companion of PSR J0514 -4002E is not detectable in the HST images.', 'System mass': "The non-detection of the companion of PSR J0514 -4002E in the HST images implies that it must be a compact object. There is no measurable excess pulse dispersion near superior conjunction, which can be produced by ionized gas emanating from either a main sequence or a  \n<!-- image -->  \nFigure 2: Stars near the position of PSR J0514 -4002E. (A) HST image in the F275W filter of the 4 '' × 4 '' region surrounding the radio timing position of PSR J0514 -4002E (green cross). The blue and red circles indicate two stars discussed in the text; the circle sizes are twice the astrometric uncertainty. (B) A 0 . 3 '' × 0 . 3 '' zoom of the green box of panel A. The radio timing position is indicated with a green ellipse with size and orientation equal to the 95% C. I. on the pulsar's position relative to the ICRS. (C) m F275W , m F275W -m F336W color-magnitude diagram of NGC 1851 derived from the HST images. The black dots represent all the cluster stars within the sampled field of view. The blue and red circles indicate the stars marked in panels A-B, which we interpret as a BSS and red giant, respectively. The photometric uncertainty for both the magnitude and color of the two marked stars is smaller than the symbol size (22).  \n<!-- image -->  \ngiant star companion(22). We therefore interpret the measured rate of advance of periastron ( ˙ ω ) as being purely of relativistic origin, with negligible contribution from the spins of the binary components (see Supplementary Text). Assuming general relativity (GR), the total mass of the system is then(30,31):  \nM ≡ m p + m c = c 3 G [ ˙ ω 3 (1 -e 2 ) ] 3 / 2 ( P b 2 π ) 5 / 2 = 3 . 887 ± 0 . 004 M ⊙ . (2)  \nThis value is 1 . 0 M ⊙ larger than the mass of the most massive double neutron star (DNS) known in the Milky Way, PSR J1913+1102, which is 2 . 8887 ± 0 . 0006 M ⊙ (32). It is also larger than the total mass of the heaviest DNS merger detected in gravitational waves, GW190425, at > 99 . 5%  \nprobability [(33), their figure 5].", 'Nature of the companion': 'Combining the measurement of the total mass with the mass function, we obtain, for an edge-on orbit ( i = 90 · ), m p ≤ 2 . 04 M ⊙ and m c ≥ 1 . 84 M ⊙ (Figure 3). This companion mass is far too high for a WD, as the upper mass limit for a rigidly rotating WD is about 1 . 47 M ⊙ (34). Smaller inclination angles ( i < 90 · ) would imply a smaller m p and larger m c . Adopting the minimum neutron star mass discussed above, m p ≥ 1 . 17 M ⊙ , we set limits of i ≥ 42 . 9 · and m c ≤ 2 . 71 M ⊙ .  \nWe do not detect additional relativistic effects that could allow the individual masses to be determined. However, our upper limits on the Shapiro delay and the Einstein delay provide additional constraints on the masses and orbital inclination. To quantify these, we have made a Bayesian estimate of the component masses based on the quality of fit to the observational data ( χ 2 ) over a grid of total mass and orbital inclination values (22). The variation of the orbital period is contaminated by the acceleration of the system in the cluster. We assume that all other relativistic effects are as predicted by GR and that m p ≥ 1 . 17 M ⊙ as above (22). We find Bayesian posteriors of M = 3 . 887 ± 0 . 004 M ⊙ , m p = 1 . 53 +0 . 18 -0 . 20 M ⊙ , m c = 2 . 35 +0 . 20 -0 . 18 M ⊙ and i = 52 +6 · -5 (median values with 68.3% confidence limits (22)). The 95% probability limits are i < 62 · , m p < 1 . 79 M ⊙ and therefore m c > 2 . 09 M ⊙ .  \nThe companion mass is therefore likely to be in the mass gap for compact objects (37), being higher than the largest precisely measured pulsar masses, m p = 2 . 08 ± 0 . 07 M ⊙ for PSR J0740 + 6620 (38) and m p = 2 . 01 ± 0 . 04 M ⊙ for PSR J0348 + 0432 (39), while simultaneously below the observed minimum mass of about 5 M ⊙ for BHs in Galactic X-ray binaries (40,41).  \nIf the companion were a massive NS, it might also be a radio pulsar. We searched for radio  \nFigure 3: Companion mass of PSR J0514 -4002E. (A) Companion mass as a function of orbital inclination for the PSR J0514 -4002E system. The solid red curve shows solutions within the 95% C. I. (see text); the dotted part indicates lower and higher masses disregarded by our Bayesian model and assuming m p ≥ 1 . 17 M ⊙ , respectively. The grey-shaded region is ruled out because of the mass function and the total mass (Eqs. 1 and 2). Depending on the (unknown) NS equation-of-state (35,36), the light-blue shaded mass band corresponds to either massive NSs or light BHs. (B) Inferred companion mass of PSR J0514 -4002E (red) compared to the largest observed masses of radio pulsars (blue), low-mass components of gravitationalwave mergers (black), and the total post-merger remnant mass of GW170817 (grey, assuming no energy and mass loss after the inspiral, i.e. an absolute upper limit). Source names followed by (C) or (M) refer to the mass of the companion star and the remnant mass of the merger product, respectively; see Table S2 for references, masses and error bars.  \n<!-- image -->  \npulsations from the companion, assuming the full allowed range of mass ratios, but did not detect any (22). We therefore cannot determine whether the companion is a massive NS or a  \nlow-mass BH.', 'Formation of the system': "The combination of the location in a dense GC (where stellar exchange encounters often occur, see above) the highly eccentric orbit, the fast spin of the pulsar and the large companion mass indicates that the PSR J0514 -4002E system is the product of a secondary exchange encounter. We propose that an earlier low-mass companion transferred mass to this pulsar, increasing the spin rate, before being replaced by the present high-mass companion in an exchange encounter. However, a more complicated evolution with multiple exchange encounters is also possible. We therefore cannot infer the nature of the companion from binary evolution models.  \nIf the mass of the primary in PSR J0514 -4002E is in the range 1.25 to 1.55 M ⊙ of the four measurements of pulsar masses in GCs (9, 11, 17, 42), then the corresponding value of m c (between 2 . 34 and 2 . 63 M ⊙ ) overlaps with the range of masses of remnants from mergers of DNSs, such as the merger product of GW170817 (Figure 3). We suggest that, prior to becoming part of the current PSR J0514 -4002E system, the companion could potentially have formed in such a merger event, regardless of whether it is a NS or BH. Although the probability of DNS mergers is generally low in GCs that have not undergone core collapse(43), NGC 1851 has a dense core (see Supplementary Text) which makes a DNS merger in the progenitor of the PSR J0514 -4002E system more probable. A DNS has been observed in a GC (PSR B2127 + 11C, in M15) with a calculated merger time of 217 Myr (11), implying that merger remnants are likely to be present in GCs.  \nOur derived companion mass overlaps with the mass estimates for the lighter components of the BH+BH or BH+NS merger candidates GW190814 (44), GW190917 and GW200210 (45). The lighter component of GW190814 has previously been interpreted as the product of an earlier merger (46) that later acquired a more massive BH companion via exchange encounters  \nin GCs (and then merged in the GW190814 event).  \nIf the companion of PSR J0514 -4002E is a light BH formed in such a merger, it would acquire a spin parameter χ c ∈ [0 . 6 , 0 . 875] during the merger (47), where χ c is the dimensionless BH spin angular momentum. A NS rotating at the maximum theoretical rate would have a similar χ c immediately after merger (48), though we expect this would decrease rapidly after formation due to electromagnetic torque. Assuming a magnetic field of 10 9 G , the spin parameter would become ≲ 0 . 3 (corresponding to the fastest known MSPs) after ∼ 30 Myr , so we regard a fast-spinning NS companion as unlikely.  \nABHcompanion with χ c ∈ [0 . 6 , 0 . 875] would induce relativistic spin-orbit coupling, causing the orbital plane to precess around the total angular momentum vector, an effect known as Lense-Thirring precession. We calculate the resulting variation of the projected semi-major axis of the pulsar's orbit ( ˙ x ) would be ≲ 1 . 7 × 10 -13 . This is slightly beyond the effect size that would be detectable in our data; the 1σ uncertainty on the measurement of ˙ x is 2 . 0 × 10 -13 (see Supplementary Text). We therefore cannot differentiate between a NS and a BH companion.", 'References and Notes': "- 1. G. W. Clark, Astrophys. J. Lett. 199 , L143 (1975).\n- 2. M. A. Alpar, A. F. Cheng, M. A. Ruderman, J. Shaham, Nature 300 , 728 (1982).\n- 3. T. M. Tauris, E. P. J. van den Heuvel, Physics of Binary Star Evolution. From Stars to X-ray Binaries and Gravitational Wave Sources (Princeton University Press, 2023).\n- 4. P. C. C. Freire, Pulsars in Globular Clusters, https://www3.mpifr-bonn.mpg. de/staff/pfreire/GCpsr.html . [Online; accessed 10-Nov-2023].  \n- 5. E. S. Phinney, Philosophical Transactions of the Royal Society of London Series A 341 , 39 (1992).\n- 6. F. Verbunt, P. C. C. Freire, Astronomy & Astrophysics 561 , A11 (2014).\n- 7. T. A. Prince, S. B. Anderson, S. R. Kulkarni, A. Wolszczan, Astrophys. J. Lett. 374 , L41 (1991).\n- 8. P. C. Freire, Y. Gupta, S. M. Ransom, C. H. Ishwara-Chandra, Astrophys. J. Lett. 606 , L53 (2004).\n- 9. R. S. Lynch, P. C. C. Freire, S. M. Ransom, B. A. Jacoby, Astrophys. J. 745 , 109 (2012).\n- 10. A. Ridolfi, et al. , Month. Not. R. Astron. Soc. 504 , 1407 (2021).\n- 11. B. A. Jacoby, et al. , Astrophys. J. Lett. 644 , L113 (2006).\n- 12. K. Kremer, et al. , Phys. Rev. Lett. 120 , 191103 (2018).\n- 13. J. Jonas, MeerKAT Team, 'The MeerKAT Radio Telescope' in Proceedings of MeerKAT Science: On the Pathway to the SKA, Stellenbosch, 25-27 May 2016 (PoS, 2016), p. 1.\n- 14. F. Camilo, Nature Astronomy 2 , 594 (2018).\n- 15. E. D. Barr, TRAPUM Discoveries, http://trapum.org/discoveries/ . [Online; accessed 10-Nov-2023].\n- 16. A. Ridolfi, et al. , Astronomy & Astrophysics 664 , A27 (2022).\n- 17. A. Ridolfi, P. C. C. Freire, Y. Gupta, S. M. Ransom, Month. Not. R. Astron. Soc. 490 , 3860 (2019).\n- 18. A. Prˇsa, et al. , Astronom. J. 152 , 41 (2016).  \n- 19. J. G. Martinez, et al. , Astrophys. J. 812 , 143 (2015).\n- 20. Y. Suwa, T. Yoshida, M. Shibata, H. Umeda, K. Takahashi, Month. Not. R. Astron. Soc. 481 , 3305 (2018).\n- 21. M. Bailes, et al. , Publications Astron. Soc. Australia 37 , e028 (2020).\n- 22. Materials and methods are available as supplementary materials.\n- 23. P. C. C. Freire, S. M. Ransom, Y. Gupta, Astrophys. J. 662 , 1177 (2007).\n- 24. I. I. Shapiro, Phys. Rev. Lett. 13 , 789 (1964).\n- 25. P. C. C. Freire, N. Wex, Month. Not. R. Astron. Soc. 409 , 199 (2010).\n- 26. M. Libralato, et al. , Astrophys. J. 934 , 150 (2022).\n- 27. S. Raso, et al. , Astrophys. J. 895 , 15 (2020).\n- 28. J. Chen, et al. , Nature Astronomy 5 , 1170 (2021).\n- 29. D. Valcin, J. L. Bernal, R. Jimenez, L. Verde, B. D. Wandelt, Journal of Cosmology and Astroparticle Physics 2020 , 002 (2020).\n- 30. H. P. Robertson, Ann. Math. 39 , 101 (1938).\n- 31. J. H. Taylor, J. M. Weisberg, Astrophys. J. 253 , 908 (1982).\n- 32. R. D. Ferdman, et al. , Nature 583 , 211 (2020).\n- 33. B. P. Abbott, et al. , Astrophys. J. Lett. 892 , L3 (2020).\n- 34. S. C. Yoon, N. Langer, Astronomy & Astrophysics 435 , 967 (2005).\n- 35. M. Shibata, E. Zhou, K. Kiuchi, S. Fujibayashi, Phys. Rev. D 100 , 023015 (2019).  \n- 36. E. R. Most, L. J. Papenfort, L. R. Weih, L. Rezzolla, Month. Not. R. Astron. Soc. 499 , L82 (2020).\n- 37. T. B. Littenberg, B. Farr, S. Coughlin, V. Kalogera, D. E. Holz, Astrophys. J. Lett. 807 , L24 (2015).\n- 38. E. Fonseca, et al. , Astrophys. J. Lett. 915 , L12 (2021).\n- 39. J. Antoniadis, et al. , Science 340 , 448 (2013).\n- 40. C. D. Bailyn, R. K. Jain, P. Coppi, J. A. Orosz, Astrophys. J. 499 , 367 (1998).\n- 41. F. Ozel, D. Psaltis, R. Narayan, J. E. McClintock, Astrophys. J. 725 , 1918 (2010).\n- 42. A. Corongiu, et al. , Astronomy & Astrophysics 671 , A72 (2023).\n- 43. C. S. Ye, et al. , Astrophys. J. Lett. 888 , L10 (2020).\n- 44. R. Abbott, et al. , Astrophys. J. Lett. 896 , L44 (2020).\n- 45. The LIGO Scientific Collaboration, et al. , arXiv e-prints p. arXiv:2111.03606 (2021).\n- 46. W. Lu, P. Beniamini, C. Bonnerot, Month. Not. R. Astron. Soc. 500 , 1817 (2021).\n- 47. S. Bernuzzi, General Relativity and Gravitation 52 , 108 (2020).\n- 48. K.-W. Lo, L.-M. Lin, Astrophys. J. 728 , 12 (2011).\n- 49. arunimad10, ewanbarr, ewanbarr/NGC1851E Additional Materials, Zenodo [dataset], DOI:10.5281/zenodo.7991220.\n- 50. adutta, d/a59d968eccf148239f1d/, Keeper (MPDL) [dataset], https://doi.org/10.17617/4.e3.  \n- 51. F. Abbate, et al. , Month. Not. R. Astron. Soc. 513 , 2292 (2022).\n- 52. W. van Straten, M. Bailes, Publications Astron. Soc. Pacific 28 , 1 (2011).\n- 53. D. L. Kaplan, et al. , Publications Astron. Soc. Pacific 117 , 643 (2005).\n- 54. W. van Straten, M. Bailes, DSPSR: Digital Signal Processing Software for Pulsar Astronomy, Astrophysics Source Code Library, record ascl:1010.006 (2010).\n- 55. A. W. Hotan, W. van Straten, R. N. Manchester, Publications of the Astronomical Society of Australia 21 , 302 (2004).\n- 56. W. van Straten, et al. , PSRCHIVE: Development Library for the Analysis of Pulsar Astronomical Data, Astrophysics Source Code Library, record ascl:1105.014 (2011).\n- 57. S. M. Ransom, New search techniques for binary pulsars, Ph.D. thesis, Harvard University, Massachusetts, https://www.cv.nrao.edu/ ˜ sransom/ransom\\_ thesis\\_2001.pdf (2001).\n- 58. S. Ransom, PRESTO: PulsaR Exploration and Search TOolkit, Astrophysics Source Code Library, record ascl:1107.017 (2011).\n- 59. D. Nice, et al. , Tempo: Pulsar timing data analysis, Astrophysics Source Code Library, record ascl:1509.002 (2015).\n- 60. P. C. C. Freire, A. Ridolfi, Month. Not. R. Astron. Soc. 476 , 4794 (2018).\n- 61. Y. J. Guo, et al. , Astronomy & Astrophysics 654 , A16 (2021).\n- 62. R. S. Park, W. M. Folkner, J. G. Williams, D. H. Boggs, Astronom. J. 161 , 105 (2021).\n- 63. T. Damour, N. Deruelle, Ann. Inst. Henri Poincar'e Phys. Th'eor 44 , 263 (1986).  \n- 64. T. Damour, J. H. Taylor, Phys. Rev. D 45 , 1840 (1992).\n- 65. H. Baumgardt, M. Hilker, Month. Not. R. Astron. Soc. 478 , 1520 (2018).\n- 66. P. C. C. Freire, et al. , Month. Not. R. Astron. Soc. 471 , 857 (2017).\n- 67. B. J. Prager, et al. , Astrophys. J. 845 , 148 (2017).\n- 68. J. M. Weisberg, Y. Huang, Astrophys. J. 829 , 55 (2016).\n- 69. M. Kramer, et al. , Phys. Rev. X 11 , 041050 (2021).\n- 70. E. M. Splaver, et al. , Astrophys. J. 581 , 509 (2002).\n- 71. P. P. Eggleton, Astrophys. J. 268 , 368 (1983).\n- 72. A. Ridolfi, PYSOLATOR: Remove orbital modulation from a binary pulsar and/or its companion, Astrophysics Source Code Library, record ascl:2003.012 (2020).\n- 73. P. B. Stetson, Publications Astron. Soc. Pacific 99 , 191 (1987).\n- 74. M. Cadelano, et al. , Astrophys. J. 875 , 25 (2019).\n- 75. M. Cadelano, et al. , Astrophys. J. 905 , 63 (2020).\n- 76. P. B. Stetson, Publications Astron. Soc. Pacific 106 , 250 (1994).\n- 77. G. Piotto, et al. , Astronom. J. 149 , 91 (2015).\n- 78. A. Bellini, J. Anderson, L. R. Bedin, Publications Astron. Soc. Pacific 123 , 622 (2011).\n- 79. Gaia Collaboration, et al. , Astronomy & Astrophysics 674 , A1 (2023).\n- 80. P. Montegriffo, CataPack, http://davide2.bo.astro.it/ ˜ paolo/Main/ CataPack.html . [Online; accessed 10-Nov-2023].  \n- 81. A. Dotter, et al. , Astrophys. J. Suppl. 178 , 89 (2008).\n- 82. W. E. Harris, Astrophys. J. 112 , 1487 (1996).\n- 83. G. Fiorentino, et al. , Astrophys. J. 783 , 34 (2014).\n- 84. S. Raso, et al. , Astrophys. J. 879 , 56 (2019).\n- 85. P. Miocchi, et al. , Astrophys. J 774 , 151 (2013).\n- 86. L. Wang, et al. , Astrophys. J. 892 , 43 (2020).\n- 87. B. Liu, D. Lai, Month. Not. R. Astron. Soc. 502 , 2049 (2021).\n- 88. C. O. Heinke, et al. , Astrophys. J. 625 , 796 (2005).\n- 89. Gaia Collaboration, et al. , Astronomy & Astrophysics 616 , A12 (2018).\n- 90. D. R. Lorimer, M. Kramer, Handbook of Pulsar Astronomy , vol. 4 (Cambridge University Press, 2004).\n- 91. S. M. Kopeikin, Astrophys. J. Lett. 467 , L93 (1996).\n- 92. P. C. C. Freire, et al. , Month. Not. R. Astron. Soc. 412 , 2763 (2011).\n- 93. K. J. Joshi, F. A. Rasio, Astrophys. J. 479 , 948 (1997).\n- 94. E. Poisson, C. M. Will, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge University Press, 2014).\n- 95. S. Chandrasekhar, The mathematical theory of black holes (Oxford University Press, 1983).\n- 96. K. S. Thorne, Astrophys. J. 191 , 507 (1974).  \n- 97. A. Sadowski, et al. , Astronomy & Astrophysics 532 , A41 (2011).\n- 98. S. K. Greif, K. Hebeler, J. M. Lattimer, C. J. Pethick, A. Schwenk, Astrophys. J. 901 , 155 (2020).\n- 99. T. Damour, G. Schafer, Nuovo Cimento B Serie 101B , 127 (1988).\n- 100. R. Blandford, S. A. Teukolsky, Astrophys. J. 205 , 580 (1976).\n- 101. N. Wex, Classical and Quantum Gravity 12 , 983 (1995).\n- 102. N. Wex, S. M. Kopeikin, Astrophys. J. 514 , 388 (1999).\n- 103. V. Venkatraman Krishnan, et al. , Science 367 , 577 (2020).\n- 104. K. Boshkayev, H. Quevedo, B. Zhami, Month. Not. R. Astron. Soc. 464 , 4349 (2017).\n- 105. K. Boshkayev, J. A. Rueda, R. Ruffini, I. Siutsou, Journal of Korean Physical Society 65 , 855 (2014).\n- 106. I. S. Shklovskii, Soviet Astronomy 13 , 562 (1970).\n- 107. P. J. McMillan, Month. Not. R. Astron. Soc. 465 , 76 (2017).\n- 108. T. Dietrich, M. Ujevic, W. Tichy, S. Bernuzzi, B. Brugmann, Phys. Rev. D 95 , 024029 (2017).", 'Acknowledgments': "Acknowledgments The MeerKAT telescope is operated by the South African Radio Astronomy Observatory (SARAO), which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and  \ncalibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time space reference systems department of the Paris Observatory. Parts of this research were conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. 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. PTUSE was developed with support from the Australian SKA Office and Swinburne University of Technology. The authors also acknowledge Max-Planck-Institut fur Radioastronomie funding to contribute to MeerTime infrastructure. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. TRAPUM observations use the FBFUSE and APSUSE computing clusters for data acquisition, storage and analysis. These clusters were funded, installed and operated by the Max-Planck-Institut fur Radioastronomie and the Max-Planck-Gesellschaft. This work has made use of data from the European Space Agency (ESA) mission Gaia ( https://www.cosmos.esa.int/gaia ), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos. esa.int/web/gaia/dpac/consortium ). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.  \nFunding AC, AP, AR, and MBu acknowledge the use of resources from the research grant 'iPeska' (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016. AP and AR also had the support from the Ministero degli  \nAffari Esteri e della Cooperazione Internazionale - Direzione Generale per la Promozione del Sistema Paese - Progetto di Grande Rilevanza ZA18GR02. SMR is a CIFAR Fellow and is supported by NSF Physics Frontiers Center awards 1430284 and 2020265. EDB, AD, PCCF, MK, VVK, NW, JB, WC, DJC, C.-H.R.C., YPM and PVP acknowledge continuing valuable support from the Max-Planck Society. BWS acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694745). C.-H.R.C. acknowledges support from the Deutsches Zentrum fur Luft- und Raumfahrt (DLR) grant NS1 under contract no. 50 OR 2214. C.P. and M.C. acknowledge funding from Italian MIUR throughout the PRIN-2017 grant awarded to the Light-on-Dark project (PI: Ferraro) through contract PRIN-2017K7REXT. MG acknowledges funding from the South African Research Chairs Initiative of the Department of Science and Technology and the National Research Foundation of South Africa.  \nAuthor contributions Conceptualization, AD, AP, BWS, MB, MK, NW, PCCF, SMR, VVK; Data curation, AD, AR, EDB, MB, MG, PVP, SMR, YPM, VVK; Formal Analysis, AD, AR, C-HRC, CP, DJC, MBu, MC, NW, PCCF, TG, VVK; Funding acquisition, AP, MB, MK; Investigation, AD, AR, C-HRC, EDB, MBu, SB, VVK; Methodology, AD, EDB, MB, PCCF, SMR; Project administration, AD, AP, AR, BWS, EDB, MB, MBu, MK, PCCF; Resources, AR, BWS, EDB, JB, MB, MBu, MK, PCCF, C-HRC; Software, AR, DJC, EDB, MB, MG, PCCF, PVP, SB, WC, SMR; Supervision, AP, BWS, DJC, EDB, MK, NW, PCCF, SMR, VVK; Validation, AD, C-HRC, NW, PCCF, TMT, VVK; Visualization, AD, EDB, MC, PCCF, TMT, VVK; Writing - original draft, AD, AR, C-HRC, CP, EDB, MC, NW, PCCF, TMT, VVK; Writing - review & editing, AC, AD, AP, BWS, C-HRC, CP, DJC, EDB, MB, MC, MG, MK, NW, PCCF, PVP, SMR, TG, TMT, VVK.  \nCompeting interests The authors have no competing interests to declare.  \nData and materials availability Our measured ToAs, timing model, and simulation code are available at https://github.com/ewanbarr/NGC1851E Additional Materials/ and archived at Zenodo (49). Folded pulsar timing archives from MeerKAT and the GBT as well as dedispersed pulsar search mode observations from MeerKAT are available at https://keeper.mpdl.mpg.de/d/a59d968eccf and archived at the Max Planck Digital Library (50). Optical and near-ultraviolet data sets can be downloaded from the MAST archive at https://mast.stsci.edu/search/ui/#/hst, using the search terms 'Target=NGC 1851', 'Proposal ID=12311,13297' and 'Filters / Gratings=F275W,F336W'.", 'Supplementary Material': 'Materials and Methods Supplementary Text Figs. S1 to S10 Tables S1 and S2 References (51-108)', 'MeerKAT observations': 'PSR J0514 -4002E was detected with the MeerKAT radio telescope on 24 occasions between 2021 and 2022, and all observations were conducted using the standard configuration for globular cluster observations (10,51). The initial three MeerKAT observations were conducted with the L-band receivers, covering the band 856 to 1712 MHz, and the following 21 observations using the UHF receivers covering the band 544 to 1088 MHz. Data acquisition used the search mode of the Pulsar Timing User Supplied Equipment (PTUSE) instrument (21). This mode recorded time and frequency resolved data that was coherently dedispersed (52) to the nominal dispersion measure of the cluster, 52.14 pc cm -3 (the DM of the only pulsar previously known in the cluster, (17)). Depending on the array configuration, MeerKAT data have time resolutions between 7 and 15 µs. Observations varied in duration between 2 and 4 hours.  \nTo maximise sensitivity to Shapiro delay, six of the UHF observations were scheduled to coincide with the orbital phases of the expected extrema of the unabsorbed Shapiro delay signal (the part of the Shapiro delay that is not easily absorbed by varying the Keplerian orbital parameters [(25), especially their eq. 31] within a single orbit). These phases were determined by synthesising low-noise ToAs using the timing model in Table 1 with the inclusion of a detectable Shapiro delay signal. The ToAs were then fit for the Keplerian orbital parameters with a timing model which assumed no Shapiro delay. The extrema of the unabsorbed Shapiro delay signal are then taken as the extrema of the model residuals from the fit. Synthesising the ToAs using Shapiro delay signals from a range of possible geometries showed negligible variation in the orbital phases corresponding to the extrema. These observations are highlighted in Table S1.', 'GBT observations': 'Using a preliminary timing model for PSR J0514 -4002E derived from only MeerKAT observations (see below), we could recover the signal of the pulsar with the GBT on 6 occasions between 2005 and 2006. Of these detections, 5 were made in observations conducted with the 800 MHz receiver covering the band 795 to 845 MHz and one in an observation conducted with the S-band receiver covering the band 1649.6 to 2249.6 MHz. For both receivers the Pulsar Spigot instrument was used for data acquisition (53). Data were recorded with a time resolution of 82 µs. Observations varied between 80 minutes and 4 hours in duration. Table S1 summarizes all the GBT and MeerKAT observations used in this work.', 'Data reduction': 'The coherently dedispersed MeerKAT search-mode data were folded with 256 bins across the pulse profile using the DSPSR software (52,54) with the timing ephemeris of PSR J0514 -4002E from (16). These folded archives were then downsampled to 256 frequency channels and manually cleaned for any residual radio frequency interference (RFI) using the routines psrzap and pazi , part of the PSRCHIVE software package (55, 56). For both the L-band and UHF observations, a standard profile template was generated by summing the pulse profiles from the corresponding observations to generate a low-noise profile. An analytic profile was then generated from this high signal-to-noise ratio (S/N) profile for each observing band by fitting the low-noise profile with two components, modelled by von Mises functions, using the paas program from PSRCHIVE. The integrated pulse profile and analytic profile for the L-Band and UHF are shown in Figures S1A and S1B, respectively. Although we use separate templates for UHF and L-band, no profile difference is observed between the bands.  \nThe pulse profiles for the individual observations were then cross-correlated with the corresponding analytic profiles for each band, and ToAs were extracted using the pat command  \nin PSRCHIVE. The ToAs were derived with 4 sub-bands and integrations of 10 to 60 minutes, depending on the length of the observation and S/N of the detection.  \nData acquired with the GBT were incoherently dedispersed and phase folded with 64 bins across the profile using the prepfold routine of the PRESTO software package (57,58) using the ephemeris of PSR J0514 -4002E derived from the MeerKAT observations. ToA estimation used a previously described method (23). For ToA estimation, a standard profile template was generated from the highest S/N GBT observation with the 800 MHz receiver (Figure S1C). Due to the unknown profile evolution between 820 MHz and 1950 MHz, only data from the 800 MHz receiver were used in the timing analyses.', 'Phase coherent solution': 'We used the TEMPO(59) timing software for the analysis of the ToAs. Initially, we extended the previous orbital solution (16) to all available MeerKAT ToAs. At this stage, we fitted arbitrary time offsets between each group of ToAs from each observation, because the rotation count between observations was not known due to insufficient precision in the previously-determined timing parameters.  \nIn a second stage, we removed these arbitrary offsets using the DRACULA algorithm (60). This allowed us to determine the rotation counts between observing epochs and, hence, a phasecoherent timing solution for all the MeerKAT data.  \nUsing this preliminary solution, we re-folded the archival data from the GBT observations (23). We detect the pulsar in six GBT observations taken with the 800 MHz and S-band receivers (Table S1). These detections are faint, so we only determine two ToAs from each detection. In the subsequent analysis, we have used only the ToAs from the five observations taken with the 800-MHz receiver.  \nGiven the > 15 -year gap between the last GBT and first MeerKAT detections of the pulsar,  \nwe investigate the validity of the MeerKAT-determined rotational model for the GBT observing epochs. To test this, we used DRACULA to determine the combination of missing rotations between the GBT observations that minimises the χ 2 of the ToA residuals. The algorithm returned 32 possible solutions with reduced χ 2 < 2 . Of these, the most favoured model is that of no missing rotations over the GBT epochs, which has reduced χ 2 = 1 . 019 . The next best solution consists of a single phase jump with reduced χ 2 = 1 . 255 . For the 458 degrees of freedom in the model, these correspond to χ 2 = 466 . 70 and χ 2 = 561 . 05 , respectively. Using equation S5 below, we find that the next best solution has 3 . 2 × 10 -21 times the Bayesian likelihood of the solution in which the GBT observations are phase connected.  \nFor all the datasets, we used pulse profile templates that have, as much as possible, consistent phase definitions. The profiles used for the GBT and the MeerKAT data were aligned by eye to a phase precision of 0.01 to ensure that no large arbitrary time offsets are introduced by TEMPO due to different physical longitudes on the neutron star to which the ToAs refer. TEMPO assumes that all ToAs refer to a consistent longitude; if there are differences in the arrival times caused by different phase definitions, the software assumes that there is a time delay between the systems, something that can bias other parameters, like T 0 (61).  \nThis consistent phase definition allows us to determine the offset between the GBT and the MeerKAT ToAs, which is 0 . 6 ± 1 . 3 ms. This is consistent with zero and smaller than the spin period, consistent with phase connection between the two data sets. The parameters in Table 1 are derived by fitting for this time offset parameter. This results in slightly larger uncertainties than in the case where we do not fit for this offset parameter, therefore a more conservative estimate of parameter precision. This small offset suggests that, within our measurement precision, both observing systems share functionally identical clocks. If we assume that the offset is zero, then the number of rotations between the last 2006 ToA and the first 2021 ToA is 82,461,896,322; this provides a χ 2 of 473, subtracting or adding 1 rotation we obtain  \nχ 2 = 484 . Using equation S5 below, we obtain for each of the two latter solutions 0.0041 times the Bayesian likelihood of the former solution; this implies the former has an overall probability of 99.19%. Additional solutions have negligible probabilities. This number of rotations corresponds to 461,453,260.433689 s, or 14.62 years. The estimated timing parameters under this assumption are consistent within 1σ with the parameters in Table 1. However, since we cannot rule out time offsets smaller than about 1.3 ms, we regard the parameters in Table 1 as more reliable, since they already take the uncertainty of this small offset into account.  \nWealso investigated including a systematic time offset between MeerKAT L-band and UHF observations. The best-fitting value of 18 ± 9 µ s is consistent with an offset arising from the mis-alignment of the profile templates between the two bands. The instrumental offset between the two bands is expected to be less than 5 ns (21).', 'Timing parameters': "We used TEMPO to analyze the ToAs, which were first converted to the Bureau International des Poids et M'esures (BIPM) 2019 terrestrial time standard (BIPM2019). Then, we use the Jet Propulsion Laboratory's DE440 Solar System ephemeris (62) to compensate for the motion of the radio telescopes relative to the Solar System barycentre (SSB). The results are reported in the Dynamic Barycentric Time (TDB) timescale. To describe the effects of the orbital motion and light propagation, we use a modification of the theory-independent DD orbital model (63, 64) that includes the orthometric parametrization of the Shapiro delay (25) (hereon the DDH model). Like the DD model, the DDH model can describe the timing of binary pulsars for a wide range of alternative theories of gravity (64). Table 1 lists the resulting timing parameters. The ToA residuals from the best fitting model as a function of observing epoch and orbital phase are shown in Figure 1. No trends are apparent, and the residuals are consistent with the ToA uncertainties. The overall best fitting model has a reduced χ 2 of 1.019, which was  \nobtained after adjusting the uncertainties of individual data sets such that the reduced χ 2 for each data set is 1.0. The adjustment factors for each dataset were: GBT 800 MHz data: 1.227 (indicating that the original uncertainties are slightly under-estimated), MeerKAT L-band data: 0.989, MeerKAT UHF data: 1.012. These adjustments result in more conservative estimates of parameter precision. They are close to unity, implying the ToA uncertainties were well estimated by the routines used to derive the ToAs. The overall reduced χ 2 is close to 1, indicating that the timing solution provides a statistically justified description of the observed ToAs.  \nWe cannot measure the proper motion of PSR J0514 -4002E because the MeerKAT data only span a period of 1.5 years. We therefore assume that the pulsar's proper motion matches that of NGC 1851 determined from an analysis of HST and Gaia astrometric data (26): µ α = 2 . 128 ± 0 . 031 mas yr -1 , µ δ = -0 . 646 ± 0 . 032 mas yr -1 . We expect that the proper motion of PSR J0514 -4002E might differ slightly from this value, but that difference should be small, given the fact that the GC has a central escape velocity of only 42 . 9 km s -1 (65), corresponding to a proper motion of 0.78 mas yr -1 . We also assume a parallax of 0 . 0858 ± 0 . 0018 mas , which corresponds to the astrometric distance estimate of 11 . 66 ± 0 . 25 kpc for NGC 1851 (26).  \nThe timing model includes a 3σ detection of the second derivative of the pulsar's spin frequency, ¨ ν = (7 . 3 ± 2 . 2) × 10 -26 Hz s -2 . This is likely to be caused by the variation of the acceleration of the system, caused by its change of position relative to the gravitational potential of the GC, or relative to nearby stars (5). Such changes have been seen for other GC pulsars with long-term timing, and are consistent with theoretical expectations (60,66,67). Our measurement of ¨ ν is of a similar magnitude to the values observed for pulsars in other GCs (66,67). Introducing this parameter results in slightly increased uncertainties of all other timing parameters, especially those that improve more with the use of the early GBT data: ˙ ν , ˙ P b and γ E .  \nIn the orbital model we adopted, relativistic effects detectable in the timing are quantified by  \npost-Keplerian (PK) parameters. We measure ˙ ω with > 1000σ significance even if we use only the MeerKAT data, and obtain an upper limit for the orthometric amplitude of the Shapiro delay, h 3 . We fit the DDH model including both Shapiro delay parameters h 3 and ς (the orthometric ratio of the Shapiro delay), but find the model does not converge because the Shapiro delay is not detected. We therefore fixed one of the parameters, ς , to a value of 0.456, which corresponds to an inclination of 49 · from our ephemeris (see below). The assumed value of ς has little effect on the detection of h 3 because the two parameters are designed to be only weakly correlated (25). Using this fixed ς , we obtain h 3 = -0 . 02 ± 0 . 91 µs, which is a non-detection with a 2σ upper limit of 1 . 80 µs.  \nIncluding the earlier GBT data provides a 3σ detection of the variation of the orbital period, ˙ P b = (18 . 1 ± 5 . 6) × 10 -12 s s -1 and a 2σ upper limit for the Einstein delay, γ E < 28 ms . The latter parameter quantifies the variation of a combination of the special relativistic time dilation and gravitational redshift with orbital phase.  \nThe ˙ ω parameter quantifies a relativistic effect which can be used to estimate the total system mass, M , if we assume GR is the correct theory of gravity. However, this is not the case for ˙ P b . For some binary pulsars, including PSR B1913+16 and PSR J0737 -3039A, the observed ˙ P b is dominated by gravitational wave damping (68,69). However, for wider binaries, this effect becomes insignificant; in the case of PSR J0514 -4002E the GR prediction for ˙ P b corresponding to the most likely masses is ˙ P b, GR = -0 . 05 × 10 -12 s s -1 , hundreds of times smaller than the observed ˙ P b . The variation of P b is instead likely to be dominated by the secular variation of the Doppler shift: this does not include the orbital variation of the Doppler shift, which is already subtracted by the timing model. This secular variation of the Doppler shift is itself dominated, in this case, by the line-of-sight component of the acceleration of the system in the gravitational field of the GC, A . Assuming this is the case,  \nA ≃ ˙ P b P b c = (8 . 4 ± 2 . 6) × 10 -9 ms -2 . (S1)  \nThis acceleration also affects the variation of the spin period derivative ˙ P . Our measured value (Table 1) implies a characteristic age of 0.46 Gyr, which is within the range of previously studied MSPs (3) but unusually low. If, instead, ˙ P is dominated by A , then  \nA ≃ ˙ P P c = 10 . 3 × 10 -9 ms -2 , (S2)  \nwhich is consistent with the estimate of A above. We subtract the estimate of A from equation S1 from the observed spin-down to estimate the intrinsic spin-down of the pulsar ( ˙ P int ), independently from the secular variation of the Doppler shift:  \n˙ P int = ˙ P -˙ P b P b P = (3 . 6 ± 4 . 8) × 10 -20 s s -1 , (S3)  \nwhich is a non-detection with 2σ upper limit of 13 . 2 × 10 -20 s s -1 . Thus, the intrinsic spindown of PSR J0514 -4002E is too small to be measured from our data.", 'Robustness of the timing parameters': 'We investigate the robustness of our measurement of the rate of advance of periastron of the pulsar by performing a self-consistent check of the observed values as a function of the number of observations. The results are shown in Figure S2; we find no substantial change in ˙ ω as more observations are added and the fluctuations are consistent with the uncertainties obtained at each stage. This stability indicates that the value of ˙ ω is robust.  \nAs an additional check, we investigated the orbital coverage of all the observations (Figure S3). We find a broad distribution of the observations in orbital phase and time of the year, which improves the determination of ˙ ω and the position and reduces the correlations between orbital and astrometric quantities.  \nThe correlation matrix for the timing parameters produced by TEMPO (Figure S4) indicates only small correlations between orbital and astrometric parameters. The largest correlations are between the three orbital parameters, x , γ E and ω 0 ; this explains their lower numerical precision  \nthan the other parameters in the timing model [(17), their equations 23 and 24]. There is also a relatively large correlation between e and h 3 (this is expected, see (25), their table 4). There are no additional large ( > 0 . 9 ) correlations between other parameters.', 'Masses from χ 2': "The observed ˙ P b is not caused by gravitational wave damping, and we have only marginal limits for the other PK parameters, so the only parameter that directly constrains M is ˙ ω . However, the upper limits on γ E and h 3 provide additional information on the masses. These limits have opposite effects: the upper limit on γ E excludes very large values of m c , while the upper limit on h 3 excludes high orbital inclinations and small values of m c .  \nTo take all these effects into account, we used the DDGR timing model, a modification of the DD model which directly fits the two masses ( M and m c ), assuming that all relativistic effects in the timing of the system are as predicted by GR (63). From the pulsar timing observations we only obtain the product Gm j ( j = c, p ). The masses in units of nominal solar masses (M ⊙ ) are given by the ratio ( Gm j ) / ( GM ) N ⊙ , where ( GM ) N ⊙ ≡ 1 . 3271244 × 10 26 cm 3 s -2 denotes the nominal solar mass parameter (18). Because the model returns this ratio, the measurement uncertainty in G cancels. In cases where absolute masses are needed in physical units we use the value G = 6 . 6743 × 10 -8 cm 3 g -1 s -2 for conversion (18).  \nAs discussed above, the variation of the orbital period ( ˙ P b ) in Table 1 is not caused by gravitational wave damping but instead by the system's acceleration in the gravitational field of NGC 1851. This is taken into account in the model by the inclusion of an extra free parameter ∆ ˙ P b , which quantifies non-relativistic contributions to ˙ P b .  \nWe obtain M = 3 . 887 ± 0 . 004 M ⊙ and m c = 2 . 44 ± 0 . 46 M ⊙ ; this implies m p = 1 . 44 ± 0 . 46 M ⊙ . The uncertainties are not Gaussian distributed; the 2σ upper limit for m p is above the physical upper limit for m p derived from M and the mass function (Figure S5). By forcing  \nTEMPO to iterate several times, we find that there is no stable numerical convergence on a particular value of m c , because the relativistic effects used to estimate m c are not detected.  \nTo estimate uncertainties and take additional relativistic effects into account in a self-consistent way, we have used a Bayesian technique described in detail elsewhere (70), with some adjustments.  \nFor each point in a two-dimensional (2-D) grid of values of M and cos i , we calculate m c using equation 1, then use M and m c as fixed parameters in a DDGR timing model. We use TEMPO to fit this model to the timing data set (a set of ToAs, { t j } ); this is done by minimizing the weighted sum of the squares of the residuals by varying all parameters in the model apart from M and m c (astrometric and spin parameters, Keplerian orbital parameters and the extra ˙ P b caused by the acceleration of the system relative to that of the Solar System projected along the line of sight). For each point in the grid, TEMPO returns a χ 2 value corresponding to the goodness of fit to { t j } , which we denote χ 2 ( M, cos i ) .  \nIf χ 2 min is the global minimum of χ 2 ( M, cos i ) , then each value of  \n∆ χ 2 ( M, cos i ) = χ 2 ( M, cos i ) -χ 2 min (S4)  \nhas a χ 2 distribution with two degrees of freedom (70). We map this to a Bayesian likelihood function p ( { t j }| M, cos i ) :  \np ( { t j }| M, cos i ) = 1 2 e -∆ χ 2 2 . (S5)  \nThe 2-D joint posterior probability density function (pdf) for M and cos i , p ( M, cos i |{ t j } ) , is given by Bayes' theorem:  \np ( M, cos i |{ t j } ) = p ( { t j }| M, cos i ) p ( { t j } ) p ( M, cos i ) , (S6)  \nwhere p ( { t j } ) is the Bayesian evidence, which is calculated from the integral of p ( M, cos i |{ t j } ) over areas of the parameter space where m p ≡ M -m c > 1 . 17 M ⊙ (19,20).  \nThe quantity p ( M, cos i ) is the Bayesian prior; for these we chose a uniform distribution in cos i because that corresponds to assuming that there is no favoured prior orientation to the orbital angular momentum. However, there is no need to sample the full cos i space: for the median value of M the constraint on m p implies that 0 < cos i < 0 . 732 , thus 90 · > i > 42 . 9 · , so we considered a uniform distribution of cos i only within this interval. We do not sample the region -1 < cos i < 0 because it cannot be distinguished from the regions between 0 and 1 given the available timing precision. We also assume a uniform prior distribution of M , due to our prior ignorance about the total mass of the system. We work with M (not the individual masses) because the measurement of ˙ ω is directly related to M , assuming GR. A different choice of prior (such as a uniform distribution of m c ) would require sampling much larger regions of parameter space, much of which would have low probability densities, greatly increasing the computational effort without substantially changing the results.  \nWe marginalize 2-D joint posterior pdfs for M , cos i , m c and m p . Because m c and m p are known for each point of the M -cos i grid, we update their marginalized posterior pdfs after calculating p ( M, cos i |{ t j } ) , as we do for M and cos i . At the end of the process, these marginalized posterior pdfs are normalized.  \nThe medians of the marginalized posterior pdfs and equivalent 1σ confidence intervals are presented in the 'Nature of the companion' section of the main text. The upper and lower limits of the confidence intervals are calculated in such a way that they include 34.13% of the total probability above and below the medians of each pdf. For a 2σ equivalent confidence level (95.45 %) around the medians, we obtain the following confidence intervals: M = 3 . 887 ± 0 . 009 M ⊙ , m p = 1 . 53 +0 . 30 -0 . 33 M ⊙ , m c = 2 . 35 +0 . 33 -0 . 30 M ⊙ and i = 52 +12 · -8 . The mass of m c in this estimate is lower than the value estimated in the DDGR model, it is more conservative because the prior (a flat distribution in cos i ) makes higher orbital inclinations more likely, resulting in lower values of m c (equation 1). The choice of sampling region is also conservative: if we  \nincluded pulsar masses below 1 . 17 M ⊙ , then we would be sampling regions with very large values of m c , which would increase the median for the latter quantity.", 'Search for plasma from the companion': "If the pulsar's binary companion is a main-sequence star, it is likely that plasma expelled from the companion fills a large fraction of the intra-binary medium. For a 2.3 M ⊙ main-sequence companion and a 1.6 M ⊙ pulsar, we calculate the Roche-lobe radius [using (71), their equation (2)] at periastron is 3.027 R ⊙ , using the nominal solar radius adopted by the International Astronomical Union (18). The corresponding Roche-lobe filling factor is 0.536. The filling of the intra-binary medium due to the plasma from the companion leads to two potentially observable effects: firstly, there could be eclipses of the pulses from the pulsar if the intervening plasma is clumpy. We do not find any evidence for any eclipses. Secondly, any intervening plasma would contain free electrons which the radio signal must pass through; this would cause changes in pulse dispersion at different points in the orbit. We investigate this in two ways: firstly, we average the timing residuals from observations recorded at similar angular orbital positions with respect to the ascending node ( ψ ), independently for the top and bottom half of the MeerKAT UHF band, to measure if there are any residual frequency-dependent trends in the ToAs. We do not find any evidence for this. Secondly, we use multi-frequency ToAs for each observation, and independently measure the DM per epoch. Figure S6 shows the averaged timing residuals and DM offset measurements ( δ DM; after subtracting the nominal DM of the pulsar, 51.93 pc cm -3 ) as a function of the orbital position ψ measured from the longitude of the ascending node. The position in the orbit when the companion is located between the observer and the pulsar along the line of sight is defined as superior conjunction, which is at an angular orbital position of 90 degrees with respect to the ascending node. We do not see any statistically significant ( > 3 σ ) deviations in the variation of DM along the orbit.  \nThe second potential indication of intra-binary plasma is variable Faraday rotation, caused by a magnetic field along the line of sight, along the orbit. This would produce an orbital phase dependent change of the rotation measure of the pulsar. Only 17 of the 24 MeerKAT observations include full polarisation information; of these only 3 detect the RM with > 3 σ significance, due to the low signal-to-noise and low polarisation fraction (16) of the pulsar. The RMs of these observations are 13 . 5 ± 2 . 5 , 10 . 7 ± 2 . 0 and 12 . 3 ± 5 . 7 rad m -2 . These are consistent with each other, and have an weighted average of 12 . 2 ± 2 . 2 rad m -2 . We find no evidence for variable RM.", 'Search for radio pulsations from the companion': "If the companion star is a NS, it could potentially produce observable pulsed radio emission under the conditions that it is a radio pulsar, its beam sweeps across the Earth, and its average pulsed flux density is above the detection limit of our data. No companion pulsar was detected in the earlier GC pulsar survey (16). We performed an additional search specifically targeted on PSR J0514 -4002E by de-dispersing each of the 24 MeerKAT observations of NGC 1851 at the DM of PSR J0514 -4002E (51.93 pc cm -3 ), using PRESTO's prepdata routine. We ran the latter in combination with radio frequency interference (RFI) masks, to filter out artificial signals from the data.  \nFor observations that covered an orbital phase interval far from periastron, we performed an acceleration search using the accelsearch routine, with a 'zmax' parameter of 50. At these orbital phases, we assume the line-of-sight acceleration undergone by the companion (and the pulsar) is constant across the 2 to 4 hours of each observation, hence we expect the acceleration search to recover the Fourier power of any pulsed signal from the companion. This assumption does not hold for the five observations that were taken around the periastron passage, for which we expect the line-of-sight acceleration to change substantially during the  \n2 to 3 hours of the observations. The expected change is so rapid that even a search for jerk (rate of change of acceleration) would not be effective. To recover the Fourier power associated with potential companion pulsations from the five observations taken around periastron, we used the PYSOLATOR(72) software, which removes the orbital modulation affecting the signal of the pulsar and companion, using the orbital parameters of the latter. In the case of the PSR J0514 -4002E system, we do not know the size of the projected semi-major axis of the companion orbit, x c , because this depends on the (known) projected semi-major axis of the pulsar orbit, x p , and the unknown mass ratio q = m p /m c , via the relation x c = q x p . Therefore, we demodulated the companion orbit assuming trial q values in the range between 0.519 to 0.784 (which reflects the 1σ uncertainties on m p and m c ) with a step size of 0.03. As the trial q value approaches the underlying true mass ratio of PSR J0514 -4002E, the effect of the orbital motion on the companion's signal is minimised, allowing searches over a range of constant accelerations to become effective. For each periastron passage observation, we performed an acceleration search with a 'zmax' parameter of 20, on all the demodulated time series. All candidate pulsar signals produced from the 24 observations were folded modulo the detected spin frequency with PRESTO's prepfold and the resulting diagnostic plots were inspected visually. We recovered several of the 15 pulsars known in NGC 1851 with DMs close to that of PSR J0514 -4002E(16). The known pulsars in the cluster are all ruled out as the companion, as they are not spatially coincident with and do not share the same orbital parameters as PSR J0514 -4002E. We found no pulsar-like signal that could be ascribed to the companion.", 'Hubble Space Telescope observations': 'We utilised archival observations of NGC 1851 with Wide Field Camera 3 (WFC3) on the Hubble Space Telescope (HST) taken in 2010 as part of the observing programme with proposal ID 12311 (Principal Investigator: Piotto) and in 2014 through proposal ID 13297 (Principal Investigator: Piotto). We downloaded from the MAST archive observations in the F275W and F336W filters, which are sensitive to hot stars such as blue stragglers and WDs. The data-set is composed of 14 images acquired with the F275W filter with exposure time of 1280 s, and 4 images acquired with the F336W filter with exposure time of 453 s.', 'Photometric analysis': 'The photometric analysis was performed using DAOPHOT (73) following methods described elsewhere (74, 75). First, we selected a few hundred bright stars to model the point spread function (PSF) of each image. These PSF models were applied to all the sources detected at more than 5 σ above the local background level. Then, we produced a reference list of stars including all the sources detected in at least half the F275W images. At the corresponding positions of these stars, PSF fitting was forced in all the other frames using the allframe routine in DAOPHOT (76). For each star we homogenized and averaged the magnitudes estimated in different images using the daomaster routine in DAOPHOT, to obtain the instrumental stellar magnitudes and uncertainties.  \nThe observed magnitudes were calibrated to the Vega photometric system by cross-correlation with a previous catalogue derived from the same observations of NGC 1851 (77).  \nThe instrumental positions were corrected for geometric distortion following previously described methods (78) and converted to the ICRS using the stars in common with the Gaia Data Release 3 (DR3) catalog (79). We used the cross-correlation software CATAXCORR(80)  \nadopting a six parameter linear transformation to convert the instrumental pixel positions to the ( α , δ ) absolute coordinate system using Gaia DR3 stars as reference. The residuals of this transformation showed a combined root mean square of 14 mas (9 mas in each coordinate), which we adopt as the 1 σ astrometric uncertainty.', 'Stars near the position of PSR J0514 -4002E': "Two stars are observed within 100 mas of the position of PSR J0514 -4002E. As discussed in the main text, we identify one as a blue straggler ( m F 275 W = 18 . 575 ± 0 . 009 , m F 336 W = 18 . 15 ± 0 . 02 ) and the other as a red giant ( m F 275 W = 19 . 047 ± 0 . 007 , m F 336 W = 17 . 60 ± 0 . 02 ); both of which have estimated masses below the lower limit on the mass for the companion of PSR J0514 -4002E ( 1 . 4 M ⊙ ) obtained without assuming that ˙ ω is relativistic. Assuming that the blue straggler star (BSS) is a H-core burning star, we estimate its mass through a theoretical mass-luminosity relation (81) assuming the cluster metallicity, distance and reddening (82) (2010 edition). We estimate that the blue straggler's magnitude is compatible with a ∼ 1 . 2 M ⊙ star, too low to be the companion of PSR J0514 -4002E. Direct estimation of the mass of the BSS is not possible with the available observations. Main sequence mass-luminosity relationships can be used to estimate BSS masses (83,84), though with a higher uncertainty on the BSS mass of 0 . 5 M ⊙ .  \nWe found no evidence of photometric variability in either of the two stars near the position of PSR J0514 -4002E. If in a binary, variability could arise from heating or eclipsing of the companion star. Although the observations sample only ∼ 25% of the binary orbital period, the lack of variability is consistent with our interpretation that neither of the two stars is associated with the pulsar.", 'Formation of massive compact binaries in NGC 1851': "Although NGC 1851 is not a core-collapsed GC, it does have a high central density relative to most other GCs [(65), their table 2]. This and its small core radius of 0 . 09 ' , which was determined using the stellar density profile (85), imply that it has a high rate of stellar interactions per binary ( γ ) [(6), their figure 1]. The high theoretical γ -value of NGC 1851 is supported by observations of three highly eccentric MSPs in binaries in this GC with degenerate companions that are much more massive than the companions associated with similar systems in the Galactic disk: PSR J0514 -4002A (8,17,23), PSR J0514 -4002D and PSR J0514 -4002E(16).  \nThe encounter rate per binary affects the size and composition of the binary population: a GCwith a large overall stellar interaction rate (like for instance, 47 Tuc) will form large numbers of LMXBs and subsequently radio MSPs, but if the γ -value is low, these binaries will likely evolve without disturbance. Most binary MSP systems found in lowγ GCs have low-mass companions (WDs) in near-circular orbits (66,86) - i.e. they resemble the MSP population in the Galaxy (3). Such GCs are expected to form tight DNS systems at rates that are much smaller than the observed rate of DNS mergers detected by gravitational waves, which is supported by simulations (43). On the other hand, in GCs with a high γ -value, several subsequent exchange encounters involve the same MSP, after the exchange that formed its parent LMXB, which could in principle form many DNS systems. This is the case for NGC 1851 (43). If the merger products of such DNSs stay in the cluster (i.e. if the merger kick is not larger than the escape velocity of the cluster), they will be available for subsequent exchange interactions.  \nWe compare the expected merger kick imparted on a DNS merger product to the escape velocity of NGC 1851. Following previous calculations [(87), their equation 10], we find a 1 . 2 + 1 . 4 M ⊙ DNS merger has an estimated merger kick of w mk = 43 . 8 kms -1 , which is  \nclose to the central escape velocity of NGC 1851 of 42 . 9 kms -1 (65). However, the closer the unknown mass ratio, q , of the two merging NSs is to unity, the smaller is the imparted merger kick ( w mk → 0 as q → 1 ), and the higher the probability is for the merger remnant to remain within the GC. Hence, we conclude that it is no problem to keep a DNS merger remnant within NGC 1851 from a kinematic point of view.  \nDense GCs experience mass segregation, in which more massive stellar populations concentrate closer to the centre, especially in the positions of pulsars and X-ray sources (88). Therefore, the most massive compact objects (NSs, black holes and merger remnants) should, when in equilibrium with the remaining stars in the cluster (when all types of objects have the same average kinetic energy), occupy a small volume near the centre of the GC. Therefore, in the centres of GCs, massive stellar remnants are expected to constitute most of the mass (65). The measured position of PSR J0514 -4002E (Table 1) is offset from the optical centre of the cluster (89) by 0 . 064 ' , which, given estimates of the cluster core radius of 0 . 09 ' and 0 . 044 ' (65,85), places this pulsar either in or close to the core. This is also true of the other massive binary pulsar known in this GC, PSR J0514 -4002A. The abundance and close proximity of massive degenerate objects in this region implies a higher probability of formation of systems like PSR J0514 -4002E.", 'Contributions to the observed post-Keplerian parameters': "Contributions to ˙ ω obs  \nThe observed rate of advance of periastron for a binary pulsar can be produced by the combination of classical contributions and relativistic effects. The observed change in the longitude of periastron ( ˙ ω obs ) is (90)  \n˙ ω obs = ˙ ω M + ˙ ω PM + ˙ ω 3rd -body + ˙ ω SO , (S7)  \nwhere ˙ ω M is the leading-order relativistic periastron advance, ˙ ω PM and ˙ ω 3rd -body are the impact of the proper motion of the system and the possible presence of a nearby third body respectively, and ˙ ω SO includes the contributions from relativistic and classical spin-orbit coupling. Below, we consider each of these effects in turn.  \nThe relativistic periastron advance can be expressed as a function of the total mass of the system and the measured Keplerian parameters, to the lowest post-Newtonian order and for negligible spin contributions, by re-arranging equation 2:  \n˙ ω M = 3 c 2 (1 -e 2 ) ( P b 2 π ) -5 / 3 ( GM ) 2 / 3 . (S8)  \nAproper motion of the binary would change the longitude of periastron due to the variations in the orientation of the orbital plane with respect to our line of sight. This effect can be quantified as (91,92)  \n˙ ω PM = 2 . 78 × 10 -7 µ T sin i cos [Θ µ -Ω] deg yr -1 , (S9)  \nwhere µ T is the magnitude of the total proper motion in mas yr -1 , Θ µ its position angle, and Ω is the longitude of the ascending node of the pulsar's orbit. The maximal contribution from ˙ ω PM is  \n˙ ω PM , max = 2 . 78 × 10 -7 µ T sin i deg yr -1 . (S10)  \nAs discussed above, we cannot measure the proper motion of the pulsar, so have assumed the cluster's proper motion value for our timing solution. However, to estimate the maximum contribution from ˙ ω PM , we use a sum of the proper motion of the cluster ( µ GC ) and the proper motion corresponding to the cluster's central escape velocity ( µ esc ). The value of sin i used in this calculation corresponds to the angle of inclination predicted by the timing model.  \nThe dynamics of the cluster (26) indicate µ esc = 0 . 78 mas yr -1 ; adding this to the proper motion of the GC ( µ GC = 2 . 22 mas yr -1 ), we derive an upper limit to the proper motion of  \nPSR J0514 -4002E as µ T ≤ 3 . 00 mas yr -1 . Using this value in equation S10, the maximum contribution to the change in the longitude of periastron from the proper motion of the system can be estimated:  \n˙ ω PM , max = 1 . 104 × 10 -6 deg yr -1 . (S11)  \nThe corresponding contribution to the observed rate of advance of periastron,  \n˙ ω PM , max = 3 . 184 × 10 -5 ˙ ω obs , (S12)  \nis negligibly small.  \nThe presence of an additional third body near the system (either bound to it, or more likely in a close flyby) could also contribute to the observed change in the longitude of periastron. However, we expect a third body that contributes to the ˙ ω obs to cause additional secular changes in other binary parameters (especially orbital inclination and eccentricity), and produce spinfrequency derivatives ( ˙ ν , ¨ ν , . . . ) that are orders of magnitude larger than we observe. We assume the geometry illustrated in Figure S7. Neglecting angular dependencies and corrections due to eccentricity, | ˙ ω 3rd -body | ∼ (3 π/P b )( m 3 /M )( a pc /r 3 ) 3 [(93), their equation (33)], where m 3 is the mass of the third body and a pc denotes the semi-major axis of the relative motion of the (inner) binary system (about 0 . 12 astronomical units (au), derived using Kepler's third law). The minimum distance r 3 for an external mass m 3 to produce a ˙ ω 3rd -body that has a nonnegligible contribution to the ˙ ω obs (more than two times its measurement uncertainty), and the gravitational acceleration consequently caused by this third body, can be determined using this equation. Using the values in Table 1 and equation (3) in (93) we find  \n∣ ∣ ∣ ∣ ˙ ν 3rd -body ˙ ν ∣ ∣ ∣ ∣ ≳ 180 ( m 3 M ⊙ ) 1 / 3 | cos Θ r | . (S13)  \nConsequently, even a small mass m 3 has to be close to the plane of the sky ( Θ r = π/ 2 ) to produce the observed ˙ ν . m 3 cannot be very small, otherwise the third body would have to be  \nvery close to the (inner) system, which would produce higher order derivatives in ν (93). We use an equivalent analysis to obtain an expression for ¨ ν [(93), their equation (3)]. We assume a (bound) circular orbit for the third body and leading order in m 3 /M , leading to  \n∣ ∣ ∣ ∣ ¨ ν 3rd -body ¨ ν ∣ ∣ ∣ ∣ ≳ 14000 ( m 3 M ⊙ ) -1 / 6 | cos Θ v | . (S14)  \nConsequently, the motion of the third body again has to be (very) closely aligned with the plane of the sky ( Θ v ≃ π/ 2 ) to explain the measurement. In our example, the combination of equations S13 and S14 requires a nearly perfect alignment of the orbital plane of the putative third body with the plane of the sky. This is a highly unlikely orientation (probability of ≲ 10 -5 , before accounting for the probability of having a third body in sufficient proximity ( ≲ (57 au)( m 3 / M ⊙ ) 1 / 3 )). We expect such a third body would produce related changes in the eccentricity and projected semi-major axis of the pulsar orbit [(93), their equations 34 and 35].  \nWe may make this result more rigorous by removing the assumption of circular motion for the third body (outer orbit) and accounting for the high eccentricity of the inner orbit. We have obtained the corresponding expressions for the tidally induced ˙ ω , ˙ e and ˙ x by (semi-analytically) integrating the equations of osculating orbital elements, using as a perturbing force the differential acceleration in the binary system caused by the external mass m 3 [see (94), their equations 3.69 and 3.77] (note, the fully analytical equations for the tidally induced ˙ ω , ˙ x and ˙ e given in (93) assume a small eccentricity for the inner orbit which is not the case in the PSR J0514 -4002E system, where we measure e ≈ 0 . 71 ). These calculations use the geometry of Figure S7, and allow the calculation of the tidally induced ˙ ω , ˙ x and ˙ e for any given distance r 3 = | ⃗r 3 | and given directional angles Φ r and Θ r . Like for the circular outer orbit above, we need the first and second frequency derivative, while this time however we allow the velocity of the third body to have any magnitude and direction (see Figure S7). Like before, we used the standard equations for the externally induced time derivatives of the pulse frequency, ˙ ν 3rd -body  \nand ¨ ν 3rd -body [(93), their equation (3)].  \nIn order to estimate an upper limit on the probability for the presence of a third body, we Monte-Carlo (MC) sample the four dimensional parameter space of Φ r , Φ v , (cos Θ r ) , and (cos Θ v ) , which all have a uniform probability distribution over their corresponding range ( [0 , 2 π ) for the first two and [ -1 , +1] for the last two). When sampling the ( Φ r , cos Θ r , Φ v , cos Θ v ) parameter space, we keep the distance r 3 as a free parameter, to maximise the tidal contribution to ˙ ω obs (i.e. find max | ˙ ω 3rd -body | ) within the constraints given by the observed frequency derivatives (we use | ˙ ν 3rd -body | < 3 | ˙ ν obs | and | ¨ ν 3rd -body | < 3 | ¨ ν obs | ) and orbital parameters (within the 3σ range of ˙ x and ˙ e ). We therefore regard our results as conservative upper limits because they are obtained without considering the low probability of having a third body in the proximity to PSR J0514 -4002E. To obtain a high accuracy in our MC sampling of the ( Φ r , cos Θ r , Φ v , cos Θ v ) parameter space, we conducted 10 10 MC realisations for each individual combination of mass m 3 and velocity v 3 , always counting those configurations that have a max | ˙ ω 3rd -body | that is larger than two times the measurement uncertainty of ˙ ω obs , while at the same time pass the spin-frequency ( ˙ ν , ¨ ν ) and orbital ( ˙ x , ˙ e ) constraints. Our procedure corresponds to calculating a (normalized) volume in the four-dimensional parameter space of ( Φ r , cos Θ r , Φ v , cos Θ v ), defined by the given boundary conditions. This integral is solved numerically using Monte Carlo integration, a standard technique for higher-dimensional integrals. The results of our MC runs are shown in Figure S8. Based on these results, we conclude that the probability for a non-negligible (more than two times the measurement uncertainty) tidal contribution to ˙ ω obs is ∼ 10 -6 , so only highly fine-tuned (directional) configurations could produce such an effect. The maximum ˙ ω 3rd -body from our MC runs that does not violate our observed values for ˙ x ( (4 ± 2) × 10 -13 ) and ˙ e ( ( -0 . 06 ± 1 . 11) × 10 -14 s -1 ) is ∼ 1% of ˙ ω obs . This is too low to influence the conclusions of this paper.  \nWe note that 10 km s -1 as the relative velocity between the third body and the inner binary  \ncorresponds to the escape velocity (lower limit for flyby) at about 100 au .  \nAsecular change in the longitude of periastron could arise from the contribution of the spins of the pulsar and its binary companion due to Lense-Thirring precession ( ˙ ω LT ). A change could also be caused by the spin-induced quadrupole moment of the rotating companion ( ˙ ω Q ). The ˙ ω SO term from equation S7 can be decomposed as  \n˙ ω SO = ˙ ω LT + ˙ ω Q . (S15)  \nWe discuss them separately and calculate the maximum contributions for different possible binary companions of PSR J0514 -4002E.  \nThe dimensionless measure of the spin angular momentum ( S c ) of a body of mass m c is (94,95)  \nχ c = c G S c m 2 c . (S16)  \nFor a NS with a mass larger than ∼ 1 M ⊙ , the maximum value of the spin parameter is χ max ∼ 0 . 7 (48). This is valid for most proposed equation of states (EoSs) and is independent of the mass of the NS. A rotating Kerr BH in general relativity can have spin parameter values going up to 1, mathematically (95) (astrophysical BHs are generally expected to have spin values below the Thorne limit, 0.998 (96); however, the maximum possible spin could also reach values closer to 1 depending on the details of the accretion process (97)). For a group of objects including any NS or BH, χ c < 1 .  \nAccording to (98), for a NS mass of 1 . 4 M ⊙ , the typical moment of inertia of is ∼ 1 . 5 × 10 45 g cm 2 ; using this, the spin angular momentum of PSR J0514 -4002E is then ≈ 1 . 7 × 10 48 g cm 2 s -1 and the spin parameter is χ ≈ 0 . 1 . This is about an order of magnitude smaller than expected for an extremely fast rotating NS or BH companion. We therefore neglect it hereafter.  \nThe contribution to ˙ ω obs from the Lense-Thirring precession of the compact companion is (99)  \n˙ ω LT = n b B c sin 2 i [ (1 -3 sin 2 i )( ˆ k · ˆ s c ) -cos i ( ˆ K 0 · ˆ s c ) ] = n b B c f c , (S17)  \nwhere n b is the orbital frequency, the unit vector ˆ k is normal to the orbital plane defined by the orbital angular momentum, ˆ s c is the unit vector along the spin of the companion, and ˆ K 0 is the unit vector along the line of sight. The explicit expressions for B c and f c are given below. Using equation S16, B c is [(99), their equation 5.17a]  \nB c = 4 + 3 m p /m c 2 c 2 M (1 -e 2 ) 3 / 2 n b S c = Gm c (3 M + m c ) 2 c 3 M (1 -e 2 ) 3 / 2 n b χ c . (S18)  \nTo estimate the maximum ˙ ω LT contribution from a NS or BH companion, we calculated the values for B c over a range of pulsar and companion masses using χ c = χ max = 1 . The factor f c can be expressed as  \nf c = (1 -3 sin 2 i )( ˆ k · ˆ s c ) -cos i ( ˆ K 0 · ˆ s c ) sin 2 i = -2 cos ϑ +sin ϑ cot i cos ϕ SO , (S19)  \nwhere ϑ (0 ≤ ϑ ≤ π ) is the angle between the spin of the companion ( S c = S c ˆ s c ) and the orbital angular momentum and ϕ SO (0 ≤ ϕ SO < 2 π ) is the angle swept by the projection of S c during its precession around ˆ k measured in the orbital plane as shown in Figure S9 (64).  \nFor our calculations, we used a uniformly distributed array of masses given by m p ∈ [1 . 1 M ⊙ , 4 . 0 M ⊙ ] and m c ∈ [1 . 4 M ⊙ , 60 . 0 M ⊙ ] , where the upper limit for the mass of the companion is set as the maximum mass allowed by an offset in the reducedχ 2 value within a 3σ deviation from the solution. A companion with a mass exceeding 60 M ⊙ would produce a large Shapiro delay, independent of the inclination, inconsistent with our measured limit for the system [(100), their equation 2.20]. For each pair of pulsar and companion mass ( m p , m c ) within this range, we evaluate the values of B c and f c to obtain the maximum ˙ ω LT . The largest absolute value of f c , for any such given pair, can be expressed as | f c | max = √ 4 + cot 2 i . The ratio ( ρ LT , 1 ) of the maximum value of ˙ ω LT with respect to ˙ ω M for every pair of mass,  \nρ LT , 1 ≡ | ˙ ω LT | max ( m p , m c , χ c = 1) ˙ ω M ( m p , m c ) , (S20)  \nis shown as a function of mass in Figure S10. The maximum value of this ratio occurs for the largest companion mass and the smallest pulsar mass, reaching  \nmax( ρ LT , 1 ) = 0 . 0072 . (S21)  \nConsequently, for the whole range of masses, i.e. m p ∈ [1 . 1 M ⊙ , 3 . 0 M ⊙ ] and m c ∈ [1 . 4 M ⊙ , 60 . 0 M ⊙ ] , the Lense-Thirring precession has negligible contribution to the total advance of periastron. The contribution to ˙ ω SO in equation S15 from a rotationally-induced quadrupole moment ˙ ω Q for a NS companion (101) or a BH counterpart (102) is  \n˙ ω Q ∼ 10 -3 ˙ ω LT . (S22)  \nThis contribution is also negligibly small.  \nIn the alternative case of a massive WD, contributions to the ˙ ω SO value from the LenseThirring precession and the quadrupole moment are of similar magnitude to each other (103). The upper mass limit for a rigidly rotating WD is about 1 . 47 M ⊙ (34), so we consider an uniformly rotating 1 . 4 M ⊙ WD for the companion of PSR J0514 -4002E. The typical spin value for such a massive WD with extreme rotation is S c ∼ 0 . 6 × 10 50 g cm 2 s -1 (104) and assuming the minimum pulsar mass of 1 . 1 M ⊙ (see above), we used equations S17, S18, and S19 above to calculate the maximum contribution from the Lense-Thirring precession:  \n| ˙ ω LT | max / ˙ ω M ( m p = 1 . 1 M ⊙ , m c = 1 . 4 M ⊙ ) = 0 . 0016 . (S23)  \nUsing Q = 0 . 6 × 10 49 g cm 2 for a rotating WD (105), we calculate  \n˙ ω Q / ˙ ω M ( m p = 1 . 1 M ⊙ , m c = 1 . 4 M ⊙ ) = 0 . 0039 . (S24)  \nThis is too low for the observed advance of periastron to be explained by a (uniformly rotating) WD, even if it is rotating at the mass-shedding limit. We conclude that ˙ ω obs is dominated by  \nthe leading order GR contribution (equation S8), and thus the inferred mass of the companion ( m c > 2 . 09 M ⊙ ) is too large for it to be a massive WD.  \nTo summarize, there is negligible contribution to the rate of advance of periastron from the proper motion of the binary or the presence of a possible third body near the system. Over the whole range of companion masses investigated, the maximum possible change in the longitude of periastron due to Lense-Thirring precession, ˙ ω LT (Figure S10), is more than a factor of 100 smaller than the observed rate of advance of periastron (Table 1). Combining this with equation S7, we conclude that the observed change in the longitude of periastron of the pulsar is relativistic, and the total mass of the system can be approximated by the leading-order periastron advance term using equation 2.  \nFigure S10 shows the constraint on the total mass of the system given by ˙ ω obs = ˙ ω M . For any combination of pulsar and companion mass ( m p , m c ), the maximum possible contribution from Lense-Thirring precession is always less than 1% of the observed rate of advance of periastron, and is therefore negligible.", 'Contributions to ˙ P b , obs': 'The observed variation in the orbital period of a binary pulsar ( ˙ P b , obs ) can be expressed as (90):  \n˙ P b , obs = ˙ P b , int + ˙ P b , GC + ˙ P b , Shk + ˙ P b , Gal , (S25)  \nwhere ˙ P b , int is the intrinsic orbital period derivative and ˙ P b , GC is the contribution from the lineof-sight acceleration of the system in the gravitational field of the cluster ( A ). The third and the fourth terms in this equation are the change in the orbital period of the binary due to the Shklovskii effect (106) and the difference between the Galactic accelerations of the binary and the Solar System.  \nThe observed orbital period derivative (Table 1) is ˙ P b = (18 . 1 ± 5 . 6) × 10 -12 s s -1 . The change in the orbital period from the Shklovskii effect, associated with the proper motion of the  \nsystem, is given by:  \n˙ P b , Shk = µ 2 d c P b = 9 . 02 × 10 -14 s s -1 , (S26)  \nwhere d is the distance between the Earth and the pulsar. Here we use the cluster distance estimate of (26), 11.66 kpc. With the available timing baseline, we cannot measure the proper motion of the binary so have used the proper motion of the cluster (see above). If the binary is gravitationally bound to the GC, the proper motion of the pulsar cannot differ much from the proper motion of the cluster; the same applies for ˙ P b , Shk .  \nWe calculate the contributions to ˙ P b , obs from the Galactic acceleration using a Milky Way potential (107), finding:  \n˙ P b , Gal = -2 . 52 × 10 -14 s s -1 . (S27)  \nThe sum of ˙ P b , Shk and ˙ P b , Gal is almost 300 times smaller than ˙ P b , obs , which implies that the latter is dominated by the effect of A :  \n˙ P b , obs ≃ ˙ P b , GC = A c P b . (S28)  \nWe conclude that ˙ P b , obs depends only on A .', 'Contributions to ˙ x obs from Lense-Thirring precession': "The observed change in the projected semi-major axis of the binary pulsar ( ˙ x obs ) could be caused by a physical change in the orbit of the pulsar, a temporal change in the angle of inclination of the system, or a combination of both. While the semi-major axis of the pulsar is unaffected by the Lense-Thirring effect, precession of the binary orbit could potentially be observable as a slow change in the angle of inclination of the system, caused by the spins of the binary components. As discussed above, the spin angular momentum of the pulsar is about an order of magnitude smaller than the maximum possible value of the spin angular momentum of the companion so we neglect it in further analysis; however, the spin of the companion could lead to a secular change in the projected semi-major axis depending on its alignment.  \nThe change in the inclination angle as an effect of the Lense-Thirring precession of the compact companion is, to the leading order (64,99):  \n( di dt ) LT = 4 + 3 m p /m c 2 c 2 M (1 -e 2 ) 3 / 2 n 2 b ( S c · ˆ i ) = Gm c (3 M + m c ) 2 c 3 M (1 -e 2 ) 3 / 2 n 2 b χ c sin ϑ sin ϕ SO . (S29)  \nSecular evolution in the angle of the inclination would occur if the spin of the companion has a non-negligible component along the direction of the ascending node. This, in turn, would change the observed projected semi-major axis:  \n˙ x LT = x cot i ( di dt ) LT . (S30)  \nTo estimate the maximal contribution to ˙ x LT , we consider a pulsar of 1 . 2 M ⊙ and a light BH companion with a mass of 2 . 7 M ⊙ . For this mass pair, the maximum observable value of the variation in the projected semi-major axis of the pulsar is ˙ x = (1 . 7 × 10 -13 ) χ c . For a Kerr BH χ c ≤ 1 . If this BH is the result of a NS binary merger (irrespective of the collapse time), the maximum value of the dimensionless spin of the companion is theoretically expected to lie in the range χ c ∈ [0 . 6 , 0 . 875] (47).  \nWe left ˙ x as a free parameter in our radio timing analysis to determine the maximum expected value. We calculated this using the DDGR model, due to the large expected correlation between γ E and ˙ x (17); this method has been used in previous measurements of ˙ x (103). Assuming a range of companion masses ( 2 . 3 M ⊙ ≤ m c ≤ 2 . 7 M ⊙ ) we obtain ˙ x ∼ (4 ± 2) × 10 -13 . This is not significant, and the 1σ uncertainty is larger than our estimated upper limit for ˙ x LT .  \nIn addition to the precession of the orbit due to spin-orbit coupling, there are several other geometric and physical effects that could influence the observed change in the projected semimajor axis: these include the proper motion of the binary, the emission of gravitational waves or mass loss from the system, a change in the aberration or Doppler shift, and the presence of a hypothetical third body around the binary. However, the contributions from any of these terms  \nare at least an order of magnitude smaller than the contribution from the Lense-Thirring effect, and are hence neglected in our analysis.  \nFigure S1: Multi-frequency pulse profiles for PSR J0514 -4002E . The summed profiles are shown in grey, and the red dashed lines show the corresponding analytic profiles. (A) Pulse profile for MeerKAT L-Band (856 to 1712 MHz) generated by summing the profiles for all observations. (B) Same as (A), but for MeerKAT UHF (544 to 1088 MHz). (C) Pulse profile of the highest S/N observation for GBT 800 MHz (795 to 845 MHz).  \n<!-- image -->  \nCumulative number of observations  \n<!-- image -->  \nFigure S2: Observed rate of advance of periastron ( ˙ ω obs ) for PSR J0514 -4002E. In both panels, the vertical errorbars indicate 1σ uncertainties. (A) ˙ ω obs as a function of the cumulative addition of observations from the MeerKAT dataset. (B) Same as panel (A), but only for the last 7 observations from the MeerKAT dataset. The last point indicates the value of ˙ ω derived with the addition of the GBT data set to the full MeerKAT dataset, and the horizontal line corresponds to this value as reported in Table 1.  \nFigure S3: Orbital and monthly coverage of PSR J0514 -4002E observations. Each point shows a single ToA estimation as a function of the time of year it was obtained and the mean anomaly of PSR J0514 -4002E's orbit to which it corresponds. A mean anomaly of 0 denotes periastron passage. Green points indicate ToAs estimated from the archival GBT dataset, while the orange and blue points indicate ToAs estimated from the MeerKAT L-Band and the UHF datasets, respectively.  \n<!-- image -->  \nFigure S4: Correlations between the timing parameters of PSR J0514 -4002E. Here we show the Pearson product-moment correlation coefficient (colour scale) between the timing parameters presented in Table 1, as reported by TEMPO. We also include J 0 , the time offset between MeerKAT L-band and MeerKAT UHF data, and J 1 , the time offset between GBT data and MeerKAT UHF data.  \n<!-- image -->  \nFigure S5: Mass and orbital inclination ( i ) constraints for the PSR J0514 -4002E system . In panel (A) , we display the constraints on cos i and the companion mass m c ; the gray shaded area is excluded by m p > 0 . In panel (B) , we display the same constraints on m p and m c , the gray area is excluded by the mass function and sin i ≤ 1 . The solid blue lines indicate the nominal and ± 1 -σ constraints on the total mass of the binary ( M ) obtained from the measurement of the rate of advance of periastron ( ˙ ω ). At the scale of the figure, these three lines appear practically superposed. The solid red and dashed magenta lines indicate upper 1, 2 and 3σ limits on h 3 and γ E respectively. The thick black lines are contours that include 95.4% of the total probability density. Panels (C), (D) and (E) present the probability density functions (pdfs) for cos i , m p and m c respectively, these are normalized to the point of maximum probability density. The black lines indicating the median (solid) and the limits of the confidence intervals with equivalent ± 1, 2 and 3σ confidence levels (dashed). These pdfs are truncated below m p = 1 . 17 M ⊙ , the smallest known NS mass.  \n<!-- image -->  \nTrue Anomaly  \n( )  \nFigure S6: Timing residuals and dispersion measure offset as a function of angular orbital position. (A) Averaged residuals between the ToAs in the MeerKAT UHF band the and timing model of PSR J0514 -4002E split into two frequency bands as a function of angular orbital position ψ with respect to the ascending node of the pulsar. (B) Same as panel A, but for DM offsets (with respect to the pulsar's DM (Table 1)). An orbital position of ψ = 90 · corresponds to superior conjunction and f denotes frequency.  \n<!-- image -->  \n/uni20D7  \nFigure S7: Geometry for third-body calculations . ( X,Y,Z ) denotes the right-handed Cartesian coordinate system, where the X-axis is aligned with the ascending node of the pulsar orbit and the Z-axis is normal to the plane of the sky. i is the angle of inclination between the pulsar orbit and the plane of the sky, and ω is the longitude of periaston. The putative third body (orange) is at distance ⃗r 3 and moves with velocity ⃗v 3 . We assume uniform probability distributions for Φ r , cos Θ r , Φ v , and cos Θ v . The angles ϕ r and θ r are defined with respect to the orbit ( ⃗ L : orbital angular momentum of the pulsar binary) [(93), their ϕ 2 and θ 2 respectively].  \n<!-- image -->  \n/uni20D7  \n/uni20D7  \nFigure S8: Upper limits on the probability of non-negligible tidal contributions to ˙ ω obs . The figure shows the upper limit on the probability that ˙ ω obs has a tidal contribution (caused by a nearby third body) that exceeds twice the measurement uncertainty. Results are shown for different masses ( m 3 ) and velocities ( v 3 ) of the third body (see text). Masses significantly below 0.1 M ⊙ are excluded, since they have to be very close to the binary system in order to meet the ˙ ω obs condition above. As a consequence, they would lead to higher order time derivatives in the spin frequency ( ... ν , .... ν , . . . ), which are not observed.  \n<!-- image -->  \nDirection of  \norbital angular  \nmomentum  \n/uni0302  \nFigure S9: Definition of the orbital orientation and geometric parameters of the system. The conventions of (64) are followed. The plane of the sky is denoted by the vectors ˆ I and ˆ J , and the vectors ˆ I ≡ ˆ i and ˆ j are a part of the orbital plane. The direction of the orbital angular momentum ˆ k is perpendicular to the orbital plane and is inclined at an angle i with respect to the line-of-sight vector ˆ K 0 . The spin angular momentum of the binary companion is defined by S c , and its direction with respect to the line-of-sight is given by the angle λ . The misalignment angle between the spin and the orbital angular momentum is denoted by ϑ , where 0 ≤ ϑ ≤ π . S c forms a precession cone around ˆ k , and the angle swept by the projection of the former in the orbital plane measured from -ˆ j is given by ϕ SO .  \n<!-- image -->  \n/uni0302  \n/uni0302  \n/uni0302  \n/uni0302  \n/uni0302  \n/uni0302  \nFigure S10: Maximum contribution to the rate of advance of periastron from LenseThirring precession . Companion masses are on a logarithmic scale; the upper limit of 60 M ⊙ is set from radio timing analysis (see text). The pulsar masses, on a linear scale, span 1.1 M ⊙ to 3.0 M ⊙ (see text). The grey shaded region is excluded by the constraint sin i ≤ 1 and the boundary of this region corresponds to i = 90 · . The red curve shows the nominal measurement of ˙ ω for PSR J0514 -4002E, which equals the relativistic rate of advance of periastron given by equation S8.  \n<!-- image -->  \nTable S1: Radio timing observations of PSR J0514 -4002E. Shapiro delay campaign observations (indicated by an asterisk) covered orbital phase ranges of 0.21-0.23, 0.61-0.63, 0.88-0.90 and 0.98-0.034.  \nTable S2: Masses and references of BH/NS components within or near the lower mass gap. For GW sources, the error bars are 90% symmetric credible intervals. For GW170817, we followed (35, 108) and estimated the upper mass limit assuming a remnant mass of 90% of sum of two merging NS components minus 0 . 04 M ⊙ of baryonic ejecta mass. For the radio pulsars, see the individual references in the table for a description of how error bars on masses are calculated in more detail."}