ture supernovae (e.g. P. Podsiadlowski et al. 2004 ) or accretioninduced collapse of a white dwarf in a binary (e.g. A. J. Ruiter et al. 2019 ). Observ ations of X-r a y binaries containing Be stars show evidence for a subpopulation of neutron stars that do indeed receive low ( 10 km s -1 ) kicks due to a yet-unknown mechanism (E. Pfahl et al. 2002b ;R. Valli et al. 2025 ). The physical processes that pr oduce low-kick neutr on stars in globular clusters likely depend on metallicity and other clust er paramet ers in a testable manner (e.g. N. Iv anov a et al. 2008 ), though observ ational evidence for specific scenarios is still scant.

The open questions do not stop at the origin of the old neutron stars recycled to form millisecond pulsars. F. Kirsten et al. ( 2022 ) reported the discovery of a fast radio burst in a metal-poor globular cluster in the nearby galaxy M81, which strongly suggests the presence of a very young, highly magnetized neutron star in this old globular cluster. The close distance to M81 and the normalcy of this cluster (e.g. it has no unusual structural properties; K. C. Dage et al. 2023 ) implies that young neutron stars are forming in globular clusters at a meaningful rate, possibly through mergers of dynamically formed white dw arf-white dw arf binaries (K. Kr emer, A. L. Pir o & D. Li 2021 ; K. Kremer et al. 2023 ). Another r ecently discover ed fast radio burst is localized with a large projected offset from a massive elliptical, most consistent with an old globular cluster host (V. Shah et al. 2025 ). These recent discover-MNRAS 547, 1-20 (2026) ies accompany existing evidence for a few young neutron stars in Galactic globular clusters (e.g. J. Boyles et al. 2011 ;K. Kremer et al. 2024 ). Overall, the ongoing uncertainty about the origin and fate of neutron stars in clusters motivates their continued study.

Terzan 5 has been a t ouchst one object in the study of neutron stars in globular clusters, hosting an early cluster X-ray burster (K. Makishima et al. 1981 ) and found to contain only the second discovered eclipsing millisecond pulsar (A. G. Lyne et al. 1990 ). Radio continuum imaging with the Very Large Arr ay (VLA) w as used to infer the presence of many y et-undiscov ered millisecond pulsars (A. S. Fruchter & W. M. Goss 1990, 2000 ), with the latter paper suggesting that Terzan 5 contained more pulsars than any other clust er, ev en though at that time only a single pulsar was known. This prediction was borne out by the discovery of 21 millisecond pulsars in Terzan 5 using the Green Bank Telescope by S. M. Ransom et al. ( 2005 ). Ther e ar e now 49 total pulsars confirmed in Terzan 5 (P. V. Padmanabh et al. 2024 ), indeed the largest known population of any cluster, including the fastestspinning millisecond pulsar known in either cluster or field (J. W. T. Hessels et al. 2006 ). As all but one pulsar in Terzan 5 is a millisecond pulsar, we use the terms interchangeably in this paper.

The abundant pulsar population in Terzan 5 is a consequence of its high mass and dense core, which leads to an e xtraor dinarily high stellar inter action r ate-the highest of any Galactic globular cluster (A. Bahramian et al. 2013 ). However, the total population of pulsars in Terzan 5 is still very uncertain. Pulsar search observations can have difficulty detecting the eclipsing 'spider' pulsars that are common in Terzan 5 as well as faint pulsars in close binaries, and are sensitivity limited even for isolated pulsars (M. Cadelano et al. 2018 ; P. V. Padmanabh et al. 2024 ).

In a previous paper (R. Urquhart et al. 2020 , hereafter U20 ) we used new S (2-4 GHz) and C (4-8 GHz) band imaging of Terzan 5 obtained as part of the MAVERIC (Milky Way ATCA and VLA Exploration of Radio sources in Clusters) survey of globular clusters (L. Shishkovsky et al. 2020 ) to discover three new candidate spider millisecond pulsars. One of these has been subsequently confirmed as the redback pulsar Ter5ar (P. V. Padmanabh et al. 2024 ). How ev er, neither the frequency band of these VLA data (split between S and C bands) nor the spatial resolution (the data were taken in the hybrid BnA configuration, which is most similar to B configuration but with impr oved r esolution for lower declination sources, like Terzan 5) was optimized for studying pulsars.

Here we remedy these shortcomings with new, very deep Sband-only imaging of Terzan 5 in the VLA's highest-resolution A configuration, reaching rms noise levels and angular resolution both more than a factor of two better than existing data. These data allow us to both discover new candidate pulsars and set the best constraints yet on the total pulsar population of Terzan 5, indicating that hundreds of pulsars potentially await discovery.

Terzan 5 was observed with the VLA across two 5.5 hr blocks, one on 2022 March 14 (hereafter 'March') and the other on 2022 April 16 (hereafter 'April'), for a total observing time of 11 hr (Project ID: 22A-396). The data were taken at S band (2-4 GHz) as the best compromise between sensitivity and spatial resolution, split into two 1.024 GHz basebands centred at 2.5 and 3.5 GHz. The VLA was in its most extended A configuration. 3C286 was used as the bandpass and #ux calibrator. J1751-2524 was used as the phase calibrator.

We used the Common Astronomy Software Application ( & & ; J. P. McMullin et al. 2007 ;CASA Team et al. 2022 ) to perform #agging, calibration and imaging. Data were processed using the VLA Calibration Pipeline, with additional manual #agging where needed. The two basebands were imaged separately, however, both 5.5 hr observing blocks were concatenated together using & & 's tclean task. We performed multi-frequency synthesis imaging; nterms = 2 and a Briggs weighting with robust = 1 were selected.

Upon inspection of the image, we found significant artefacts as a result of the strong #ux density variations of the bright eclipsing r edback sour ce Ter5A (see Section 3.2 for the full details). We made time-resolved light curves, which were used to e x cise this sour ce fr om both baseband images using the uvsub command within & & . Finally, each image was primary-beam corrected. We note that Terzan 5 is not at the exact centre of the image: the phase centre is locat ed ∼ 1 w est of the core ( ≈ 90 per cent of the primary beam response at 2.5 GHz and ≈ 80 per cent at 3.5 GHz).

The final images have a restoring beam of 1 . 41 × 0 . 79 (PA = 0.3 • ) and local rms of ∼ 3 µJy beam -1 for a central frequency of 2.5 GHz and restoring beam of 1 . 05 × 0 . 58 (PA = 0.2 • ) and local rms ∼ 2 µJy beam -1 for 3.5 GHz. They are shown in Fig. 1 .

For consistency with U20 , we assume the same values for the Terzan 5 core radius ( 9 . 6 ; W. E. Harris 1996 ) and distance (5.9 kpc; E. V alenti, F . R. Ferraro & L. Origlia 2007 ). Our analysis of the radio data is conducted using observed #ux densities, so our inferences about the radio pulsar population do not depend on the specific distance assumed.

To clarify the nature of the published diffuse #ux measurement of the core of Terzan 5 from A. S. Fruchter & W. M. Goss ( 1990 ), we reimaged their archiv al C-configur ation r adio data from 1989 Jun 12, adding in additional data taken with the same setup and configuration obtained on 1989 Sep 10 (project code TEST), BnCconfiguration data from 1990 Oct 5 (project code AF206), Dconfiguration data from 1991 Apr 29 (project code TEST), and Bconfiguration data from 1990 Sep 7 (project code AF206). These data were all obtained in the continuum mode of the historic VLA, with two spectral windows each of 50 MHz bandwidth, centred at 1.465 and 1.515 GHz, for a mean frequency of 1.490 GHz. We reduced these data using standard routines in AIPS (E. W. Greisen 2003 ), and concat enat ed all data t ogether using DBCON . Imaging with a Briggs robust value = 0 yields an image with 10 . 2 × 7 . 3 (PA = -1 . 7 • ) synthesized beam and an rms noise of 34.9 µJy beam -1 . The resulting multi-configuration image has enhanced sensitivity to lower surface brightness emission and allows the bright pulsars Ter5A and Ter5C to be distinguished from the diffuse radio emission in the cluster core. We discuss the interpretation of this image in Section 4.2 .

For the new VLA A configuration images, we used the Python Blob Detection and Source Finder (PyBDSF; N. Mohan & D. Rafferty 2015 ) software package, version 1.9.0rc, to perform source finding. For consist ency, w e follow ed the same basic procedure described in U20 . The primary-beam corrected 2.5 GHz and 3.5 GHz images were searched independently. Within the halflight radius of Terzan 5 ( r < 0 . 73 ), we searched for, and selected, sources with a signal-to-noise (S/N) ratio of > 3 . Beyond the halflight radius we restrict our selection to more significant sources with S/N > 5 . As for previously published MAVERIC observations, we do not analyse sources beyond 3 . 7 from the cluster centre, equivalent to ∼ 4 half-light radii; objects beyond this radius ar e overwhelmingly back gr ound sour ces. All sour ces ar e listed in Table A1 ; source ID numbers 1-43 match those of U20 , with sources 44-89 being new sources not present in U20 .

After this independent source finding, we force-fit for the #uxes at the timing positions of pulsars that were not yet detected, and retained the #uxes in our catalogue only if the resulting detection was at least 2 σ . This led to #ux density measurements for five additional known pulsars (Ter5S, Ter5T, Ter5aa, Ter5ac, Ter5ah), which are listed in Table A1 .

