© The Authors 2024

1. Introduction

Radio pulsar searches in globular clusters (GCs) have yielded 325 discoveries in 42 different GCs¹, demonstrating that GCs are unusually efficient pulsar factories (see e.g. Camilo & Rasio 2005; Ransom 2008; Freire 2013). The high core densities in GCs (up to ∼10⁶ stars pc⁻³) encourage the formation of binary systems and exchange interactions between binaries. In many of these exchanges, those neutron stars (NSs) that have crossed the “death line” (that is their radio emission turned off and hence have become undetectable), become members of binary systems. The evolution of their companions leads to mass transfer onto these NSs, spinning them up to millisecond spin periods. During this stage, these systems are observed as low mass X-ray binaries (LMXBs). The latter are ∼10³ times more abundant per unit stellar mass in GCs compared to the Galaxy (Clark 1975). Once accretion stops, these NSs become detectable as radio millisecond pulsars (MSPs; see e.g. Alpar et al. 1982; Radhakrishnan & Srinivasan 1982; Ivanova 2013, and references therein).

The remarkable rotational stability of MSPs (Ṗ ∼ 10⁻²⁰ s s⁻¹) has allowed for precise measurements of astrometric, spin and binary parameters of various systems in GCs. For MSPs in compact binary systems, when this precision is coupled with long timing baselines (>10 yr), it enables the measurement of general relativistic effects via measurement of Post-Keplerian (PK) parameters like the rate of advance of periastron (), the Shapiro delay and the Einstein delay, which have allowed measurements of a few NS masses (e.g. Lynch et al. 2012; Ridolfi et al. 2019; Corongiu et al. 2023) and in one case tests of gravity theories (Jacoby et al. 2006). Apart from individual systems, an ensemble of precisely timed pulsars can help constrain the structural properties of GCs and place upper limits on the mass of a potential intermediate mass black hole in their centres (e.g. Prager et al. 2017; Perera et al. 2017; Freire et al. 2017; Abbate et al. 2018, 2019). Furthermore, these pulsars can also probe the ionised gas in the intra-cluster medium (Freire et al. 2001a; Abbate et al. 2018) and also along the line-of-sight (Martsen et al. 2022).

The high core densities in GCs also promote the formation of a wide range of unique binary pulsars whose properties stand out from the pulsars typically found in the Galactic field. For example, most known binary pulsars in the Galactic field with substantial eccentricity are double neutron star systems where a supernova explosion of the companion induces an eccentricity (Tauris et al. 2017). However, binary pulsars in GCs can undergo multiple close encounters with neighbouring stars inducing a significant eccentricity in these objects (Phinney 1992; Heggie & Rasio 1996). These encounters can further lead to an exchange of pulsar companions where usually a higher mass object replaces the lighter companion that spun up the pulsar (e.g. Verbunt & Freire 2014), also termed a secondary exchange encounter. This is evidenced by the GCs with the largest rate of such encounters having several highly eccentric systems with fast spinning pulsars and a companion mass larger than expected from standard recycling scenarios (e.g. DeCesar et al. 2015; Ridolfi et al. 2019), which are therefore very likely to be secondary exchange products. Apart from the possibility of finding exotic and unique pulsars, the enhanced sensitivity of current generation telescopes like MeerKAT and FAST, provides further motivation for continuing GC pulsar searches. One such globular cluster rich with pulsars is Terzan 5.

Terzan 5 (Ter5 hereafter) has historically been one of the most widely studied GCs and has yielded remarkable results across multiple wavelengths. Located in the inner bulge of the Galaxy at a distance of D = 6.62 ± 0.15 kpc from Earth (Baumgardt & Vasiliev 2021), it is believed to be a remnant of primordial structures that are integral to Galaxy formation (Ferraro et al. 2009). It is massive ((1.09 ± 0.08) × 10⁶ M_⊙; Baumgardt & Vasiliev 2021) and is among the GCs with the highest stellar encounter rate (Γ ∼ 6800; Bahramian et al. 2013). This is commensurate with Ter5 holding the record for hosting the largest number of confirmed MSPs in a GC (39; prior to this work) of which more than 50% are in binary systems (20; prior to this work). In contrast, NGC 104 (also known as 47 Tuc) which hosts the second highest number of confirmed MSPs, has a stellar encounter rate nearly a factor of seven lower than Ter5 (Γ ∼ 1000; Bahramian et al. 2013). Owing to the excess reddening from the Galactic bulge at optical and infra-red wavelengths, Ter5 has been observed and analysed more in the high-energy end of the spectrum. Apart from evidence for gamma-ray emission in GeV (Abdo et al. 2010) and TeV ranges (H. E. S. S. Collaboration 2011), at least 50 sources have been identified in X-rays (Heinke et al. 2006). Recent deep X-ray studies have established cross-matches to several known pulsars (see Bogdanov et al. 2021, and references therein). Conversely, sources with significant X-ray variability with no known pulsar counterparts have been attributed as potential candidates for “spider” type pulsars (e.g. Urquhart et al. 2020). These are pulsars where the companion material is being ablated away by the pulsar wind and this outflow material can obscure the pulsed radio emission leading to dynamic eclipses across an orbital period cycle. Depending on the companion mass, the “spider” systems can be classified as redbacks (M_c ∼ 0.2 − 0.4 M_⊙) or black-widow systems (M_c ≪ 0.1 M_⊙; Roberts 2013). The X-ray emission in these spider pulsars is usually caused by an intra-binary shock (see e.g. van der Merwe et al. 2020). These results have established the need for synergistic studies across the electromagnetic spectrum for an enhanced understanding not only of individual pulsars but also to trace the origins and formation of Ter5.

The radio pulsar searches in Ter5 have benefited from a combination of large telescopes combined with advanced instrumentation and unique search techniques. The first Ter5 discovery, PSR J1748−2446A, was found using the Very Large Array using a standard periodicity search and is still the most compact binary known in Ter5 (P_b = 0.07 d; Lyne et al. 1990). Later, a Fourier domain acceleration search technique was developed which provides improved sensitivity to binary pulsars (Ransom et al. 2002). This technique, along with the S-band receiver and the SPIGOT pulsar backend (Kaplan et al. 2005) on the Robert C. Byrd Green Bank Telescope (GBT) yielded a flurry of discoveries (Ransom et al. 2005). Soon after, this setup also yielded the current record holder for the fastest spinning pulsar PSR J1748−2446ad (P = 1.39 ms; Hessels et al. 2006). More recent discoveries have benefited from the application of different search techniques. For example, application of an additional acceleration derivative or “jerk” dimension to the binary pulsar search space helped discover PSR J1748−2446am (Andersen & Ransom 2018). PSR J1748−2446ae was found using a dynamic power spectrum search technique where the orbital period is roughly of the order of the observation time span searched (Ransom et al., in prep.). Three more pulsars were found via stacking of Fourier power spectra from hundreds of hours of archival GBT data (Cadelano et al. 2018).

The latest pulsar discovery Ter5an, benefited significantly from the enhanced sensitivity of the MeerKAT radio telescope despite using just 42 out of the 64 antennas (Ridolfi et al. 2021). Apart from demonstrating the value of MeerKAT for pulsar surveys, this study demonstrated the benefits of using archival data from the GBT. Based on an initial orbital solution from MeerKAT, the archival data allowed the quick determination of a long baseline (16 yr) phase connected solution. This large baseline even allowed the detection of PK effects like orbital period derivative (Ṗ_b) and rate of advance of periastron (), which yielded measurements of the systemic acceleration of the system in the potential of the cluster and its total mass (M_tot = 2.97 ± 0.52 M_⊙). Archival data has since proved valuable for timing several other MeerKAT GC pulsar discoveries (Douglas et al. 2022; Ridolfi et al. 2022; Gautam et al. 2022a; Abbate et al. 2022; Vleeschower et al. 2022).

Recently, Martsen et al. (2022) published spectral indices and flux densities for most of the known pulsars in Ter5. The corresponding pseudo-luminosity showed a turnover in the “logN–logS” relation which describes the cumulative number of sources detectable at a given telescope sensitivity. This suggested that all the searches conducted so far are incomplete (see Fig. 3 in Martsen et al. 2022). It also agrees with previous population synthesis simulations suggesting that more than 100 detectable pulsars (allowing for beaming) reside in Ter5 (Bagchi et al. 2011; Chennamangalam et al. 2013), thus motivating further searches.

In this paper, we describe ten new millisecond pulsars discovered in Ter5 with the MeerKAT radio telescope. The paper is structured as follows. We present the observations with MeerKAT and GBT in Sect. 2. The corresponding data reduction and techniques applied for searching and timing the new pulsars are presented in Sect. 3. We discuss the physical properties of individual discoveries in Sect. 4 and also provide phase connected timing solutions for nine of these pulsars. Section 5 discusses the individual mass constraints of the pulsar and companion in some selected binary pulsar discoveries. In Sect. 6, we provide a discussion on the scientific outcomes of the new discoveries individually as well as collectively and the prospects for future studies involving Ter5. Finally, we state our conclusions in Sect. 7.

2. Observations

2.1. MeerKAT

Ter5 was observed with the MeerKAT radio telescope using the L-Band (856–1712 MHz) receiver on 5 September 2020 (hereafter Epoch 1) and on 6 January 2021 (hereafter Epoch 2). This was done as part of the Globular Cluster pulsar survey within the umbrella of the Transients and Pulsars with MeerKAT (TRAPUM) large survey project (Stappers & Kramer 2016). Each observation spanned 4 h and used 288 synthesised beams spanning roughly a 3′ radius. Of the 288 beams, 38 were placed on positions of known pulsars². The rest of the 250 beams were placed in a hexagonal packed tiling centered at the position of PSR J1748−2446N (right ascension and declination ; Ransom et al., in prep.) since this pulsar is located within the central core of Ter5.

The point spread function (PSF) of the individual beam and the corresponding tiling pattern was generated using the MOSAIC software suite³ (Chen et al. 2021) assuming an overlap factor of 70%. This indicates that neighbouring beams shared a common boundary at a distance where the power level drops to 70% from the centre of the PSF. Figure 1 gives a detailed visualisation of the beam tiling pattern as well as the positions of known pulsars and the discoveries from this work overlaid across both observations.

Fig. 1. TRAPUM beam tilings of the two search observations (Epoch 1 and Epoch 2; see text) of Ter5, performed with MeerKAT on 05 Sep 2020 (left) and on 06 Jan 2021 (right) at L-Band. There are 288 beams in total per observation. Each individual beam has been overlaid with the corresponding point spread function. The vertical bar shows the reference gray scale for the fractional sensitivity with reference to the boresight position of the individual beam. The different concentric circles indicate the core radius (red), half-light radius (yellow-dashed) and half-mass radius (blue dotted). The known pulsars are marked in each tiling as black crosses. The discoveries made in each epoch are denoted as red stars along with their individual names. The green crosses indicate the positions of the discoveries but from the alternate epoch where the discovery was first made. The positions of the discoveries comes from timing.

