© ESO 2022

1. Introduction

Globular clusters (GCs) have historically been popular targets for radio pulsar surveys (e.g., Ransom 2008). This is due to their nature, which makes them extremely fecund pulsar seedbeds. GCs are ancient (ages of several billions of years) stellar systems that are kept together by their own gravity, resulting in a spherical distribution of stars. In their cores, the stellar number density can reach up to 10³⁻⁶ pc⁻³, that is up to ∼10⁶ times the typical values found in the Solar neighborhood.

In such crowded environments, it is fairly common that an old, “dead” neutron star (NS) gravitationally interacts with binary systems within the cluster. If the encounter is close enough, this can cause the ejection of the lighter member of the previous binary and the formation of a new binary system, consisting of a NS and a main sequence (MS) star (Hills 1975). The MS star then evolves, fills its Roche lobe, and starts transferring matter to the NS, which can thus be “recycled”. During this phase, the system is seen as a low-mass X-ray binary (LMXB). Many of these systems then evolve to form millisecond pulsars (MSPs, e.g., Alpar et al. 1982; Radhakrishnan & Srinivasan 1982; Papitto et al. 2013). The formation of LMXBs by dynamical means is the reason why there are, in unit of stellar mass, about three orders of magnitude more LMXBs and MSPs in GCs than in the Galactic disk (Clark 1975; van den Berg 2020).

Searches for pulsars in GCs are first and foremost limited by sensitivity: owing to the large distances (in most cases > 5 kpc) of these stellar systems, only the brightest pulsars can be detected. After a decade characterized by a slow discovery rate, mainly due to a sensitivity limit being reached at the existing radio telescopes, a new surge in the total number of pulsars in GCs has recently occurred. In the last three years, the total GC pulsar population has increased by more than 50%, reaching, at the time of writing, 235 known pulsars¹ in 36 different clusters. While a few of the new discoveries resulted from instrumentation upgrades made at the Parkes (Dai et al. 2020) and GMRT radio telescopes (Gautam et al. 2022), the vast majority were found by the two newly built major radio facilities: the Chinese Five-hundred-meter Aperture Spherical Telescope (FAST, Nan et al. 2011) and the South African MeerKAT radio telescope (Jonas & MeerKAT Team 2016; Camilo 2018). Despite their very different designs and total collecting areas, both FAST and MeerKAT are providing a major leap in the available raw sensitivity for pulsar searching experiments in their respective hemispheres. As a result, since its first GC observation taken in 2018, FAST has discovered over 30 new pulsars in GCs (Wang et al. 2020; Pan et al. 2020, 2021a,b; Qian & Pan 2021; Yan et al. 2021), some of which are in clusters with no previously known pulsars. Similarly, MeerKAT began observing GCs in 2019, with a first “census” campaign carried out in the context of the MeerTime² (Bailes et al. 2020) Large Survey Project (LSP). During this, a sample of nine promising GCs hosting previously known pulsars were observed. These observations used the Pulsar Timing User Supplied Equipment (PTUSE) backend, which, at the time, could record a single tied array beam on the sky. Hence those observations used only the central 44 antennas of the array, so as to trade off some sensitivity for an increased field of view. These data served a variety of purposes, among which were to: (a) test the telescope performance; (b) time and characterize previously known pulsars (Abbate et al. 2020a); and (c) perform a first search for new pulsars. The latter turned out to be successful, with eight new MSPs found across six different clusters (Ridolfi et al. 2021).

Besides being scientifically very fruitful, the MeerKAT GC census was also key in planning future observations and paved the way for a more thorough GC pulsar survey, which started in mid 2020 as part of the TRansients And PUlsars with MeerKAT (TRAPUM³) LSP (Stappers & Kramer 2016). Unlike the census observations, the TRAPUM GC pulsar survey uses the Filterbanking Beamformer User Supplied Equipment (FBFUSE) and the Accelerated Pulsar Search User Supplied Equipment (APSUSE) computing clusters as backend to form up to 288 tied-array beams on the sky, tiled radially out from the GC center using an optimal hexagonal tessellation. This means that we can use up to all of the 64 MeerKAT antennas, substantially increasing the sensitivity compared to the initial GC census. This survey also extends the sample of target GCs to a total of 28, and includes observations using different frequency bands, depending on the distance to the cluster. One of these is the relatively unexplored Ultra High Frequency (UHF) band. The latter, which at MeerKAT operates between 544 and 1088 MHz, is proving to be particularly useful, as shown in this work.

One of the GCs selected for the TRAPUM GC survey is NGC 1851, which is the main focus of this paper. This is a fairly compact cluster, having a core radius of r_c = 0.09 arcmin and a half-mass radius of r_hm = 0.85 arcmin (Miocchi et al. 2013). With center coordinates being at right ascension α_J2000 = 05^h14^m06.76^s and declination δ_J2000 = −40° 02′47.6″ (Harris 2010), it is located in the southern constellation of Columba at a distance of 11.95 ± 0.13 kpc (Baumgardt & Vasiliev 2021), far below the Galactic plane (at Galactic coordinates , ). Until recently, NGC 1851 was known to be the host of a single, but very peculiar binary pulsar, PSR J0514−4002A (also known as NGC 1851A). Discovered by Freire et al. (2004) at 350 MHz with the GMRT, and subsequently timed with the Green Bank Telescope (GBT) and the GMRT at similar frequencies (Freire et al. 2007; Ridolfi et al. 2019), NGC 1851A is a 4.99-ms pulsar in a wide and extremely eccentric (e = 0.89) orbit. This characteristic made it possible to precisely measure two post-Keplerian parameters, which in turn led to a precise determination of the system component masses, something that is rarely possible. NGC 1851A is most likely the result of a secondary exchange interaction, where the original low-mass star that recycled the pulsar was replaced by the current, massive (M_c = 1.22 M_⊙) companion. The existence of such a peculiar binary motivated our searches in this cluster, as it informed our suspicions that other such systems may be present therein.