All #ux densities are determined assuming objects are point sources: det ect ed sources are fit with Gaussians fixed to the dimensions of the synthesized beam. This is appropriate for true Terzan 5 sources, all of which are expected to be point sources in these data, but will underestimate the #ux of extended backgr ound sour ces, which ar e not our focus.

Owing to the back gr ound #ux of unr esolved pulsars in the cor e (Fig . 2 ), the rms noise in this region is slightly higher than outside the core, leading to a bright er det ection limit. We find 3 σ detection limits of > 11 . 2 µJy at 2.5 GHz and > 7 . 8 µJy at 3.5 GHz in the core, which we take as our pulsar completeness limits.

To test the astr ometry, we compar e our new VLA positions of 18 bright, isolated pulsars to their precisely known positions determined by previous pulsar timing surveys. We identify a small but significant shift in declination (Dec. = -0.14 arcsec) and a negligible shift in right ascension (RA = 0.008 arcsec), both of which are applied to the VLA coordinates. The rms residuals betw een the tw o sets of positions ar e RA = 0.10 ar csec and Dec.

MNRAS 547, 1-20 (2026) = 0.06 ar csec, r eassuringly consistent with the 1 σ positional uncertainties estimated from the synthesized beam.

The initial positional uncertainties are determined by & & imfit . How ev er, as suggest ed by the VLA Observational Status Summary 1 , we use 10 per cent of the full width at half-maximum (FWHM) of the synthesized beam as a minimum value for the positional uncertainties.

We cross-match the complete list of sources found in the two sub-band images. Sources are considered to match if their 3 σ error ellipses overlap. The final catalogue includes 84 sources, with 37 associated with 38 known pulsars (Table A1 ). We also recover the transitional millisecond pulsar candidate CX1 (A. Bahramian et al. 2018 ). The final International Celestial Reference System (ICRS) RA and Dec. of each source are calculated using the v ariance-w eight ed mean position. If a source is only det ect ed in one of the subband images, we report the 3 σ upper limit, calculated using the local noise measur ed fr om the PyBDSF rms images.