Total intensity SIGPROC format filterbank data for each individual synthesised beam were generated using Max-Planck-Institut für Radioastronomie’s (MPIfR’s) Filterbanking Beamformer User Supplied Equipment (FBFUSE) and was written to disk on MPIfR’s Accelerated Pulsar Search User Supplied Equipment (APSUSE) high performance compute cluster. A detailed description of FBFUSE and APSUSE can be found in Padmanabh et al. (2023), and references therein. The filterbank data were recorded at a sampling time of 76 μs and 4096 frequency channels with a channel bandwidth of 0.208 MHz. Additionally, we also recorded a single beam via the Pulsar Timing User Supplied Equipment (PTUSE) backend (Bailes et al. 2020) placed on the position of PSR J1748−2446N in Epoch 1 and 2. The PTUSE data were recorded with a 9.57 μs time resolution with 1024 channels spanning 856 MHz after applying coherent dedispersion at a dispersion measure (DM) of 237 pc cm⁻³, which is roughly the average of the DMs of all the previously known pulsars. Due to the limited capacity for long-term offline storage, the PTUSE data were later decimated in time from 9.57 μs to 75.29 μs and to 512 frequency channels. This ensured that despite the reduction in size, the data fidelity was minimally affected (since the average dispersion measure for the pulsars in Ter5 is known, and the decimation in frequency accounts for that).

These observations offered multiple advantages compared to previous Ter5 observations conducted with MeerKAT as described in Ridolfi et al. (2021). Firstly, both observations were done with 56 antennas each thus providing a 33% boost in telescope gain compared to the observations made by Ridolfi et al. (2021). Secondly, the semi-major axis of an individual tied-array beam is ∼15 arcsec compared to 0.5 arcmin from before, owing to the usage of outer antennas that in turn increased the maximum baseline of the MeerKAT interferometer. Thus, any new discovery would already have a reasonable constraint on the position. This, however, comes with the increased computational and data storage cost of having to process a large number of beams. Finally, the large number of beams tiled across a given patch of sky provides a nearly uniform sensitivity across a wide area compared to observations with single dish telescopes like Parkes or GBT where the sensitivity is reduced away from the boresight position.

Apart from the above two observations, we also conducted five observations between 26 June 2023 and 7 July 2023 as part of a Shapiro delay campaign and to also derive orbital solutions for four of the ten discoveries made in Ter5 (further details are given later in Sect. 4). Finally, we also utilised archival PTUSE data recorded as part of the Globular cluster theme under the MeerTIME project (Bailes et al. 2020) between May 2019 and February 2020 (reported in Ridolfi et al. 2021). Table 1 summarises all the MeerKAT observations and the respective configuration parameters used.

2.2. GBT

Apart from MeerKAT data, we also used archival data from GBT observations for deriving long-term phase connected solutions for all our discoveries. GBT data from three backends, namely, SPIGOT (Kaplan et al. 2005), GUPPI (DuPlain et al. 2008), and VEGAS (Prestage et al. 2015) were used. Pulsar/GC observations with these back-ends are described in more detail by Ransom et al. (2005), Cadelano et al. (2018) and Martsen et al. (2022), respectively. This represents a total of ∼130 observations across a timeline spanning 2004–2021.

3. Data analysis

3.1. Cleaning and subbanding

Before searching, the TRAPUM filterbank data from all beams of all epochs were put through a crucial preprocessing step. The data were first cleaned using the Inter-Quartile Range Mitigation algorithm (Morello et al. 2022), thus significantly reducing the impact of radio frequency interference (RFI). Following this, groups of 16 channels were dedispersed at the nominal cluster DM of 237 pc cm⁻³ reducing the total number of channels to 256. This subbanding process not only reduced data volume but also sped up subsequent processing steps. Post subbanding, the raw data recorded with 4096 channels were deleted to free up storage space on APSUSE.

3.2. Search strategy

We used PULSAR_MINER (Ridolfi et al. 2021), a Python-based wrapper for PRESTO (Ransom 2011), for implementing acceleration searches on the subbanded filterbank data. PULSAR_MINER has already proven to be successful in discovering multiple pulsars across various globular clusters (Ridolfi et al. 2021, 2022; Vleeschower et al. 2022; Abbate et al. 2022; Chen et al. 2023). First, we dedispersed the time series across 500 different DM trials between 230 and 250 pc cm⁻³ with a step size of 0.05 pc cm⁻³. The upper and lower DM trial limits were chosen to extend the range slightly beyond the lowest (Ter5Q; DM = 234.50 pc cm⁻³) and highest DM (Ter5D; DM = 243.83 pc cm⁻³) among the known pulsars in Ter5. Our search used incoherent harmonic summing with 8 harmonics and a threshold of z_max = 200 where z_max is the maximum spectral drift (in terms of frequency bins) due to linear acceleration. The limit placed on the value of z was to strike a balance between sensitivity to binary pulsars and the computational expense of running the searches. In general, z depends on the duration of the observation (T_obs) and the line-of-sight component of the orbital acceleration of the pulsar (a_l) as (see e.g. Andersen & Ransom 2018)

(1)

where P is the pulsar’s spin period, h is the harmonic number and c denotes the speed of light. Using Eq. (1), we can derive the maximum acceleration a_l, max we are sensitive to based on the value of z_max. For example, a pulsar spinning at P = 2 ms and z = 200 would experience a_l, max = 9.2 ms⁻² for a 1 h data segment when searching only the fundamental frequency (h = 1).

In order to be sensitive to compact binary pulsars, we searched non-overlapping data segments of 0.5 h, 1 h, 2 h, and the full 4 h time span. This segmented search approach enabled the best sensitivity to pulsars with an orbital period above 5 h. We also applied a sifting algorithm to retain candidates which are detectable in at least 3 neighbouring DM trials and are above a 4σ threshold in Fourier significance. Candidates that made this cut were folded and visually inspected. These criteria led to an average of 1000 candidates per beam. The parameters of promising candidates were used to fold data from neighbouring beams from the same epoch and the closest beams in the alternate epoch. A detection in more than one epoch gave a strong indication that the pulsar candidate was real.

The PULSAR_MINER pipeline using the above-mentioned configuration parameters was run on the ATLAS⁴ supercomputer operated by the Max Planck Institute for Gravitational Physics in Hanover, Germany. In order to speed up processing, PULSAR_MINER was run with PRESTO-ON-GPU⁵ enabling the acceleration search routine to be run on GPUs and speed up this step by a factor of 20–30.

3.3. Deriving orbital solutions

Discoveries that showed a non-zero acceleration or a change in barycentric period across epochs strongly indicated that the pulsars were in binary systems. As a first step, we plotted the observed spin period and acceleration parameters from multiple epochs on a period-acceleration diagram (Freire et al. 2001b). If all the points tracked out an ellipse, this gave an immediate indication that the pulsar is most likely in a near-circular orbit. We then fit a parabola (as explained in Freire et al. 2001b) derive the spin period of the pulsar P and two binary parameters, the orbital period (P_b) and the semi-major axis of the pulsar’s orbit (a) projected along line of sight expressed in time units (x ≡ a sin i/c, where i is the generally unknown orbital inclination). This ephemeris was used in turn as an input to fit_circular_orbit.py from PRESTO, which also fits for the time of passage through periastron (T₀) by assuming that a single sinusoidal modulation will fit all detections. If the initial spin period vs. acceleration pattern was not an ellipse, we used fitorb.py from PRESTO to also fit for eccentricity (e) and longitude of periastron (ω).

In order to obtain more data points for better constraining the orbital solution, we first used the initial orbital solution (obtained with relatively few epochs) to predict the expected spin period and acceleration in other epochs. We then folded data from other epochs using the predicted spin and acceleration values as starting input parameters. This method was successful in increasing the overall number of detections and led to a more robust orbital solution when put through another iteration of fit_circular_orbit.py or fitorb.py. The routine to calculate the expected spin period and acceleration for a given epoch is now available as a Python routine termed binary_info.py in PRESTO.

Depending on the segment length in which the discovery was made, a different strategy was applied to derive a rough orbital solution. For example, if discoveries were made in the shortest data segments of 30 min, we also attempted to get detections in other 30-min data segments across the 4-h observation span by restricting the input DM, spin and spin derivative parameters to search. Multiple closely spaced detections helped to easily break the degeneracy in obtaining a unique orbital solution. However, if the discovery was made in the full 4-h data span, the pulsar was searched with a refined search space but in archival GBT data. If enough detections could be obtained (particularly, containing closely spaced detections spanning up to a day), we could be confident of fitting for an orbital solution. For those pulsars where we could not obtain enough detections with GBT data, we used the five follow-up observations spanning 10 days (explained in Sect. 2.1 earlier) as a means to obtain closely spaced detections and ease the orbit solving process.

3.4. Timing

In order to obtain a precise estimate of the astrometric, spin and binary parameters, we attempted to extract times-of-arrival (TOAs) from all available data (from MeerKAT and GBT). Firstly, we folded the data using prepfold from PRESTO with the best ephemeris derived from the orbital solution described earlier. For isolated pulsars, the best barycentric spin period and DM served as starting points to create an ephemeris to build on. Folding the data quickly revealed more detections across several epochs from which more TOAs could be extracted and in turn help in obtaining long-term phase connected solutions. Topocentric TOAs were extracted from the folded PRESTO archives files (pfd format) using get_TOAs.py after cross-correlating an analytical template to all detected pulse profiles. This analytical template was obtained using pygaussfit.py from PRESTO by fitting multiple Gaussian profiles with varying centroid positions, widths and heights. Depending on the brightness of each detection, the number of TOAs extracted varied from observation to observation.

Finally, all the TOAs were fit for a timing model using TEMPO (Nice et al. 2015). We applied a scaling factor (known as EFAC) for each of the different subsets of TOAs derived from different observing back-ends. This factor increases the uncertainty values of the TOAs and in turn makes the overall reduced χ² value close to one. Some discoveries had a large number of bright detections which allowed for manually phase connecting the pulsar across the entire timing baseline. However, there were some pulsars where manual attempts failed. In these cases, we applied an automated timing procedure termed as the Algorithmic Pulsar Timer for Binaries⁶ (APTB; Taylor et al. 2024) that could help determine the exact number of rotation counts between the sparsely sampled data points. APTB expands on the Algorithmic Pulsar Timer (Phillips & Ransom 2022) algorithm by implementing techniques required for phase connecting binary systems in an automated manner. Furthermore, APTB uses robust statistical tests including the dependence on the reduced χ² (similar to DRACULA, developed by Freire & Ridolfi 2018) for decision making while building phase connection between data points. APTB also adds new parameters to fit on-the-fly, depending on lengths of data that are already phase connected from previous steps. A detailed description of this automated routine is described in Taylor et al. (2024). One pulsar (Ter5ap; see Sect. 4.2) could not be solved by APTB (or DRACULA), but was eventually solved by a new method currently in development (Clark & van Haasteren, in prep.), which greatly speeds up the phase-connection procedure by exploiting covariances between the uncertain rotation counts between observations to avoid costly model re-fitting.