Here we present the results of our search for new pulsars in NGC 1851 using MeerKAT observations made in the context of the TRAPUM GC pulsar survey, as well as follow-up observations taken within the MeerTime program. In Sect. 2 we describe the observations, their instrumental set-ups, and the data analysis methods. We present the results in Sect. 3 and discuss their implications in Sect. 4. Finally, in Sect. 5, we summarize our findings and examine possible scientific developments.

2. Observations and data analysis

We have observed NGC 1851 with MeerKAT on 13 occasions, from January to November 2021 (see Table 1 for the full list of observations), each time using at least 56 antennas. The first seven observations were carried out as part of the TRAPUM GC pulsar survey, hence, we recorded only the total intensity with the APSUSE computing cluster as the main backend. The first three of them were made using the L-band receivers, at a central frequency f_c = 1284 MHz, and with a bandwidth Δf = 856 MHz. The latter was split into 4096 channels and sampled every 76.56 μs. The other four TRAPUM observations were made with the UHF receivers, at a central frequency f_c = 816 MHz, and with Δf = 544 MHz, also split into 4096 channels but sampled every 60.24 μs. In all of the TRAPUM observations, the computing power of FBFUSE was used to synthesize either 277 or 288 tied-array beams (Chen et al. 2021). These allowed us to cover an area of several square arcminutes around the nominal center of NGC 1851, which extended a few times beyond the half-mass radius of the cluster (see Fig. 1).

Fig. 1. TRAPUM tilings of two observations of NGC 1851, performed with MeerKAT on 2021 May 8 at L-band (left) and on 2021 May 21 at UHF-band (right), with similar telescope elevations. Both panels show the same portion of the sky, for an easier comparison of the area covered by the tilings in the two different bands. The ellipses show the tied-array beams up to 84% of their boresight power at their respective central frequencies, that is, 816 MHz for the UHF band, 1284 MHz for the L band. The blue ellipses show the beams that have been searched for this work. The dashed circle shows the half-mass radius of NGC 1851, located at 0.85 arcmin from the cluster center. The position of the previously known pulsar, NGC 1851A, is shown by the white cross.

The other six observations were carried out as part of the MeerTime GC pulsar timing program. Their set-up was chosen based on the results obtained from the TRAPUM data previously taken: we used the UHF receivers and PTUSE as the backend. The latter produced a single tied-array beam, synthesized at the center of NGC 1851. This was enough to observe all the new pulsars at once with maximum sensitivity (see Sect. 3). The observing band was divided into 1024 channels, coherently de-dispersed at a dispersion measure (DM) of 52.1 pc cm⁻³ (i.e., the DM of the only pulsar previously known in the cluster), and recorded in full-Stokes mode, every 9.41 μs. Since the PTUSE data were also stored as search-mode files, they were used to confirm and follow-up the pulsars found in the TRAPUM data, and also to search for new ones.

For some of our analyses, we also took advantage of the large set of archival observations of NGC 1851 taken with the GBT and the GMRT telescopes. The GBT data were taken with central frequencies of 350, 820, and 1950 MHz in the years 2004–2006 and are described by Freire et al. (2007). The GMRT data were acquired with a central frequency of 400 MHz from 2017 to 2020, with observational set-ups that are described by Ridolfi et al. (2019).

2.1. Sensitivity

Given the observing parameters described above, we can use the radiometer equation (Dewey et al. 1985) to estimate our typical sensitivity: we find that we are sensitive to pulsars with minimum mean flux densities of S_L = 10 μJy and S_UHF = 14 μJy in the L band and UHF band, respectively. These are estimated considering, for both observing bands: 60 antennas, for a total telescope gain of 2.63 K Jy⁻¹; 2 h of integration time; two polarizations summed; a correction factor of 1.1 to account for possible digitization and other losses of various nature; a pulse duty cycle of 8 percent; a DM of 52 pc cm⁻³; a detection signal-to-noise ratio (S/N) of 10. For the L-band observations, we assume an effective bandwidth (i.e., usable bandwidth after the removal of radio frequency interference, henceforth RFI) of 650 MHz, and a total system temperature of 26 K⁴. The latter is the sum of the receiver temperature (T_rec = 18 K), and contributions from the sky (T_sky = 3.5 K), atmosphere (T_atm = 2 K), and ground spillover (T_spill = 2.5 K). For the UHF observations, we consider an effective bandwidth of 500 MHz, and a total system temperature of 31.5 K (with the contributions being T_rec = 20 K, T_sky = 5.5 K, T_atm = 1.5 K, T_spill = 4.5 K). For extremely fast-spinning pulsars, with spin periods ≲2 ms, the minimum detectable mean flux densities are higher because of pulse smearing effects, in turn due to the sampling time and finite channel width becoming important.

2.2. Data reduction