We calculate the spectral index α (defined as S ν ∝ ν α , where ν is the frequency and S ν is the #ux density at a given frequency) of each radio source. For non-detections, 3 σ upper limits are used. The spectral indices are modelled using the Bayesian Markov chain Monte Carlo software JAGS (M. Plummer 2012 ) assuming a power law model and a uniform prior on α between -3.5 and 3.5. The modelling self-consistently takes into account both measurements (including uncertainties) and 3 σ upper limits. The median of the posterior distribution of the spectral indices and 1 σ uncertainties are reported in Table A1 . For the sources only det ect ed in one subband, a 3 σ upper or lower limit on α is report ed inst ead; these are prior-dependent.

We follow the same pr ocedur e for finding X-ray counterparts to our radio sources that is outlined in U20 . Where available, we use the X-ray properties from U20 , otherwise we use the catalogue of A. Bahramian et al. ( 2020 ), which uses ∼ 750 ks of ar chival Chandr a /ACIS observations. All X-ray matches were fit with an absorbed power -la w model. Two absorption components were used; one fixed to the line-of-sight cluster absorbing column ( n H = 2 . 07 × 10 22 cm 2 ), while the second is left free to account for intrinsic absorption. Sources with X-ray matches have their 0.5-10 keV X-ray luminosities listed in Table A1 .

The shallower catalogue of U20 contained 43 sources, close to half of the current catalogue. 38 of these are found in our new catalogue. In Fig. 3 , we compare the 2.5 GHz #ux densities of the sources det ect ed at both epochs at this frequency as a function 1 ht tps://science.nr ao.edu/facilities/vla/docs/manuals/oss/performance/ positional-accuracy Figure 3. Top panel: Flux ratio between the 2.5 GHz 2022 March 14 and the 2.6 GHz 2012 observations for sources det ect ed at both epochs. Unfilled circles are confirmed pulsars. There is some evidence that the 2022 #ux densities for central sources are lower, likely due to scintillation. Bottom panel: Ratio of #uxes between March and April epochs. In both panels, red dotted lines represent the core and half-light radii. Sources within the half-light radius that are not confirmed pulsars are identified with their IDs. of radius. An immediate impression is that many of the sources within the half-light radius of the cluster are fainter in the new data compared to in U20 , though this is only of modest statistical significance, with a median #ux ratio of 0 . 74 ± 0 . 23 .

To the extent this is real, a possible explanation is refractive scintillation, expect ed t o hav e an rms value of 33 per cent at 2.5 GHz (using the methodology of P. J. Hancock et al. 2019 ). Hence the amplitude of the observed #ux variations is consistent with those expected from scintillation. There is existing evidence that the P. J. Hancock et al. ( 2019 ) model does a good job of r epr oducing the Ter 5 scintillation pr operties: it pr edicts an rms variation of 26 per cent at 1.75 GHz, compared to an observed value of 25 ± 4 per cent (A. R. Martsen et al. 2022 ). The sources at large projected radii do not share the apparent behaviour of the central sources, but this is not str aightforw ard t o int erpr et, as these ar e likely back gr ound sour ces that may have larger angular sizes that would not be expect ed t o show scintillation.

MNRAS 547, 1-20 (2026)

A cav eat t o the interpr etation of r efractive scintillation for the central #ux variations is that A. R. Martsen et al. ( 2022 ) did not find clear spatial correlations among the pulsar #ux variations, the opposite of the appar ently corr elated variations w e observ e. Unfortunat ely, there is no ov erlap in time betw een their published observations and our new continuum data that would allow a direct comparison of the #ux densities at a similar time.

The time-scale of refractive scintillation is harder to predict as it depends on the relative screen velocity, which is a priori unknown. Assuming a value of 100 km s -1 , it is 4.6 d at 2.5 GHz. Hence, in principle, #ux variations could also be det ect ed between the March and April epochs. We compare these in the bottom panel of Fig. 3 . There is no clear evidence for variation on this ∼ month-long time-scale, perhaps broadly consistent with the idea that the time differ ence corr esponds to sever al refr active time-scales, but these observ ations are not very constraining.

The consistency of the March and April #ux densities at all radii, as well as the large-radius consistency between 2012 and 2022, tog ether sugg est that there is no issue with the #ux density scale of the 2022 March observations. Future time series observations of Terzan 5 would be valuable to compare to predictions of scintillation models and see whether there is additional evidence of spatially correlated #ux variations.

Five of the sour ces fr om U20 ar e not present in the new catalogue. The non-detection of Ter5-VLA42 is expected-it is the transient low -mass X-r ay binary EXO 1745-248, which w as only det ect ed at the higher frequencies in the previous data. Ter5-VLA21 and Ter5-VLA25 are located far outside the half-light radius and had r adio spectr al indices α = -0 . 7 to -0.9: they are almost certainly back gr ound sour ces unassociat ed with the clust er.

The remaining two sources, Ter5-VLA34 and Ter5-VLA40, were both located inside the core radius and hence are very likely to be associated with the cluster. On the basis of its Xray luminosity ( L X ∼ 8 × 10 31 erg s -1 ; 0.5-10 keV), hard X-ray photon index, and poorly characterized radio spectral slope, U20 classified Ter5-VLA34 as a likely redback millisecond pulsar. While it is formally not det ect ed in our new data, there is some #ux present at this position in the 2.5 GHz image, with an appr o ximate value of ∼ 8 µJy; there is no #ux at 3.5 GHz. In the newer X-ray analysis of G. Kumawat et al. ( 2025 ) the source has a similar X-ray luminosity and is still found t o hav e a hard power law, with = 1 . 1 ± 0 . 5 . It seems plausible that the identification of Ter5-VLA34 as a spider is accurate but that it is potentially eclipsed in the new data, leading to it being fainter in the radio. In this case it may be discovered by future pulsar search observations. Alt ernativ ely, the X-ray source could be an accreting white dwarf, and the radio source an unrelated faint pulsar.

The case for Ter5-VLA40 is somewhat similar to Ter5-VLA34, in that U20 identified it as a likely redback pulsar and it is not present in the new catalogue. But there are some differences: its radio spectral index was well-measured in U20 and consistent with a pulsar ( α = -2 . 7 +0 . 8 -0 . 5 ), while its X-ray luminosity is lower ( L X ∼ [1 -2] × 10 31 erg s -1 ). Hence there is undoubtedly a yet-untimed pulsar at this location. If it is a spider, its non-detection in the new data could be explained as an eclipse, though as for Ter5-VLA34, it could also be a chance alignment between a non-spider pulsar and an unrelated X-ray source.

The U20 locations of T er5-VLA34, T er5-VLA40, and T er5-VLA42 are plotted in Fig. 1 .

Ter5-VLA31 was previously identified in U20 as an accreting stellar-mass black hole candidate due to its #at radio spectrum ( α = +0 . 2 ± 0 . 2 ), r adio/X-r ay #ux ratio, with an X-ray counterpart of L x = (2 . 1 ± 0 . 5) × 10 31 erg s -1 , and its location within the core.

In our new radio data, we find some evidence for a steeper r adio spectr al index ( α = -1 . 4 +0 . 6 -0 . 7 ), though with much larger uncertainties than before, owing to the lack of higher frequency data. The actual #ux density estimates are not formally inconsistent with the 2012 measurements: at 2.5 GHz we now find 27 ± 4 µJy versus < 20 µJy before, and at 3.5 GHz the measurements agree: 17 ± 3 µJy compared to 22 ± 4 µJy before. Hence the evidence that there is a change at all is not strong but only suggestive.

To the extent that the measured properties of the source have changed, one possibility is that a subset of the measurements are affected by source confusion in the core, which could vary between the epochs due to scintillation. We note that there is a new timed pulsar (Ter5av; P. V. Padmanabh et al. 2024 ) located near the position of Ter5-VLA31. A priori we do not expect to detect Ter5av in our data: its expected #ux density is ∼ 5 µJy at 2.5 GHz.

Examining the 2012 images at all fr equencies fr om U20 , VLA31 does not necessarily have a consistent position among them. At 2.6 GHz it cannot be confidently separated from the much brighter Ter5P. At 3.4 GHz there is a high-confidence source potentially consistent with Ter5av (separated by ∼ 0 . 25 -0 . 3 arcsec ), while the 5.0 GHz image shows a w ell-det ect ed ( 32 ± 3 µJy) source that definitely does not match Ter5av (separated by 0 . 96 arcsec ). P. V. Padmanabh et al. ( 2024 ) also suggest the possible association of the X-ray source that w e associat e with Ter5-VLA31 with Ter5av, but this pulsar is not in a short-period orbit and shows no eclipses, so it does not appear to be a spider pulsar that would be expected to have X-ray emission at the level observed.

One plausible explanation of these data is that there are two nearby sources, Ter5av and a #atter-spectrum source, and our measured #ux values for Ter5-VLA31 are a mix of these that vary with fr equency, r esolution, and time. It could instead be true that Ter5av is too faint to contribute to most of our data. In either case, the identification of the #atter-spectrum source as a candidate black hole is still plausible.

We can rule out an alt ernativ e explanation of Ter5-VLA31 as a transitional millisecond pulsar as its X-ray luminosity is much lower than the expected accretion state value of L X ∼ 10 33 -10 34 erg s -1 (e.g . M. Linar es 2014 ), and ther e is no evidence for X-ray variability over a large number of X-ray observations spanning 13 yr (G. Kumawat et al. 2025 ), including annual observations from 2011 to 2014.

Future deep C or X-band data would be ideal to reduce source confusion around this source by minimizing pulsar emission, revealing the #at-spectrum emission if present.

Figur e 4. Mar ch (top) and April (bottom) #ux density and spectral index measurements for the eclipsing redback Ter5A, with filled circles detections and unfilled circles 3 σ upper limits. The dashed lines represent the expected conjunction times for the companion to be in front of the pulsar. The length, depth, and timing of eclipses vary substantially in these data. Data used to create each light curve are reported in Table A2 .

Ter5A is the brightest radio continuum source in the cluster (it is at a projected radius of 0.6 from the centre and hence outside the field of view of Fig. 1 ). First discovered by A. G. Lyne et al. ( 1990 ), it is an eclipsing redback millisecond pulsar with orbital period of 1.8 hr. Previous studies of the eclipses have led to the conclusion that the predominant eclipse mechanism is fr ee-fr ee absorption from the ionized wind from the low-mass companion star (F. A. Rasio, S. L. Shapiro & S. A. Teukolsky 1991 ; S. E. Thorsett & D. J. Nice 1991 ).

Each of the two 5.5-hr e x ecution blocks of our VLA observations covers close to 3 orbital periods of Ter5A. We imaged Ter5A in each scan in both blocks, leading to a light curve with typical sampling of ∼ 6 -7 min for each block in both 2.5 and 3.5 GHz. These light curves are listed in Table A2 and shown in Fig. 4 . In each panel the inferior conjunctions of the companion are also marked, using the Ter5A ephemeris from A. C. Rosenthal et al. ( 2025 ). These are the epochs where the absorbing material might naively be expected to be thickest.

On 2022 March 14, the light curve captures the egress of one eclipse, two full eclipses, and most of a fourth. In each case the pulsar is not det ect ed at 2.5 GHz for a varying length of time, ranging from 40 min for the first full eclipse to only 20 min for the second full eclipse. We have taken these eclipse lengths as the time from the start of the first fully eclipsed scan to the end of the last fully eclipsed scan; they are necessarily quantized by the scan lengths. The timing of the eclipses is also not stable: the first full eclipse is centred at conjunction, while the second and especially the last eclipses start 'early', suggesting the absorbing material is concentrated ahead of the companion in its orbit, rather than trailing it.

In addition to this variation in eclipse lengths, the behaviour at 3.5 GHz also varies among the eclipses. In the first full eclipse, it is det ect ed in all but one scan, at an average #ux value ∼ 35 per cent of the preceding out-of-eclipse #ux. But it is undet ect ed for the last two scans of the shorter second full eclipse, and only det ect ed in a single scan in the last eclipse.

Variations in the eclipse properties are also seen in the 2022 April 16 light curve (Fig. 4 ,bottom). This second light curve also shows a 'bonus' eclipse occurring around superior conjunction of the companion, only 30 min after the previous eclipse ends, as well as an unusually long complete eclipse at both frequencies, lasting at least 1.1 hr, at the end of the data set.

In March the spectral indices show that absorption is present at the majority of epochs even outside of the obvious eclipses: as the 2.5 GHz #ux is more absorbed than the 3.5 GHz #ux, the spectral index becomes less steep or in a few cases is even inverted such that the pulsar is brighter at higher frequency. In April there is less evidence for absorption outside of eclipses.

Eclipses of varying length and sometimes at unexpected phases have been previously observed for Ter5A in continuum at lower fr equencies (e.g . S. E. Thorsett & D . J . Nice 1991 ; O . M. Smirnov et al. 2025 ). These new data show multiple eclipses of Ter5A in continuum at the relatively high frequency of 3.5 GHz, though we note that another redback, PSR J1740-5340 in the globular cluster NGC 6397, has shown continuum eclipses at an even higher frequency of 5.5 GHz (Y. Zhao et al. 2020 ). Our high-S/N multiband measurements of eclipses in Ter5A should prove useful in constraining the physical parameters of the eclipsing material in a future work.

Ther e ar e curr ently 49 pulsars in Terzan 5 that have been confirmed via pulsar timing (P. V. Padmanabh et al. 2024 ), of which 48 (all but Ter5al) have timing positions. Of these 48 pulsars, 38 ar e pr esent in our point sour ce catalogue as 37 separate sour ces: the pulsars Ter5Z and Ter5ae are too close to be separated in our data, and are det ect ed t ogether as a single source.

While few of the timed pulsars have been observed at 2.5 GHz, we can predict their 2.5 GHz #uxes from timing #uxes calculated from the radiometer equation. The most reliable of these are pr obably the long -term 2.0 GHz #ux es and spectral indices for 31 pulsars listed in A. R. Martsen et al. ( 2022 ), but for the 17 pulsars not included in that paper, we take the 1.3, 1.4, or 2.0 GHz #ux density from the relevant papers (S. M. Ransom et al. 2005 This comparison is given in Fig. 5 . Similar to the comparison between 2012 and 2022 continuum #uxes in Fig. 3 , we see that the 2022 continuum pulsars are mostly fainter than predicted by their timing #ux. For T er5A, T er5O, and T er5P this could be due 2022) for the 31 pulsars where this information is available. For the remaining pulsars that do not have timing spectral information, we use an assumed α = -1 . 8 . Upper limits are denoted with unfilled circles. Red and black points indicate eclipsing redback (Ter5A, Ter5P, Ter5ad, Ter5ar) and black widow (Ter5O, T er5aq, T er5at) sour ces, r espectively. The diagonal dashed line r epr esents unity. The horizontal dotted line indicates our new candidate pulsar Ter5-VLA62. We label notable pulsars (e.g . r edbacks, bright sour ces and outliers).

to eclipses, but that is not the case for the non-spider pulsars. Since there is not a similar offset between the U20 2012 continuum #uxes and the timing -pr edicted #ux es, scintillation in our data is likely the main culprit (Section 3.1.1 ). How ev er, it may not be the only mechanism at play. A. R. Martsen et al. ( 2022 ) find that many of these pulsars hav e #att er spectra than typical, with a mean spectral index of α = -1 . 35 . They attribute this to a selection effect from observations at 2.0 GHz rather than the lower frequencies typical of pulsar sear ches. Ther e ar e at least two pulsars in our data (T er5M and T er5Y) that we find to have much steeper spectral indices than in A. R. Martsen et al. ( 2022 ) at high confidence. It will be worth seeing whether such differences persist in future data sets, or whether they are due to scintillationinduced variations in our measured spectral indices.

Of the t en undet ect ed pulsars, nine (T er5U, T er5ak, T er5am, T er5ap, T er5aq, T er5as, T er5at, T er5av, T er5aw) have predicted 2.5 GHz #ux densities from ∼ 3 -10 µJy, so their non-detections ar e e xpected.

The last undet ect ed pulsar is the redback Ter5ad, which was also not det ect ed in U20 . Giv en its report ed 1.95 GHz #ux density of 80 µJy (J. W. T. Hessels et al. 2006 ), for a typical spectral index this would imply a #ux density of ∼ 50 µJy at 2.5 GHz and ∼ 30 µJy at 3.5 GHz, where it should be easily det ect ed in our new data. While we cannot rule out the possibility that it is eclipsed in all our data sets (both old and new), it is plausible that its intrinsic #ux density is fainter than the value reported in J. W. T. Hessels et al. ( 2006 ), and a combination of intrinsic faintness and eclipses is preventing its detection.

All the det ect ed pulsars in the centre of the clust er are labelled in Fig. 1 .

In Table A1 , there are seven continuum sources in the core without a clear association. Of these seven, the candidate spider pulsars Ter5-VLA34 and Ter5-VLA40 from U20 were discussed above, as was the candidate black hole Ter5-VLA31.

Of the remaining four sources, Ter5-VLA62 is the only unassociat ed st eep -spectrum source within the core and is a new candidate millisecond pulsar. Ter5-VLA62 is significantly det ect ed in our lower subband (2.5 GHz; 18 ± 4 µJy) and is not coincident with a known pulsar. This source was not det ect ed in U20 , and does not have an X-ray counterpart.

T er5-VLA86, T er5-VLA88, and T er5-VLA89 are all det ect ed only at 3.5 GHz, and were not detected by U20 . Each source lies in a crowded region of the core and at 2.5 GHz appears to be contaminated by emission from neighbouring bright pulsars, which could hide real 2.5 GHz #ux. Hence we cannot rule out the possibility that these sources are indeed pulsars. Of these, the only one with a potential X-ray match is Ter5-VLA89, located 0 . 7 ar csec fr om the X-ray sour ce CX9. This X-ray sour ce is classified as a quiescent neutron star low-mass X-ray binary (G. Kumawat et al. 2025 ). Such binaries are not expect ed t o show radio emission and the offset is large enough that these sources are probably unrelated.

Given the large pulsar population of Terzan 5 and the detection of radio continuum sources at our detection limit in the 2012 data, the central scientific goal of the new observations was t o bett er charact erize the population of pulsars by increasing both the depth and spatial resolution of the VLA imaging.

We construct an empirical pulsar luminosity function from our continuum data, starting from the 37 continuum sources present in our point source catalogue. We split the combined #ux density of the overlapping pulsars Ter5ae and Ter5Z based on their relative 1.4 GHz #ux densities (A. R. Martsen et al. 2022 ), assigning a #ux density of 36 . 6 ± 3 . 8 µJy to Ter5ae and 19 . 6 ± 3 . 8 µJy to Ter5Z. This gives 38 confirmed pulsars with 2.5 GHz continuum #ux estimates.

As highlighted above in Section 3.4 , Ter5-VLA62 has not yet been det ect ed as a timed pulsar, but has a steep radio spectrum and is located in the cluster core. We consider Ter5-VLA62 to be high-probability pulsar and include it in our luminosity function, giving 39 pulsars. Two additional pulsar candidates identified in U20 (Ter5-VLA34 and VLA40) are not det ect ed in the new imaging; these are plausibly real pulsars that are fainter in these data due t o scintillation, but w e do not include them here.

Fig. 6 shows the 2.5 GHz continuum luminosity function of the 39 high-probability pulsars, with a power -la w fit to the complete ( > 11 . 2 µJy) sample ov erplott ed. The superficial turnov er in the data at fainter #uxes is due t o incomplet eness, though there is a hint of a change in slope around ∼ 20 -30 µJy, above the completeness limit of the sample.

To further inv estigat e the pulsar population, we next calculate the total #ux of pulsars in the cluster, including unresolved #ux. ). The vertical grey shaded region indicates where we are incomplete ( < 11 . 2 µJy). The obvious turnover in the luminosity function is due to this incompleteness, though there is a hint of a change in slope in the more complete region.

Fig. 2 shows that there is w ell-det ect ed 'residual' #ux from the pulsars that are fainter than the detection limit of our new A configur ation observ ations. This #ux is apparent because in many beams ther e ar e a sufficient number of sub-threshold pulsars whose summed #ux is detectable. We can constrain the total pulsar population in Terzan 5 by estimating the int egrat ed amount of #ux in these subthreshold pulsars.

Owing to the limited sensitivity of our A configuration observations to low surface brightness emission, we are not sensitive to the #ux from subthreshold pulsars located at large distances fr om the centr e. Ther efor e, t o calculat e the total diffuse #ux, w e hav e t o choose a radius in which w e hav e a high-S/N measurement of the diffuse #ux and then extrapolate to get an integrated measurement over all radii.

To do this, we assume the spatial distribution of the faint pulsars producing the diffuse #ux follows the brighter resolved pulsars. To be self-consist ent, w e only consider the radial distribution of the detected pulsars in our data. Of the 39 pulsars, we find that about half (20/39; 51.2 per cent) are found within a radius of 9.6 arcsec (the core radius), and about three-quarters (29/39; 74.4 per cent) are found within a radius of 13.3 arcsec. We int egrat e the unresolved #ux in our data out to these radii, and extrapolate based on the above fractions to estimate the total unresolv ed #ux int egrat ed ov er all radii. The uncertainties in these quantities include both thermal noise and the uncertainty in the extr apolation. Separ at ely, w e add up the resolved #ux from pulsars out to these radii, but when calculating the total integrated resolv ed #ux ov er all radii, w e use the actual (not extrapolated) measurements, since pulsars above our detection limit should be detectable at all radii within our field of view. Table 1 shows the int egrat ed resolv ed (from det ect ed sources list ed in our catalogue) and unresolved #ux within these radii, as well as an extrapolated v alue integr ated over all radii. The unresolved #ux is bright: it is comparable to or larger than the summed #ux of resolved pulsars (e x cluding Ter5A and Ter5C). We get consistent values from both radius choices, in the range 1.8-2.0 mJy, consistent to within the listed uncertainties in the quantities.

Our new measurements are the first that allow simultaneous measurements of the resolved and diffuse #ux at the same frequency and hence are w ell-suit ed for modelling the pulsar population.

We acknowledge several additional sources of uncertainty in our new measurements. First, our 2.5 GHz residual image (Fig. 2 ) shows imperfect subtraction of the point sources (visible for Ter5C, right outside the core). To assess how this might affect our diffuse #ux measurement, we compared the summed diffuse + r esolved sour ce #ux measur ement to the int egrat ed #ux out t o 9.6 and 13.3 arcsec in the original unsubtracted 2.5 GHz image. We find agreement within 4 per cent for 9.6 arcsec and 1 per cent for 13.3 arcsec. The differences are likely a combination of the imperfect subtraction of point sources as well as edge effects: several sources cross these radial limits and are only partially counted in a dir ect pix el integration. These small differences are subdominant to the other uncertainties in our measurement so they do not meaningfully affect our analysis.

A second source of uncertainty is that A-configuration images lack sensitivity to low surface brightness emission, so it is difficult to be sure that we are sensitive to 'all' the emission within our chosen radii. To address this pot ential concern, w e ne xt compar e our results with those from archival lower angular resolution data. To anticipate these results, we find no evidence that we are missing substantial diffuse #ux, at least within the core.

The most direct comparison is to the results of A. S. Fruchter & W. M. Goss ( 2000 ), who reported that their 1.49 GHz A-configuration images had ∼ 2 mJy 'central' #ux (e x cluding Ter5A and Ter5C) that was diffuse compared to their beam of 5 . 9 × 4 . 1 arcsec. It is difficult to precisely define the outer edge of the A. S. Fruchter & W. M. Goss ( 2000 ) diffuse #ux measurement, but it does not appear t o ext end bey ond the 10 . 7 ar csec pr ojected radius of Ter5C. When combined with the spectral index they measured for the diffuse #ux of α = -1 . 65 ± 0 . 16 , this corresponds to a predicted MNRAS 547, 1-20 ( 2026) 'central' #ux of ∼ 850 µJy at the 2.5 GHz central frequency of our new observ ations. This v alue is about 25 per cent smaller than, but generally consistent with, the measurement of 1.1 mJy that we find within the cluster core, especially given that the outer boundary of their measurement is not well-defined.

The earliest paper with radio continuum observations of Terzan 5 (A. S. Fruchter & W. M. Goss 1990 ) reported an int egrat ed #ux density of 1 . 9 ± 0 . 2 mJy at 1.49 GHz (beam of 18 × 10 arcsec ) from C-configuration observations. While this value agrees with their 1.49 GHz A-configuration measurement of diffuse #ux, it was not clear to us how far out from the centre this measur ement e xtended, nor whether it successfully separated Ter5C from the diffuse #ux. This was our motivation for reimagining these older data, together with other old data taken in the same setup but a r ange of configur ations, to enable sensitivity to both diffuse #ux and the two bright pulsars (see description of this additional imaging in Section 2.1 ).

In this multi-configuration image (Fig. 7 ), which has a large beam ( 10 . 2 × 7 . 3 ) that is comparable to the core radius, the diffuse central emission is separated from Ter5C 2 , and has a #ux density consistent with the value of 1 . 9 ± 0 . 2 mJy reported by A. S. Fruchter & W. M. Goss ( 1990 ). The diffuse emission is not clearly det ect ed bey ond the core.

These comparisons of measurements made at different angular resolutions show that there is no evidence our new Aconfiguration data are missing #ux in the core of the cluster, corresponding to the 9.6 arcsec r adius of our centr al measurement.

2 In this image, Ter5C has a 1.49 GHz #ux density of about 1.8 mJy, much more luminous than the value of 1.1 mJy reported by A. R. Martsen et al. ( 2022 ), likely due to scintillation.

We cannot make a confident statement about whether this is also true of our int egrat ed measurement out t o 13.3 arcsec, owing t o the noise levels in these old data as well as the effects of Ter5C. Since the extrapolated total measur ements fr om differ ent radii in our data agree to within ∼ 10 per cent , to the extent we trust the inner measurement, we can have confidence that the large radius measurement is not too much in error.

We can now proceed to use these unresolved #ux measurements in our modelling of the pulsar luminosity function.

We first fit a hierarchical Bayesian power -la w model with free parameter α p (where d N/d L = L -αp ) t o the observ ed (in the VLA data) individual pulsar #ux densities and uncertainties and the observed number of pulsars, assuming the luminosity function is censored below the 3 σ limit of our 2.5 GHz catalogue ( 11 . 2 µJy).

These model fits are implemented using the Bayesian Markov Chain Monte Carlo software JAGS (M. Plummer 2012 ). We directly fit the #ux densities (rather than conv erting t o pseudoluminosities), but for ease of reading still refer to a 'luminosity function'. No binning is used.

Fitting the 33 pulsars in our continuum catalogue above this limit, we find α p = 2 . 49 +0 . 29 -0 . 25 if pulsars A and C are excluded, and, as expected, a #atter slope ( α p = 2 . 16 +0 . 22 -0 . 19 ) if they are included. Note that there is a large gap between Ter5C (319 µJy) and the third-brightest pulsar in our data (Ter5P; 51 µJy).

N e xt we extend the model to also fit the measured unresolved #ux below our continuum sensitivity limit. This r equir es an assumption of the minimum #ux density, which we take to be either 0.1 or 1 µJy. These values correspond to 2.5 GHz pseudoluminosities of 3.5 or 35 µJy kpc for a distance of 5.9 kpc, equivalent to 1.4 GHz pseudo-luminosities of about 0.01 or 0.1 mJy kpc for an individual pulsar spectral index of α = 1 . 8 .

We perform model fits within each of the radii defined above (9.6 and 13.3 arcsec) for which we also have measurements of the unresolved #ux listed in Table 1 . Since we also censor the observed pulsars at 11 . 2 µJy, the number of pulsars fit is slightly smaller than discussed above within these radii, 19 and 25, respectiv ely. In these fits, w e also obtain a normalization that corresponds to the total pulsar population N, including faint unobserved pulsars.

Assuming a minimum #ux density of 1 µJy, we find a powerlaw slope α p = 2 . 08 +0 . 11 -0 . 10 ( α p = 2 . 05 +0 . 09 -0 . 09 ) and an all-cluster pulsar population of N = 489 +61 -57 ( N = 408 +46 -47 ) inferr ed fr om the inner 50 per cent (75 per cent) samples. For a minimum #ux density of 0.1 µJy, the best-fitting slopes are slightly #atter, with α p = 1 . 87 ± 0 . 07 ( α p = 1 . 85 ± 0 . 07 ), but as expected the total pulsar population is much larger, N = 2. Results of lognormal fits to the VLA data. Each r adial r ange has two rows, one for each minimum pulsar #ux density assumed (1 or 0.1 µJy). The listed uncertainties on the lognormal parameters µ and σ and the total number of pulsars N are equal-tailed 68 per cent intervals. For µ, because there is support to the low edge of its prior (-2.5), the listed results are prior-dependent: see main text.

Min. µ σ N (arcsec) ( µJy) (log µJy) (log µJy) < 9 . 6 1 -1 . 24 +0 . 99 -0 . 85 0 . 97 +0 . 17 -0 . 24 419 +61 -53 0.1 -1 . 03 +0 . 90 -0 . 96 0 . 98 +0 . 22 -0 . 26 1028 +338 -357 < 13 . 3 1 -1 . 36 +0 . 91 -0 . 79 1 . 02 +0 . 16 -0 . 21 349 +48 -40 0.1 -1 . 16 +0 . 85 -0 . 87 1 . 04 +0 . 20 -0 . 23 875 +246 -256 the distance of Terzan 5. Adopting this more realistic minimum #ux density would imply ∼ 1500 -2000 pulsars, approaching the total maximum population of neutron stars that might plausibly be present in a cluster like Terzan 5 (C. S. Ye et al. 2020 ), most of which are unlikely to have been recycled as observable pulsars.

Hence, while we cannot formally reject a power -la w model on the basis of available data, it seems much more likely that the pulsars fainter than those found in current searches have a more complex luminosity function.

An alt ernativ e model for the luminosity function of pulsars is a lognormal distribution. A low-N, censored lognormal distribution can look similar to a power -la w, so they are difficult to distinguish when only part of the distribution is observed.

As for the power -la w model, we first fit a lognormal model to the individual pulsars without consideration of the unresolved #ux. In a lognormal model the log of the #ux densities follows a normal distribution with mean µ and standard deviation σ . Our fits use minimally informative priors on µ ([-2.5,2.5], corresponding to mean #ux densities from about 0.003 µJy to 316 µJy) and σ ([0,4]). These priors ar e br oader than used in other recent papers modelling Terzan 5 pulsars (e.g. J. Chennamangalam et al. 2013 ;J. Berteaud et al. 2024 ). Nonetheless, µ and σ strongly covary and there is still support for µ to the low end of the prior, even though the prior range is large. How ev er, w e not e that high values of µ that might superficially appear to fit the luminosity function in Fig. 6 (e.g. µ ∼ 1 . 3 , σ ∼ 0 . 5 ) turn over too quickly at the faint end: they produce only ∼ 100 µJy of unresolved #ux, about a factor of 5-7 lower than observed. Hence if the lognormal model is correct, the mean is at lower values.

We next expand the models to include the unresolved #ux, again assuming minimum #ux densities of either 1 or 0.1 µJy and fitting the unresolved #ux within the radius containing 50 per cent or 75 per cent of the pulsars. The parameters obtained from the two different radii are consistent within the uncertainties. The results of these fits are given in Table 2 , and the results for the larger radius are plotted in Fig. 8 ; the plots for the smaller radius are similar.

In all cases, because there is at least some support to the low edge of the µ prior, the posteriors for µ and σ are priordependent. For this reason, the exact values of µ and σ listed in Table 2 should not be taken too seriously. For all the models there is posterior support for values of µ that are less than the minimum #ux density; in these cases the pulsars would solely be drawn from the right tail of the distribution.

The inferred number of pulsars is less prior-dependent than µ or σ , though it does depend on the assumed minimum #ux density. In all cases the number of inferred pulsars is large: the lowest median value for any of the fits (75 per cent radius, 1 µJy minimum) is N = 349 +48 -40 , while the median N for the more realistic minimum #ux density of 0 . 1 µJy is at least double this. Even the 1 per cent quantile of the lowest-N fit is N = 264 , implying that 200 pulsars remain to be discovered in Terzan 5. This is a str aightforw ard consequence of the large amount of stillunresolved #ux: since it is primarily produced by pulsars fainter than the current detection limit, the number of pulsars needed is large.

In addition to fitting these models with a single censored #ux density limit, we also considered more realistic models where the detection probability (in our data) is given by a logistic distribution. This distribution monotonically increases from a low probability of detecting a very faint source to a high probability of detecting a bright pulsar through a transition region of int ermediat e det ection probability. While the detection probability is not known a priori as a function of #ux density, we considered a reasonable range of parameters that seemed broadly consistent with our data. For both lognormal and power law models, we found pulsar population numbers slightly lower than but consistent with the values obtained from the simpler models. Since a more complex censoring model has poorly constr ained par amet ers but giv es qualitativ ely similar results, w e keep the results from the simpler censoring model fits as those to report.

In the previous subsection we found that a lognormal model was consistent with the available data, but that the model parameters wer e r elatively poorly constrained. Here we consider whether the pixel distribution of the residual #ux in the present data contains additional information about the luminosity function of subthr eshold sour ces.

It is definitely the case that the properties of subthreshold pulsars will affect the spatial variations of the residual #ux. While a given #ux could in principle be produced by a smaller number of brighter pulsars or a larger number of fainter pulsars, these two cases will show differing Poissonian #uctuations, with a larger signal for the smaller number of brighter pulsars; higher angular resolution will also increase the signal, in a manner akin to surface brightness #uctuations in galaxy light (J. P. Blakeslee, E. A. Ajhar & J. L. Tonry 1999 ).

In practice, the shot noise signal will be strongly affected by the thermal noise #uctuations in the image, which are larger than the #ux densities of most of the faint pulsars in which we are int erest ed. Hence w e proceed by a more elaborate route: simulating first the intrinsic pulsar population and then a radio image of this population, which we finally analyse in a manner akin to the original data as described in Section 2 .

We have discussed the radii containing an observed fraction of about 50 per cent or 75 per cent of the known pulsars, but have not otherwise modelled their radial distribution. To do this, we MNRAS 547, 1-20 ( 2026) assume the generalized King model previously used for the distribution of X-ray sources in Terzan 5 (C. O. Heinke et al. 2006 ). This model assumes a surface density (r) ∝ [1 + (r/r c ) 2 ] (1 -3 q ) / 2 , where r is the projected radius, r c is the core radius (for Terzan 5, 9.6 arcsec, identical to observed radius containing half the pulsars in our data), and q is a free parameter that r epr esents the mean mass ratio between the pulsars and the visible cluster stars. Using the 48 millisecond pulsars with precise positions, we find a bestfit ting v alue of q = 1 . 74 ± 0 . 16 . Using this model, the projected radius that contains 50 per cent (75 per cent) of the millisecond pulsars is 8.3 arcsec (13.3 arcsec), v alues toler ably close to those observed.

Pulsars are expect ed t o be point sources superposed on the thermal noise. Thus we first simulate a noise map with rms equal to the observed 2.5 GHz rms ( 3 µJy beam -1 ) at the location of Terzan 5. We use the & & task simobserve to simulate our underlying noise map; we use a single noise map for all simulations.

All pulsars are point sources and thus are represented with 2-D Gaussians matching the synthesized beam ( 1 . 4 × 0 . 8 ). As in Section 4.3.2 , we assume a lognormal analytic model to assign #uxes, but do so over a uniform grid in µ and σ : we do not use any of the r esults fr om the pr evious analysis in Section 4.3.2 for these simulations. We perform 100 simulations per set of parameters, and assume a minimum #ux of 0 . 1 µJy in all cases.

Sour ce positions ar e drawn fr om the pr eviously discussed King model (Section 4.4.1 ). While w e experiment ed with allowing the value of the relative concentration parameter q to vary, we found that this did not significantly affect any results for reasonable variations in q . This is perhaps expected given our focus on the cor e. Hence we fix ed q to the previously derived best-fitting value of q = 1 . 74 for all the results reported below.

To construct the simulat ed images, simulat ed sour ces ar e added until the cumulative source #ux of the cluster within the core ( r < 9 . 6 ) matches the observed total (residual plus source) #ux within the same radius. This essentially sets the normalization of our analytic models. N e xt, all significant ( > 3 σ ) sources ar e e xtracted, ag ain in line with our pr ocedur e with the real data. We emphasize that to make this process realistic, no prior information on the added sources is used. Hence, due to source blending, the sources extracted will not necessarily match the sources added. After all the significant sources are remov ed, w e hav e a simulated residual image, whose properties can be compared to MNRAS 547, 1-20 (2026) the observed residual image. We repeat this process to produce 100 simulated images for each set of parameters.

To compare the simulated and real images, we calculate the onedimensional Kolmorgorov-Smirnov (KS) statistic on the pixel #uxes: this is the maximum difference of the cumulative luminosity functions of the pixel values of two images. Since our focus is on r epr oducing the observed r esidual emission, we only consider the pixel #ux distributions within the core radius ( ≤ 9 . 6 ), though the radial distribution of generated sour ces e xtends to larger radii in our simulations. Since there are 100 simulated images per parameter pair, there are 100 corresponding KS statistics per parameter pair, as each simulated image is compared to the real image separately.

In Fig. 9 we show the median KS statistic from comparing the simulated to real data as a function of µ and σ . The lowest median KS statistics are found in a broad swath of lognormal model parameter values that appear broadly consistent with the r esults fr om the r esolved + unr esolved luminosity function fitting in Section 4.3.2 . In detail, these new image simulations prefer slightly lower σ at fixed µ (or, alt ernativ ely, low er µ at fixed σ ) than the direct luminosity function fitting, but the overall conclusions are similar: a large number of undiscovered pulsars are needed to explain the residual #ux.

We have plotted the KS statistics and not the results from the KS test , and hence the values in Fig. 9 are relative rather than absolute. This is because there are multiple pixels across a synthesized beam, making the individual simulated pixels correlated. Thus a standard one-dimensional KS test (where the samples are independent) is not applicable: that is, our calculated KS statistics cannot be immediately used to accept or reject a null hypothesis that the two pixel distributions are drawn from the same underlying distribution at a particular confidence level.

Inst ead, t o help int erpret these KS statistics, w e also calculat ed the KS statistics in a pairwise fashion among the 100 simulations themselves for each µ/ σ set, making these measurements between the 4950 100 2 unique pairs. The typical median KS statistics for these parameter pairs were ∼ 0 . 03 -0.04. Since these simulations r epr esent identical model parameters, the implication is that when comparing the real image to the simulated images, median KS statistics 0 . 03 -0 . 04 imply the observed and simulated residual #ux distributions are not inconsistent. This condition is met by the central trough of dark colours in Fig. 9 .

As examples to illustrate the simulations, in Fig. 10 we show three pairs of images, with the left side the original or simulated image, and the right side the residual image. The top row is the real data, while the middle and bottom rows represent examples of the range of acceptable simulations. In each case it is clear that the simulated residual images appear to be reasonable matches for the observed residual #ux distribution in the core.

Fig . 11 compar es the r eal observed luminosity function to the simulated luminosity function for the same parameter pairs shown in Fig. 10 . While each model does a reasonable job of representing the brighter detected pulsars, the differing behaviour among the faintest pulsars is why the number of pulsars inferred is so different ( N ∼ 2000 versus N ∼ 600 ) in the two cases.

Unfortunately, at least here it does not appear that the image simulations offer much stronger constraints on the pulsar luminosity function than modelling the resolved pulsars + unresolved #ux: deeper and/or higher resolution data are needed to detect the necessary shot noise signal. We emphasize that these simulations do not prove that a lognormal model is the correct one, just that it can produce residual #ux distributions not inconsistent with that observed.

One disadvantage of using radio data to constrain the pulsar population in Terzan 5 is that there is no str aightforw ar d r elationship between the radio luminosity and the fundamental physical parameters of the pulsar such as the spindown luminosity.

In principle, the γ -ray emission from Terzan 5 can offer complementary constraints on its pulsars, as the γ -ray luminosity of a pulsar is essentially uncorrelated with the radio luminosity (D. A. Smith et al. 2023 ). Instead, the γ -ray luminosity is closely linked to the pulsar spindown luminosity (e.g. A. A. Abdo et al. 2013 ). This can be both an advantage and a disadvantage: while γrays give more direct physical information about the pulsars, the int egrat ed γ -ray #ux from a globular cluster can be dominated by even a single source if a young pulsar is present, as in NGC 6624 (P. C. C. Fr eir e et al. 2011 ). Terzan 5 has no known young pulsars.

A. A. Abdo et al. ( 2010 ) used the measured γ -ray luminosity of Terzan 5, an assumed mean spindown luminosity ˙ E = 1 . 8 × 10 erg s -1 , γ -ray conversion efficiency of 8 per cent to estimate a total population of N = 180 +100 -90 pulsars. The most updated measurement of the 0.1-100 GeV #ux of Terzan 5 is 9 . 2 ± 0 . 2 × 10 -11 erg s -1 cm -2 (A. Amerio, D. Hooper & T. Linden 2024 ). At a distance of 5.9 kpc, this is equivalent to L γ = 3 . 8 × 10 35 erg s -1 . Notably, this is the largest γ -ray luminosity of any cluster in the A. Amerio et al. ( 2024 ) sample by ∼ 50 per cent, consistent with the fact that Terzan 5 has the largest known population of pulsars of any Galactic cluster. Simply applying the mean L γ = 1 . 44 × 10 33 erg s -1 assumed in A.

MNRAS 547, 1-20 ( 2026) MNRAS 547, 1-20 (2026) Figure 11. Comparison of observed (red) and simulated (blue) pulsar luminosity functions for the lognormal models shown in Fig. 10 : the top model has ( µ = -2 . 2 , σ = 1 . 1 ) while the bottom model has ( µ = 0 . 3 , σ = 0 . 4 ). The observed pulsars are those within 13 . 3 arcsec of the centre. We overlay 25 samples of each model t o demonstrat e #uctuations, and the grey shaded region indicates the 11 . 2 µJy incompleteness limit. Both models r equir e a large population of faint undiscovered pulsars, and the e xtr emely larger number of ultra-faint pulsars in the top model is evident.

A. Abdo et al. ( 2010 ) to this updated total luminosity gives an estimated N ∼ 260 .

Using a different methodology, A. Amerio et al. ( 2024) find a best-fitting mean L γ of 1 . 8 × 10 33 erg s -1 for globular clusters by modelling the γ -ray emission together with the number of known pulsars and the interaction rates. This would imply N ∼ 180 . How ev er, giv en that some clusters modelled (including Terzan 5) likely have many more pulsars than currently known, this reported mean luminosity will be biased higher than the true mean. It is also the case that quite a few field radio millisecond pulsars are not det ect ed as γ -ray pulsars, suggesting the typical observed γ -ray efficiency could be higher than the true mean efficiency (D. A. Smith et al. 2023 ). Both factors would t end t o reduce the mean L γ and lead t o a great er inferred pulsar population.

As a limiting case, consider the lognormal fits to the radio luminosities above, which allow a population of ∼ 1000 pulsars depending on the assumptions. This would r equir e a mean L γ of 4 × 10 32 erg s -1 in Terzan 5. While lower than found by A. Amerio et al. ( 2024 ), this value is close to that inferred from modelling disc millisecond pulsars by I. Holst & D. Hooper ( 2025 ), which is L γ ∼ 6 × 10 32 erg s -1 .

Another approach is to compare Terzan 5 to a cluster whose millisecond pulsar census is likely more complete: 47 Tuc. This w ell-studied clust er has a minimum of 42 millisecond pulsarsfoot_2 and a γ -ray luminosity that is a factor of ∼ 5 . 8 lower than Terzan 5 (A. Amerio et al. 2024 ). Making the questionable assumption that the pulsars in the two clusters have the same mean L γ would then imply N 240 in Terzan 5.

Due to the uncertainties on the L γ luminosity function, it does not appear that the γ -ray properties offer strong independent constraints on the pulsar population in Terzan 5 at present. How ev er, the γ -ray data do not contradict our radio finding that a large number of pulsars are yet undetected.

Based on #ux densities measur ed fr om pulsar timing over many epochs, A. R. Martsen et al. ( 2022 ) found that their cumulative 1.4 GHz luminosity function showed an apparent change in slope near the faint end of their luminosity function (at a 1.4 GHz #ux density of ∼ 23 µJy), which they attributed to observational incompleteness. The most recent Terzan 5 pulsar timing paper, P. V. Padmanabh et al. ( 2024 ), brie#y revisits the luminosity function of the timed pulsars, again finding a potential slope change in the cumulative luminosity function at a 1.28 GHz #ux density ∼ 25 µJy, which they argue is 'suggesting that we are approaching a minimum luminosity cutoff'.

We disagree: we think A. R. Martsen et al. ( 2022 ) were correct in attributing a slope change to observational incompleteness. From the radiometer equation, P. V. Padmanabh et al. ( 2024) estimat e their det ection limit t o be in the range of ∼ 11 -17 µJy at 1.28 GHz, depending on the overlap between synthesized beams (they tiled Terzan 5 with 288 separate beams out to a radius of 3 ). Uncertainties in the pulsar duty cycle, scintillation, and scattering at L band imply that the detection efficiency at this listed limit is below 100 per cent. Indeed, given how close this detection limit is to the apparent slope change in the cumulative luminosity function at ∼ 25 µJy, observational incompleteness is a more str aightforw ard explanation, r ather than an actual change in the luminosity function.

We emphasize it is possible that the luminosity function is indeed changing its slope around the current pulsar detection limits, but that the presently available data are too uncertain to clearly establish this. This is primarily due to the modest #ux range of existing data and the small numbers of pulsars: e x cluding pulsars A and C, the remaining sources cover at most 1 dex in #ux density to the detection limit. Our continuum 2.5 GHz luminosity function (Fig. 6 ) allows a change in the slope of the luminosity function near our 2.5 GHz detection limit, but does not demand it.

What is not uncertain is whether current pulsar searches are near a cutoff in the luminosity function. As we have shown, the amount of residual #ux is still so large that many pulsars must be undet ect ed-and also much fainter-than known pulsars to selfconsistently explain the observations. This does not con#ict with MNRAS 547, 1-20 ( 2026) other physical constraints or observations: as mentioned above, there is at least one field millisecond pulsar that would have a 2.5 GHz #ux density around 0 . 1 µJy at the distance of Terzan 5 (P. Wang et al. 2021 ). In addition, since the radio luminosity of a pulsar is an insignificant fraction of its spindown luminosity, there is no particular physical reason to expect current radio observations ar e r eaching such a luminosity cut off. Inst ead, w e need to look to observations themselves to constrain models.

M. Bagchi, D. R. Lorimer & J. Chennamangalam ( 2011 ) found that the observed luminosity functions of pulsars in globular clusters could be fit by either power -la w or lognormal functions and that the total number of pulsars in Terzan 5 strongly depended on the assumed model and its parameters (including the minimum pulsar pseudo-luminosity).

In a follow-up paper, J. Chennamangalam et al. ( 2013 ) fit a lognormal model to the extrapolated 1.4 GHz #ux densities of 25 pulsars and the t otal int egrat ed #ux from A. S. Fruchter & W. M. Goss ( 2000 ), alt ernativ ely assuming br oad or narr ow priors on the parameters of the lognormal distribution. For their 'broad' priors there is still support to the lower edge of the prior distribution in µ, which appears to truncate the posterior. In any case, for these assumptions they find a median and 95 per cent interval of N = 142 +310 -110 , only marginally consistent with our findings. While we are not entirely sure of the reason for this disagreement, it is at least partially due to their modelling of the total (diffuse + point source #ux). They assume a total 1.4 GHz #ux of 5.2 mJy, of which they attribute 1.9 mJy (37 per cent) to Ter5A, and overall that known pulsars sum to 4.2 mJy (81 per cent) of the total assumed #ux. Hence there is little #ux that needs to be attribut ed t o faint er, y et-undet ect ed pulsars. The value of 5.2 mJy appears to be an underestimate: as we discussed above in Section 4.2 , the diffuse 1.4 GHz #ux within 10 from A. S. Fruchter & W. M. Goss ( 2000 ) is ∼ 2 mJy. The comparison to the lower resolution imaging of A. S. Fruchter & W. M. Goss ( 1990 ), and our reimaging of this lower resolution data (Section 4.2 ), both show that these diffuse #ux measurements only pertain to the cluster core. As this radius contains only about half the known pulsars, the diffuse #ux int egrat ed ov er all radii should be double this value, ∼ 4 mJy. Adding in pulsars Ter5A and Ter5C, visible as point sources in the 1.4 GHz image, gives a total #ux of 6.7-7.8 mJy depending on whether the #ux densities from this image or A. R. Martsen et al. ( 2022 ) are assumed for the two bright pulsars. This would imply a much larger amount of 'residual' diffuse #ux t o be attribut ed t o faint pulsars, leading to a larger N. As our independent measurement of the diffuse #ux is consistent with that of A. S. Fruchter & W. M. Goss ( 2000 ) (see Section 4.2 ), we believe this larger value of the diffuse #ux is the correct one, implying a need for more pulsars than inferred by J. Chennamangalam et al. ( 2013 ).

J. Berteaud et al. ( 2024 ) compared Bayesian and simulationbased inference for inferring the pulsar population of a globular cluster, using Terzan 5 as a case study. While they used #ux densities for a larger sample of pulsars, they assumed the same diffuse #ux as in J. Chennamangalam et al. ( 2013 ), and found broadly similar results, with a total pulsar population of N = 158 +294 -104 . They emphasized that their results were extremely sensitive to the diffuse #ux and that new measurements of this quantity, as well as deeper surveys, were needed to improve constraints on the pulsar population. This is borne out by our results.

The specific number of pulsars inferred from our fits depends on the models assumed. It may be that the luminosity function is mor e comple x than a lognormal; for e xample, a br oken power law would add another free parameter, and of course yet more complex functions could be considered. But these are not yet demanded by the data.

Ov erall, w e conclude our #ux measurements of individual pulsars alone have not substantially improved our understanding of the luminosity function of pulsars in Terzan 5: the main advance is the simultaneous measurement of these with the residual diffuse #ux. For any of our model fits, the mean #ux is fainter than current detection limits, implying that the specific luminosity function parameters are not well-constrained, and the total number of pulsars is still very uncertain. Nevertheless the pulsar population appears much larger than inferred from previous work, with a 'minimum' reasonable estimate of pulsars of ∼ 250 and a true value that is likely higher.

We hav e present ed new deep high-resolution S -band (2-4 GHz) imaging of Terzan 5, detecting 38 of its 49 timed pulsars as continuum sources. Our central result is that ev en aft er subtracting these pulsars from the images, there is still a substantial diffuse residual #ux associated with a large population of pulsars below our detection limit. Our data alone are consistent with a pulsar population from ∼ 250 up to 1000 pulsars, though the values at the high end may be unlikely or impossible from a pulsar formation standpoint.

While with the present data the image simulations did not lead to meaningfully stronger conclusions about the pulsar population in Terzan 5, that should change with the ultra-deep, long baseline continuum imaging enabled by the next-generation VLA (E. J. Murphy et al. 2018 ) and SKA1-Mid (G. P . Swart, P . E. Dewdney & A. Cremonini 2022 ). We will e xplor e simulated ngVLA observations of Terzan 5 in a future work, also adding physical effects such as scintillation not yet included in our simulations.

How ev er, w e need not wait for next-generation facilities to make pr ogr ess: the continued detection of timed pulsars at the sensitivity limits of existing data (P. V. Padmanabh et al. 2024 ) already indicates the utility of ongoing searches for new pulsars in Terzan 5. Many pulsars just slightly fainter than those known are ripe for discovery.

We acknowledge the comments of an anonymous r efer ee, which helped improve the paper .

The VLA data are available via the NRAO data archive (Project code 22A-396) https://data.nrao.edu/portal/ .

Table A1. Catalogue of Terzan 5 sources. The columns are: (1) source ID number, (2-5) position and positional uncertainty, (6) alt ernativ e name for previously known sources [Ter5 CX1 is the transitional millisecond pulsar candidate (A. Bahramian et al. 2018 ) and EXO 1745-248 is a transient lowmass X-ray binary (A. J. Tetarenko et al. 2016 )], (7) distance from the cluster core, (8-9) source #ux densities or 3 σ upper limit for non-detection, (10) spectral index, (11) 2012 2.6 GHz #ux density from U20 , (12) 1-10 keV X-ray luminosity from A. Bahramian et al. ( 2020 ). # RA (ICRS) σ RA Dec. (ICRS) σ Dec . Alt ID Rad. S 2 . 5 GHz S 3 . 5 GHz α S 2 . 6 GHz L X (h:m:s) (arcsec) ( • : : ) (arcsec) (arcmin) ( µ Jy) ( µJy) ( µJy) ( 10 31 erg s -1 ) 1 17:48:02.247 0.06 -24:46:37.82 0.11 Ter5A 0.60 714 ±3 234 ±2 -3 . 31 +0 . 03 -0 . 03 698 ±7 0 . 05 +1 . 29 -0 . 05 2 17:47:50.917 0.06 -24:44:51.39 0.11 -3.68 363 ±3 299 ±2 -0 . 58 +0 . 03 -0 . 03 327 ±7 -3 17:48:06.404 0.06 -24:45:04.95 0.11 -1.70 243 ±3 200 ±3 -0 . 59 +0 . 05 -0 . 06 287 ±6 -4 17:48:10.921 0.06 -24:45:40.13 0.11 -1.75 146 ±3 145 ±3 -0 . 03 +0 . 09 -0 . 09 192 ±6 -5 17:48:05.041 0.06 -24:46:41.30 0.11 Ter5P 0.07 51 ±4 29 ±3 -1 . 7 +0 . 3 -0 . 3 224 ±7 43 . 8 +1 . 7 -1 . 6 6 17:47:55.637 0.06 -24:44:58.97 0.11 -2.73 126 ±3 95 ±2 -0 . 8 +0 . 1 -0 . 1 173 ±6 -7 17:48:04.534 0.06 -24:46:34.83 0.11 Ter5C 0.18 319 ±4 116 ±3 -3 . 01 +0 . 07 -0 . 08 275 ±7 -8 17:48:09.899 0.06 -24:48:57.66 0.11 -2.50 94 ±3 87 ±3 -0 . 2 +0 . 1 -0 . 1 108 ±6 -9 17:48:05.106 0.06 -24:46:34.51 0.11 Ter5V 0.18 29 ±4 21 ±3 -1 . 0 +0 . 6 -0 . 5 93 ±7 2 . 6 +3 . 1 -1 . 6 10 17:48:16.813 0.06 -24:46:41.59 0.11 -2.72 82 ±3 71 ±4 -0 . 4 +0 . 2 -0 . 2 60 ±6 -11 17:47:53.154 0.06 -24:46:52.01 0.11 -2.66 67 ±3 57 ±2 -0 . 5 +0 . 2 -0 . 2 71 ±6 -12 17:47:57.125 0.06 -24:45:21.08 0.11 -2.24 67 ±3 42 ±2 -1 . 4 +0 . 2 -0 . 2 76 ±6 -13 17:48:08.874 0.06 -24:47:34.14 0.11 -1.23 40 ±3 32 ±3 -0 . 6 +0 . 3 -0 . 3 54 ±6 -14 17:48:01.734 0.06 -24:46:28.13 0.11 -0.76 58 ±3 53 ±2 -0 . 2 +0 . 2 -0 . 2 41 ±7 -15 17:48:08.064 0.06 -24:46:00.98 0.11 -1.03 43 ±3 31 ±2 -1 . 0 +0 . 3 -0 . 3 21 ±7 -16 17:47:49.281 0.06 -24:46:19.72 0.11 -3.56 44 ±3 37 ±2 -0 . 5 +0 . 3 -0 . 3 38 ±6 -17 17:48:04.623 0.06 -24:46:40.80 0.11 Ter5M 0.08 31 ±4 10 ±3 -3 . 0 +0 . 5 -0 . 4 58 ±7 19 . 7 +2 . 9 -2 . 7 18 17:48:04.956 0.06 -24:46:45.66 0.11 Ter5ae/Z 0.03 56 ±4 24 ±3 -2 . 5 +0 . 4 -0 . 4 66 ±7 1 . 9 +2 . 7 -1 . 0 19 17:47:52.460 0.06 -24:48:58.42 0.11 -3.59 36 ±3 28 ±2 -0 . 7 +0 . 4 -0 . 4 38 ±7 -20 17:47:50.438 0.06 -24:46:34.73 0.11 -3.28 45 ±3 39 ±2 -0 . 4 +0 . 3 -0 . 2 39 ±6 -21 17:48:17.922 0.16 -24:46:31.894 0.14 -2.97 < 11 < 11 -26 ±7 -22 17:48:05.100 0.06 -24:46:44.40 0.11 Ter5Y 0.06 35 ±4 13 ±3 -2 . 7 +0 . 6 -0 . 5 46 ±7 2 . 3 +0 . 5 -0 . 4 23 17:48:08.566 0.08 -24:44:15.75 0.14 -2.62 28 ±3 < 8 ---24 17:48:14.271 0.08 -24:46:50.93 0.14 -2.14 28 ±3 < 9 -38 ±6 -25 17:48:15.084 0.16 -24:47:47.332 0.14 -2.55 < 10 < 10 -35 ±6 -26 17:48:04.869 0.06 -24:46:46.55 0.11 Ter5I 0.03 45 ±4 20 ±3 -2 . 5 +0 . 4 -0 . 4 --27 17:48:04.915 0.06 -24:46:53.87 0.11 Ter5N 0.16 36 ±4 24 ±3 -1 . 2 +0 . 4 -0 . 4 62 ±7 0 . 5 +2 . 1 -0 . 5 28 17:48:03.407 0.06 -24:46:35.57 0.11 Ter5E 0.36 28 ±3 19 ±2 -1 . 0 +0 . 5 -0 . 5 46 ±7 0 . 02 +0 . 13 -0 . 01 29 17:48:05.116 0.06 -24:46:38.13 0.11 Ter5F 0.12 31 ±4 20 ±3 -1 . 3 +0 . 5 -0 . 5 41 ±7 0 . 02 +0 . 17 -0 . 02 30 17:48:04.677 0.06 -24:46:51.45 0.11 Ter5O 0.12 43 ±4 13 ±3 -3 . 1 +0 . 4 -0 . 3 53 ±7 2 . 3 +0 . 4 -0 . 4 31 17:48:05.023 0.06 -24:46:43.52 0.11 -0.04 27 ±4 17 ±3 -1 . 4 +0 . 6 -0 . 7 -2 . 1 +0 . 5 -0 . 5 32 17:48:04.734 0.06 -24:46:35.90 0.11 Ter5L 0.15 27 ±4 15 ±3 -1 . 6 +0 . 7 -0 . 7 30 ±7 0 . 03 +0 . 19 -0 . 02 33 17:48:03.904 0.06 -24:46:47.69 0.11 Ter5K 0.22 20 ±4 9 ±3 -2 . 2 +0 . 9 -0 . 8 24 ±7 0 . 02 +0 . 19 -0 . 02 34 17:48:04.677 0.16 -24:46:48.743 0.14 -0.08 < 11 < 8 --8 . 0 +0 . 8 -0 . 8 35 17:48:04.572 0.06 -24:46:42.22 0.11 CX1 0.07 30 ±4 13 ±3 -2 . 5 +0 . 6 -0 . 6 50 ±7 43 . 9 +1 . 9 -1 . 6 36 17:48:04.754 0.06 -24:46:43.06 0.11 Ter5ab 0.03 28 ±4 10 ±3 -2 . 7 +0 . 7 -0 . 5 53 ±7 10 . 2 +1 . 7 -1 . 8 37 17:48:10.329 0.06 -24:47:59.08 0.11 -1.76 35 ±3 26 ±3 -0 . 9 +0 . 4 -0 . 4 43 ±6 -38 17:48:04.618 0.06 -24:46:45.71 0.11 Ter5ar 0.06 30 ±4 18 ±3 -1 . 6 +0 . 6 -0 . 6 38 ±7 25 . 5 +1 . 3 -1 . 3 39 17:48:05.924 0.06 -24:46:05.71 0.11 Ter5D 0.69 28 ±3 14 ±2 -2 . 0 +0 . 6 -0 . 6 37 ±7 0 . 01 +0 . 07 -0 . 01 40 17:48:04.418 0.16 -24:46:48.780 0.14 -0.12 < 11 < 8 -36 ±7 1 . 5 +0 . 9 -0 . 4 41 17:48:03.726 0.06 -24:47:23.45 0.11 -0.70 12 ±3 9 ±2 -0 . 6 +1 . 1 -1 . 2 31 ±7 -42 17:48:05.225 0.08 -24:46:47.661 0.14 EXO1745 0.10 < 11 < 8 --179 +3 -3 43 17:48:04.835 0.08 -24:46:42.17 0.14 Ter5W 0.04 28 ±4 < 8 ---44 17:48:12.548 0.08 -24:43:55.59 0.14 -3.31 28 ±3 < 10 ---45 17:47:54.477 0.06 -24:45:45.17 0.11 -2.55 27 ±3 23 ±2 -0 . 4 +0 . 4 -0 . 4 --46 17:48:04.335 0.06 -24:47:05.14 0.11 Ter5Q 0.36 26 ±3 13 ±2 -2 . 0 +0 . 7 -0 . 7 -0 . 06 +0 . 16