To describe the orbits, we have used several orbital models available in TEMPO based on the theory-independent Damour & Deruelle (“DD”) orbital model (Damour & Deruelle 1986), which is used to estimate the Keplerian and PK parameters of the pulsar’s orbit described in the next subsection. The first, known as “DDFWHE”, is nearly identical to the DD model, except that it uses the orthometric parameterisation of the Shapiro delay (Freire & Wex 2010). A second variation of the DD model, known as “ELL1”, is especially suited for very low-eccentricity orbits (Lange et al. 2001). A third model, known as “BTX” (Shaifullah et al. 2016), is derived from the Blandford-Teukolsky (“BT”) model (Blandford & Teukolsky 1976) and allows a non-predictive description of the random orbital period variations observed for several eclipsing systems. The special “DDGR” model (Damour & Deruelle 1986) assumes the validity of general relativity (GR) to estimate masses directly and self-consistently from all observed relativistic effects. Details on these models and their parameters are given later in Sect. 4 as and when they are used to describe different pulsars.

3.5. Derived parameters from Keplerian and Post-Keplerian parameter measurements

While fitting for different timing models with the measured TOAs, we also derived constraints on the properties of some of the discovered binaries based on estimates from Keplerian and PK parameters. The Keplerian parameters correspond to the five aforementioned orbital parameters, namely: P_b, x, e, ω, and T₀. When dealing with pulsars in compact binary systems, higher order relativistic effects become more prominent in the timing analysis, which can be quantified by additional PK parameters in a theory-independent way. These PK parameters depend on the Keplerian orbital parameters and the individual masses as prescribed by a specific gravity theory. For all cases below, we make the assumption that GR describes the strong-field gravity regime, unless stated otherwise. We give details of the derived and Keplerian and PK parameters below.

(i) Mass function: From Kepler’s third law, we can derive the mass function f of the system from two of the Keplerian parameters of the pulsar’s orbit

(2)

where i is the orbital inclination M_p is the pulsar mass and M_c is the companion mass divided by a solar mass (M_⊙) and is an exact constant, the solar mass parameter (, Prša et al. 2016, where G is Newton’s gravitational constant) in time units. This high precision is warranted by the fact that the product G M_⊙ is known much more precisely than either G or M_⊙.

Owing to the sin i ≤ 1 limit, one can derive a lower limit for M_c assuming a certain fixed value for M_p. In this paper, we assume M_p = 1.35 M_⊙ for deriving the minimum companion mass unless stated otherwise.

(ii) Rate of advance of periastron (): This PK parameter quantifies the rate of change in the longitude of the periastron, which determines the orientation of the orbit with respect to the observer’s line of sight. If the effect is purely relativistic, then the effect is always positive. Assuming the validity of GR we can express as a function of the total mass of the system (M_tot = M_p + M_c) as (e.g. Taylor & Weisberg 1982):

(3)

where n_b = 2π/P_b.

(iii) Shapiro delay: In cases where the orbit of a binary pulsar is viewed nearly edge-on, there is also a possibility of measuring a time delay in the pulsar signal owing to the gravitational influence of the companion. This is known as the Shapiro delay (Shapiro 1964). There are several ways to parameterise it. When there is a weak detection of the Shapiro delay, especially if the orbital inclination is not very close to 90 deg, then the use of the “orthometric” parameters (Freire & Wex 2010) avoids strong correlations. Assuming the validity of GR, the orthometric ratio (ς) and amplitude (h₃) are given by:

(4)

(5)

Even a faint detection of this effect can, when combined with a measurement of , result in very precise mass measurements (e.g., Martinez et al. 2015; Stovall et al. 2019; McKee et al. 2020). As described later in Sect. 5, the same technique has been applied to derive precise mass measurements for three binaries in Ter5.

(iv) Variation of the orbital period (Ṗ_b): In some systems, this variation is caused by the emission of gravitational waves, which provide a precise test of gravity theories (e.g., Taylor & Weisberg 1982, 1989; Kramer et al. 2021 and references therein). However, for the binaries described in this paper, this effect is very small compared to the effect on Ṗ_b,obs caused by the change in the Doppler factor of the binary system due to its motion. Differentiating the expression for the Doppler factor as a function of time, and assuming a negligible intrinsic variation of the orbital period, we obtain (Damour & Taylor 1991):

(6)

where a_cluster represents the acceleration contribution from the pulsar in the gravitational potential of the cluster, a_gal is the contribution from the Galactic field potential and the last term is the contribution from the Shklovskii effect (Shklovskii 1970) owing to transverse motion of the pulsar. This term depends on the total proper motion (where μ_α is the proper motion contribution along right ascension and μ_δ along declination) as well as the distance of Ter5 from the Earth D. Of these effects, a_cluster is usually dominant. Thus, by estimating a_gal and a_shk, we can determine a_cluster from precise measurements of Ṗ_b,obs. Those estimates can be used to constrain the mass model of the cluster (see e.g. Prager et al. 2017).

In addition, it allows us to uniquely solve for the true value of the intrinsic spin period derivative of the pulsar (Ṗ_true) as:

(7)

Subtracting Eq. (6) from (7) means we can estimate Ṗ_true as

(8)

Once Ṗ_true is obtained, other derived parameters like the characteristic surface magnetic field (B = 3.2 × 10⁻¹⁹ PṖ_true) and characteristic age (τ_c = P/2Ṗ_true)) can also be obtained.

(v) Variation of the projected semi-major axis (ẋ): A change in the projected semi-major axis can be a consequence of changes in the physical size of the orbit or a change in i from the changing viewing angle due to the motion of the binary. There are multiple factors that contribute to it including PK effects (see e.g. Lorimer & Kramer 2004). Here we assume that the contribution from effects like aberration, Doppler modulation, gravitational wave damping and spin orbit coupling are negligible and the primary contribution is from proper motion. The observed ẋ can be written as (see Arzoumanian et al. 1996; Kopeikin 1996):

(9)

where μ_α and μ_δ are the proper motion terms in right ascension and declination (expressed, in this equation, in mas yr⁻¹) respectively; i is the orbital inclination and Ω is the longitude of ascending node. We can place an upper limit on the maximum contribution from proper motion and constrain i:

(10)

where μ is also expressed in mas yr⁻¹. Placing this limit on the inclination angle can in turn place limits on the individual masses of the system.

(vi) Einstein delay (γ): This PK parameter quantifies a delay caused by the variation of the gravitational redshift and special relativistic time delay with orbital phase. In GR, it can be related to the binary masses as:

(11)

However, this effect is hard to measure for wide, eccentric binaries with a small change of ω within the timing baseline. This has only been achieved successfully for one wide binary, PSR J0514−4002A (with P_b = 18.8 d, Ridolfi et al. 2019). This was done by showing first that the effect of γ_E on the timing of these systems is a linear variation of the projected semi-major axis of the pulsar’s orbit (ẋ_γ, see their Eq. (25)). Therefore, this term superposes inevitably with the ẋ from other causes, like the effect of the proper motion discussed above, ẋ_μ. Only if the expected values of ẋ_μ ≪ ẋ_γ can measure γ, as was done for PSR J0514−4002A.

3.6. Localisation

The narrow width of the synthesised tied-array beams allows for an instantaneous localisation of new discoveries with an uncertainty of tens of arcseconds (at L-Band). However, if a pulsar is detected in multiple neighbouring beams with different signal-to-noise (S/N), one can better constrain the position. The SeeKAT software (Bezuidenhout et al. 2023) implements such an algorithm by measuring the likelihood of the true pulsar position weighted by the point spread function (PSF) and the corresponding different S/N values in neighbouring beam positions. This method can provide a sub-arcsec localisation and has demonstrably eased the phase connection process when fitting for the position of the pulsar (see Bezuidenhout et al. 2023, and references therein). As explained later in Sect. 4, we used SeeKAT to localise two of the faint discoveries and this quickly led to a unique phase connected timing solution for one of these.

3.7. Radiometer flux density estimate

In order to obtain the limiting flux density of the survey (S_min), we used the radiometer equation (Morello et al. 2020) given as:

(12)

where n_pol is the number of polarisations which is 2. The S/N threshold was chosen to be 10. β accounts for digitisation losses and was chosen as 1.0 due to very minimal loss (< 0.01%) in 8-bit data (Kouwenhoven & Voûte 2001). The gain of the telescope G was chosen as 2.5 K Jy⁻¹, corresponding to 56 dishes of MeerKAT representative of most of our observations. The system temperature was chosen as 26 K after taking into account the receiver temperature and the sky temperature at L-Band. t_obs is the observation time which was chosen as 4 h. About 25% of the band is unusable leading to the effective bandwidth BW_eff of 642 MHz. The duty cycle of the pulsar δ was chosen as 10%. Finally, we assumed a search efficiency factor of ϵ = 0.7 based on the work by Morello et al. (2020) when accounting for FFT based searches up to 8 incoherent harmonic sums. Finally, this is the best case scenario assuming maximum sensitivity per synthesised beam. Given the overlap factor of 0.7 between beams, the worst case scenario is a limiting flux density of 16.78 μJy.

Additionally, we also used Eq. (12) without ϵ to obtain flux density estimates and corresponding pseudo-luminosity values (L = S_minD² where D is the distance to Ter5) for all the discovered pulsars. These values are given in Tables 2–4. It is important to note that the flux values are subject to significant changes owing to refractive scintillation as well as different spectral indices per pulsar (e.g. Martsen et al. 2022).

3.8. X-ray source cross matching

As mentioned in Sect. 1, many millisecond pulsars in Ter5 are known to have associations with X-ray sources and such a multi-wavelength study could help better understand the properties of each individual source. We used the catalogue provided by Bahramian et al. (2020) of all X-ray sources in globular clusters to check for associations with any discoveries made. Figure 2 depicts the region within 1.5 times the core radius of the Ter5 field, with X-ray sources, known pulsars, and new pulsar discoveries plotted. We used a conservative simple cut of 0.5 arcsec as the maximum separation between the X-ray source and the discovered pulsars to suggest any association. This cut was chosen to mitigate the impact of any errors from previously applied source extraction methods. We comment on potential (not certain) cross-matches between the position obtained from radio timing and those reported from X-ray imaging in Bahramian et al. (2020) in the next section.

Fig. 2. Ter5 field shown with each axis approximately spanning 1.5 times the core radius of the cluster. The positions of the X-ray sources (black rings), known pulsars (blue dots) as well new discoveries (red dots) are also overlaid. Particular emphasis is made on the VLA-38 source (brown star) obtained from radio imaging and whose coordinates are reported in Urquhart et al. (2020). The timing position of Ter5ar (green star) is shown to overlap significantly with VLA-38 as well as with an X-ray source (CXOU J174804.63−244645.2) in the backdrop. The dashed lines surrounding the X-ray sources and VLA-38 are the 3-σ positional uncertainties.

4. Discoveries and their properties

We have so far completed the searches and candidate viewing for all beams that lie within the core radius of Ter5 as well as beams that were placed on the positions of known pulsars. This totals to roughly 45 beams each for Epoch 1 and 2 (Table 1). Based on the search strategy described in Sect. 3.2, we have confirmed ten new pulsars in Ter5. Nine of these discoveries are in binary systems. Figure 3 shows a collage of pulse profiles for each of these discoveries. We now give a detailed description of each of the ten Ter5 discoveries.

Fig. 3. Integrated pulse profiles of all the newly discovered pulsars in Ter5. These plots are obtained after summing together individual profiles from different epochs and aligning the profiles with respect to a reference template profile.

4.1. Ter5ao