Before being analyzed, all the observations were preprocessed. In particular, all the TRAPUM beams, which are recorded as SIGPROC’s⁵ “filterbank” search-mode files, were first incoherently de-dispersed at the nominal cluster DM of 52.1 pc cm⁻³. Then, their frequency channels were summed in groups of 16, bringing their total number from 4096 to 256, substantially reducing the total data volume. While doing so, the Inter-Quartile Range Mitigation algorithm⁶ (Morello et al. 2022) was applied to significantly reduce the presence of RFI in the data. This algorithm is very effective with the MeerKAT data in both the L and UHF bands. The resulting files were, in the vast majority of cases, almost completely free of RFI, greatly easing their subsequent analysis. Once the processed files were produced and checked, the original 4096-channel data were deleted. The MeerTime PTUSE observations, which are recorded as PSRFITS (Hotan et al. 2004) search-mode files, were decimated by a factor of two in the frequency domain (preemptively applying appropriate shifts to the channels, to account for the cluster DM of 52.1 pc cm⁻³), a factor of eight in the time domain, and the two polarizations were summed. This was done with psrfits_subband, one of the tools of the PSRFITS_UTILS⁷ package. The resulting files, which therefore have 512 channels, a time resolution of 75.29 μs, and only total-intensity information, were used to search for new pulsars. The original full-Stokes, native-resolution PTUSE data were instead used to carry out polarimetric measurements, and they will also be preserved for future high time resolution timing and more in-depth polarization studies.

2.3. Search

All the 13 NGC 1851 observations were searched for new pulsars using PULSAR_MINER⁸, an automated pipeline based on the PRESTO⁹ pulsar searching package (Ransom et al. 2002; Ransom 2011). Each search-mode file was cleaned using PRESTO’s rfifind and then searched in the Fourier domain with PRESTO’s accelsearch. The latter can recover a pulsar’s periodic signal in a situation where its Doppler shift is varying rapidly because of its orbital motion¹⁰, by considering signals whose frequency changes by up to a number z_max of Fourier bins around a given frequency. Therefore, we ran accelsearch with a z_max = 0 to target isolated pulsars, and with a z_max = 200 to target binary pulsars, in both cases summing up to 8 harmonics. This was done on the full-length data, as well as on segments of 120 (for the 4-h observation), 60 and 30 min, with DM trials spanning the range 49.5−54.5 pc cm⁻³. We refer to Sect. 3 of Ridolfi et al. (2021) for further details on all the steps of the searching pipeline. Due to the huge data volume and consequent computational time required, in the case of the TRAPUM observations, we restricted our initial search to only those beams that covered the half-mass radius of NGC 1851 (Fig. 1). If the same candidate was found in more than one epoch and in beams with compatible sky positions, it was marked as a “confirmed” pulsar and further inspected.

Finally, in addition to the periodicity search, we carried out a single-pulse search. This was done by running PRESTO’s single_pulse_search.py with default parameters on all the full-length de-dispersed time series produced by PULSAR_MINER and by visually inspecting the resulting diagnostic plots.

2.4. Binary orbits and timing

A high fraction of the confirmed pulsars exhibited some nonzero acceleration, a likely sign of being part of a binary system. For these pulsars, we first tried to determine their rough orbital parameters using the following methods. The first attempt was made using the Period–Acceleration diagram technique described in Freire et al. (2001a). After that, we fit the observed spin period as a function of time, P_obs(t), to an orbital model, with initial parameters derived from the previous method. The P_obs(t) fitting greatly benefitted from the six closely spaced observations taken in May 2021 and was implemented using PRESTO’s fit_circular_orbit.py for the circular systems, and either PRESTO’s fitorb.py or the FORTRAN code fitorbit¹¹ for the eccentric systems. The derived orbital parameters were later further refined through timing (see below).

To better characterize all the newly found pulsars, we extracted times-of-arrival (ToAs) from the available detections and did a timing analysis. We did so by first folding the search-mode data with the DSPSR¹² pulsar software package (van Straten & Bailes 2011). For the isolated pulsars, we constructed an ephemeris that contained the best barycentric spin period, P, and DM as determined by the search pipeline. For the binary pulsars, the ephemeris also included the best-fitting orbital model as previously derived. Using this to fold all the observations often resulted in additional faint detections of the pulsar that were not found by the search pipeline.

The resulting folded archives were processed with PSRCHIVE¹³ (Hotan et al. 2004; van Straten et al. 2012) to extract topocentric ToAs from each detection. This was done by cross-correlating the observed pulse profiles with an analytical pulse template of the pulsar. These ToAs were then fit to an appropriate timing model with the TEMPO¹⁴ pulsar timing software. Because of the faintness and the sparsity of the detections of several of the new pulsars, we used DRACULA¹⁵ (Freire & Ridolfi 2018) to determine the exact number of rotations of the pulsar between groups of ToAs. The use of this algorithm showed that, for most pulsars, many possible solutions can still be obtained using the existing ToA sets. This means that additional observations are required for the determination of the correct timing solutions. For this reason, the timing analysis of these pulsars will be presented elsewhere.

2.5. DM and flux densities

For each pulsar, we ran PSRCHIVE’s pdmp on all the available UHF detections to obtain, among other things, a best-fitting DM value for each of them. The average of the best-fitting DMs was calculated to give the best estimate of the DM associated with the pulsar. The standard error on the average of the same measurements was adopted as the associated uncertainty.

As described in Sect. 3.2, the pulsars in NGC 1851 show significant variations in their apparent brightness from epoch to epoch and even within the same observation. Therefore, a single measurement of the mean flux density S for a pulsar, is not representative of its true intrinsic value. Rather, to obtain a more realistic estimate of S, we measured the latter from each observation and considered the average value ⟨S⟩, and its standard deviation as associated uncertainty. This was done using only the UHF-band observations, since they are the most numerous ones and have the highest detection rate for all the discoveries. For each observation we took the S/N value estimated by pdmp, and used the radiometer equation to calculate the corresponding mean flux density in the UHF band, S_UHF, for all the detections. The nondetections were also taken into account by arbitrarily assigning S_UHF = 7 μJy to them, as this is exactly half of the estimated minimum detectable S for the UHF band.