Ter5ao (PSR J1748−2446ao) was the first Ter5 pulsar to be discovered using data from Epoch 1 of the TRAPUM GC pulsar survey. It was found with no hint of acceleration. On conducting a refined search at a single DM trial and lowering the threshold for the Fourier significance by a factor of 2, the pulsar was redetected in Epoch 2. The change in barycentric period confirmed the binary nature of the system. Based on these two detections, the GBT data were searched for more detections after refining the DM obtained from the MeerKAT detections.

The majority of the detections showed a negative acceleration similar to Epoch 1 indicating that the pulsar possibly spends more time behind the companion than in front of it. This suggested that the orbit could be significantly eccentric. Using multiple detections from the GBT data, an orbital solution was derived using fitorb.py. The initial estimated orbital parameters were P_b ≃ 57.55 d, x ≃ 62.3 lt-s and e = 0.32. Based on these orbital parameters, the minimum companion mass was estimated (using Eq. (2)) to be M_c, min ≃ 0.7 M_⊙. Using this ephemeris, we were able to obtain many more detections over 19 yr of GBT data and generated more than 300 topocentric TOAs (including MeerKAT and GBT data). We manually phase connected these TOAs and the final timing solution is provided in Table 2.

The long timing baseline of 19 yr enabled significant measurements of several important astrophysical parameters. Firstly, we were able to measure proper motion in right ascension (μ_α = −1.6(2) mas yr⁻¹). Given that Ter5 is located close to the ecliptic (b = 1.68 deg), constraining the proper motion in declination via timing is difficult. We were also able to get significant detections of spin frequency derivatives up to the fourth order (see Table 2). Using our measurement of and Eq. (3), we derive M_tot = 3.154(17) M_⊙. Combining this with the mass function in Eq. (2), we can additionally constrain M_p < 2.23 M_⊙ and M_c > 0.92 M_⊙.

If the system has a low orbital inclination, the masses of the pulsar and the companion would be closer to each other and the system would be a double neutron star (DNS) system. However, if the system is nearly edge-on, the pulsar mass could be high (> 2 M_⊙). This would be evident if a Shapiro delay signature were to be detected in the data. We thus carried out an observing campaign with 5 epochs between 26 Jun 2023 and 07 July 2023 since the time of superior conjunction was expected on 30 Jun 2023. We did not detect a Shapiro delay, so no precise estimates of the individual masses can be made at the moment. We defer the discussion for the Shapiro delay constraints and individual mass measurements to Sect. 5.

We also obtained a 2.7-σ detection of ẋ =(4.4 ± 1.6) × 10⁻¹⁴ implying a maximum inclination of 50(9) deg (based on Eq. (10)). If this is confirmed more precisely with continued timing, it will rule out the largest pulsar masses. We were also able to measure the orbital period derivative Ṗ_b = (164.7 ± 7.1) × 10⁻¹² s s⁻¹ which is consistent with the predicted value from the cluster acceleration alone (166.38 × 10⁻¹² s s⁻¹). Using Eq. (8), we determined the intrinsic spin period derivative to be Ṗ = 1.13 × 10⁻²⁰ s s⁻¹ which is consistent with other MSPs. The post-fit residuals for Ter5ao and other pulsars discussed below are depicted in Figs. 4 and 5.

Fig. 4. Post-fit residuals after fitting a timing model to the times of arrival (TOAs) obtained from GBT and MeerKAT data spanning about 19 yr. Each colour represents TOAs obtained from different back-ends of both the telescopes (as shown in the legend). The secondary y-axis also shows the corresponding residuals as a fraction of the rotational spin phase. These residuals are shown for Ter5ao, Ter5ap, Ter5aq, and Ter5as, respectively. The residuals for Ter5ar is depicted here for only 2 yr with phase connection.

Fig. 5. Post-fit residuals as shown in Fig. 5 but for the remaining four pulsars with long-term timing solutions. These are Ter5at, Ter5au, Ter5av, Ter5aw and Ter5ax.

4.2. Ter5ap

Ter5ap (PSR J1748−2446ap) was first discovered in Epoch 2. Folding the neighbouring beams revealed multiple detections, suggesting that the candidate was real. Furthermore, it was detected with a significant acceleration of −0.280(5) ms⁻². A corresponding refined search of Epoch 1 data revealed no detection. We then attempted to detect the pulsar in GBT data, but only managed to obtain a detection in one epoch. More data were thus needed to obtain an orbital solution. Using data from the campaign between 26 Jun and 07 July 2023, we made 3 more detections. Using multiple iterations of fitorb.py and also aided by the binary_info.py routine (see Sect. 3.3) we were able to get multiple detections in GBT data that eventually helped in getting an orbital solution with an orbital period of P_b ≃ 21.3 d and x ≃ 13.10 lt-s along with an eccentricity of e = 0.90. These parameters imply a highly eccentric orbit where the pulsar has a positive acceleration for roughly 1 day out of the 21.3 day orbit. Figure 6 shows the dependence of the observed spin period on the mean anomaly and acceleration as predicted by the orbital solution.

Fig. 6. Initial orbital solution for Ter5ap depicted. The left plot shows the line of sight acceleration as a function of orbital phase (which is in turn a function of the mean anomaly). The red line indicates the orbital solution model and the blue dots are data points. Note that the error bars on the data points are smaller than the marking symbol size. All points have a negative acceleration given that the pulsar spends very little time near periastron (roughly a day out of a 21.3 day orbital period). There are no positive acceleration data points which would better constrain the orbital solution. This is further demonstrated by showing the plot on the right depicting acceleration v/s the observed spin period. One data point clearly does not obey this initial solution.

The sparsity of detections from data spanning 19 yr proved difficult for long-term phase connection initially despite multiple manual attempts. However, we were able to phase connect the data with a new technique as mentioned in Sect. 3.4 (Clark & van Haasteren, in prep.). We initially used the GUPPI and VEGAS data to find phase-connected timing solutions using this technique. We found multiple distinct solutions that resulted in similar reduced χ² values, all of which extrapolated well to the MeerKAT and SPIGOT TOAs. Upon further inspection, we found that these were related to each other, with spin frequencies separated by n/P_b for small integers |n| < 5 relative to the solution with the least χ² value. These were caused by a lack of detections near periastron, around which integer pulsar rotations could be gained or lost without deteriorating the residual χ² value. Folding GBT observations closest to periastron with these solutions led to new detections, including in an observation spanning the periastron on MJD 57573, which eventually led to one solution being unambiguously preferred over the others. The final timing solution is presented in Table 2.

Similar to Ter5ao, we were also able to measure astrometric as well as Post-Keplerian parameters. We detected proper motion only along the right ascension direction (μ_α = −2.8(5) mas yr⁻¹). The timing solution also yielded implying a total mass of M_tot = 1.991 ± 0.004 M_⊙. Besides this, we also obtained a significant detection of the orbital period derivative Ṗ_b = (130.356 ± 21.691) × 10⁻¹² s s⁻¹ which is consistent with the predicted value (144.856  ×  10⁻¹² s s⁻¹). We used Eq. (8) again to obtain the intrinsic spin period derivative (Ṗ = 4.284 × 10⁻²⁰ s s⁻¹) and in turn derived the surface magnetic field (B = 4.05 × 10⁸ G) and characteristic age (τ_c = 1.38 Gyr) which are all consistent with the MSP population. In Sect. 5, we discuss the implications of the non detection of γ for Ter5ap on the companion mass.

4.3. Ter5aq

Ter5aq (PSR J1748−2446aq) was discovered in Epoch 2 and was found in a 1 h segment at an acceleration of −2.471(75) ms⁻² and a spin period of 12.52 ms. It was confirmed in Epoch 1 after dedispersing and searching at the DM value from Epoch 2. The Epoch 1 detection showed a different barycentric period and a residual drift in pulse phase indicative of an unaccounted jerk term. These observations suggested a compact binary pulsar. We then tried to obtain detections in other 1 h segments. However, it was detected in just one other segment suggesting that the pulsar could be eclipsing.

Using these two detections we were able to obtain a preliminary orbital solution with P_b ≃ 0.12 d and x ≃ 0.03 lt-s, implying M_c, min ∼ 0.013 M_⊙. This solution worked well for Epoch 2 and confirmed that Ter5aq is a black-widow system. The orbital solution yielded multiple detections in data from the corresponding closest GBT observation epochs.

However, this solution resulted in significant phase drift for Epoch 1 suggesting that the orbital period needed improving. Using more detections from GBT, we searched for an improved estimate of T₀ using the SPIDER_TWISTER⁷ routine (Ridolfi et al. 2016). It is useful for detecting pulsars in spider-type systems which experience significant orbital period variability (for well-studied examples, see Shaifullah et al. 2016; Ridolfi et al. 2016), but also for short-period binaries where the orbital period is not known precisely at first. The routine was used to search for the best T₀ value within a limited orbital phase range for the observation epoch of interest. After detecting the pulsar and measuring the local T₀, we could further constrain the orbital period by fitting for an integer number of orbits between the different T₀ values. This improved the S/N of the detections in Epoch 1. After generating TOAs for Epoch 1 and Epoch 2 with this ephemeris, we were able to better fit for the orbital parameters. Multiple iterations of this process after including more data points eventually yielded detections from most of the GBT observations across 19 yr.

Attempts to manually phase connect the TOAs failed and hence we used the APTB algorithm (Taylor et al. 2024). Initially, 13 yr of data (excluding early GBT data from SPIGOT) were given as an input to APTB. The intricacies involved in these iterations are explained in detail in Sect. 6 of Taylor et al. (2024). Using this solution, we were able to manually extend the timing baseline to 19 yr after inclusion of TOAs from SPIGOT data. The full timing solution is provided in Table 3.

The final timing solution yielded a significant detection of proper motion only along the right ascension direction (μ_α = −1.7(5) mas yr⁻¹) and up to third order spin frequency derivatives (see Table 3). Owing to the circular nature of the system, we were not able to obtain any significant PK parameter measurements despite the long 19-yr baseline. In particular, we do not detect the random orbital variability seen in several other black-widows systems. For this reason we can describe the orbital motion using the ELL1 orbital model.

The position of Ter5aq is in close proximity to an X-ray source (CXOU J174804.75−244642.5) with a separation of ∼0.37 arcsec. This source is reported to be associated with Ter5ab by Urquhart et al. (2020). Although Ter5ab is closer to the X-ray source (0.09 arcsec), it is isolated. However, the prominence of X-ray luminosity arising from eclipsing systems like Ter5aq suggests otherwise and thus warrants further examination.

4.4. Ter5ar

Ter5ar (PSR J1748−2446ar) was found in a 30 min segment in Epoch 1 at a spin period of 1.95 ms. Inspecting the other segments quickly revealed that the pulsar was found in 7 out of 8 segments in Epoch 1. Using these detections, we were able to derive a robust orbital solution with P_b ≃ 0.51 d and x ≃ 1.5 lt-s implying a M_c, min ∼ 0.34 M_⊙. Using the derived ephemeris, we were able to get detections with GBT data. We then applied the same strategy as for Ter5aq, by constraining the orbital period using integer orbits between T₀ values (provided by SPIDER_TWISTER) to better constrain the orbital period. Using the new orbital period we were able to obtain a detection in Epoch 2. We then refined the orbital solution further by generating TOAs for both MeerKAT epochs and fitting via TEMPO. This solution revealed an eclipse in Epoch 2 confirming that Ter5ar is an eclipsing redback system. We then extended the data baseline to obtain TOAs from all the archival GBT data.

When attempting manual phase connection, it was quickly clear that the timing of this system is not trivial. This was not surprising owing to significant long-term changes in the measurement of T₀. To demonstrate this, we estimated the T₀ value using SPIDER_TWISTER for every epoch across a 19 yr baseline. The variation of T₀ across the entire baseline of TOAs obtained is shown in Fig. 7. The ΔT₀ that is the difference between the T₀ expected using a constant orbital period and the measured T₀ can vary by tens of seconds.

Fig. 7. Significant orbital variability in Ter5ar demonstrated. The difference between expected and observed measurement of epoch of periastron (T0) is shown as a function of time in MJD. A stable orbit would show a flat horizontal line but the plot shows the variations ranging from –20 s to 60 s, thus demonstrating the difficulty in timing such systems.

A common way to describe such orbits is by using the aforementioned BTX orbital model. We applied this model to 2 yr of GBT data obtained with SPIGOT and were able to manually get a phase connected solution only after applying six orbital frequency derivatives and up to fourth order spin frequency derivatives. It was difficult to extrapolate data points beyond 2 yr due to the scale of the orbital variability (as shown in Fig. 7). Moreover, the cadence of observations after these 2 yr also reduced from approximately once a month to once in 3 months. This two year timing solution is given in Table 3. Interestingly, APTB was independently able to phase connect data spanning a year, but the high orbital variability prevented convergence with longer datasets (Taylor et al. 2024). A deeper analysis of the timing variability across the entire data span and in turn obtaining a long-term timing solution is work in progress.

Urquhart et al. (2020) conducted deep radio continuum imaging of Ter5 at 2–8 GHz yielding 24 sources where 19 of them could be associated with previously known pulsars and X-ray binary systems. They observed 3 such sources to have a steep spectral index in radio data and a hard X-ray photon index, suggesting a spider-type system with an intra-binary shock. One of these systems, namely VLA-38 was shown to have an X-ray light curve with a periodicity of 12.32 h similar to Ter5ar. The position obtained from radio timing (; ) is consistent with the position quoted for VLA-38 in radio imaging and CXOU J174804.63−244645.2 from X-ray imaging (as reported in Urquhart et al. 2020) all to within 0.6 arcsec (see Fig. 2). On fixing the timing position to the VLA-38 source coordinates we were still able to hold phase connection. Considering this along with the matching orbital period from X-ray as well as radio observations, we can unambiguously link Ter5ar with VLA-38.

4.5. Ter5as

Ter5as (PSR J1748−2446as) was discovered in Epoch 1 at a spin period of 2.32 ms. Although it could not be detected in Epoch 2 after searching in a restricted DM range, it was visible in multiple epochs with GBT data. All the barycentric spin periods were consistent with the pulsar being an isolated system. Thus, obtaining more detections across all the epochs was relatively simple compared to the binary pulsars discussed earlier. Consequently, we were able to manually obtain a phase coherent timing solution spanning 19 yr.

The timing solution revealed significant proper motion in right ascension (μ_α = −2(1) mas yr⁻¹) as well as in declination (μ_δ = −12(6) mas yr⁻¹). The timing solution also yielded a high first order spin frequency derivative value that is ḟ = −4.7295(1) × 10⁻¹⁴ Hz s⁻¹ as well as a significant term. This could explain a significant drift in spin frequency with time thus reducing the usefulness of stacking across long time span datasets. The full timing solution is given in Table 4 and the post-fit timing residuals are plotted in Fig. 4.

4.6. Ter5at

Ter5at (PSR J1748−2446at) was discovered in Epoch 1 in a 30 min segment at a spin period of 2.188 ms and at an acceleration of 0.43(11) ms⁻². Examining the other 30 min segments revealed three more detections. Two segments showed the pulse to be fading indicating a possible ingress and egress of a radio eclipse and suggesting another spider-type system. Using these detections we were able to derive an initial orbital solution of P_b ≃ 0.22 d and x ≃ 0.1 lt-s implying a M_c, min ∼ 0.03 M_⊙. Folding the entire Epoch 1 time span with this ephemeris revealed a radio eclipse, confirming that Ter5at is an eclipsing black widow system.

We then extracted TOAs from the detections in Epoch 1 and Epoch 2 to get a better constrained timing solution. Using similar strategies to those applied to Ter5aq, we obtained multiple detections of Ter5at in all the datasets spanning 19 yr. We again used APTB to obtain a fully phase connected timing solution for Ter5at across a 13 yr time span (excluding SPIGOT data). We then extracted TOAs for the remaining 6 yr to extend the timing baseline to 19 yr. Owing to noticeable changes in the residuals from orbital variability, we switched to the BTX timing model and included two orbital frequency derivatives (see Table 3).

The timing solution yielded a significant value of proper motion (μ_α = −1.93(9) mas yr⁻¹; μ_δ = −7(3) mas yr⁻¹) and needed the first three spin frequency derivatives. Assuming that Eqs. (6) and (8) hold, we observe that the expected Ṗ_b value from the cluster acceleration is a factor of 3 lower than the measured value suggesting other effects may be at play. Additionally, we observed that Ter5at is located just 0.29 arcsec from the X-ray source CXOU J174805.37−244646.7 suggesting a possible association, which is supported by the eclipsing nature of the system.

4.7. Ter5au

Ter5au (PSR J1748−2446au) was discovered in Epoch 1 in the full 4 h observation at a spin period of 4.54 ms and at an acceleration of 0.113(3) ms⁻² and was also independently detected in Epoch 2.

Similar to Ter5ap, Ter5au was initially detected in just a few GBT epochs and required more detections to obtain a reasonable orbital solution. We thus used the 5 follow-up campaign epochs (that is Obs ID 08L-orb to 12L-orb in Table 1) to search for Ter5au and obtained three detections. Using the MeerKAT detections and the initial GBT detections, we were able to obtain a preliminary orbital solution using fitorb.py yielding P_b ≃ 5.97 d and x ≃ 6.55 lt-s with a hint of eccentricity (e ∼ 0.02) implying M_c, min ∼ 0.28 M_⊙. Using this orbital solution as an input for binary_info.py, we were able to get several detections across the GBT data span. Most detections showed the pulses were not drifting in rotational phase. We then generated TOAs for all the data and were able to manually phase connect all the available data. The full timing solution is given in Table 2.

We were able to detect proper motion in right ascension alone (μ_α = −2.4(2) mas yr⁻¹), and also a significant Ṗ_b,obs = (−16.357 ± 1.186) × 10⁻¹² s s⁻¹. Using Eq. (8), we obtained the intrinsic spin period derivative Ṗ = 3.72 × 10⁻²⁰ s s⁻¹ which is a typical value of the MSP population.

Furthermore, we obtained a significant detection of implying a total mass of M_tot = 1.83 ± 0.07 M_⊙. The measured mass function and the nominal value of M_tot would imply that M_p < 1.53 M_⊙ and M_c > 0.30 M_⊙. This indicates that the companion is not a He WD: for the orbital period of this system, Tauris & Savonije (1999) predict a WD mass of ∼0.24 M_⊙. This means that the companion is more likely to be a relatively light carbon-oxygen (CO) WD. The relatively large mass of the WD companion, the spin of the pulsar and the orbital period make the system look very similar to PSR J1614−2230, which is thought to have evolved via case A Roche lobe overflow (RLO, Tauris et al. 2011) when the donor star (that is the progenitor of the companion) was still in the main sequence stage.

4.8. Ter5av

Ter5av (PSR J1748−2446av) was discovered in Epoch 1 at a spin period of 1.85 ms and an acceleration of −0.1642(3) ms⁻² and was independently detected in Epoch 2. We were then able to get a few detections in GBT data using the DM as a constraint. However, these detections were not enough to get a well constrained orbital solution. We therefore searched the follow-up campaign epochs (similar to the strategy applied for Ter5ap and Ter5au) and obtained three more detections.

Using these GBT and MeerKAT detections, we determined a preliminary orbital solution with fitorb.py of P_b ≃ 3.38 d and x ≃ 1.25 lt-s implying M_c, min ∼ 0.07 M_⊙. The binary_info.py script was used with the orbital solution to help yield many more GBT detections. Another iteration with fitorb.py, after including the new detections, yielded a low non-zero eccentricity (e ∼ 0.0002). This pulsar had a broad profile with faint detections across most of the GBT data, resulting in low timing precision and difficulty in obtaining phase connection.

To aid the timing process, we used SeeKAT to get a better position for Ter5av. During one of the five follow-up campaign epochs, we placed seven beams that were hexagonally packed with the central beam centered on the position of the best beam detection in Epoch 1. This beam tiling pattern resulted in a robust position from SeeKAT with sub-arcsecond precision. We generated TOAs for all the data and used the newly obtained position as a starting point for attempting phase connection. This step was integral to successfully phase connect the data manually. The final timing solution is presented in Table 4.

We were able to detect proper motion in right ascension alone (μ_α = −1.7(2) mas yr⁻¹). We were also able to detect Ṗ_b = (−10.122 ± 1.909) × 10⁻¹² s s⁻¹ and using Eq. (8), we obtain the intrinsic spin period derivative value of Ṗ = 2.3 × 10⁻²⁰ s s⁻¹ with a derived surface magnetic field of B = 2.05 × 10⁸ G and characteristic age of τ_c = 1.3 Gyr. We also identified a potential association with X-ray source CXOU J174805.05−244643.1 which is within 0.43 arcsec of the position of Ter5av.

4.9. Ter5aw

Ter5aw (PSR J1748−2446aw) was discovered in Epoch 1 in a 1-h segment with a spin period of 13.04 ms, an acceleration of 5.77(2) ms⁻² and confirmed in Epoch 2. Searching 1-h long data segments across a few epochs of GBT data revealed several detections. These detections helped obtain a starting orbital solution with P_b ≃ 0.73 d and x ≃ 3.72 lt-s implying M_c, min ∼ 0.8 M_⊙. With no significant detection of eccentricity, we suspected that the companion is likely a high mass white dwarf. Using this starting orbital solution, we were able to generate TOAs for all GBT data. We then proceeded with manual phase connection and got a fully phase connected timing solution for all the GBT and MeerKAT data. The final timing solution is provided in Table 4.

The 19-yr timing solution yielded a significant detection of proper motion (μ_α = −2.2(1) mas yr⁻¹; μ_δ = −14(3) mas yr⁻¹). Furthermore, we also detected spin frequency derivatives up to the fourth order, similar to Ter5ao (see Table 4). The timing solution also yielded the detection of implying a total mass of M_tot = 2.12 ± 0.61 M_⊙.

We also detected the lowest order orbital period derivative Ṗ_b,obs = (5.924 ± 0.011) × 10⁻¹² s s⁻¹ which is slightly lower than the upper limit derived from Ṗ (6.139 × 10⁻¹² s s⁻¹). Using the measured Ṗ_b value we obtained estimates for the intrinsic spin period derivative (Ṗ = 8.341 × 10⁻²⁰ s s⁻¹), surface magnetic field (B = 1.05 × 10⁹ G) and characteristic age (τ_c = 2.48 Gyr) which are consistent with the population of pulsars which are mildly recycled.