2.6. Polarimetric properties

We used the six most recent UHF observations listed in Table 1, which were taken with PTUSE in full-Stokes mode, to measure the polarimetric properties of the pulsars in NGC 1851. After folding the native-resolution data for each pulsar with DSPSR, we used PSRCHIVE’s pac routine on the folded archives to correct for the time-varying parallactic angle. This was sufficient to obtain fully calibrated folded archives for the latest four observations (obs. ids 10U, 11U, 12U, 13U), as these were calibrated, correcting for the instrumental response, directly while being recorder. The other two observations (obs. ids 08U and 09U) were not calibrated during the recording, hence we applied a calibration solution accounting for the system’s response a posteriori on the folded archives (Serylak et al. 2021). For each pulsar, we then selected the best detections of the sample and summed them together with PSRCHIVE’s psradd, in order to maximize the S/N of the resulting archive. We then ran PSRCHIVE’s rmfit on the latter to measure the associated rotation measure (RM): the code brute-force searched for the best RM within a specified range (which we chose to be between −500 and 500 rad m⁻²) that maximizes the measured fraction of linearly polarized flux density. If the pulsar was sufficiently polarized and the S/N high enough, rmfit returned the best RM and its uncertainty, evaluated by fitting a Gaussian function to the polarized flux as a function of trial RM. The so obtained RM value was used to de-Faraday the high-S/N summed archive of the pulsar, and to plot the resulting linearly and circularly polarized integrated profiles.

2.7. Localization

In the absence of timing solutions (which generally result in very precise positions for the pulsars), we exploited the hundreds of beams synthesized in each of the seven TRAPUM observations of NGC 1851 to obtain improved localizations for all the newly found pulsars. To do so, we used a custom version of the SeeKAT¹⁶ tied-array beam localization software (Bezuidenhout et al., in prep.).

For a given observation, SeeKAT takes the S/N values with which a pulsar is detected in the different beams. Using the ratios of the S/N in pairs, it then performs a maximum-likelihood analysis to find the source location that minimizes the differences in the ratios, using the beam positions and accounting for their point spread functions. The latter are calculated with the Mosaic¹⁷ software (Chen et al. 2021), using the observation mid-point date and time, as well as the central observing frequency. The custom version of SeeKAT used for this work is also capable of combining information from multiple observations (up to seven, in our case), returning probability contours for the pulsar position that are significantly more accurate than those obtained using a single observation.

3. Results

Our searches resulted in the discovery of 13 new pulsars in NGC 1851, which we designate from now on with the cluster’s name and a letter, which reflects the order of discovery. No significant single pulses were found in any observation or beam, within the DM range explored (49.5−54.5 pc cm⁻³). All the discoveries are MSPs with spin periods in a very narrow range, 2.81–6.63 ms, the only exception being NGC 1851I, which is a 32.65-ms mildly recycled pulsar. In Table 2 we list the basic parameters derived for all the discoveries, while in Fig. 2 we show their high-S/N, UHF-band integrated pulse profiles. In Table 3, we list the detections of each of these pulsars in the MeerKAT observations. In the remainder of this section, we describe the characteristics of the new pulsars, and then outline some of their collective properties.

Fig. 2. High-S/N integrated pulse profiles of the 14 pulsars in NGC 1851, each obtained summing the brightest detections available in the MeerKAT UHF band. Bottom panels: total intensity (black line) and, for the pulsars where some significant polarization was measurable, linear (red line) and circular (blue line) polarization profiles. Top panels: position angle (PA) of the linear polarization, showing measurements that are at least 3-σ significant. We follow the PSR/IEEE convention (van Straten et al. 2010), where the left-hand circular polarization is positive and the PA angle increases counterclockwise on the sky. The horizontal section of the short bars on the bottom-right of each panel shows the effective time resolution, which accounts for the sampling time (9.41 μs) and the intra-channel smearing; in all cases it is < 14 μs.

3.1. Characteristics of the new pulsars

About half of the new discoveries (NGC 1851B, C, J, K, M and N) are, based on our current understanding of their properties, isolated pulsars. This fraction is similar to that of Terzan 5, M28, NGC 6440 and NGC 6441. Pulsars C, J, K and M were also detected on a few occasions in archival GMRT data taken at 400 MHz, while both B and C were also seen a few times in the old GBT 350-MHz observations. However, all of these detections had low S/N.

One of these isolated pulsars, NGC 1851M, is extremely faint even in MeerKAT data: it was blindly detected by the search pipeline in only four observations. The direct folding of all the MeerKAT observations yielded an additional four detections, all of which were extremely faint. Because of this, we could not derive an improved position for this pulsar using SeeKAT. Nonetheless, we can assert that it must be located within ∼20−30 arcsec from the cluster center, as it is always detected with highest S/N in the central UHF beam of the TRAPUM tiling. We also note that NGC 1851M is the only pulsar, among the discoveries, that was not detected in the L band (see Table 3), suggesting that it has a particularly steep spectrum.

An additional five pulsars (NGC 1851F, G, H, I and L) are in binary systems with low-mass companions (as estimated from the mass function) and low-eccentricity orbits. Pulsars F, G and H have orbital periods in the range ∼2−16.94 d, which are typical of pulsar–white dwarf (WD) binaries. Pulsars I and L, on the other hand, have orbital periods of ∼1 d, like those of a few “redbacks”. These are a subclass of the “spider” binary pulsars (see Roberts 2013), which are compact systems where the companion shows signs of being irradiated by the pulsar, and often exhibit other peculiar features. One of these is the presence of eclipses, where the radio pulses are either delayed or are undetectable for a good fraction of the orbit. Their absence in NGC 1851I and L, however, is a strong indication that the companion stars are not extended like the companions of redbacks, but instead are also low-mass WDs. Another characteristic of redbacks is the presence of long-term, erratic variations of the orbital period. We do not yet have enough data to decide on the stability of the orbit of any of these systems. Ultimately, optical observations of the companion stars are required to determine the actual nature of these systems.

At least two of these binaries, NGC 1851H and I, seem to have significant nonzero eccentricities. Despite being small, they are several orders of magnitude larger than similar pulsar–WD systems in the Galaxy. This is commonly observed in GCs and is due to perturbations caused by close stellar encounters, which are likely to happen in the crowded environment. Of this set of pulsars, G is the only pulsar detected in archival data, specifically in a few GMRT 400-MHz observations.

The most interesting new pulsars are the two eccentric binary MSPs, NGC 1851D and E. As for the other pulsars, no phase-connected timing solutions are yet available. However, their orbital parameters have been determined to high precision (see Fig. 3). These systems resemble NGC 1851A, the previously known binary pulsar in the cluster: they have large orbital eccentricities and massive companions which suggests they are very likely to be the result of secondary exchange encounters, in which the original low-mass He WD that recycled these pulsars was replaced by the current, more massive companion. NGC 1851D is in a 5.69-d orbit with an eccentricity of e = 0.86. From its mass function, for a pulsar mass M_p = 1.4 M_⊙, we derive a companion mass M_c > 0.48 M_⊙. Hence, we can preliminarily conclude that its companion is likely to be a massive carbon–oxygen WD. The pulsar has also been detected often in the existing GMRT and GBT data. Using these data we were able to extend its orbital solution to May 2005, when the GBT started observing NGC 1851A regularly (Freire et al. 2007). NGC 1851E is in a 7.44-d, less eccentric (e = 0.71) orbit, but the companion is more massive: from the mass function, for M_p = 1.4 M_⊙, we obtain M_c > 1.53 M_⊙, which indicates that the companion is most likely a NS. Unlike D, it is not detectable in the bulk of the archival observations (which were made at lower frequencies, in order to time NGC 1851A, a steep spectrum pulsar), but it is detected in the few GBT observations made with the 820 MHz and S-band receivers (as described in Freire et al. 2007).

Fig. 3. Spin period as a function of time as observed in May 2021 for the two eccentric binary pulsars NGC 1851D (left) and NGC 1851E (right). The former is in a 5.69-d orbit with an eccentricity of 0.86. The latter is in a 7.44-d orbit with an eccentricity of 0.71.

All the newly found NGC 1851 pulsars were detected with the highest S/N in the central beam of the tiling of each TRAPUM observation. This is confirmed by the localizations obtained with SeeKAT as described in Sect. 2.7: all the maximum-likelihood pulsar positions are at ≲2 arcsec from the nominal cluster center (Fig. 4). However, given the uncertainties, some objects could eventually display a larger offset. Much more accurate radio timing positions will soon be available once phase-connected solutions are achieved with future follow-up observations.

Fig. 4. Positions of the new pulsars determined by SeeKAT, relative to the center of the cluster (the intersection of the horizontal and vertical dashed lines) and the position of NGC 1851A, indicated by the red star. The gray scale background indicates probability density, while the purple, blue and green ellipses represent the 68%, 95%, and 99.7% probability contours, respectively. The black dashed circle indicates the core radius of NGC 1851, located at 0.09 arcmin from its nominal center (Harris 2010). The detections of NGC 1851M were too faint to obtain a reasonable position estimate with SeeKAT. All pulsars are detected with highest S/N in the central beams of all TRAPUM observations.

All the discoveries have DMs in the range ∼50.96−52.58 pc cm⁻³. Including the previously known binary pulsar PSR J0514−4002A, the DMs of the pulsars in NGC 1851 have an average of 51.82 pc cm⁻³ and a standard deviation of 0.47 pc cm⁻³. The overall spread of 1.62 pc cm⁻³ is just 3% of the average value, in line with the DM distribution of all GCs known to host two or more pulsars.

For ten out of the 14 pulsars in NGC 1851, we were able to measure their basic polarimetric properties as described in Sect. 2.6. All the pulsars have best RM values between −6 and −14 rad m⁻² (see Table 2), although these are affected by large uncertainties due to the faintness of their polarized flux. While a few pulsars show high fractions of linear polarization (> 10 percent), pulsar K stands out for having an almost 100% linearly polarized leading profile component (Fig. 2). We defer a more thorough analysis of the polarization properties of the NGC 1851 pulsars and their possible applications to a future publication.

3.2. Flux density variations

In many clusters located well outside the Galactic plane (such as 47 Tuc and NGC 6752) the low average DM is often accompanied by large (i.e., more than 50% in magnitude) and rapid (i.e., over time scales of < 1 h) variations in the apparent brightness of the pulsars, caused by diffractive scintillation (Camilo et al. 2000; D’Amico et al. 2002). This is also what we see in NGC 1851: all the new discoveries show very large flux density variations (sometimes by a factor of 3 or more) that can occur over time scales as short as a few tens of minutes (Fig. 5). When not averaged over long integration times, such variations manifest themselves as very small-bandwidth “scintles” in the frequency band (Fig. 6), something that is particularly evident in the UHF-band observations. This is exactly what is expected in the case of diffractive scintillation, where the scintle size increases in both time and bandwidth in proportion to the fourth power of the frequency (Rickett 1977).