Finally, we obtained a hint of Shapiro delay in this system: the orthometric amplitude in our DDFWHE solution (h₃ = 2.13 ± 0.88 μs) has a low (2.4-σ) significance. For this reason, and also because of the relatively low precision of , no precise masses can be derived yet (see discussion in Sect. 5). However, with its relatively slow spin period and massive WD companion, the system strongly resembles PSR J1952+2630, a pulsar–massive WD system (Gautam et al. 2022b).

4.10. Ter5ax

Ter5ax (PSR J1748−2446ax) was discovered in Epoch 2 in the full 4 h data span at a spin period of 1.94 ms and an acceleration of 0.005(1) ms⁻². A search in Epoch 1 revealed no detection. However, a blind search in GBT data revealed a few more detections. But the total number of detections and the wide spacing between the observation epochs was insufficient to solve for a unique orbital solution. We thus searched the data from the follow-up campaign and obtained detections in 4 out of 5 epochs.

Using all the MeerKAT detections and the few GBT detections, we were able to obtain a preliminary orbital solution with P_b ≃ 30.20 d and x ≃ 14.31 lt-s implying a M_c, min ∼ 0.20 M_⊙, thus suggesting a He WD companion. This preliminary orbital solution served as an input to binary_info.py and helped obtain many more detections in GBT data. Refitting for the orbital solution yielded a small non-zero eccentricity (e = 0.0091(2)). We then used this solution to re-fold all the data back till 2009 and were able manually to phase connect all the data. This solution was then used to generate TOAs from SPIGOT data back till 2004 and we were eventually able to obtain a 19 yr phase connected timing solution (Table 2).

We were able to detect significant proper motion along both directions (μ_α = −2.7(2) mas yr⁻¹; μ_δ = −18(6) mas yr⁻¹). We also measured spin frequency derivatives up to the fourth order (see Table 2). Without including the SPIGOT data, there was a significant detection of implying a very uncertain total mass of M_tot = 2.30 ±  1.06 M_⊙. However, on including the SPIGOT data, suggesting a total mass of M_tot ∼ 6.68 M_⊙ if all the contribution to is considered to be relativistic. Unless the system has a very low orbital inclination, the contribution to most likely comes from additional components apart from GR. Given that the detection after including SPIGOT data is at least 5 sigma significant, it points to an unresolved discrepancy.

To further understand this, we write down the relativisitic and classical effects that can cause secular changes in the observed rate of periastron advance (). This is given by (see e.g. Lorimer & Kramer 2004):

(13)

where corresponds to the contribution from the relativistic effects within the orbit (Eq. (3)), is the contribution from proper motion, from spin-orbit coupling and is the contribution from a secondary outer companion.

However, the pulsar is faint in the less sensitive SPIGOT data and RFI could be be causing unexpected issues. Additionally, the measurement of before and after adding the SPIGOT data are barely compatible at a 2-σ level. The timing data for this pulsar is currently being further investigated to understand this better.

5. Mass estimates

5.1. Bayesian map

In order to investigate the mass constraints on the systems where we measure (that is Ter5ao, ap, au, and aw), we made a map of the quality of fit (the χ²) for a grid of M_tot and cos i. This uniform grid amounts to a uniform prior on the total mass, which is constrained by the detection of omega-dot via Eq. (3), and an isotropic prior for the orbital axis. For each point in this grid, we calculate i and then, using Eq. (2), we calculate M_c.

We used the DDGR timing model, which is similar to the DD model but assumes that GR is the correct theory of gravity (Damour & Deruelle 1986). For each point in the M_tot–cos i grid, the M_tot and M_c values are fixed and all other parameters are allowed to vary. The resulting χ² describes how well the values of M_tot and M_c can describe the timing data.

After this stage, we use the Bayesian technique described by Splaver et al. (2002). The likelihood is calculated from the χ² using

(14)

where X represents the set of all TOAs with uncertainties from all epochs and is the smallest value of χ² in the whole grid (corresponding, therefore, to the best fit, which will therefore have the highest probability density). Using Bayes’ theorem, we then derive 2-dimensional joint posterior probability distribution functions (pdfs) for M_c and cos i (see contours in left main panels in Fig. 8), for M_c and M_p (right main panels in Fig. 8), and 1-D pdfs for M_c, M_p, cos i, (see lateral panels) and M_tot.

Fig. 8. Mass-mass diagrams for Ter5ao (top), Ter5ap (middle), and Ter5aw (bottom). The main panels show Mc as a function of cos i (left) and Mp (right). In the left panels, the regions marked in gray are excluded because they imply a negative pulsar mass. For the panels on the right, the gray regions are constrained by the mass function and sin i < 1. The black contours include 68.3, 95.4 and 99.7% of all probability in the 2-D joint posterior probability distribution functions (pdfs) of each panel. The red lines indicate constraints derived from the nominal values of the total mass derived from the nominal value of and the ±1-σ uncertainties of its measurement. The minimum companion masses for Ter5ao and Ter5ap are constrained by and the limit of sin i < 1 giving Mc > 0.8 M⊙ and Mc > 0.278 M⊙ respectively, this results in Mp < 2.23 M⊙ and Mp < 1.72 M⊙ respectively. The blue lines represent the constraints derived from the nominal and ±1-σ limits of the orthometric amplitude of the Shapiro delay, h3. The orange lines designate the 2 and 3-σ upper limits derived for Ter5ap. The top panels depict the normalised 1-D pdfs for cos i, Mp and, on the side, Mc. The medians of the pdfs are depicted by the solid black lines, the dotted lines indicate successive ±1 and 2-σ equivalent percentiles around the median.

5.2. Results

The χ² maps confirm the fact mentioned earlier that there is still no significant detection of the Shapiro delay for any of the systems: indeed, in the top left side panels of the distributions obtained for both Ter5ao and Ter5aw, we see significant probability tails for the lower inclinations. For Ter5aw, all mass values are at this stage too imprecise to be astrophysically useful: the companion mass is consistent with what is expected for a CO WD (see e.g. Shamohammadi et al. 2023), but the mass is too uncertain for any firm conclusions.

Our χ² map for Ter5au (not displayed in Fig. 8) confirms that, apart from the well-measured M_tot, there are no additional constraints on cos i, which has a nearly flat pdf (the median for i is 60 deg, as we would expect from the assumed prior).

However, for Ter5ap, there is a significant decrease of the probability for lower inclinations. This is caused not by a detection of the Shapiro delay, by the small value of the Einstein delay measured for this system, γ = −4.8 ± 4.4 ms, which excludes large companion masses, as we can see in the middle plot of Fig. 8. This small measurement of γ cannot be due to a measurement of ẋ_μ that is being interpreted as γ. Evaluating Eqs. (9) and (25) of Ridolfi et al. (2019), we find that, for Ter5ap; ẋ_γ ≫ ẋ_μ, which implies that the effect of γ should be dominant.

The opposite is true for Ter5ao, where the measurements of ẋ do constrain the orbital orientation, as described in Sect. 4.1. However, for Ter5ao, the constraints on ẋ_μ were not taken into account in the Bayesian mass estimates, as they would require a 3-D map in M_c, cos i and orbital orientation (Ω). The results of such a map would be qualitatively similar, with the further exclusion of a small range of orbital inclinations close to 90 deg, but no additional constraints for lower inclinations. Quantitatively they would also be similar, since ẋ is not measured with high significance. This means that for this pulsar there are no solid constraints on the individual masses. The deduced constraints on the individual masses and orbital inclination of Ter5ao, Ter5ap, Ter5au, and Ter5aw are summarised in Table 5.

6. Discussion and future prospects

The ten new discoveries made in Ter5 not only show a wide range of properties but also raise open questions, particularly pertaining to formation and stellar evolution channels.

6.1. Eccentric binaries

From the large measured , we know that Ter5ao possesses a large binary mass (3.166 ± 0.024 M_⊙), a large minimum companion mass (M_c, min = 0.93 M_⊙) and significant eccentricity (e = 0.32). In the Galactic field, the most likely possibility would be that the system is a double neutron star binary, where a second supernova explosion from the progenitor of the companion has induced the observed eccentricity (Tauris et al. 2017). However, given the pulsar’s spin period (P = 2.27 ms), small intrinsic Ṗ (1.13 × 10⁻²⁰ s s⁻¹) and small magnetic field (1.6 × 10⁸ G), a prolonged episode of mass accretion from a low-mass companion via Case B RLO (Tauris 2011) would be required, or potentially Case A (Tauris et al. 2011). This implies that the companion star was relatively light and gradually evolved into a low-mass WD in a circular orbit, as seen for most such systems in the Galaxy (see e.g. review by Tauris 2011; Tauris et al. 2011; Tauris & van den Heuvel 2023).

However, because of the large stellar encounter rate in Ter5 – in particular the large stellar encounter rate per binary (γ) that Verbunt & Freire (2014) predicted for this GC – it is possible that this low-mass WD companion was replaced by a more massive degenerate object in a secondary exchange encounter. Given the chaotic nature of this process, it almost invariably results in highly eccentric orbits (see Sect. 5 in Phinney 1992). Many such eccentric and high mass companion systems already exist in Ter5 like J, Q, U and ai (Ransom et al., in prep.) as well as globular clusters with a high γ (see e.g. Balakrishnan et al. 2023, and references therein).

Given this possible origin of the companion of Ter5ao, we cannot infer anything about its nature based on considerations of stellar evolution. This degenerate companion could be another NS; if confirmed, it would make Ter5ao simultaneously the fastest spinning pulsar in any known DNS system and the widest orbit for any known DNS. This could be confirmed either with the measurement of a large mass, or the detection of radio pulsations from the companion of Ter5ao, the latter of which is currently being investigated. The hypothesis of an equal-mass DNS implies a relatively low orbital inclination of 35 deg.

If the system is closer to edge-on (i ∼ 90 deg), then the mass of the pulsar could be up to 2.23 M_⊙, and in this case the companion would have a mass of only 0.93 M_⊙ – very likely a massive WD. This pulsar mass would be larger than the largest well-measured mass of a neutron star previously obtained from PSR J0740+6620 (2.08 ± 0.07 M_⊙; Fonseca et al. 2021) and would serve as an excellent test bed for constraining the equation of state of super-dense matter (Özel & Freire 2016). However, such high orbital inclinations would result in a detectable Shapiro delay. Our current measurement is not precise enough to claim such a detection. Currently we do not have a significant detection of ẋ for Ter5ao. We note, however, that a detection of a large value of ẋ from future observations could exclude large orbital inclinations.

For Ter5ap, we can introduce tighter constraints on the individual masses, from the non-detection of the Einstein delay γ_E. The pulsar mass (95.4% C.L.) is the largest measured for a pulsar in a globular cluster, the previous one being PSR J1910−5959A (; Corongiu et al. 2023). An exchange encounter is also a possibility for the formation of this highly eccentric (e = 0.905) system. Among recycled pulsars, the eccentricity of this system is only second to NGC 6652A (DeCesar et al. 2015). The latter authors suggested that the most suitable explanation for the high eccentricity and massive companion (M_c, min = 0.73 M_⊙) of NGC 6652A is the origin of the system in an exchange encounter, as discussed in the case of Ter5ao. However, for Ter5ap, the companion is much lighter (M_c = 0.294 M_⊙). This means that the latter system could either have formed in an exchange encounter, or it could have maintained its original He WD companion, with the eccentricity being later raised by gravitational perturbations from encounters with nearby stars. The time required to induce a certain eccentricity (assuming e > 0.01) can be quantified as (see Eq. (5) in Rasio & Heggie 1995)