Fig. 5. Short-term flux density variations as seen in the time domain in the isolated pulsar NGC 1851C in the UHF band (left panel) and in the binary pulsar NGC 1851E in the L band (right panel).

Fig. 6. Example of flux density variations as seen in the frequency domain in NGC 1851M in the UHF band over ∼30 min of integration time (left panel) and in NGC 1851B in the L band over ∼1 h of integration time (right panel).

What we see in the detections of the new pulsars can also be used to address a few hypotheses proposed by Freire et al. (2007) to explain the large flux density variations that they saw in NGC 1851A. One of their argument was that those variations could be “eclipses” caused by gas outflowing from the companion star. This hypothesis can now be excluded, because these changes are seen in all of our new discoveries, even those that are isolated pulsars. The other possibility that they proposed was refractive scintillation caused by plasma structures in NGC 1851. However, refractive scintillation generally has long timescales on the order of weeks to years, not hours, and in addition, typically only modulates the measured flux densities by up to a factor of ∼2. It is for these reasons that Freire et al. (2007) speculated on the presence of plasma “lenses” within NGC 1851 itself as a way of producing these changes on a much faster scale. However, in the case of refractive scintillation, the flux variations should be achromatic, that is, independent of the radio frequency, which is not what we see in our MeerKAT data. A detailed discussion of the scintillation phenomenon seen in the NGC 1851 pulsars is beyond the scope of this paper and will be addressed in a future study.

4. Discussion

The 13 new discoveries presented here, represent, thus far, the most successful search of a single GC with MeerKAT. They provide us with a far clearer view of the pulsar population of NGC 1851, and give us the opportunity to relate the pulsar characteristics with the cluster properties.

This population is unique in having three highly eccentric binary MSPs with massive companions, which are, in all likelihood, the result of secondary exchange encounters. Although several other clusters have a significant number of MSPs in eccentric binaries, like Terzan 5 and M28, those binaries have generally lower companion masses; this means that they are not as clearly identifiable as secondary exchange products.

As discussed by Verbunt & Hut (1987), the number of LMXBs and MSPs in a globular cluster depends primarily on the overall exchange encounter rate, Γ, which depends on the cluster’s structural parameters. However, apart from their sizes, the pulsar populations in different clusters present strikingly different characteristics, like the fraction of binary pulsars, the presence of apparently young pulsars, and the number of systems formed in secondary exchange encounters. As described by Verbunt & Freire (2014), the so-called “encounter rate for a single binary”, γ, explains these differences, but not others, like the distribution of MSP spin periods. This parameter, which they normalize to the value of the cluster M 4 (γ_{M 4}), is a measure of the probability for an already-formed binary system to undergo a close encounter with a third body, which could perturb the orbit, replace one of the stars, or even disrupt the whole system. In GCs with a large γ – especially those with collapsed cores – we are much more likely to find slow, mildly recycled pulsars and products of secondary exchange interactions. This is because a NS that accretes matter in a LMXB phase has a smaller probability of being fully recycled and becoming an MSP before the system is disrupted due to a dynamical interaction, hence a significant fraction of pulsars in such clusters are partially recycled. Even if the system completes recycling, the resulting binary millisecond pulsar will also have a much higher probability of disruption. Because of such disruptions, the fraction of isolated pulsars should be larger for clusters with larger γ.

This general picture is confirmed by the observations: in several core-collapsed, high-γ clusters (NGC 6624, NGC 6752, NGC 6517, Terzan 1 and M 15) the vast majority of the pulsars are isolated, and a good fraction are slow or mildly recycled; also, of the few known binaries in these clusters, a good fraction (PSR B2127+11C in M 15, Jacoby et al. 2006; PSR J1807−2500B in NGC 6544, Lynch et al. 2012; and PSR J1823−3021G in NGC 6624, Ridolfi et al. 2021) are eccentric binary pulsars with massive companions that are clearly the result of secondary exchange encounters.

On the contrary, GCs with a smaller γ provide an environment where a binary can survive, on average, for much longer. A LMXB will evolve undisturbed, and will almost inevitably evolve like the LMXBs in the Galaxy. The MSP population in such clusters will thus resemble that of the Galactic disk, with most of the pulsars being fully recycled MSPs, either isolated or in circular binaries with a He–WD, or with a more main sequence-like star, as in the case of the aforementioned spider systems. Striking examples of clusters of this type are 47 Tucanae (Freire et al. 2017) and M 62 (Lynch et al. 2012); these examples are even more striking because of their large total interaction rate, Γ, which is the reason why the absolute number of MSPs is so large.

Between these two groups are the clusters with intermediate values of γ. This is reflected in the observed pulsar populations, which have mixed characteristics. A major example is Terzan 5: of 39 pulsars that the latter is known to host, 19 are isolated and 20 are in binaries; almost all of them are MSPs, but a few are mildly recycled; there are many circular binaries in pulsar–WD, spiders (i.e., black widows and redbacks), but also several wide, eccentric, and sometimes massive, binaries. Almost identical considerations apply to M28, as well as NGC 6440 and NGC 6441.

With an encounter rate per binary of γ = 12.44 γ_{M 4}, NGC 1851 is also an intermediate-γ GC. This is reflected in the characteristics of its pulsar population: (a) the number of isolated pulsars (6) is comparable to the number of binary pulsars (8); (b) all the pulsars are MSPs, with the exception of NGC 1851I, which is a mildly recycled pulsar; (c) among the binaries, three (pulsars A, D and E) are extremely eccentric and clearly result from secondary exchange encounters, while all the others look to be “standard” pulsar–WD systems.

Notwithstanding, NGC 1851 also shows some marked differences with some of the other intermediate-γ clusters. Firstly, similarly to NGC 6441 but at odds with Terzan 5, NGC 6440 and M28, none of the eight binary pulsars are definite spiders (although I and L could be wide-orbit, noneclipsing redbacks). However, this could be due to the small number statistics and to the fact that spider pulsars are notoriously elusive: the presence of long eclipses and their extremely compact orbits make them challenging targets for search algorithms. Nevertheless, this is unlikely to represent the full explanation: M28 has the same number of known pulsars (14), but fully half of its pulsar population consists of spiders: five black widows (M28G, J, L, M, N) and two redbacks, M28H and I, the latter being a transitional MSP system (Papitto et al. 2013). More observations, especially with the S-band receivers of MeerKAT, and the use of more sophisticated search techniques may eventually reveal some spiders in NGC 1851.

Another interesting aspect is that all the 14 pulsars in NGC 1851 are within or just outside the core. This situation is similar to that of M28, where, with the exception of M28F, all pulsars are also very close to the center (Bégin 2006). Although some concentration is to be expected from mass segregation, which makes heavier objects (such as NSs) sink toward the center of the cluster, it is certainly much more centrally condensed than the pulsar populations of Terzan 5 and 47 Tucanae, which are generally understood to be in dynamical equilibrium with the remaining stellar populations in the cluster (see e.g., Grindlay et al. 2002 and references therein). The reasons for this peculiarity of NGC 1851 and M28 are not clear: either their core radii have not been accurately measured, or the pulsars are confined by a deeper, narrower gravitational potential caused by a central concentration of massive objects that do not affect the overall distribution of luminous stars. A detailed investigation of this issue will be carried out once timing solutions are available for all pulsars in NGC 1851.

5. Conclusions and prospects

In this paper, we have presented 13 new pulsars found by TRAPUM in the GC NGC 1851 using the MeerKAT radio telescope. Six of these are solitary pulsars, while the other seven are part of binary systems. Among these discoveries, there are a few noteworthy objects. NGC 1851D is a MSP in a wide orbit with a heavy WD companion. The system has an eccentricity of e = 0.86, which makes it the 5th most eccentric binary known. NGC 1851E is also a MSP in a very eccentric (e = 0.71) binary with a NS as companion. In all likelihood, both binaries formed in secondary exchange encounters where the original star that recycled the pulsar was replaced by the current, heavier companion. NGC 1851 is a cluster with intermediate interaction rate per single binary, γ. This is reflected in its pulsar population, which shows properties that are in between those found in core-collapsed and dynamically relaxed clusters.

The large ensemble of pulsars in NGC 1851 opens up interesting scientific possibilities for future studies. One of the obvious extensions of the work presented here using the current TRAPUM data, will be the search for additional pulsars in all the beams beyond the half-mass radius, which have so far been excluded due to computational and time constraints. This could reveal, for instance, systems that may have been ejected from the core during secondary exchange interactions, and that have not had time yet to sink back in the core (these types of objects are more common in GCs with a large γ, see Verbunt & Freire 2014). To complement this, all the beams, starting from the inner ones, will also be searched using more sophisticated (but often more computationally demanding) methods, such as jerk searches (Andersen & Ransom 2018), phase-modulation searches (Ransom et al. 2003), template-banking algorithms over 3 or 5 Keplerian parameters (Balakrishnan et al. 2022), as well as fast-folding algorithms (Morello et al. 2020); the first three have a higher potential of unveiling fast pulsars in extremely compact binaries, while the latter is more sensitive to very slow, isolated pulsars. With such searches, we will know whether the intriguing lack of spiders is real, or a result of current limitations in our search techniques. Given the intermediate structural characteristics of NGC 1851, the likelihood of finding both types of pulsars is not negligible.

The polarization properties of the pulsars presented in this work will be further studied as more follow-up observations are taken in full-Stokes mode. If more precise RM measurements are possible, these could potentially be exploited to probe the Galactic magnetic field at arcsecond scales in a different direction than that explored by Abbate et al. (2020b), who used the pulsars in 47 Tucanae.

With a relatively small number of additional observations, we can expect to derive phase-connected timing solutions for all the newly found pulsars. These will immediately give access to a much more precise determination of their projected positions, which will allow us to look for possible X-ray counterparts in Chandra data, as well for the optical counterparts of the binary pulsars in archival HST and future JWST data.

Continued timing in the coming years will enable the measurement of the proper motions and of the higher-order spin period derivatives for all pulsars. This will allow us to constrain the cluster gravitational potential and probe the presence of nonluminous matter (such as a possible intermediate-mass back hole) at its center (e.g., Freire et al. 2017; Abbate et al. 2018, 2019). This and a more in-depth study of the scintillation patterns observed in the pulsars could shed light on the possible presence of an intracluster medium, as already seen in 47 Tuc (Freire et al. 2001b; Abbate et al. 2018).