(15)

where n is the number density of stars near the pulsar, v is the one dimensional velocity dispersion in the core, P_b is the orbital period and e is the observed eccentricity. We take v = 15.6 km s⁻¹ for Ter5 from Baumgardt & Vasiliev (2021). In order to derive n we make assumptions similar to Lian et al. (2023). They used n ∝ ρ_c where ρ_c is the core luminosity density of the globular cluster. Using the values for NGC 5024 as a reference (as given in Lian et al. 2023 for convenience) and ρ_c = 1.38 × 10⁵ L_⊙ pc⁻³ for Ter5 as given in Harris (2010), we get n = 3.3 × 10⁵ pc⁻³. Applying these values to Ter5ap, we can estimate the time it would take for close encounters to raise the eccentricity to this value: t_{> e, ap} ∼ 0.94 Gyr, which fits well within the age of the cluster (12 Gyr; Ferraro et al. 2016).

6.2. Less eccentric binaries

Apart from Ter5ao and ap, the other Ter5 discoveries show a range of eccentricities from essentially zero for the eclipsing systems to 0.025 for Ter5au. These small eccentricities indicate that none of these companions were exchanged, and because of this we can make some inferences about the nature of these systems from basic considerations of stellar evolution.

For instance, from the mass functions, and total mass measurements, we can say that Ter5au and Ter5aw likely have CO WD companions. As discussed above, the orbital and spin periods make the Ter5au system remarkably similar to PSR J1614−2230, which likely evolved through case A RLO (Tauris et al. 2011). Other systems like PSRs J1125−6014 (Shamohammadi et al. 2023) and J1933−6211 (Geyer et al. 2023) are very similar and likely had a similar origin. The only difference about Ter5au is the much larger eccentricity, likely caused by its location in a dense GC. A measurement of the individual masses of Ter5au would be important to confirm the nature of the companion mass. For Ter5aw, the companion mass is consistent with other CO WD companion pulsar binaries known in the Galactic field (McKee et al. 2020; Shamohammadi et al. 2023). Moreover, this system is remarkably similar, in its companion mass and orbital and spin periods, to PSR J1952+2630 (Gautam et al. 2022b).

Given their lower masses, Ter5av and ax likely have He WD companions, with Ter5av likely seen at a low orbital inclination. None of these systems shows eclipses. Their eccentricities (in particular, that of Ter5ax, e ∼ 0.009) are still significantly larger (in the case of Ter5av, by three orders of magnitude) than those predicted by the relation between eccentricity and orbital period for He WD companions (Eq. (7.35) from Phinney 1992 also known as the Phinney relation). This eccentricity could be induced from exchange encounters or from close flybys of nearby stars as discussed in Sect. 6.1. For the latter case and assuming low eccentricities (e < 0.01), Eq. (15) is slightly modified as (see e.g. Rasio & Heggie 1995; Lian et al. 2023)

(16)

For Ter5ax, we use Eq. (16) to derive t_{> e, ax} ∼ 0.3 Gyr. This is consistent with the estimated age of Ter5 (12 Gyr, Ferraro et al. 2016). For the other He WD system, Ter5av, we get a similar value t_{> e, av} ∼ 0.3 Gyr. A similar argument can also be applied to Ter5au (which also has a large e compared to pulsar – CO WDs observed in the disk), this yields t_{> e, au} ∼ 1.31 Gyr, also consistent with the age of the cluster. Ter5av has a relatively lower eccentricity but still higher than expected from the Phinney relation (Phinney 1992).

Among the discoveries showing similar properties to spider systems, Ter5aq (P = 12.52 ms) seems to be an outlier. Figure 9 shows the spin period vs minimum companion mass for all known black-widow pulsars. The spin period of Ter5aq is much slower than the rest of the black widow population including the Galactic as well as globular cluster pulsars. Although there are relatively slow spinning redback systems like Ter5A (P = 11.56 ms), the relatively lower companion mass of black-widow systems suggests that the mass loss is higher and the pulsar could be more recycled and hence spin faster in black-widow systems. However, given that the companion material could undergo evaporation, this may not always be true. Given that there is no constraint on the intrinsic spin period derivative, it is currently difficult to comment further on the recycled nature of Ter5aq.

Fig. 9. Spin period plotted as a function of the minimum companion mass for all known black widow systems. No clear distinction is seen between these pulsars whether they are in the Galactic field or within globular clusters. Ter5aq stands out given its relatively slow spin period amongst the rest of the population.

While Ter5ar and Ter5at are most likely spider systems especially given their eclipsing nature, the nature of the companion of Ter5av is debatable. Firstly, we do not observe eclipses with Ter5av. Despite this, if assumed to be a spider with a high orbital inclination, the minimum companion mass (M_c, min ∼ 0.07 M_⊙) places it in between the black widow and redback populations. However, the relatively high orbital period of Ter5av suggests it is not a spider system but instead has a He WD companion (see Fig. 10) with a low inclination angle.

Fig. 10. Orbital period plotted as a function of minimum companion mass for all binary pulsars. The information was obtained from the ATNF catalogue (Manchester et al. 2005). The different colours represent the different types of companions. All Ter5 binaries discovered in this work are marked in blue and their respective letters are also indicated.

Ter5at and Ter5ar are eclipsing black widow and redback systems respectively whose spin period, orbital period and derived minimum companion masses are consistent with the known spider population (see Fig. 10). The discovery of Ter5ar and its corresponding association with VLA-38 from Urquhart et al. (2020) as well as the matching orbital period between radio and X-ray data demonstrates the importance of continued synergies between imaging and time-domain searching through a multiwavelength approach. Furthermore, identifying steep spectrum sources in radio imaging similar to Urquhart et al. (2020) are particularly useful with interferometers like MeerKAT where narrow synthesised beams could be placed on these positions. Besides this, apart from Ter5ar, three other discoveries show likely associations with X-ray sources. All except Ter5av show eclipsing properties further cementing the well known association of spider systems with bright X-ray sources, as mentioned earlier. It also provides a platform to better understand the properties of these systems from a multi-wavelength approach.

6.3. The pulsar population of Ter5

As a collective, the ten new discoveries offer scope on understanding several aspects of globular cluster properties. Nine out of the ten are binary systems making the overall fraction of binaries in Ter5 as 59%. This is the sixth highest percentage among all globular clusters with more than 10 known pulsars following M62, NGC 362, M28, 47 Tuc, and NGC 1851. M62 (which has 10 known pulsars) is a unique case where all known pulsars are in binaries. Moreover, the evolution history and core-collapsed nature of M62 is still a matter of open debate (Vleeschower et al. 2024). NGC 362 and NGC 1851 are known to have compact cores thus making the high fraction of binaries less surprising. On the other hand, Ter5 is established as a bulge cluster similar to M28 (e.g. Ferraro et al. 2009). Ferraro et al. (2016) suggested that Ter5 likely underwent multiple episodes of star formation in its history given the trimodal distribution of stars based on iron content. This in turn could explain why no other globular cluster currently possesses more confirmed pulsars than Ter5. Moreover, the high mass and density of Ter5 contribute to a high stellar encounter rate which in turn encourages a large fraction of binary pulsars to be formed. It is also worth noting that selection effects from the types of searches conducted so far play a major role in the binary fraction. For example, at least five isolated MSPs in Ter5 have been discovered by specifically applying stack searches (Ransom et al., in prep.). Such a method applied to other globular clusters could also boost their respective isolated pulsar numbers.

Prager et al. (2017) already undertook a study using 34 pulsars to help place constraints on a possible black-hole in the core of Ter5 and better constrain the structural properties of the cluster. Their results favoured the argument that Ter5 is a fragment of the Galactic bulge rather than a remnant of a dwarf galaxy. Given the addition of 15 more pulsars since, it is worth conducting a similar study to improve the constraints and confirm these findings. Finally, a large collection of pulsars in Ter5 could also potentially probe the presence of ionised gas in the intra-cluster medium, similar to the work done on 47 Tuc by Freire et al. (2001b). However, a similar study for Ter5 in the past has yielded ambiguous results owing to high DM in the foreground unlike 47 Tuc. Nevertheless, multiple pulsars at multiple DMs and positions can better ascertain the presence or absence of ionised gas within the cluster. Martsen et al. (2022) already obtained precise rotation measure (RM) values from 28 pulsars to deduce constraints on the parallel component of the magnetic field component along the line of sight to Ter5. Polarisation studies of the latest discoveries could help further update these constraints.

6.4. Prospects

There are various reasons to continue searching for pulsars in Ter5. Firstly, multiple simulations using different methods have predicted a larger population of pulsars yet to be discovered. Bagchi et al. (2011) conducted Monte-Carlo simulations and modelled the observed population of pulsars as the tail end of a broader intrinsic luminosity distribution function. They used a log-normal distribution and estimated ∼150 pulsars to be detectable in Ter5. Chennamangalam et al. (2013) built on this work and used Bayesian statistics to constrain the luminosity function parameters. In the best case scenario, they predicted more than 200 detectable pulsars. Adding to this, Fig. 11 is adapted from Fig. 3 in Martsen et al. (2022) and shows the number of pulsars brighter than a certain luminosity value as a function of the pseudo-luminosity. The figure has been updated with the luminosity values of the ten newly discovered pulsars and it still shows a turnover suggesting that we are approaching a minimum luminosity cutoff and there are potentially some more pulsars to be found. Its noteworthy that half of these discoveries are fainter than any previously known pulsars.

Fig. 11. Cumulative number of pulsars with pseudo-luminosity greater than a given pseudo-luminosity (L1284) as a function of the pseudo-luminosity. This figure is adapted from Fig. 3 in Martsen et al. (2022). We have used the known pulsar values as quoted in Martsen et al. (2022) and added the radiometer flux based estimates (see Sect. 3.7) for the ten new pulsars. It is evident that the linear trend tapers when going from right to left.

Secondly, the searches conducted so far have multiple limitations. All the discoveries made with MeerKAT so far have resulted from searches limited to the beams within the core of the cluster alone. Searches in the outer beams especially within the half-mass and half-light radius still have the potential to boost discovery numbers. These pulsars in the outer regions of the cluster are expected to have lower mass companions if in binaries (owing to mass segregation) or isolated pulsars that have been flung out from interactions. It is also worth noting that previous searches have used several datasets to stack Fourier power spectra to boost the S/N of possible pulsar candidates (Cadelano et al. 2018). While being effective for isolated pulsars, binary pulsars however would undergo significant drift in spin frequency leading to leakage of Fourier power over several bins and in turn reducing the effective S/N. It was thus assumed that any isolated pulsars hidden in Ter5 data would be recovered from the stack searches. However, the discovery of Ter5as proved otherwise. This is due to the significant drift in spin frequency with time (in turn reducing the stacked Fourier power) caused by the spin frequency derivatives. A significant fraction of discoveries with eclipsing properties also showed that there are possibly other pulsars yet to be seen owing to observing at unsuitable orbital phases.