For the binary pulsars (especially for the very eccentric systems, NGC 1851D and E) the long-term timing will yield precise measurements of the rate of advance of periastron (). Additional post-Keplerian parameters are desirable, particularly for determining the masses of the individual components. A detection of the Shapiro delay, even at a very low level, can, when combined with a measurement of , yield precise masses (e.g., NGC 6544B in Lynch et al. 2012); however, this is especially difficult in this case given how faint the pulsars are, and how low is their timing precision. For NGC 1851E, it will be possible to measure the Einstein delay (γ_E), but, as in the case of NGC 1851A, this will take at least a decade (Ridolfi et al. 2019). This is not likely to happen for NGC 1851D: because of its longitude of periastron (∼358deg, close to 0deg), other geometric effects caused by the proper motion of the system will make γ_E undetectable (see detailed discussion in Sect. 3.7 of Ridolfi et al. 2019). For all the binaries, the measurements of the orbital period derivative would give access to reliable estimation of the intrinsic spin-down rate, surface magnetic field and characteristic age of the pulsars.

Many of these goals will greatly benefit from archival data taken at other radio facilities (particularly GBT and GMRT) if the pulsars are also frequently detected in them once their phase-connected solutions are determined.

The discoveries presented in this paper are among the first made by the TRAPUM GC pulsar survey. The latter exploits the outstanding sensitivity of all the 64 antennas available at MeerKAT, as well as advanced beamforming techniques to search for pulsars in a large portion of 28 selected GCs. As more and more clusters are observed, we expect to significantly increase the known population of pulsars in several other GCs, allowing new studies and a much deeper understanding of the conditions and the physical processes that drive the evolution of these stellar systems.

Acknowledgments

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time space reference systems department of the Paris Observatory. MeerTime data is housed on the OzSTAR supercomputer at Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government. 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. We thank the GMRT staff for help with the observations that were used in this work. The GMRT is operated by the National Centre for Radio Astrophysics (NCRA) of the Tata Institute of Fundamental Research (TIFR), India. This work was partly performed on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government. 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. APo, AR and MBu 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. APo, AR, MBu also acknowledge support from the Ministero degli Affari Esteri e della Cooperazione Internazionale – Direzione Generale per la Promozione del Sistema Paese – Progetto di Grande Rilevanza ZA18GR02. SMR is a CIFAR Fellow and is supported by NSFP hysics Frontiers Center awards 1430284 and 2020265. SMR is a CIFAR Fellow and is supported by the NSF Physics Frontiers Center award 1430284. FA, MK, and PVP gratefully acknowledges support from ERC Synergy Grant “BlackHoleCam” Grant Agreement Number 610058. AR, TG, PCCF, VVK, MK, FA, EDB, PVP, DJC, WC and APa acknowledge continuing valuable support from the Max-Planck Society. BWS and MCB acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694745). LV acknowledges financial support from the Dean’s Doctoral Scholar Award from the University of Manchester. RPB acknowledges support from the ERC under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 715051; Spiders). MBa acknowledges support through ARC grant CE170100004. We are grateful to the anonymous referee for useful comments, which helped improve the manuscript.

References

All Tables

All Figures

Fig. 1. TRAPUM tilings of two observations of NGC 1851, performed with MeerKAT on 2021 May 8 at L-band (left) and on 2021 May 21 at UHF-band (right), with similar telescope elevations. Both panels show the same portion of the sky, for an easier comparison of the area covered by the tilings in the two different bands. The ellipses show the tied-array beams up to 84% of their boresight power at their respective central frequencies, that is, 816 MHz for the UHF band, 1284 MHz for the L band. The blue ellipses show the beams that have been searched for this work. The dashed circle shows the half-mass radius of NGC 1851, located at 0.85 arcmin from the cluster center. The position of the previously known pulsar, NGC 1851A, is shown by the white cross.

In the text

Fig. 2. High-S/N integrated pulse profiles of the 14 pulsars in NGC 1851, each obtained summing the brightest detections available in the MeerKAT UHF band. Bottom panels: total intensity (black line) and, for the pulsars where some significant polarization was measurable, linear (red line) and circular (blue line) polarization profiles. Top panels: position angle (PA) of the linear polarization, showing measurements that are at least 3-σ significant. We follow the PSR/IEEE convention (van Straten et al. 2010), where the left-hand circular polarization is positive and the PA angle increases counterclockwise on the sky. The horizontal section of the short bars on the bottom-right of each panel shows the effective time resolution, which accounts for the sampling time (9.41 μs) and the intra-channel smearing; in all cases it is < 14 μs.

In the text

	Fig. 3. Spin period as a function of time as observed in May 2021 for the two eccentric binary pulsars NGC 1851D (left) and NGC 1851E (right). The former is in a 5.69-d orbit with an eccentricity of 0.86. The latter is in a 7.44-d orbit with an eccentricity of 0.71.
In the text

Fig. 4. Positions of the new pulsars determined by SeeKAT, relative to the center of the cluster (the intersection of the horizontal and vertical dashed lines) and the position of NGC 1851A, indicated by the red star. The gray scale background indicates probability density, while the purple, blue and green ellipses represent the 68%, 95%, and 99.7% probability contours, respectively. The black dashed circle indicates the core radius of NGC 1851, located at 0.09 arcmin from its nominal center (Harris 2010). The detections of NGC 1851M were too faint to obtain a reasonable position estimate with SeeKAT. All pulsars are detected with highest S/N in the central beams of all TRAPUM observations.

In the text

	Fig. 5. Short-term flux density variations as seen in the time domain in the isolated pulsar NGC 1851C in the UHF band (left panel) and in the binary pulsar NGC 1851E in the L band (right panel).
In the text

	Fig. 6. Example of flux density variations as seen in the frequency domain in NGC 1851M in the UHF band over ∼30 min of integration time (left panel) and in NGC 1851B in the L band over ∼1 h of integration time (right panel).
In the text