The searches conducted on MeerKAT data so far have also been limited to acceleration searches down to 30 min segments. This implies that the searches are reasonably sensitive to binary pulsars with orbital periods of the order of hours rather than minutes. Applying jerk searches (Andersen & Ransom 2018) or coherent template bank based searches (Allen et al. 2013; Balakrishnan et al. 2022) can thus provide a possibility to find these highly compact binaries whose orbital period is of the order of few minutes. Such binaries can provide an excellent platform for testing GR in stronger gravitational fields than those of known binary pulsars. We are currently working with MeerKAT data to search for such systems using Einstein@Home (Anderson et al. 2006), a volunteer distributed computing project that has already been successful in the past for discovering several radio pulsars (see Knispel et al. 2015, and references therein).

Ter5 will also be observed in the near future with the recently installed S-Band (1.7–3.5 GHz) receiver system (Kramer et al. 2016; Padmanabh et al. 2023) at MeerKAT. This frequency band has historically been successful with the GBT primarily owing to the lower impact of dispersion compared to L-Band. The enhanced sensitivity of MeerKAT would only boost the possibility of finding several more pulsars in Ter5.

Thirdly, the detection of systems like Ter5ao, which are clearly the products of secondary exchange interactions, implies that other such systems might be found in Ter5. These could include MSP–MSP binaries, or even MSP–stellar mass black hole binaries, which would represent new laboratories for tests of gravity theories.

Finally, the long timing baselines of the ten new pulsars has only been possible from the use of archival data from the GBT. Thus any further discoveries made would also benefit from the already existing rich archival dataset for obtaining quick and robust timing solutions in the future.

7. Conclusions

We have discovered ten new pulsars using the MeerKAT telescope in the Ter5 globular cluster as part of the TRAPUM globular cluster survey. This has brought the total number of known pulsars in this cluster to 49, the highest for any globular cluster known. We also presented long-term timing solutions nearing two decades for nine of these pulsars mainly using archival GBT data. These include astrometric parameters like proper motion as well as PK parameter measurements including Ṗ_b that helped derive intrinsic spin period derivative values for several pulsars.

Ter5ao is an eccentric, wide-orbit pulsar with a large minimum companion mass and a rapid spin period, suggesting that it is the result of a secondary exchange encounter. We were able to detect for this system, which yields a total mass of 3.17 ± 0.04 M_⊙ to 95% C.L. The system is either a DNS, or has a massive WD companion. The pulsar could be quite massive (> 2 M_⊙), but at the moment we have no constraints on the individual component masses. Ter5ap has the second highest eccentricity for any recycled pulsar known (e = 0.905) and the highest mass measured for a pulsar in a globular cluster, but unlike other similar systems, its companion is relatively light. For this reason, we cannot conclude that it is the result of an exchange encounter. Two systems, Ter5aq and Ter5at, are black widow systems confirmed by the detection of radio eclipses. Ter5ar is an eclipsing redback system which is associated with the radio counterpart to the source named as VLA-38 by Urquhart et al. (2020). Ter5as is the lone isolated pulsar discovery among the ten pulsars. Ter5au, Ter5ap, Ter5av, and Ter5ax likely have WD companions, but their orbital eccentricities have most likely been acquired from gravitational perturbations from the dense surrounding environment within the core of Ter5. We have a potential Shapiro delay in Ter5aw (currently with 2.4-σ significance) leading to individual mass measurements of and to 68.3% C.L., suggesting a rare CO WD companion. The wide variety of pulsars found in Ter5 demonstrates the richness of this cluster, not only to our understanding of stellar evolution but also for probing the intricacies of globular cluster dynamics and environments.

Acknowledgments

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time space reference systems department of the Paris Observatory. MeerTime data is housed on the OzSTAR supercomputer at Swinburne University of Technology. PTUSE was developed with support from the Australian SKA Office and Swinburne University of Technology. 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. The authors also acknowledge MPIfR funding to contribute to MeerTime infrastructure. TRAPUM observations used the FBFUSE and APSUSE computing clusters for data acquisition, storage and analysis. These clusters were funded and installed by the Max-Planck-Institut für Radioastronomie and the Max-Planck-Gesellschaft. P.V.P., P.C.C.F., A.R., F.A., E.D.B., W.C., D.J.C., C.J.C., A.D., M.K., V.V.K. acknowledge continuing valuable support from the Max-Planck Society. S.M.R. is a CIFAR Fellow and is supported by the NSF Physics Frontiers Center award 1430284. M.E.D. acknowledges funding from the National Science Foundation Physics Frontier Center award No. 2020265. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. L.Z. acknowledges financial support from ACAMAR Postdoctoral Fellowship and the National Natural Science Foundation of China (Grant No. 12103069). The research activities described in this paper were carried out with contribution of the NextGenerationEU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures, Project IR0000034 – “STILES – Strengthening the Italian Leadership in ELT and SKA”. V.V.K. acknowledges financial support from the European Research Council (ERC) starting grant, “COMPACT” (Grant agreement number 101078094), under the European Union’s Horizon Europe research and innovation programme. M.Ba. acknowledges support through ARC grant CE170100004. B.W.S. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694745). A.P., A.R. and M.Bu. gratefully acknowledge financial support by the research grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016. A.P. also acknowledges the support from the “Italian Ministry of Foreign Affairs and International Cooperation”, grant number ZA23GR03. L.V. acknowledges financial support from the Dean’s Doctoral Scholar Award from the University of Manchester. Finally, we thank the anonymous referee for their suggestions to bring clarity to the text.

References

All Tables

All Figures

Fig. 1. TRAPUM beam tilings of the two search observations (Epoch 1 and Epoch 2; see text) of Ter5, performed with MeerKAT on 05 Sep 2020 (left) and on 06 Jan 2021 (right) at L-Band. There are 288 beams in total per observation. Each individual beam has been overlaid with the corresponding point spread function. The vertical bar shows the reference gray scale for the fractional sensitivity with reference to the boresight position of the individual beam. The different concentric circles indicate the core radius (red), half-light radius (yellow-dashed) and half-mass radius (blue dotted). The known pulsars are marked in each tiling as black crosses. The discoveries made in each epoch are denoted as red stars along with their individual names. The green crosses indicate the positions of the discoveries but from the alternate epoch where the discovery was first made. The positions of the discoveries comes from timing.

In the text

Fig. 2. Ter5 field shown with each axis approximately spanning 1.5 times the core radius of the cluster. The positions of the X-ray sources (black rings), known pulsars (blue dots) as well new discoveries (red dots) are also overlaid. Particular emphasis is made on the VLA-38 source (brown star) obtained from radio imaging and whose coordinates are reported in Urquhart et al. (2020). The timing position of Ter5ar (green star) is shown to overlap significantly with VLA-38 as well as with an X-ray source (CXOU J174804.63−244645.2) in the backdrop. The dashed lines surrounding the X-ray sources and VLA-38 are the 3-σ positional uncertainties.

In the text

	Fig. 3. Integrated pulse profiles of all the newly discovered pulsars in Ter5. These plots are obtained after summing together individual profiles from different epochs and aligning the profiles with respect to a reference template profile.
In the text

Fig. 4. Post-fit residuals after fitting a timing model to the times of arrival (TOAs) obtained from GBT and MeerKAT data spanning about 19 yr. Each colour represents TOAs obtained from different back-ends of both the telescopes (as shown in the legend). The secondary y-axis also shows the corresponding residuals as a fraction of the rotational spin phase. These residuals are shown for Ter5ao, Ter5ap, Ter5aq, and Ter5as, respectively. The residuals for Ter5ar is depicted here for only 2 yr with phase connection.

In the text

	Fig. 5. Post-fit residuals as shown in Fig. 5 but for the remaining four pulsars with long-term timing solutions. These are Ter5at, Ter5au, Ter5av, Ter5aw and Ter5ax.
In the text

Fig. 6. Initial orbital solution for Ter5ap depicted. The left plot shows the line of sight acceleration as a function of orbital phase (which is in turn a function of the mean anomaly). The red line indicates the orbital solution model and the blue dots are data points. Note that the error bars on the data points are smaller than the marking symbol size. All points have a negative acceleration given that the pulsar spends very little time near periastron (roughly a day out of a 21.3 day orbital period). There are no positive acceleration data points which would better constrain the orbital solution. This is further demonstrated by showing the plot on the right depicting acceleration v/s the observed spin period. One data point clearly does not obey this initial solution.

In the text

	Fig. 7. Significant orbital variability in Ter5ar demonstrated. The difference between expected and observed measurement of epoch of periastron (T0) is shown as a function of time in MJD. A stable orbit would show a flat horizontal line but the plot shows the variations ranging from –20 s to 60 s, thus demonstrating the difficulty in timing such systems.
In the text

Fig. 8. Mass-mass diagrams for Ter5ao (top), Ter5ap (middle), and Ter5aw (bottom). The main panels show Mc as a function of cos i (left) and Mp (right). In the left panels, the regions marked in gray are excluded because they imply a negative pulsar mass. For the panels on the right, the gray regions are constrained by the mass function and sin i < 1. The black contours include 68.3, 95.4 and 99.7% of all probability in the 2-D joint posterior probability distribution functions (pdfs) of each panel. The red lines indicate constraints derived from the nominal values of the total mass derived from the nominal value of and the ±1-σ uncertainties of its measurement. The minimum companion masses for Ter5ao and Ter5ap are constrained by and the limit of sin i < 1 giving Mc > 0.8 M⊙ and Mc > 0.278 M⊙ respectively, this results in Mp < 2.23 M⊙ and Mp < 1.72 M⊙ respectively. The blue lines represent the constraints derived from the nominal and ±1-σ limits of the orthometric amplitude of the Shapiro delay, h3. The orange lines designate the 2 and 3-σ upper limits derived for Ter5ap. The top panels depict the normalised 1-D pdfs for cos i, Mp and, on the side, Mc. The medians of the pdfs are depicted by the solid black lines, the dotted lines indicate successive ±1 and 2-σ equivalent percentiles around the median.

In the text

	Fig. 9. Spin period plotted as a function of the minimum companion mass for all known black widow systems. No clear distinction is seen between these pulsars whether they are in the Galactic field or within globular clusters. Ter5aq stands out given its relatively slow spin period amongst the rest of the population.
In the text

	Fig. 10. Orbital period plotted as a function of minimum companion mass for all binary pulsars. The information was obtained from the ATNF catalogue (Manchester et al. 2005). The different colours represent the different types of companions. All Ter5 binaries discovered in this work are marked in blue and their respective letters are also indicated.
In the text

Fig. 11. Cumulative number of pulsars with pseudo-luminosity greater than a given pseudo-luminosity (L1284) as a function of the pseudo-luminosity. This figure is adapted from Fig. 3 in Martsen et al. (2022). We have used the known pulsar values as quoted in Martsen et al. (2022) and added the radiometer flux based estimates (see Sect. 3.7) for the ten new pulsars. It is evident that the linear trend tapers when going from right to left.

In the text