<figure><img src="image_1.png" /></figure>

HAL
open science

# Characterizing the Gamma-ray Emission Properties of the Globular Cluster M5 with the Fermi-LAT

X Hou, W Zhang, P.C.C Freire, D.F Torres, J Ballet, D.A Smith, T.J

Johnson, M Kerr, C.C Cheung, L Guillemot, et al.

► To cite this version:

X Hou, W Zhang, P.C.C Freire, D.F Torres, J Ballet, et al. Characterizing the Gamma-ray Emission Properties of the Globular Cluster M5 with the Fermi-LAT. Astrophys.J., 2024, 964 (2), pp.118. 10.3847/1538-4357/ad3210. hal-04446870

HAL Id: hal-04446870

https://hal.science/hal-04446870

Submitted on 16 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Characterizing the Gamma-ray Emission Properties of the Globular Cluster M5 with the Fermi-LAT

X Hou, W Zhang, P.C.C Freire, D.F Torres, J Ballet, D.A Smith, T.J

Johnson, M Kerr, C.C Cheung, Lucas Guillemot, et al.

► To cite this version:

X Hou, W Zhang, P.C.C Freire, D.F Torres, J Ballet, et al. Characterizing the Gamma-ray Emission Properties of the Globular Cluster M5 with the Fermi-LAT. The Astrophysical Journal, 2024, 964 (2), pp.118. 10.3847/1538-4357/ad3210. hal-04446870

HAL Id: hal-04446870

https://hal.science/hal-04446870v1

Submitted on 16 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Characterizing the Gamma-Ray Emission Properties of the Globular Cluster M5 with the Fermi-LAT

X. Hou<sup>1,2,3</sup>, W. Zhang<sup>4,5</sup>, P. C. C. Freire<sup>6</sup>, D. F. Torres<sup>4,5,7</sup>, J. Ballet<sup>8</sup>, D. A. Smith<sup>9</sup>, T. J. Johnson<sup>10</sup>, M. Kerr<sup>11</sup>, C. C. Cheung<sup>11</sup>, L. Guillemot<sup>12,13</sup>, J. Li<sup>14,15</sup>, L. Zhang<sup>16,17</sup>, A. Ridolfi<sup>6,18</sup>, P. Wang<sup>16,19</sup>, D. Li<sup>16,20,21</sup>, J. Yuan<sup>22</sup>, and N. Wang<sup>22</sup>

<sup>1</sup> Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216, People's Republic of China; xhou@ynao.ac.cn

<sup>2</sup> Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming 650216, People's Republic of China

<sup>3</sup> Center for Astronomical Mega-Science, Chinese Academy of Sciences, Beijing 100012, People's Republic of China

<sup>4</sup> Institute of Space Sciences (ICE, CSIC), Campus UAB, 08193 Barcelona, Spain; dtorres@ice.csic.es

<sup>5</sup> Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

<sup>6</sup> Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>7</sup> Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain

<sup>8</sup> Université Paris Saclay and Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France

<sup>9</sup> Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France

<sup>10</sup> College of Science, George Mason University, Fairfax, VA 22030, resident at Naval Research Laboratory, Washington, DC 20375, USA

<sup>11</sup> Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

<sup>12</sup> Laboratoire de Physique et Chimie de l'Environnement et de l'Espace—Université d'Orléans/CNRS, F-45071 Orléans Cedex 02, France

<sup>13</sup> Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France

<sup>14</sup> CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei 230026, People's Republic of China

<sup>15</sup> School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, People's Republic of China

<sup>16</sup> National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, People's Republic of China

<sup>17</sup> Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia

<sup>18</sup> INAF—Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius (CA), Italy

<sup>19</sup> Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal University, Beijing 102206, People's Republic of China

<sup>20</sup> University of Chinese Academy of Sciences, Chinese Academy of Sciences, Beijing 100101, People's Republic of China

<sup>21</sup> Zhijiang Lab, Hangzhou 311121, People's Republic of China

<sup>22</sup> Xinjiang Astronomical Observatory, Chinese Academy of Sciences, Urumqi, Xinjiang 830011, People's Republic of China

Received 2024 January 15; revised 2024 March 7; accepted 2024 March 7; published 2024 March 22

## Abstract

We analyzed the globular cluster M5 (NGC 5904) using 15 yr of gamma-ray data from the Fermi Large Area Telescope (LAT). Using rotation ephemerides generated from Arecibo and FAST radio telescope observations, we searched for gamma-ray pulsations from the seven millisecond pulsars (MSPs) identified in M5. We detected no significant pulsations from any of the individual pulsars. In addition, we searched for possible variations of the gamma-ray emission as a function of orbital phase for all six MSPs in binary systems, but we did not detect any significant modulations. The gamma-ray emission from the direction of M5 is well described by an exponentially cutoff power-law spectral model, although other models cannot be excluded. The phase-averaged emission is consistent with being steady on a timescale of a few months. We estimate the number of MSPs in M5 to be between 1 and 10, using the gamma-ray conversion efficiencies for well-characterized gamma-ray MSPs in the Third Fermi-LAT Catalog of Gamma-ray Pulsars, suggesting that the sample of known MSPs in M5 is (nearly) complete, even if it is not currently possible to rule out a diffuse component of the observed gamma rays from the cluster.

Unified Astronomy Thesaurus concepts: Globular star clusters (656); Millisecond pulsars (1062); Gamma-ray sources (633)

## 1. Introduction

Globular clusters (GCs) are the oldest and densest stellar systems bound by gravity. Due to the high stellar density (>1000 pc⁻³; Sollima & Baumgardt 2017) and frequent dynamical interactions between stars in GCs, the formation rate per unit mass of low-mass X-ray binaries (LMXBs) is orders of magnitude higher in GCs than in the Galactic field (Clark 1975; Katz 1975). LMXBs are more abundant in GCs as a natural consequence of such a dynamical formation scenario, and a linear correlation between the number of LMXBs in GCs

and the stellar encounter rate Γ_c has been expected and confirmed by observations of GCs (Gendre et al. 2003; Pooley et al. 2003; de Menezes et al. 2023).

Millisecond pulsars (MSPs; usually defined as those having spin period *P* ≤ 30 ms) are generally believed to be descendants of LMXBs (e.g., Alpar et al. 1982; Bhattacharya & van den Heuvel 1991). Of the 305 pulsars detected in radio in 40 GCs in the Milky Way (MW) halo,<fnref n="23" /> 80% are MSPs. In comparison, 427 known MSPs are not associated with GCs,<fnref n="24" /> only 10% of the known Galactic pulsar population. A positive correlation between the MSP population in GCs and Γ_c has also been reported (e.g., Hui et al. 2010; Bahramian et al. 2013),

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="23"><sup>23</sup> https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html</footnote>

<footnote n="24"><sup>24</sup> http://astro.phys.wvu.edu/GalacticMSPs/GalacticMSPs.txt</footnote>
<figure><img src="image_2.png" /></figure>

# Characterizing the Gamma-Ray Emission Properties of the Globular Cluster M5 with the Fermi-LAT

X. Hou<sup>1,2,3</sup>, W. Zhang<sup>4,5</sup>, P. C. C. Freire<sup>6</sup>, D. F. Torres<sup>4,5,7</sup>, J. Ballet<sup>8</sup>, D. A. Smith<sup>9</sup>, T. J. Johnson<sup>10</sup>, M. Kerr<sup>11</sup>, C. C. Cheung<sup>11</sup>, L. Guillemot<sup>12,13</sup>, J. Li<sup>14,15</sup>, L. Zhang<sup>16,17</sup>, A. Ridolfi<sup>6,18</sup>, P. Wang<sup>16,19</sup>, D. Li<sup>16,20,21</sup>, J. Yuan<sup>22</sup>, and N. Wang<sup>22</sup>

<sup>1</sup> Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216, People's Republic of China; xhou@ynao.ac.cn

<sup>2</sup> Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming 650216, People's Republic of China

<sup>3</sup> Center for Astronomical Mega-Science, Chinese Academy of Sciences, Beijing 100012, People's Republic of China

<sup>4</sup> Institute of Space Sciences (ICE, CSIC), Campus UAB, 08193 Barcelona, Spain; dtorres@ice.csic.es

<sup>5</sup> Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

<sup>6</sup> Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>7</sup> Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain

<sup>8</sup> Université Paris Saclay and Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France

<sup>9</sup> Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France

<sup>10</sup> College of Science, George Mason University, Fairfax, VA 22030, resident at Naval Research Laboratory, Washington, DC 20375, USA

<sup>11</sup> Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

<sup>12</sup> Laboratoire de Physique et Chimie de l'Environnement et de l'Espace—Université d'Orléans/CNRS, F-45071 Orléans Cedex 02, France

<sup>13</sup> Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France

<sup>14</sup> CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei 230026, People's Republic of China

<sup>15</sup> School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, People's Republic of China

<sup>16</sup> National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, People's Republic of China

<sup>17</sup> Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia

<sup>18</sup> INAF—Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius (CA), Italy

<sup>19</sup> Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal University, Beijing 102206, People's Republic of China

<sup>20</sup> University of Chinese Academy of Sciences, Chinese Academy of Sciences, Beijing 100101, People's Republic of China

<sup>21</sup> Zhijiang Lab, Hangzhou 311121, People's Republic of China

<sup>22</sup> Xinjiang Astronomical Observatory, Chinese Academy of Sciences, Urumqi, Xinjiang 830011, People's Republic of China

Received 2024 January 15; revised 2024 March 7; accepted 2024 March 7; published 2024 March 22

## Abstract

We analyzed the globular cluster M5 (NGC 5904) using 15 yr of gamma-ray data from the Fermi Large Area Telescope (LAT). Using rotation ephemerides generated from Arecibo and FAST radio telescope observations, we searched for gamma-ray pulsations from the seven millisecond pulsars (MSPs) identified in M5. We detected no significant pulsations from any of the individual pulsars. In addition, we searched for possible variations of the gamma-ray emission as a function of orbital phase for all six MSPs in binary systems, but we did not detect any significant modulations. The gamma-ray emission from the direction of M5 is well described by an exponentially cutoff power-law spectral model, although other models cannot be excluded. The phase-averaged emission is consistent with being steady on a timescale of a few months. We estimate the number of MSPs in M5 to be between 1 and 10, using the gamma-ray conversion efficiencies for well-characterized gamma-ray MSPs in the Third Fermi-LAT Catalog of Gamma-ray Pulsars, suggesting that the sample of known MSPs in M5 is (nearly) complete, even if it is not currently possible to rule out a diffuse component of the observed gamma rays from the cluster.

Unified Astronomy Thesaurus concepts: Globular star clusters (656); Millisecond pulsars (1062); Gamma-ray sources (633)

## 1. Introduction

Globular clusters (GCs) are the oldest and densest stellar systems bound by gravity. Due to the high stellar density (>1000 pc⁻³; Sollima & Baumgardt 2017) and frequent dynamical interactions between stars in GCs, the formation rate per unit mass of low-mass X-ray binaries (LMXBs) is orders of magnitude higher in GCs than in the Galactic field (Clark 1975; Katz 1975). LMXBs are more abundant in GCs as a natural consequence of such a dynamical formation scenario, and a linear correlation between the number of LMXBs in GCs

and the stellar encounter rate Γ_c has been expected and confirmed by observations of GCs (Gendre et al. 2003; Pooley et al. 2003; de Menezes et al. 2023).

Millisecond pulsars (MSPs; usually defined as those having spin period *P* ≤ 30 ms) are generally believed to be descendants of LMXBs (e.g., Alpar et al. 1982; Bhattacharya & van den Heuvel 1991). Of the 305 pulsars detected in radio in 40 GCs in the Milky Way (MW) halo,<fnref n="23" /> 80% are MSPs. In comparison, 427 known MSPs are not associated with GCs,<fnref n="24" /> only 10% of the known Galactic pulsar population. A positive correlation between the MSP population in GCs and Γ_c has also been reported (e.g., Hui et al. 2010; Bahramian et al. 2013),

<figure><img src="image_3.png" /></figure>

<footnote n="23"><sup>23</sup> https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html</footnote>

<footnote n="24"><sup>24</sup> http://astro.phys.wvu.edu/GalacticMSPs/GalacticMSPs.txt</footnote>
which provided evidence for the dynamical origin of MSPs as had long been predicted, given the close relation between MSPs and LMXBs.

GCs have been established as a class of gamma-ray emitters using data from the Large Area Telescope (LAT; Atwood et al. 2009) on board the Fermi Gamma-ray Space Telescope<fnref n="25" /> launched in 2008. Up to now, gamma rays coincident with the directions of about 39 GCs have been reported (Abdo et al. 2010; Kong et al. 2010; Tam et al. 2011; Zhou et al. 2015; Zhang et al. 2016, 2022; Lloyd et al. 2018; de Menezes et al. 2019; Abdollahi et al. 2020; Yuan et al. 2022a, 2022b; de Menezes et al. 2023), all of which are listed in the Harris (1996) catalog<fnref n="26" /> of MW GCs. As a main class of LAT gamma-ray sources, MSPs are reasonably thought to be responsible for the collective gamma-ray emission from GCs, as first predicted by Chen (1991) and suggested by the spectral similarities between the observed gamma-ray MSPs and GCs. Indeed, evidence of correlation between the gamma-ray luminosity $L_γ$ and $\Gamma_c$ established by various studies has provided support to the dynamical formation of MSPs and the MSP origin of gamma rays in GCs (Abdo et al. 2010; Hui et al. 2011; Bahramian et al. 2013; Hooper & Linden 2016; Zhang et al. 2016; de Menezes et al. 2019, 2023; Feng et al. 2024). Recent analyses show that all GCs are point-like sources in gamma rays, implying that MSPs are mostly concentrated in their cores (de Menezes et al. 2019, 2023).

Exceptionally, gamma-ray pulsations from individual MSPs have also been reported in three GCs: PSR J1823−3021A in NGC 6624 (Freire et al. 2011), PSR B1821−24 in NGC 6626 (M28, Johnson et al. 2013; Wu et al. 2013), and PSR J1835−3259B in NGC 6652 (Zhang et al. 2022). The first two are isolated MSPs, while the third one is in a binary system (Gautam et al. 2022). They are all energetic pulsars with relatively high spin-down power and are bright in gamma rays, indicating that they probably dominate the gamma-ray emission from the host GCs. We have been unable to confirm the gamma-ray pulsations from PSR J1717+4308A in NGC 6341 (M92) reported by Zhang et al. (2023b).

M5 (NGC 5904) is a bright GC (visual magnitude $V \approx 5.6$) with a small $\Gamma_c$, which is also small when we divide it by the number of stars in the cluster, i.e., the encounter rate per formed binary, $\gamma_b$ (Verbunt & Freire 2014). The cluster is at a distance of $7.48 \pm 0.60$ kpc with a half-mass radius of 5.6 pc.<fnref n="27" /> Radio observations using Arecibo and FAST have resulted in the detection of seven MSPs in M5 (Anderson et al. 1997; Mott & Freire 2003; Hessels et al. 2007; Pan et al. 2021; Zhang et al. 2023a). The first discoveries were then known as B1516+02A and B1516+02B; they are now known as PSRs J1518+0204A and J1518+0204B. The following discoveries are known as J1518+0204C–J1518+0204G, and we will refer to these pulsars as M5A, M5B, M5C, M5D, M5E, M5F, and M5G, respectively. Improved timing solutions for all these pulsars have been derived based on observations using the two telescopes (Zhang et al. 2023a).

The basic properties of the seven MSPs in M5 are presented in Table 1. M5A is an isolated MSP, while the other six are in binary systems. Based on optical observations, the companions of M5D, M5E, and M5F are very likely low-mass He white dwarfs (WDs). M5C, M5D, M5E, and M5G have X-ray

counterparts detected by Chandra, from which thermal X-ray emission is observed. M5C (eclipsing) and M5G (noneclipsing) are black widow (BW) systems that show no or little nonthermal X-ray emissions, indicating that the intrabinary shock produces weak synchrotron radiation. The improved measurement of the periastron advance rate $\dot{\omega} = 0.01361 \text{ yr}^{-1}$ of M5B is compatible with a heavy neutron star (NS), consistent with previous studies. Although its inclination is not well constrained, M5B is probably not edge-on (Zhang et al. 2023a).

M5 was initially proposed as a gamma-ray emitter by Zhou et al. (2015) with a marginal detection ($3.2\sigma$) using 6 yr of LAT Pass 7 Reprocessed data and later confirmed by Zhang et al. (2016) with $4.4\sigma$ using 7 yr of Pass 8 data (Atwood et al. 2013; Bruel et al. 2018). The Fermi-LAT Fourth Source Catalog (4FGL; Abdollahi et al. 2020) first associated M5 with the LAT source 4FGL J1518.8+0203 ($6.7\sigma$) using 8 yr of Pass 8 data.

In this paper, we performed Fermi-LAT analysis of M5 to investigate its gamma-ray emission properties using an updated data set and the most recent LAT catalog. From Table 1, the angular separations between the seven MSPs are ~0.01–0.02; thus, the LAT angular resolution (~0.1°; Figure 1 in 4FGL) is insufficient to spatially separate the individual MSPs. Motivated by this, we have performed timing analysis in addition to spectral analysis taking advantage of the most updated timing solutions for the seven MSPs in M5. Either detection or nondetection of individual pulsars usefully constrains the number of MSPs in M5 and the GC gamma-ray emission models and serves as a useful test of the dynamical formation scenario of GC MSPs. We describe the data set and analysis methods in Section 2, and then we present the analysis results in Section 3. Finally, we discuss the implications of the results and conclude in Section 4.

## 2. Data Set and Analysis Methods

According to the latest 4FGL-DR4 source list (glc_psc_v32. fit;<fnref n="28" /> Ballet et al. 2023), M5 is associated with 4FGL J1518.8+0203, which is 4.1 away from the M5 center and has a gamma-ray significance of $6.8\sigma$. We used 15 yr (2008 August 4–2023 August 4) of Pass 8 data from the Fermi-LAT within 10° around the M5 center ($\alpha, \delta$) = (229.6384, 2.0810) in the J2000 frame.<fnref n="29" /> SOURCE class events in the energy range of 0.1–500 GeV have been selected and the standard event filter “DATA QUAL > 0 & LAT CONFIG==1” has been applied to get data of good quality.<fnref n="30" /> To avoid contamination from solar flares and gamma-ray bursts (GRBs), we have excluded time intervals when solar flares and GRBs occurred (Abdollahi et al. 2022; Ballet et al. 2023).

To reduce the contamination from the low-energy Earth limb emission, we followed similar point-spread function (PSF) and zenith-angle cuts adopted in the LAT 4FGL catalog (Abdollahi et al. 2020). Specifically, in the 0.1–0.3 GeV band, only PSF2 and PSF3 events were selected with zenith angles <90°; in the 0.3–1 GeV band, PSF1, PSF2, and PSF3 events were selected with zenith angles <100°; and in the 1–500 GeV band, all events (PSF0, PSF1, PSF2, and PSF3) with zenith angles

<footnote n="25">²⁵ http://fermi.gsfc.nasa.gov/ssc</footnote>

<footnote n="26">²⁶ 2010 edition: https://physics.mcmaster.ca/~harris/mwgc.dat.</footnote>

<footnote n="27">²⁷ https://people.smp.uq.edu.au/HolgerBaumgardt/globular/</footnote>

<footnote n="28">²⁸ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/14yropr 書</footnote>

<footnote n="29">²⁹ http://simbad.u-strasbg.fr/simbad/</footnote>

<footnote n="30">³⁰ https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data_preparation.html</footnote>
which provided evidence for the dynamical origin of MSPs as had long been predicted, given the close relation between MSPs and LMXBs.

GCs have been established as a class of gamma-ray emitters using data from the Large Area Telescope (LAT; Atwood et al. 2009) on board the Fermi Gamma-ray Space Telescope<fnref n="25" /> launched in 2008. Up to now, gamma rays coincident with the directions of about 39 GCs have been reported (Abdo et al. 2010; Kong et al. 2010; Tam et al. 2011; Zhou et al. 2015; Zhang et al. 2016, 2022; Lloyd et al. 2018; de Menezes et al. 2019; Abdollahi et al. 2020; Yuan et al. 2022a, 2022b; de Menezes et al. 2023), all of which are listed in the Harris (1996) catalog<fnref n="26" /> of MW GCs. As a main class of LAT gamma-ray sources, MSPs are reasonably thought to be responsible for the collective gamma-ray emission from GCs, as first predicted by Chen (1991) and suggested by the spectral similarities between the observed gamma-ray MSPs and GCs. Indeed, evidence of correlation between the gamma-ray luminosity $L_γ$ and $\Gamma_c$ established by various studies has provided support to the dynamical formation of MSPs and the MSP origin of gamma rays in GCs (Abdo et al. 2010; Hui et al. 2011; Bahramian et al. 2013; Hooper & Linden 2016; Zhang et al. 2016; de Menezes et al. 2019, 2023; Feng et al. 2024). Recent analyses show that all GCs are point-like sources in gamma rays, implying that MSPs are mostly concentrated in their cores (de Menezes et al. 2019, 2023).

Exceptionally, gamma-ray pulsations from individual MSPs have also been reported in three GCs: PSR J1823−3021A in NGC 6624 (Freire et al. 2011), PSR B1821−24 in NGC 6626 (M28, Johnson et al. 2013; Wu et al. 2013), and PSR J1835−3259B in NGC 6652 (Zhang et al. 2022). The first two are isolated MSPs, while the third one is in a binary system (Gautam et al. 2022). They are all energetic pulsars with relatively high spin-down power and are bright in gamma rays, indicating that they probably dominate the gamma-ray emission from the host GCs. We have been unable to confirm the gamma-ray pulsations from PSR J1717+4308A in NGC 6341 (M92) reported by Zhang et al. (2023b).

M5 (NGC 5904) is a bright GC (visual magnitude $V \approx 5.6$) with a small $\Gamma_c$, which is also small when we divide it by the number of stars in the cluster, i.e., the encounter rate per formed binary, $\gamma_b$ (Verbunt & Freire 2014). The cluster is at a distance of $7.48 \pm 0.60$ kpc with a half-mass radius of 5.6 pc.<fnref n="27" /> Radio observations using Arecibo and FAST have resulted in the detection of seven MSPs in M5 (Anderson et al. 1997; Mott & Freire 2003; Hessels et al. 2007; Pan et al. 2021; Zhang et al. 2023a). The first discoveries were then known as B1516+02A and B1516+02B; they are now known as PSRs J1518+0204A and J1518+0204B. The following discoveries are known as J1518+0204C–J1518+0204G, and we will refer to these pulsars as M5A, M5B, M5C, M5D, M5E, M5F, and M5G, respectively. Improved timing solutions for all these pulsars have been derived based on observations using the two telescopes (Zhang et al. 2023a).

The basic properties of the seven MSPs in M5 are presented in Table 1. M5A is an isolated MSP, while the other six are in binary systems. Based on optical observations, the companions of M5D, M5E, and M5F are very likely low-mass He white dwarfs (WDs). M5C, M5D, M5E, and M5G have X-ray

counterparts detected by Chandra, from which thermal X-ray emission is observed. M5C (eclipsing) and M5G (noneclipsing) are black widow (BW) systems that show no or little nonthermal X-ray emissions, indicating that the intrabinary shock produces weak synchrotron radiation. The improved measurement of the periastron advance rate $\dot{\omega} = 0^{\circ} 01361 \text{ yr}^{-1}$ of M5B is compatible with a heavy neutron star (NS), consistent with previous studies. Although its inclination is not well constrained, M5B is probably not edge-on (Zhang et al. 2023a).

M5 was initially proposed as a gamma-ray emitter by Zhou et al. (2015) with a marginal detection ($3.2\sigma$) using 6 yr of LAT Pass 7 Reprocessed data and later confirmed by Zhang et al. (2016) with $4.4\sigma$ using 7 yr of Pass 8 data (Atwood et al. 2013; Bruel et al. 2018). The Fermi-LAT Fourth Source Catalog (4FGL; Abdollahi et al. 2020) first associated M5 with the LAT source 4FGL J1518.8+0203 ($6.7\sigma$) using 8 yr of Pass 8 data.

In this paper, we performed Fermi-LAT analysis of M5 to investigate its gamma-ray emission properties using an updated data set and the most recent LAT catalog. From Table 1, the angular separations between the seven MSPs are $\sim 0^{\circ} 01-0^{\circ} 02$; thus, the LAT angular resolution ($\sim 0^{\circ} 1$; Figure 1 in 4FGL) is insufficient to spatially separate the individual MSPs. Motivated by this, we have performed timing analysis in addition to spectral analysis taking advantage of the most updated timing solutions for the seven MSPs in M5. Either detection or nondetection of individual pulsars usefully constrains the number of MSPs in M5 and the GC gamma-ray emission models and serves as a useful test of the dynamical formation scenario of GC MSPs. We describe the data set and analysis methods in Section 2, and then we present the analysis results in Section 3. Finally, we discuss the implications of the results and conclude in Section 4.

## 2. Data Set and Analysis Methods

According to the latest 4FGL-DR4 source list (glc_psc_v32. fit;<fnref n="28" /> Ballet et al. 2023), M5 is associated with 4FGL J1518.8+0203, which is 4.1 away from the M5 center and has a gamma-ray significance of $6.8\sigma$. We used 15 yr (2008 August 4–2023 August 4) of Pass 8 data from the Fermi-LAT within $10^{\circ}$ around the M5 center ($\alpha, \delta) = (229^{\circ} 6384, 2^{\circ} 0810$) in the J2000 frame.<fnref n="29" /> SOURCE class events in the energy range of 0.1–500 GeV have been selected and the standard event filter “DATA QUAL > 0 & LAT CONFIG == 1” has been applied to get data of good quality.<fnref n="30" /> To avoid contamination from solar flares and gamma-ray bursts (GRBs), we have excluded time intervals when solar flares and GRBs occurred (Abdollahi et al. 2022; Ballet et al. 2023).

To reduce the contamination from the low-energy Earth limb emission, we followed similar point-spread function (PSF) and zenith-angle cuts adopted in the LAT 4FGL catalog (Abdollahi et al. 2020). Specifically, in the $0.1–0.3$ GeV band, only PSF2 and PSF3 events were selected with zenith angles $<90^{\circ}$; in the $0.3–1$ GeV band, PSF1, PSF2, and PSF3 events were selected with zenith angles $<100^{\circ}$; and in the $1–500$ GeV band, all events (PSF0, PSF1, PSF2, and PSF3) with zenith angles

<footnote n="25">²⁵ http://fermi.gsfc.nasa.gov/ssc</footnote>

<footnote n="26">²⁶ 2010 edition: https://physics.mcmaster.ca/~harris/mwgc.dat.</footnote>

<footnote n="27">²⁷ https://people.smp.uq.edu.au/HolgerBaumgardt/globular/</footnote>

<footnote n="28">²⁸ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/14yropr 書</footnote>

<footnote n="29">²⁹ http://simbad.u-strasbg.fr/simbad/</footnote>

<footnote n="30">³⁰ https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data_preparation.html</footnote>
<table><caption>**Table 1**
Basic Properties of the Seven MSPs in M5</caption><thead><tr><th>Name</th><th>R.A.<br>(deg)</th><th>Decl.<br>(deg)</th><th>P<br>(ms)</th><th>P<sub>orb</sub><br>(day)</th><th>Eccent.<sup>a</sup></th><th>E&#x0307;<sup>b</sup><br>(10<sup>34</sup> erg s<sup>-1</sup>)</th><th>X-Ray<br>Counterpart</th><th>Optical<br>Counterpart</th><th>Comment</th></tr></thead><tbody><tr><td>M5A</td><td>229.6388</td><td>2.0910</td><td>5.55</td><td>...</td><td>...</td><td>&lt;1.56</td><td>No</td><td>No</td><td>Isolated</td></tr><tr><td>M5B</td><td>229.6311</td><td>2.0876</td><td>7.95</td><td>6.858</td><td>0.138</td><td>&lt;0.25</td><td>No</td><td>No</td><td>Heavy, not likely edge-on</td></tr><tr><td>M5C</td><td>229.6366</td><td>2.0799</td><td>2.48</td><td>0.687</td><td>0</td><td>&lt;12.36</td><td>Yes</td><td>Yes</td><td>Eclipsing BW</td></tr><tr><td>M5D</td><td>229.6268</td><td>2.0833</td><td>2.99</td><td>1.222</td><td>lower than M5B</td><td>0.28–4.98</td><td>Yes</td><td>Yes</td><td>He WD companion</td></tr><tr><td>M5E</td><td>229.6388</td><td>2.0772</td><td>3.18</td><td>1.097</td><td>lower than M5B</td><td>&lt;5.27</td><td>Yes</td><td>Yes</td><td>He WD companion</td></tr><tr><td>M5F</td><td>229.6350</td><td>2.0867</td><td>2.65</td><td>1.610</td><td>lower than M5B</td><td>&lt;8.30</td><td>No</td><td>Yes</td><td>He WD companion</td></tr><tr><td>M5G</td><td>229.6197</td><td>2.0875</td><td>2.75</td><td>0.114</td><td>0</td><td>0.04–3.34</td><td>Yes</td><td>No</td><td>Noneclipsing BW</td></tr></tbody></table>

**Notes.**

<105° were used. This reduces the contribution of the Earth limb contamination to the total background to less than 10%.

We built a spatial-spectral model for M5 by including the 4FGL-DR4 sources within 20° around the M5 center. The Galactic interstellar emission model (“gll_iem_v07.fits”) and the isotropic emission spectrum (“iso_P8R3_SOUR-CE_V3_v1.txt”), which takes into account the extragalactic emission and the residual instrumental background,<fnref n="31" /> were also included. Energy dispersion has been taken into account by adding two extra energy bins except for the isotropic component. We have performed a summed likelihood analysis in a 14° × 14° region of interest (ROI) around M5. The significance of a given source in the model is characterized by the test statistic (TS).<fnref n="32" /> In this work, we used the fermipy<fnref n="33" /> package (v1.2.0; Wood et al. 2017), in which Fermitools<fnref n="34" /> (v2.2.0) is integrated.

## 3. Analysis Results

### 3.1. Localization

Since the LAT cannot resolve M5, we modeled it as a point source in our analysis. We first localized it in the higher energy range of 1–500 GeV to take advantage of LAT’s better spatial resolution. Data were binned using a 0.05 × 0.05 pixel size and 10 logarithmic energy bins per decade. The ROI was optimized by fitting all sources in the ROI to ensure that all parameters are close to their global likelihood maxima. We tested three spectral models for M5 during the localization: the default LogParabola (LP) model in 4FGL-DR4,

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-(\alpha + \beta \ln(E/E_0))}, \quad (1) $$

the simple power-law (PL) model,

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\Gamma_0}, \quad (2) $$

and the PLSuperExpCutoff4 (PLEC4) model that is used for pulsar emission in the Third Fermi-LAT Catalog of Gamma-ray Pulsars (3PC; Smith et al. 2023),

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\Gamma+d/b} \exp \left[ \frac{d}{b^2} \left( 1 - \left( \frac{E}{E_0} \right)^b \right) \right], \quad (3) $$

where $N_0$ is the normalization; $\alpha$ and $\Gamma_0$ are the spectral indices in the LP and PL models, respectively; $\beta$ is the curvature index; $E_0$ is the reference energy and is fixed to the catalog value for the LP model and to 1 GeV for the PL and PLEC4 models; $\Gamma$ is the local spectral index at $E_0$ in the PLEC4 model; $d$ is the local curvature at $E_0$; and $b$ is fixed to 2/3.

Only the normalizations of the Galactic/isotropic diffuse components were set free to vary during the localization. We note that changing the spectral model and leaving the position free simultaneously in each fit will lead to “un-nested models.” However, the localizations determined from the three models are completely consistent (Table 2). Both PLEC4 and LP models have the same number of degrees of freedom and fit the data quite well. Since the TS values obtained with both models are basically the same, we chose to continue our analysis with the PLEC4 model simply because it is physically motivated as a superposition of curvature radiation spectra for a range of electron energies (Bednarek & Sitarek 2007; Venter & de Jager 2008; Venter et al. 2009; Cheng et al. 2010).

Figure 1 (left panel) shows the residual map for the full ROI for the localization generated by removing the contribution of all 4FGL sources in the ROI. No significant residuals (>4σ) were found, indicating a good modeling of the ROI. The TS excess map (Figure 2, left panel) presents the localization result in a 1° × 1° region centered at M5 by removing the contribution of all 4FGL sources except M5, i.e., all sources except M5 are included in the model.

### 3.2. Spectral Fit

We performed a broadband spectral fit for M5 in the energy range of 0.1–500 GeV using the best-fit localization obtained with the PLEC4 model. Data were binned using a 0.1° × 0.1° pixel size and 10 logarithmic energy bins per decade. This time, all the spectral parameters of sources within 5° and those of the Galactic/isotropic diffuse components were set free to vary, along with the normalizations of significantly variable sources within 5°–8° around M5.

<footnote><sup>a</sup> The eccentricities of M5C and M5G have been set to 0 based on the assumption that the orbits of BW pulsars are circular owing to tidal dissipation (see Table 1 in Zhang et al. 2023a).</footnote>

<footnote><sup>b</sup> Spin-down powers $\dot{E}$ have been calculated based on the intrinsic spin-down rate $\dot{P}_{\text{int}}$ upper limit corrected for accelerations caused by the gravitational field of the GC for the line of sight of each pulsar (see Table 2 in Zhang et al. 2023a).</footnote>

<footnote n="31">³¹ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html</footnote>

<footnote n="32">³² TS = 2(log $\mathcal{L}$ – log $\mathcal{L}_0$), where log $\mathcal{L}$ and log $\mathcal{L}_0$ are the logarithms of the maximum likelihood of the complete source model and of the model without the target source included, respectively (Mattox et al. 1996).</footnote>

<footnote n="33">³³ http://fermipy.readthedocs.io/en/latest/index.html</footnote>

<footnote n="34">³⁴ https://fermi.gsfc.nasa.gov/ssc/data/analysis/software</footnote>
<table><caption>**Table 1**
Basic Properties of the Seven MSPs in M5</caption><thead><tr><th>Name</th><th>R.A.<br>(deg)</th><th>Decl.<br>(deg)</th><th>P<br>(ms)</th><th>P<sub>orb</sub><br>(day)</th><th>Eccent.<sup>a</sup></th><th>E&#x0307;<sup>b</sup><br>(10<sup>34</sup> erg s<sup>-1</sup>)</th><th>X-Ray<br>Counterpart</th><th>Optical<br>Counterpart</th><th>Comment</th></tr></thead><tbody><tr><td>M5A</td><td>229.6388</td><td>2.0910</td><td>5.55</td><td>...</td><td>...</td><td>&lt;1.56</td><td>No</td><td>No</td><td>Isolated</td></tr><tr><td>M5B</td><td>229.6311</td><td>2.0876</td><td>7.95</td><td>6.858</td><td>0.138</td><td>&lt;0.25</td><td>No</td><td>No</td><td>Heavy, not likely edge-on</td></tr><tr><td>M5C</td><td>229.6366</td><td>2.0799</td><td>2.48</td><td>0.687</td><td>0</td><td>&lt;12.36</td><td>Yes</td><td>Yes</td><td>Eclipsing BW</td></tr><tr><td>M5D</td><td>229.6268</td><td>2.0833</td><td>2.99</td><td>1.222</td><td>lower than M5B</td><td>0.28–4.98</td><td>Yes</td><td>Yes</td><td>He WD companion</td></tr><tr><td>M5E</td><td>229.6388</td><td>2.0772</td><td>3.18</td><td>1.097</td><td>lower than M5B</td><td>&lt;5.27</td><td>Yes</td><td>Yes</td><td>He WD companion</td></tr><tr><td>M5F</td><td>229.6350</td><td>2.0867</td><td>2.65</td><td>1.610</td><td>lower than M5B</td><td>&lt;8.30</td><td>No</td><td>Yes</td><td>He WD companion</td></tr><tr><td>M5G</td><td>229.6197</td><td>2.0875</td><td>2.75</td><td>0.114</td><td>0</td><td>0.04–3.34</td><td>Yes</td><td>No</td><td>Noneclipsing BW</td></tr></tbody></table>

**Notes.**

<105° were used. This reduces the contribution of the Earth limb contamination to the total background to less than 10%.

We built a spatial-spectral model for M5 by including the 4FGL-DR4 sources within 20° around the M5 center. The Galactic interstellar emission model (“gll_iem_v07.fits”) and the isotropic emission spectrum (“iso_P8R3_SOUR-CE_V3_v1.txt”), which takes into account the extragalactic emission and the residual instrumental background,<fnref n="31" /> were also included. Energy dispersion has been taken into account by adding two extra energy bins except for the isotropic component. We have performed a summed likelihood analysis in a 14° × 14° region of interest (ROI) around M5. The significance of a given source in the model is characterized by the test statistic (TS).<fnref n="32" /> In this work, we used the fermipy<fnref n="33" /> package (v1.2.0; Wood et al. 2017), in which Fermitools<fnref n="34" /> (v2.2.0) is integrated.

## 3. Analysis Results

### 3.1. Localization

Since the LAT cannot resolve M5, we modeled it as a point source in our analysis. We first localized it in the higher energy range of 1–500 GeV to take advantage of LAT’s better spatial resolution. Data were binned using a 0.05 × 0.05 pixel size and 10 logarithmic energy bins per decade. The ROI was optimized by fitting all sources in the ROI to ensure that all parameters are close to their global likelihood maxima. We tested three spectral models for M5 during the localization: the default LogParabola (LP) model in 4FGL-DR4,

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-(\alpha + \beta \ln(E/E_0))}, \qquad (1) $$

the simple power-law (PL) model,

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\Gamma_0}, \qquad (2) $$

and the PLSuperExpCutoff4 (PLEC4) model that is used for pulsar emission in the Third Fermi-LAT Catalog of Gamma-ray Pulsars (3PC; Smith et al. 2023),

$$ \frac{dN}{dE} = N_0 \left(\frac{E}{E_0}\right)^{-\Gamma+d/b} \exp\left[\frac{d}{b^2}\left(1 - \left(\frac{E}{E_0}\right)^b\right)\right], \qquad (3) $$

where $N_0$ is the normalization; $\alpha$ and $\Gamma_0$ are the spectral indices in the LP and PL models, respectively; $\beta$ is the curvature index; $E_0$ is the reference energy and is fixed to the catalog value for the LP model and to 1 GeV for the PL and PLEC4 models; $\Gamma$ is the local spectral index at $E_0$ in the PLEC4 model; $d$ is the local curvature at $E_0$; and $b$ is fixed to 2/3.

Only the normalizations of the Galactic/isotropic diffuse components were set free to vary during the localization. We note that changing the spectral model and leaving the position free simultaneously in each fit will lead to “un-nested models.” However, the localizations determined from the three models are completely consistent (Table 2). Both PLEC4 and LP models have the same number of degrees of freedom and fit the data quite well. Since the TS values obtained with both models are basically the same, we chose to continue our analysis with the PLEC4 model simply because it is physically motivated as a superposition of curvature radiation spectra for a range of electron energies (Bednarek & Sitarek 2007; Venter & de Jager 2008; Venter et al. 2009; Cheng et al. 2010).

Figure 1 (left panel) shows the residual map for the full ROI for the localization generated by removing the contribution of all 4FGL sources in the ROI. No significant residuals (>4σ) were found, indicating a good modeling of the ROI. The TS excess map (Figure 2, left panel) presents the localization result in a 1° × 1° region centered at M5 by removing the contribution of all 4FGL sources except M5, i.e., all sources except M5 are included in the model.

### 3.2. Spectral Fit

We performed a broadband spectral fit for M5 in the energy range of 0.1–500 GeV using the best-fit localization obtained with the PLEC4 model. Data were binned using a 0.1° × 0.1° pixel size and 10 logarithmic energy bins per decade. This time, all the spectral parameters of sources within 5° and those of the Galactic/isotropic diffuse components were set free to vary, along with the normalizations of significantly variable sources within 5°–8° around M5.

<footnote><sup>a</sup> The eccentricities of M5C and M5G have been set to 0 based on the assumption that the orbits of BW pulsars are circular owing to tidal dissipation (see Table 1 in Zhang et al. 2023a).</footnote>

<footnote><sup>b</sup> Spin-down powers $\dot{E}$ have been calculated based on the intrinsic spin-down rate $\dot{P}_{\text{int}}$ upper limit corrected for accelerations caused by the gravitational field of the GC for the line of sight of each pulsar (see Table 2 in Zhang et al. 2023a).</footnote>

<footnote n="31">³¹ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html</footnote>

<footnote n="32">³² TS = 2(log $\mathcal{L}$ – log $\mathcal{L}_0$), where log $\mathcal{L}$ and log $\mathcal{L}_0$ are the logarithms of the maximum likelihood of the complete source model and of the model without the target source included, respectively (Mattox et al. 1996).</footnote>

<footnote n="33">³³ http://fermipy.readthedocs.io/en/latest/index.html</footnote>

<footnote n="34">³⁴ https://fermi.gsfc.nasa.gov/ssc/data/analysis/software</footnote>
<table><caption>**Table 2**
Fermi-LAT Analysis Results for M5</caption><thead><tr><th colspan="5">1–500 GeV Localization</th></tr><tr><th>Model</th><th>TS</th><th>GLON<br>(deg)</th><th>GLAT<br>(deg)</th><th>95% Localization<br>Uncertainty<br>(deg)</th></tr></thead><tbody><tr><td>LP</td><td>73.3</td><td>3.8458 ± 0.0253</td><td>46.7606 ± 0.0282</td><td>0.07</td></tr><tr><td>PL</td><td>71.1</td><td>3.8409 ± 0.0261</td><td>46.7598 ± 0.0286</td><td>0.07</td></tr><tr><td>PLEC4</td><td>73.8</td><td>3.8453 ± 0.0254</td><td>46.7622 ± 0.0283</td><td>0.07</td></tr></tbody></table><table><thead><tr><th colspan="7">0.1–500 GeV Spectral Fits</th></tr><tr><th></th><th>TS</th><th>α</th><th>β</th><th>Γ (Γ<sub>0</sub>)</th><th>d</th><th>Photon Flux<br>(10<sup>-9</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Energy Flux<br>(10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>LP</td><td>90.9</td><td>2.11 ± 0.20</td><td>0.44 ± 0.18</td><td>...</td><td>...</td><td>1.41 ± 0.60</td><td>1.75 ± 0.32</td></tr><tr><td>PL</td><td>74.7</td><td>...</td><td>...</td><td>2.24 ± 0.09</td><td>...</td><td>4.02 ± 0.89</td><td>2.90 ± 0.41</td></tr><tr><td>PLEC4</td><td>91.2</td><td>...</td><td>...</td><td>1.70 ± 0.25</td><td>0.63 ± 0.24</td><td>1.65 ± 0.67</td><td>1.78 ± 0.32</td></tr></tbody></table>

Note. LP stands for LogParabola, PL for power law, and PLEC4 for PLSuperExpCutoff4. Γ₀ is the spectral index for the PL model, and Γ is the local spectral index for the PLEC4 model.

Similarly, we tested the three aforementioned spectral models, and we summarize the best-fit results in Table 2. The PLEC4 and LP models are indistinguishable, while the PL model is the worst. We again took PLEC4 as the best-fit model based on its physical motivation. Figure 1 (right panel) shows the residual map for the full ROI for the broadband fit in which no significant residuals (>$4\sigma$) are present.

Taking the best-fit PLEC4 model in 0.1–500 GeV, we computed the spectral energy distribution (SED) by performing a maximum likelihood analysis in 10 logarithmically spaced energy bins over 0.1–500 GeV. Background sources that were free in the broadband fit were kept free in the fit of each bin. The normalization of M5 in each bin is fit using a PL spectral parameterization with a fixed index. At lower energy bands (the first six bins), the index was fixed to the local slope of the broadband PLEC4 model. At higher energy bands (the last four bins), the index was fixed to 4 in order to be consistent with the local slope used at lower energies, given that the PLEC4 model is both very steep and very low in normalization at higher energies. Upper limits on the flux at the 95% confidence level (CL) were computed for bins when the source had TS < 4. Figure 3 shows the SED, along with the PLEC4, LP, and PL models.

### 3.3. Gamma-Ray Pulsation Search

A gamma-ray pulsation search for each pulsar has been performed by selecting photons within 2° of the pulsar position. Spin phases of gamma-ray photons were calculated using the Fermi plug-in (Ray et al. 2011) for TEMPO2 and the radio ephemerides for each pulsar (Zhang et al. 2023a), respectively. The ephemerides for M5A, M5B, M5C, M5D, and M5E are valid from before Fermi was launched to 2022 November, while those for M5F and M5G are valid for the time range of 2020 November 16–2022 December 14, since these two MSPs were the newest in M5 and only detected by FAST. We used the weighted H-test (Kerr 2011) statistic to quantify the pulsation significance, with weights computed by employing the Simple Weights method as described in Bruel (2019) and Smith et al. (2019) for both the full LAT data set and that in the time range of the ephemerides validity. No significant pulsations have been detected from any of the seven MSPs. The largest H-test value is 13.8 found for M5A

corresponding to around 2.9σ, which decreases slightly when considering the six trials used in the search.

### 3.4. Gamma-Ray Variability

#### 3.4.1. Long-term Light Curves

To investigate the long-term gamma-ray flux variability of M5, we computed a light curve with a 90-day binning over 0.1–500 GeV (Figure 4, left panel). The best fit of the PLEC4 model in the 0.1–500 GeV band obtained previously was used as the starting point for each time bin, but this time only the normalizations were kept free. An independent binned likelihood analysis was performed for each bin to get the flux of M5. Upper limits at the 95% CL were calculated when M5 had TS < 4. We followed the same method as presented in Acero et al. (2015) to quantify the variability significance and obtained TS<sub>var</sub> = 62.97. In a $\chi^2$ distribution with 60 degrees of freedom, the 99% confidence TS threshold above which the variability would be considered probable is 88.38. Thus, the gamma-ray emission from M5 is consistent with being steady on a timescale of a few months.

However, there is a TS peak of ~28 in the time bin 21 spanning MJD 56572–56662 (2013 October 7–2014 January 5). To further investigate the peak, we computed a TS map for this bin with M5 removed from the source model (Figure 2, right panel). A potential excess appeared near M5. We then localized this excess to $(l, b) = (3^\circ, 71, 47^\circ, 58)$ with a 95% localization uncertainty of 0.29, which, despite being large, does not enclose M5 (0.78 offset from the M5 center). Therefore, this excess is not related to M5. We verified that the excess does not appear in the bins before and after bin 21. A likelihood fit in this bin with the excess added resulted in a TS of around 20 and 10 for M5 and the excess, respectively. It is interesting to note that a TeV source that is offset by 4' from the center of the GC Terzan 5 and extended well beyond the HESS PSF was detected with unconfirmed origin (H. E. S. S. Collaboration 2011).

We have searched multiple catalogs for possible counterparts within the 95% localization uncertainty of the excess. We first checked the Fermi-LAT Long-Term Transient Source Catalog (1FLT; Baldini et al. 2021),<fnref n="35" /> which contains transient sources

<footnote n="35">³⁵ https://heasarc.gsfc.nasa.gov/W3Browse/all/fermirltrns.html</footnote>
<table><caption>**Table 2**
Fermi-LAT Analysis Results for M5</caption><thead><tr><th colspan="5">1-500 GeV Localization</th></tr><tr><th>Model</th><th>TS</th><th>GLON<br>(deg)</th><th>GLAT<br>(deg)</th><th>95% Localization<br>Uncertainty<br>(deg)</th></tr></thead><tbody><tr><td>LP</td><td>73.3</td><td>3.8458 ± 0.0253</td><td>46.7606 ± 0.0282</td><td>0.07</td></tr><tr><td>PL</td><td>71.1</td><td>3.8409 ± 0.0261</td><td>46.7598 ± 0.0286</td><td>0.07</td></tr><tr><td>PLEC4</td><td>73.8</td><td>3.8453 ± 0.0254</td><td>46.7622 ± 0.0283</td><td>0.07</td></tr></tbody></table><table><thead><tr><th colspan="7">0.1-500 GeV Spectral Fits</th></tr><tr><th></th><th>TS</th><th>α</th><th>β</th><th>Γ (Γ<sub>0</sub>)</th><th>d</th><th>Photon Flux<br>(10<sup>-9</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Energy Flux<br>(10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>LP</td><td>90.9</td><td>2.11 ± 0.20</td><td>0.44 ± 0.18</td><td>...</td><td>...</td><td>1.41 ± 0.60</td><td>1.75 ± 0.32</td></tr><tr><td>PL</td><td>74.7</td><td>...</td><td>...</td><td>2.24 ± 0.09</td><td>...</td><td>4.02 ± 0.89</td><td>2.90 ± 0.41</td></tr><tr><td>PLEC4</td><td>91.2</td><td>...</td><td>...</td><td>1.70 ± 0.25</td><td>0.63 ± 0.24</td><td>1.65 ± 0.67</td><td>1.78 ± 0.32</td></tr></tbody></table>

Note. LP stands for LogParabola, PL for power law, and PLEC4 for PLSuperExpCutoff4. Γ₀ is the spectral index for the PL model, and Γ is the local spectral index for the PLEC4 model.

Similarly, we tested the three aforementioned spectral models, and we summarize the best-fit results in Table 2. The PLEC4 and LP models are indistinguishable, while the PL model is the worst. We again took PLEC4 as the best-fit model based on its physical motivation. Figure 1 (right panel) shows the residual map for the full ROI for the broadband fit in which no significant residuals (>4σ) are present.

Taking the best-fit PLEC4 model in 0.1–500 GeV, we computed the spectral energy distribution (SED) by performing a maximum likelihood analysis in 10 logarithmically spaced energy bins over 0.1–500 GeV. Background sources that were free in the broadband fit were kept free in the fit of each bin. The normalization of M5 in each bin is fit using a PL spectral parameterization with a fixed index. At lower energy bands (the first six bins), the index was fixed to the local slope of the broadband PLEC4 model. At higher energy bands (the last four bins), the index was fixed to 4 in order to be consistent with the local slope used at lower energies, given that the PLEC4 model is both very steep and very low in normalization at higher energies. Upper limits on the flux at the 95% confidence level (CL) were computed for bins when the source had TS < 4. Figure 3 shows the SED, along with the PLEC4, LP, and PL models.

### 3.3. Gamma-Ray Pulsation Search

A gamma-ray pulsation search for each pulsar has been performed by selecting photons within 2° of the pulsar position. Spin phases of gamma-ray photons were calculated using the Fermi plug-in (Ray et al. 2011) for TEMPO2 and the radio ephemerides for each pulsar (Zhang et al. 2023a), respectively. The ephemerides for M5A, M5B, M5C, M5D, and M5E are valid from before Fermi was launched to 2022 November, while those for M5F and M5G are valid for the time range of 2020 November 16–2022 December 14, since these two MSPs were the newest in M5 and only detected by FAST. We used the weighted H-test (Kerr 2011) statistic to quantify the pulsation significance, with weights computed by employing the Simple Weights method as described in Bruel (2019) and Smith et al. (2019) for both the full LAT data set and that in the time range of the ephemerides validity. No significant pulsations have been detected from any of the seven MSPs. The largest H-test value is 13.8 found for M5A

corresponding to around 2.9σ, which decreases slightly when considering the six trials used in the search.

### 3.4. Gamma-Ray Variability

#### 3.4.1. Long-term Light Curves

To investigate the long-term gamma-ray flux variability of M5, we computed a light curve with a 90-day binning over 0.1–500 GeV (Figure 4, left panel). The best fit of the PLEC4 model in the 0.1–500 GeV band obtained previously was used as the starting point for each time bin, but this time only the normalizations were kept free. An independent binned likelihood analysis was performed for each bin to get the flux of M5. Upper limits at the 95% CL were calculated when M5 had TS < 4. We followed the same method as presented in Acero et al. (2015) to quantify the variability significance and obtained TS<sub>var</sub> = 62.97. In a χ² distribution with 60 degrees of freedom, the 99% confidence TS threshold above which the variability would be considered probable is 88.38. Thus, the gamma-ray emission from M5 is consistent with being steady on a timescale of a few months.

However, there is a TS peak of ~28 in the time bin 21 spanning MJD 56572–56662 (2013 October 7–2014 January 5). To further investigate the peak, we computed a TS map for this bin with M5 removed from the source model (Figure 2, right panel). A potential excess appeared near M5. We then localized this excess to (l, b) = (3°, 71, 47°, 58) with a 95% localization uncertainty of 0°, 29, which, despite being large, does not enclose M5 (0°, 78 offset from the M5 center). Therefore, this excess is not related to M5. We verified that the excess does not appear in the bins before and after bin 21. A likelihood fit in this bin with the excess added resulted in a TS of around 20 and 10 for M5 and the excess, respectively. It is interesting to note that a TeV source that is offset by 4' from the center of the GC Terzan 5 and extended well beyond the HESS PSF was detected with unconfirmed origin (H. E. S. S. Collaboration 2011).

We have searched multiple catalogs for possible counterparts within the 95% localization uncertainty of the excess. We first checked the Fermi-LAT Long-Term Transient Source Catalog (1FLT; Baldini et al. 2021),<fnref n="35" /> which contains transient sources

<footnote n="35">³⁵ https://heasarc.gsfc.nasa.gov/W3Browse/all/fermiltrns.html</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Residual maps of M5 for the full ROI in $\sigma$ units (in Galactic coordinates). Black plus signs are 4FGL-DR4 sources included in the ROI. Left: in the 1–500 GeV band for the localization of M5. Right: in the 0.1–500 GeV band for the spectral fit for M5.</figcaption></figure>

detected above $4\sigma$ on monthly timescales using 10 yr of Fermi-LAT data. Given the low significance of the excess, we expected to find no counterpart in the 1FLT catalog, and indeed there is none. Then, through HEASARC,<fnref n="36" /> we searched for blazar candidates in the FERMILBLAZ (Kovačević et al. 2019), BZCat (Massaro et al. 2015), CGRaBS (Healey et al. 2008), CRATES (Healey et al. 2007), WIBRaLs2, and KDEBLLACs catalogs (D’Abrusco et al. 2019; de Menezes et al. 2019), but no matches were found for the excess. We note as well that, given the high Galactic latitude of M5, it is highly unlikely to find novae, and indeed no novae were reported in the relevant time period.<fnref n="37" />

After adding the potential excess to the source model, the $TS_{var}$ of the 90-day light curve (Figure 4, left panel) decreased to 57.9. Nevertheless, the $TS$ of M5 in bin 21 after taking the excess into account remains an outlier, although $TS$ is not a good measure of variability (it can vary owing to nearby sources or just exposure). To quantify whether the apparent high $TS$ of M5 in bin 21 is significant or not, we compared the best-fit $log L$ in this bin with the $log L$ of the fit when fixing the flux of M5 to the average, both with the excess added in this bin. We then computed $\sqrt{2\Delta log L}$, i.e., $Sqrt_TS_History$ in the 4FGL catalog, which gives the significance in $\sigma$ units for 1 degree of freedom. We obtained $2.8\sigma$, moderately significant. The apparent $TS$ peak in bin 21 is thus compatible with being a statistical fluctuation. Besides, we zoomed in on bin 21 by computing a 10-day light curve (Figure 4, right panel) for it. The $TS_{var}$ is 15.2, while the 99% confidence $TS$ threshold with 8 degrees of freedom is 20.1. Thus, the gamma-ray emission inside bin 21 is compatible with a constant signal and is also compatible with bin 21 being a positive statistical fluctuation as concluded above.

### 3.4.2. Orbital Modulation

Binary MSPs can exhibit orbitally modulated emission in multiple wavelengths such as optical, X-ray, and gamma rays. Detecting modulation constrains the binary system properties and theoretical models for the multiwavelength emission of

such systems. Orbital phases were assigned to each event using the Fermi plug-in (Ray et al. 2011) for TEMPO2 and the radio timing solution for each pulsar (Zhang et al. 2023a), respectively. We adopted two approaches to study the orbital modulation of the six binary MSPs in M5, similar to what was done in Johnson et al. (2015). We first created a counts light curve of 5° around the M5 center with a 30 s time bin and folded it with the orbital period of each pulsar.

In the first approach, we calculated the LAT exposure for each time bin in the counts light curve to account for any possible orbital variations of the exposure (Ackermann et al. 2012). By binning the exposure in 1000 bins of orbital phase and normalizing, we built a null distribution of what the orbital modulation should look like if the exposure variation versus orbital phases is the only “modulation” present. We then used this null distribution to get the exposure-corrected orbital phase for each event in a region of 2° around M5. Finally, we employed the weighted H-test (Kerr 2011) to quantify the orbital modulation for each pulsar. Weights were calculated in a similar way to that in the gamma-ray pulsation search (Section 3.3). No strong evidence of modulation was found for any of the six MSPs. The largest H-test value found was 10.3 for M5B corresponding to $2.4\sigma$, without trial corrections.

In the second approach, we computed the orbital flux for each MSP with 10 bins per orbit following the same methodology as for the long-term light curve presented above. One extra step before performing the binned likelihood fit in each orbital bin is to correct the potential LAT exposure variations across the orbit by creating orbital-phase-selected good time intervals from the orbital-period-folded 30 s count light curve that was generated previously. Similar to the long-term light-curve analysis, we computed $TS_{var}$ to quantify the flux variability. The orbital flux modulations using this approach are shown in Figure 5. The $\chi^2$ distribution with 9 degrees of freedom corresponds to a 99% confidence $TS$ threshold of 21.7. None of the six MSPs has $TS_{var}$ larger than this value, with the largest being 14.7 found for M5B. Thus, no significant variability was found for them. Despite M5B’s apparently high $TS$ and flux in the orbital bin of 0.5–0.6 (Table 3 and Figure 5), the variability over the full orbit is not significant.

<footnote n="36">³⁶ https://heasarc.gsfc.nasa.gov/W3Browse/all/</footnote>

<footnote n="37">³⁷ https://asd.gsfc.nasa.gov/Koji.Mukai/novae/novae.html</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Residual maps of M5 for the full ROI in $\sigma$ units (in Galactic coordinates). Black plus signs are 4FGL-DR4 sources included in the ROI. Left: in the 1–500 GeV band for the localization of M5. Right: in the 0.1–500 GeV band for the spectral fit for M5.</figcaption></figure>

detected above $4\sigma$ on monthly timescales using 10 yr of Fermi-LAT data. Given the low significance of the excess, we expected to find no counterpart in the 1FLT catalog, and indeed there is none. Then, through HEASARC,<fnref n="36" /> we searched for blazar candidates in the FERMILBLAZ (Kovačević et al. 2019), BZCat (Massaro et al. 2015), CGRaBS (Healey et al. 2008), CRATES (Healey et al. 2007), WIBRaLs2, and KDEBLLACs catalogs (D’Abrusco et al. 2019; de Menezes et al. 2019), but no matches were found for the excess. We note as well that, given the high Galactic latitude of M5, it is highly unlikely to find novae, and indeed no novae were reported in the relevant time period.<fnref n="37" />

After adding the potential excess to the source model, the $TS_{var}$ of the 90-day light curve (Figure 4, left panel) decreased to 57.9. Nevertheless, the $TS$ of M5 in bin 21 after taking the excess into account remains an outlier, although $TS$ is not a good measure of variability (it can vary owing to nearby sources or just exposure). To quantify whether the apparent high $TS$ of M5 in bin 21 is significant or not, we compared the best-fit $\log L$ in this bin with the $\log L$ of the fit when fixing the flux of M5 to the average, both with the excess added in this bin. We then computed $\sqrt{2\Delta\log L}$, i.e., $\text{Sqrt\_TS\_History}$ in the 4FGL catalog, which gives the significance in $\sigma$ units for 1 degree of freedom. We obtained $2.8\sigma$, moderately significant. The apparent $TS$ peak in bin 21 is thus compatible with being a statistical fluctuation. Besides, we zoomed in on bin 21 by computing a 10-day light curve (Figure 4, right panel) for it. The $TS_{var}$ is 15.2, while the 99% confidence $TS$ threshold with 8 degrees of freedom is 20.1. Thus, the gamma-ray emission inside bin 21 is compatible with a constant signal and is also compatible with bin 21 being a positive statistical fluctuation as concluded above.

### 3.4.2. Orbital Modulation

Binary MSPs can exhibit orbitally modulated emission in multiple wavelengths such as optical, X-ray, and gamma rays. Detecting modulation constrains the binary system properties and theoretical models for the multiwavelength emission of

such systems. Orbital phases were assigned to each event using the Fermi plug-in (Ray et al. 2011) for TEMPO2 and the radio timing solution for each pulsar (Zhang et al. 2023a), respectively. We adopted two approaches to study the orbital modulation of the six binary MSPs in M5, similar to what was done in Johnson et al. (2015). We first created a counts light curve of 5° around the M5 center with a 30 s time bin and folded it with the orbital period of each pulsar.

In the first approach, we calculated the LAT exposure for each time bin in the counts light curve to account for any possible orbital variations of the exposure (Ackermann et al. 2012). By binning the exposure in 1000 bins of orbital phase and normalizing, we built a null distribution of what the orbital modulation should look like if the exposure variation versus orbital phases is the only “modulation” present. We then used this null distribution to get the exposure-corrected orbital phase for each event in a region of 2° around M5. Finally, we employed the weighted H-test (Kerr 2011) to quantify the orbital modulation for each pulsar. Weights were calculated in a similar way to that in the gamma-ray pulsation search (Section 3.3). No strong evidence of modulation was found for any of the six MSPs. The largest H-test value found was 10.3 for M5B corresponding to $2.4\sigma$, without trial corrections.

In the second approach, we computed the orbital flux for each MSP with 10 bins per orbit following the same methodology as for the long-term light curve presented above. One extra step before performing the binned likelihood fit in each orbital bin is to correct the potential LAT exposure variations across the orbit by creating orbital-phase-selected good time intervals from the orbital-period-folded 30 s count light curve that was generated previously. Similar to the long-term light-curve analysis, we computed $TS_{var}$ to quantify the flux variability. The orbital flux modulations using this approach are shown in Figure 5. The $\chi^2$ distribution with 9 degrees of freedom corresponds to a 99% confidence $TS$ threshold of 21.7. None of the six MSPs has $TS_{var}$ larger than this value, with the largest being 14.7 found for M5B. Thus, no significant variability was found for them. Despite M5B’s apparently high $TS$ and flux in the orbital bin of 0.5–0.6 (Table 3 and Figure 5), the variability over the full orbit is not significant.

<footnote n="36">³⁶ https://heasarc.gsfc.nasa.gov/W3Browse/all/</footnote>

<footnote n="37">³⁷ https://asd.gsfc.nasa.gov/Koji.Mukai/novae/novae.html</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Left: 1° × 1° TS excess map of M5 (in Galactic coordinates) in the 1–500 GeV band with a bin size of 0.05. Overlaid are the LAT best localization and 95% localization uncertainty in this work (red plus sign and circle), 4FGL-DR4 position and 95% localization uncertainty (green cross and ellipse), and M5 center (blue cross). Right: 2.5 × 2.5 TS excess map of M5 in the 0.1–500 GeV band with a bin size of 0.1. In addition to the markers for M5, the green plus sign and circle show the best localization and 95% localization uncertainty for the excess in bin 21 of the 90-day light curve (see Section 3.4.1 and Figure 4).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Gamma-ray spectrum of M5 (red circles), along with the broadband models of PLEC4 with uncertainty (black line and shaded region), LP with uncertainty (blue line and shaded region), and PL with uncertainty (green line and shaded region). Flux upper limits at the 95% CL are shown as red arrows for bins when the source had TS < 4.</figcaption></figure>

Similar to the long-term light-curve analysis, we computed $\sqrt{2\Delta \log L}$ to quantify whether the high TS and flux in the bin of 0.5–0.6 of M5B are significant compared to the phase-averaged flux (Table 2). We obtained $3.4\sigma$. Taking into account the 60 trials (number of orbital bins 10 multiplied by the number of pulsars 6), the probability to get a $3.4\sigma$ excess is 4%, corresponding to $2\sigma$ after trial corrections. The flux and TS variation of M5B are thus not significant and are probably due to statistical fluctuations. To further investigate the orbital modulation, we then compared the spectral characteristics in the orbital bin of 0.5–0.6 and 0.6–0.5 by reperforming likelihood fits in the two bins. This is done similarly to the standard orbital modulation analysis described above, but this time the spectral shape parameters of M5 were also set free to vary in order to test any shape changes in addition to the overall flux variation. The results are presented in Table 3. The flux difference between 0.5–0.6 and 0.6–0.5, as well as the

phase-averaged fit (Table 2), is less than ~3σ, consistent with being a statistical fluctuation as stated above. The spectral shape parameters are, on the other hand, consistent within uncertainties.

## 4. Discussion and Concluding Remarks

From the spectral parameter values we obtained in the fit, we can evaluate the energy at which the SED peaks as

$$E_p = E_0 \left[ \left( 1 + \frac{b}{d} (2 - \Gamma) \right)^{\frac{1}{b}} \right] \quad (4)$$

and the curvature at the SED peak $d_p$ as

$$d_p = d + b(2 - \Gamma), \quad (5)$$

as outlined in 3PC (Smith et al. 2023), with $d_p$ reaching a maximum of 4/3 for synchrotron or curvature radiation from monoenergetic electrons. $E_p$ and $d_p$ are correlated. The SED peak width is inversely proportional to $d_p$ such that high curvature indicates a narrow spectrum corresponding to a narrow range of electron energies and low curvature indicates a broad spectrum with contributions from a broader range of electron energies. We obtained $E_p = 1.5$ GeV and $d_p = 0.83$ for M5, putting it near the upper right corner of the $d_p$ versus $E_p$ plot (see Figure 20 in the 3PC). Following the 3PC and taking into account the values we obtained, if the emission detected from the GC comes mostly from a single source, or if it comes from many sources behaving similarly, the electron population producing the emission is rather broad.

### 4.1. Nondetection of Individual Pulsars and Implications

While 305 radio pulsars in GCs are currently known, with most of them being MSPs, only three individual GC MSPs have been detected in gamma rays, or about 1%. This is partly due to the large distances of GCs, which imply that any MSPs detected must be exceptionally energetic (especially in the earlier phases of the Fermi mission: the first two GC MSPs detected in gamma rays are by far the most powerful known,
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: 1° × 1° TS excess map of M5 (in Galactic coordinates) in the 1–500 GeV band with a bin size of 0.05. Overlaid are the LAT best localization and 95% localization uncertainty in this work (red plus sign and circle), 4FGL-DR4 position and 95% localization uncertainty (green cross and ellipse), and M5 center (blue cross). Right: 2.5 × 2.5 TS excess map of M5 in the 0.1–500 GeV band with a bin size of 0.1. In addition to the markers for M5, the green plus sign and circle show the best localization and 95% localization uncertainty for the excess in bin 21 of the 90-day light curve (see Section 3.4.1 and Figure 4).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Gamma-ray spectrum of M5 (red circles), along with the broadband models of PLEC4 with uncertainty (black line and shaded region), LP with uncertainty (blue line and shaded region), and PL with uncertainty (green line and shaded region). Flux upper limits at the 95% CL are shown as red arrows for bins when the source had TS < 4.</figcaption></figure>

Similar to the long-term light-curve analysis, we computed $\sqrt{2\Delta \log L}$ to quantify whether the high TS and flux in the bin of 0.5–0.6 of M5B are significant compared to the phase-averaged flux (Table 2). We obtained $3.4\sigma$. Taking into account the 60 trials (number of orbital bins 10 multiplied by the number of pulsars 6), the probability to get a $3.4\sigma$ excess is 4%, corresponding to $2\sigma$ after trial corrections. The flux and TS variation of M5B are thus not significant and are probably due to statistical fluctuations. To further investigate the orbital modulation, we then compared the spectral characteristics in the orbital bin of 0.5–0.6 and 0.6–0.5 by reperforming likelihood fits in the two bins. This is done similarly to the standard orbital modulation analysis described above, but this time the spectral shape parameters of M5 were also set free to vary in order to test any shape changes in addition to the overall flux variation. The results are presented in Table 3. The flux difference between 0.5–0.6 and 0.6–0.5, as well as the

phase-averaged fit (Table 2), is less than ~3σ, consistent with being a statistical fluctuation as stated above. The spectral shape parameters are, on the other hand, consistent within uncertainties.

## 4. Discussion and Concluding Remarks

From the spectral parameter values we obtained in the fit, we can evaluate the energy at which the SED peaks as

$$E_p = E_0 \left[ \left( 1 + \frac{b}{d} (2 - \Gamma) \right)^{\frac{1}{b}} \right] \quad (4)$$

and the curvature at the SED peak $d_p$ as

$$d_p = d + b(2 - \Gamma), \quad (5)$$

as outlined in 3PC (Smith et al. 2023), with $d_p$ reaching a maximum of 4/3 for synchrotron or curvature radiation from monoenergetic electrons. $E_p$ and $d_p$ are correlated. The SED peak width is inversely proportional to $d_p$ such that high curvature indicates a narrow spectrum corresponding to a narrow range of electron energies and low curvature indicates a broad spectrum with contributions from a broader range of electron energies. We obtained $E_p = 1.5$ GeV and $d_p = 0.83$ for M5, putting it near the upper right corner of the $d_p$ versus $E_p$ plot (see Figure 20 in the 3PC). Following the 3PC and taking into account the values we obtained, if the emission detected from the GC comes mostly from a single source, or if it comes from many sources behaving similarly, the electron population producing the emission is rather broad.

### 4.1. Nondetection of Individual Pulsars and Implications

While 305 radio pulsars in GCs are currently known, with most of them being MSPs, only three individual GC MSPs have been detected in gamma rays, or about 1%. This is partly due to the large distances of GCs, which imply that any MSPs detected must be exceptionally energetic (especially in the earlier phases of the Fermi mission: the first two GC MSPs detected in gamma rays are by far the most powerful known,
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Long-term light curve and TS evolution for M5 in the energy range of 0.1–500 GeV. Left: 90-day binning. Right: a zoom-in around bin 21 of the 90-day light curve with a binning of 10 days (from MJD 56572 to MJD 56662). Flux upper limits at the 95% CL are shown as arrows for bins when the source had TS < 4. Blue dotted line: average flux in the broad band. Red dotted lines: $1\sigma$ uncertainties on the average flux. The magenta marker indicates the flux and TS in bin 21 after taking into account the excess near M5.</figcaption></figure>

with spin-down power $\dot{E} > 8 \times 10^{35}$ erg s$^{-1}$). However, this might also be due to the fact that not many GC pulsars have ephemerides that are accurate for a large fraction of the Fermi mission's 15 yr. Thus, whenever such ephemerides become available, it is important to verify whether the pulsars can be detected in gamma rays: any such detection would automatically imply an exceptionally energetic MSP.

From Table 1, we can see that some of the MSPs in M5 have $\dot{E}$ upper limits that are comparable with not only those of gamma-ray MSPs (Smith et al. 2023) but even those of gamma-ray MSPs detected in GCs: for instance, M5C has an upper $\dot{E}$ limit of $1.2 \times 10^{35}$ erg s$^{-1}$, while for J1835–3259B (where we can estimate the intrinsic spin-down from the measured spin-down and the variation of the orbital period reported by Gautam et al. 2022) $\dot{E} = 1.8 \pm 1.2 \times 10^{35}$ erg s$^{-1}$. Furthermore, the detectability of the pulsars in M5 is enhanced by the fact that the photon background toward M5 is significantly smaller than that toward NGC 6652, whose gamma-ray emission is probably dominated by a single pulsar.

The nondetection suggests that the true values of $\dot{E}$ for the M5 pulsars are well below the upper limits derived by Zhang et al. (2023a). This is not surprising given previously observed trends among the GC pulsars. The three pulsars detected in gamma rays are all located in GCs with a very high $\gamma_b$. Two of them (J1823–3021A and J1835–3259B) are located in core-collapsed GCs, NGC 6624 and NGC 6652, which generally have the highest values of $\gamma_b$ (Verbunt & Freire 2014).

The latter authors have remarked that most of the relatively slow (thus higher magnetic field *B*) pulsars are located in high-$\gamma_b$ clusters. They ascribed this to the high stellar encounter rate itself: this will disrupt many pulsar binaries, leading to a high rate of isolated pulsars in core-collapsed GCs, which has been repeatedly confirmed (e.g., Abbate et al. 2022, 2023). This high encounter rate could also be responsible for the prevalence of slow pulsars in these clusters if LMXBs are also being disrupted, leaving behind partially recycled NSs, which will appear as slower radio pulsars with larger *B*-fields than Galactic MSPs. It is possible that B1821–24A and J1823–3021A were formed in this way: although they were spun up, their *B*-fields had not been fully ablated, resulting in the very large magnetic braking torque and the unusually large $\dot{E}$.

Nevertheless, the disruption of compact binary systems in GCs by close stellar encounters is still an open topic of research and was recently contested by de Menezes et al. (2023). These authors used the Heggie–Hills law (Heggie 1975; Hills 1975) for binary encounters in combination with Fermi-LAT and Chandra data to argue that compact binary ionization would happen only in the unrealistic scenario where the dispersion velocity of stars in the cores of GCs is greater than the GCs' escape velocity.

In GCs like M5, with much smaller values of $\gamma_b$, any LMXBs, once formed (in exchange interactions), are not likely to be disturbed again, recycling their respective NSs right through to the end. Thus, we would expect all pulsars in these clusters to have not only fast spins but also small *B*-fields, as generally observed for MSPs in the Galactic disk. This is confirmed in 47 Tuc, where we can say that all MSPs in binaries (for which the cluster acceleration can be estimated from precise measurements of variation of the orbital period) have small values of $\dot{P}$, similar to those of Galactic MSPs (Freire et al. 2017). This is confirmed by the fact that none of the 23 pulsars with long-term timing solutions in 47 Tuc, or the fewer known in $\omega$ Centauri, are individually detectable in gamma rays (Dai et al. 2023). For the same reasons, we do not really expect the occurrence of pulsars with high *B*-fields in low-$\gamma_b$ GCs like M5.

Nevertheless, attempting to detect powerful gamma-ray MSPs in low-$\gamma_b$ GCs like M5 serves as a useful test of these ideas. If we find a single bright gamma-ray MSP in these clusters, we will know that their formation is not necessarily linked to processes that occur almost exclusively in high-$\gamma_b$ GCs, like LMXB disruption.

## 4.2. Gamma-Ray Emission of the Cluster as a Whole

Following the above discussion, the gamma-ray emission from M5 is therefore the resulting collective emission of individual pulsars in the cluster, similar to the findings in the study of 47 Tuc and $\omega$ Centauri. We can thus follow the method presented in Johnson et al. (2013) to estimate the total number
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Long-term light curve and TS evolution for M5 in the energy range of 0.1–500 GeV. Left: 90-day binning. Right: a zoom-in around bin 21 of the 90-day light curve with a binning of 10 days (from MJD 56572 to MJD 56662). Flux upper limits at the 95% CL are shown as arrows for bins when the source had TS < 4. Blue dotted line: average flux in the broad band. Red dotted lines: $1\sigma$ uncertainties on the average flux. The magenta marker indicates the flux and TS in bin 21 after taking into account the excess near M5.</figcaption></figure>

with spin-down power $\dot{E} > 8 \times 10^{35} \text{ erg s}^{-1}$). However, this might also be due to the fact that not many GC pulsars have ephemerides that are accurate for a large fraction of the Fermi mission's 15 yr. Thus, whenever such ephemerides become available, it is important to verify whether the pulsars can be detected in gamma rays: any such detection would automatically imply an exceptionally energetic MSP.

From Table 1, we can see that some of the MSPs in M5 have $\dot{E}$ upper limits that are comparable with not only those of gamma-ray MSPs (Smith et al. 2023) but even those of gamma-ray MSPs detected in GCs: for instance, M5C has an upper $\dot{E}$ limit of $1.2 \times 10^{35} \text{ erg s}^{-1}$, while for J1835–3259B (where we can estimate the intrinsic spin-down from the measured spin-down and the variation of the orbital period reported by Gautam et al. 2022) $\dot{E} = 1.8 \pm 1.2 \times 10^{35} \text{ erg s}^{-1}$. Furthermore, the detectability of the pulsars in M5 is enhanced by the fact that the photon background toward M5 is significantly smaller than that toward NGC 6652, whose gamma-ray emission is probably dominated by a single pulsar.

The nondetection suggests that the true values of $\dot{E}$ for the M5 pulsars are well below the upper limits derived by Zhang et al. (2023a). This is not surprising given previously observed trends among the GC pulsars. The three pulsars detected in gamma rays are all located in GCs with a very high $\gamma_b$. Two of them (J1823–3021A and J1835–3259B) are located in core-collapsed GCs, NGC 6624 and NGC 6652, which generally have the highest values of $\gamma_b$ (Verbunt & Freire 2014).

The latter authors have remarked that most of the relatively slow (thus higher magnetic field *B*) pulsars are located in high-$\gamma_b$ clusters. They ascribed this to the high stellar encounter rate itself: this will disrupt many pulsar binaries, leading to a high rate of isolated pulsars in core-collapsed GCs, which has been repeatedly confirmed (e.g., Abbate et al. 2022, 2023). This high encounter rate could also be responsible for the prevalence of slow pulsars in these clusters if LMXBs are also being disrupted, leaving behind partially recycled NSs, which will appear as slower radio pulsars with larger *B*-fields than Galactic MSPs. It is possible that B1821–24A and J1823–3021A were formed in this way: although they were spun up, their *B*-fields had not been fully ablated, resulting in the very large magnetic braking torque and the unusually large $\dot{E}$.

Nevertheless, the disruption of compact binary systems in GCs by close stellar encounters is still an open topic of research and was recently contested by de Menezes et al. (2023). These authors used the Heggie–Hills law (Heggie 1975; Hills 1975) for binary encounters in combination with Fermi-LAT and Chandra data to argue that compact binary ionization would happen only in the unrealistic scenario where the dispersion velocity of stars in the cores of GCs is greater than the GCs' escape velocity.

In GCs like M5, with much smaller values of $\gamma_b$, any LMXBs, once formed (in exchange interactions), are not likely to be disturbed again, recycling their respective NSs right through to the end. Thus, we would expect all pulsars in these clusters to have not only fast spins but also small *B*-fields, as generally observed for MSPs in the Galactic disk. This is confirmed in 47 Tuc, where we can say that all MSPs in binaries (for which the cluster acceleration can be estimated from precise measurements of variation of the orbital period) have small values of $\dot{P}$, similar to those of Galactic MSPs (Freire et al. 2017). This is confirmed by the fact that none of the 23 pulsars with long-term timing solutions in 47 Tuc, or the fewer known in $\omega$ Centauri, are individually detectable in gamma rays (Dai et al. 2023). For the same reasons, we do not really expect the occurrence of pulsars with high *B*-fields in low-$\gamma_b$ GCs like M5.

Nevertheless, attempting to detect powerful gamma-ray MSPs in low-$\gamma_b$ GCs like M5 serves as a useful test of these ideas. If we find a single bright gamma-ray MSP in these clusters, we will know that their formation is not necessarily linked to processes that occur almost exclusively in high-$\gamma_b$ GCs, like LMXB disruption.

## 4.2. Gamma-Ray Emission of the Cluster as a Whole

Following the above discussion, the gamma-ray emission from M5 is therefore the resulting collective emission of individual pulsars in the cluster, similar to the findings in the study of 47 Tuc and $\omega$ Centauri. We can thus follow the method presented in Johnson et al. (2013) to estimate the total number
<figure><img src="image_9.png" /><figcaption>Figure 5. Orbital flux modulation for the six binary MSPs in M5. Upper limits at the 95% CL are shown as arrows for bins when the source had TS < 4. Blue dotted line: average flux in the broadband fit with PLEC4 model. Red dotted lines: $1\sigma$ uncertainties on the average flux.</figcaption></figure>

of MSPs in M5 as

$$N_{\text{MSP}} = \frac{L_{\gamma}}{\langle \dot{E} \rangle \langle \eta_{\gamma} \rangle}, \qquad (6)$$

where $\langle \dot{E} \rangle = (1.8 \pm 0.7) \times 10^{34}$ erg s$^{-1}$ is the average spin-down power in GCs (Abdo et al. 2009), $\langle \eta_{\gamma} \rangle$ is the average gamma-ray efficiency of MSPs, and $L_{\gamma} = (1.2 \pm 0.2 \pm 0.2) \times 10^{34}$ erg s$^{-1}$ is

the gamma-ray luminosity$^{38}$ of M5 based on the spectral fit reported in Table 2. For $\langle \eta_{\gamma} \rangle$, 3PC includes 20 MSPs with proper-motion and distance measurements yielding systematic uncertainties on $\eta_{\gamma}$, after Shlovskii corrections, smaller than 50%, similar to

<footnote>$^{38}$ The first uncertainty on $L_{\gamma}$ comes from the statistical uncertainty in the spectral fit, and the second one is systematic induced by the distance uncertainty.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 5. Orbital flux modulation for the six binary MSPs in M5. Upper limits at the 95% CL are shown as arrows for bins when the source had TS < 4. Blue dotted line: average flux in the broadband fit with PLEC4 model. Red dotted lines: $1\sigma$ uncertainties on the average flux.</figcaption></figure>

of MSPs in M5 as

$$N_{\text{MSP}} = \frac{L_{\gamma}}{\langle \dot{E} \rangle \langle \eta_{\gamma} \rangle}, \qquad (6)$$

where $\langle \dot{E} \rangle = (1.8 \pm 0.7) \times 10^{34} \text{ erg s}^{-1}$ is the average spin-down power in GCs (Abdo et al. 2009), $\langle \eta_{\gamma} \rangle$ is the average gamma-ray efficiency of MSPs, and $L_{\gamma} = (1.2 \pm 0.2 \pm 0.2) \times 10^{34} \text{ erg s}^{-1}$ is

the gamma-ray luminosity$^{38}$ of M5 based on the spectral fit reported in Table 2. For $\langle \eta_{\gamma} \rangle$, 3PC includes 20 MSPs with proper-motion and distance measurements yielding systematic uncertainties on $\eta_{\gamma}$, after Shlovskii corrections, smaller than 50%, similar to

<footnote>$^{38}$ The first uncertainty on $L_{\gamma}$ comes from the statistical uncertainty in the spectral fit, and the second one is systematic induced by the distance uncertainty.</footnote>
<table><caption>**Table 3**
Orbital-phase-resolved Spectral Fits for M5B</caption><thead><tr><th>Bin</th><th>TS</th><th>Γ</th><th>d</th><th>Photon Flux<br>(10<sup>-9</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Energy Flux<br>(10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0.5–0.6</td><td>42</td><td>2.08 ± 0.33</td><td>0.61 ± 0.38</td><td>7.70 ± 4.67</td><td>5.70 ± 1.70</td></tr><tr><td>0.6–0.5</td><td>61</td><td>1.72 ± 0.33</td><td>0.56 ± 0.28</td><td>1.50 ± 0.84</td><td>1.58 ± 0.39</td></tr></tbody></table>

the criteria in Johnson et al. (2013). Instead of $\langle\eta_{\gamma}\rangle$, we used the FWHM range $0.07 < \eta_{\gamma} < 0.4$ of the $\eta_{\gamma}$ distribution, as in Figure 24 of the 3PC, but with Shlovskii corrections for only the 20 well-characterized MSPs. This gives $1.7 < N_{\text{MSP}} < 9.5$. The current number of seven known MSPs in M5 is in this range. There may be a few more (but not a large number of) MSPs to discover in M5.

However, such estimates are to be taken with care, as we may see (or not) a pulsar just based on geometry. Measured efficiency can be severely affected by orientation. 3PC assumes a beaming factor $f_{\Omega}=1$, while recent theoretical studies of LAT pulsars found that in general $f_{\Omega} < 1$ is expected (Kalapotharakos et al. 2022), adding another uncertainty on the estimates.

Although it is generally believed that the gamma-ray emissions of GCs are from MSPs, there are two main distinct models. The pulsar magnetosphere model proposes that gamma rays are produced via curvature radiation of relativistic electrons/positrons in the pulsar magnetosphere. The gamma rays from GCs are thus expected to originate from the cumulative contribution of all MSPs in the cluster (e.g., Venter & de Jager 2008; Venter et al. 2009). The inverse Compton (IC) model, on the other hand, suggests that gamma rays are generated by the IC scattering between relativistic electrons/ positrons in the pulsar wind of MSPs in GCs and background soft photons (Bednarek & Sitarek 2007; Venter et al. 2009; Cheng et al. 2010), which will lead to intrinsic unpulsed emission of GCs.

Currently, both models can explain the GeV gamma-ray spectra of GCs equally well (e.g., Abdo et al. 2010; Cheng et al. 2010). The IC model, on the other hand, also predicts TeV gamma rays. However, observations with CANGAROO III, VERITAS, H.E.S.S., and MAGIC of GCs have not been successful (see, e.g., Kabuki et al. 2007; Aharonian et al. 2009; Anderhub et al. 2009; McCutcheon 2009; H. E. S. S. Collaboration 2013; MAGIC Collaboration 2019), with Terzan 5 being the only one to be claimed to shine in the TeV band (H. E. S. S. Collaboration 2011). Diffuse radio and X-ray emissions from GCs can also be produced by synchrotron radiation and IC scattering, as predicted by the IC model (Cheng et al. 2010). Observational support for such a scenario has been provided by the discovery of extended radio and X-ray emissions around Terzan 5 (Eger et al. 2010; Clapson et al. 2011) and 47 Tuc (Wu et al. 2014) with possibly nonthermal origin. See Tam et al. (2016) for a review of the observations and modelings of gamma-ray emission from GCs.

Song et al. (2021) claimed to have found evidence of a power-law high-energy tail in the gamma-ray spectra of GCs beyond the exponential cutoff power-law component by analyzing Fermi-LAT data for 157 MW GCs, which they interpreted in terms of the IC model, although more data are needed to assure their findings. They also claimed that the very soft high-energy tail is the reason behind the difficulty in detecting GCs with current TeV telescopes, and that it would be possible to detect GCs with more sensitive

TeV telescopes such as CTA (Actis et al. 2011) and LHAASO (Cao et al. 2019).

In conclusion, based on current observations, the most obvious possibility to explain the GeV gamma-ray emission of GCs is that it is produced by the collective emission of individual (but yet-undetected) gamma-ray pulsars in the cluster. However, further multiwavelength observations with more sensitive telescopes such as CTA (Actis et al. 2011) and LHAASO (Cao et al. 2019) in TeV, SKA (Tan et al. 2015) in radio, and EP in X-rays (Yuan et al. 2022c) are expected to provide better constraints on the two aforementioned emission models.

Acknowledgments

We thank the anonymous referee for the very useful comments and suggestions that helped to improve the paper. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work was performed in part under DOE contract DE-AC02-76SF00515.

The authors acknowledge support from the National Natural Science Foundation of China under grant Nos. 12373051, 12041303, 12273038, U2031117, and 11988101. X.H. is also supported by the National Key R&D Program of China (2023YFE0101200). D.F.T is supported by grant PID2021-1245810B-I00 funded by MCIN/AEI/10.13039/501100011033 and 2021SGR00426 of the Generalitat de Catalunya, by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and by MCIN with funding from European Union NextGeneration EU (PRTR-C17.11). W.Z. is supported by the National Scholarship Council (PhD fellowship from the China Scholarship Council (CSC; No. 202107030003)). P.W. is also supported by the Youth Innovation Promotion Association CAS (id. 2021055), the Cultivation Project for FAST Scientific Payoff and Research Achievement of CAMS-CAS, and the National SKA Program of China (No. 2020SKA0120200). Work at NRL is supported by NASA. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.
<table><caption>**Table 3**
Orbital-phase-resolved Spectral Fits for M5B</caption><thead><tr><th>Bin</th><th>TS</th><th>Γ</th><th>d</th><th>Photon Flux<br>(10<sup>-9</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Energy Flux<br>(10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0.5–0.6</td><td>42</td><td>2.08 ± 0.33</td><td>0.61 ± 0.38</td><td>7.70 ± 4.67</td><td>5.70 ± 1.70</td></tr><tr><td>0.6–0.5</td><td>61</td><td>1.72 ± 0.33</td><td>0.56 ± 0.28</td><td>1.50 ± 0.84</td><td>1.58 ± 0.39</td></tr></tbody></table>

the criteria in Johnson et al. (2013). Instead of $\langle\eta_{\gamma}\rangle$, we used the FWHM range $0.07 < \eta_{\gamma} < 0.4$ of the $\eta_{\gamma}$ distribution, as in Figure 24 of the 3PC, but with Shlovskii corrections for only the 20 well-characterized MSPs. This gives $1.7 < N_{\text{MSP}} < 9.5$. The current number of seven known MSPs in M5 is in this range. There may be a few more (but not a large number of) MSPs to discover in M5.

However, such estimates are to be taken with care, as we may see (or not) a pulsar just based on geometry. Measured efficiency can be severely affected by orientation. 3PC assumes a beaming factor $f_{\Omega}=1$, while recent theoretical studies of LAT pulsars found that in general $f_{\Omega} < 1$ is expected (Kalapotharakos et al. 2022), adding another uncertainty on the estimates.

Although it is generally believed that the gamma-ray emissions of GCs are from MSPs, there are two main distinct models. The pulsar magnetosphere model proposes that gamma rays are produced via curvature radiation of relativistic electrons/positrons in the pulsar magnetosphere. The gamma rays from GCs are thus expected to originate from the cumulative contribution of all MSPs in the cluster (e.g., Venter & de Jager 2008; Venter et al. 2009). The inverse Compton (IC) model, on the other hand, suggests that gamma rays are generated by the IC scattering between relativistic electrons/ positrons in the pulsar wind of MSPs in GCs and background soft photons (Bednarek & Sitarek 2007; Venter et al. 2009; Cheng et al. 2010), which will lead to intrinsic unpulsed emission of GCs.

Currently, both models can explain the GeV gamma-ray spectra of GCs equally well (e.g., Abdo et al. 2010; Cheng et al. 2010). The IC model, on the other hand, also predicts TeV gamma rays. However, observations with CANGAROO III, VERITAS, H.E.S.S., and MAGIC of GCs have not been successful (see, e.g., Kabuki et al. 2007; Aharonian et al. 2009; Anderhub et al. 2009; McCutcheon 2009; H. E. S. S. Collaboration 2013; MAGIC Collaboration 2019), with Terzan 5 being the only one to be claimed to shine in the TeV band (H. E. S. S. Collaboration 2011). Diffuse radio and X-ray emissions from GCs can also be produced by synchrotron radiation and IC scattering, as predicted by the IC model (Cheng et al. 2010). Observational support for such a scenario has been provided by the discovery of extended radio and X-ray emissions around Terzan 5 (Eger et al. 2010; Clapson et al. 2011) and 47 Tuc (Wu et al. 2014) with possibly nonthermal origin. See Tam et al. (2016) for a review of the observations and modelings of gamma-ray emission from GCs.

Song et al. (2021) claimed to have found evidence of a power-law high-energy tail in the gamma-ray spectra of GCs beyond the exponential cutoff power-law component by analyzing Fermi-LAT data for 157 MW GCs, which they interpreted in terms of the IC model, although more data are needed to assure their findings. They also claimed that the very soft high-energy tail is the reason behind the difficulty in detecting GCs with current TeV telescopes, and that it would be possible to detect GCs with more sensitive

TeV telescopes such as CTA (Actis et al. 2011) and LHAASO (Cao et al. 2019).

In conclusion, based on current observations, the most obvious possibility to explain the GeV gamma-ray emission of GCs is that it is produced by the collective emission of individual (but yet-undetected) gamma-ray pulsars in the cluster. However, further multiwavelength observations with more sensitive telescopes such as CTA (Actis et al. 2011) and LHAASO (Cao et al. 2019) in TeV, SKA (Tan et al. 2015) in radio, and EP in X-rays (Yuan et al. 2022c) are expected to provide better constraints on the two aforementioned emission models.

Acknowledgments

We thank the anonymous referee for the very useful comments and suggestions that helped to improve the paper. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work was performed in part under DOE contract DE-AC02-76SF00515.

The authors acknowledge support from the National Natural Science Foundation of China under grant Nos. 12373051, 12041303, 12273038, U2031117, and 11988101. X.H. is also supported by the National Key R&D Program of China (2023YFE0101200). D.F.T is supported by grant PID2021-124581OB-I00 funded by MCIN/AEI/10.13039/501100011033 and 2021SGR00426 of the Generalitat de Catalunya, by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and by MCIN with funding from European Union NextGeneration EU (PRTR-C17.11). W.Z. is supported by the National Scholarship Council (PhD fellowship from the China Scholarship Council (CSC; No. 202107030003)). P.W. is also supported by the Youth Innovation Promotion Association CAS (id. 2021055), the Cultivation Project for FAST Scientific Payoff and Research Achievement of CAMS-CAS, and the National SKA Program of China (No. 2020SKA0120200). Work at NRL is supported by NASA. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.
ORCID iDs

X. Hou https://orcid.org/0000-0003-0933-6101

W. Zhang https://orcid.org/0000-0003-2839-1325

P. C. C. Freire https://orcid.org/0000-0003-1307-9435

D. F. Torres https://orcid.org/0000-0002-1522-9065

J. Ballet https://orcid.org/0000-0002-8784-2977

D. A. Smith https://orcid.org/0000-0002-7833-0275

T. J. Johnson https://orcid.org/0000-0002-2771-472X

M. Kerr https://orcid.org/0000-0002-0893-4073

C. C. Cheung https://orcid.org/0000-0002-4377-0174

L. Guillemot https://orcid.org/0000-0002-9049-8716

J. Li https://orcid.org/0000-0003-1720-9727

L. Zhang https://orcid.org/0000-0001-8539-4237

A. Ridolfi https://orcid.org/0000-0001-6762-2638

P. Wang https://orcid.org/0000-0002-3386-7159

D. Li https://orcid.org/0000-0003-3010-7661

J. Yuan https://orcid.org/0000-0002-5381-6498

N. Wang https://orcid.org/0000-0002-9786-8548

References

Abbate, F., Ridolfi, A., Barr, E. D., et al. 2022, MNRAS, 513, 2292

Abbate, F., Ridolfi, A., Freire, P. C. C., et al. 2023, A&A, 680, A47

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Sci, 325, 845

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, A&A, 524, A75

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53

Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23

Ackermann, M., Ajello, M., Ballet, J., et al. 2012, Sci, 335, 189

Actis, M., Agnetta, G., Aharonian, F., et al. 2011, ExA, 32, 193

Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 499, 273

Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Natur, 300, 728

Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2009, ApJ, 702, 266

Anderson, S. B., Wolsczan, A., Kulkarni, S. R., & Prince, T. A. 1997, ApJ, 482, 870

Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Bahramian, A., Heinke, C. O., Sivakoff, G. R., & Gladstone, J. C. 2013, ApJ, 766, 136

Baldini, L., Ballet, J., Bastieri, D., et al. 2021, ApJS, 256, 13

Ballet, J., Bruel, P., Burnett, T. H., Lott, B. & The Fermi-LAT collaboration 2023, arXiv:2307.12546

Bednarek, W., & Sitarek, J. 2007, MNRAS, 377, 920

Bhattacharya, D., & van den Heuvel, E. P. J. 1991, PhR, 203, 1

Bruel, P. 2019, A&A, 622, A108

Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, arXiv:1810.11394

Cao, Z., della Volpe, D., Liu, S., et al. 2019, arXiv:1905.02773

Chen, K. 1991, Natur, 352, 695

Cheng, K. S., Chernyshov, D. O., Dogiel, V. A., Hui, C. Y., & Kong, A. K. H. 2010, ApJ, 723, 1219

Clapson, A. C., Domainko, W., Jamrozy, M., Dyrda, M., & Eger, P. 2011, A&A, 532, A47

Clark, G. W. 1975, ApJL, 199, L143

D'Abrusco, R., Álvarez Crespo, N., Massaro, F., et al. 2019, ApJS, 242, 4

Dai, S., Johnston, S., Kerr, M., et al. 2023, MNRAS, 521, 2616

de Menezes, R., Cafardo, F., & Nemmen, R. 2019, MNRAS, 486, 851

de Menezes, R., Di Pierro, F., & Chiavassa, A. 2023, MNRAS, 523, 4455

Eger, P., Domainko, W., & Clapson, A. C. 2010, A&A, 513, A66

Feng, L., Cheng, Z., Wang, W., Li, Z., & Chen, Y. 2024, RAA, 24, 025001

Freire, P. C. C., Abdo, A. A., Ajello, M., et al. 2011, Sci, 334, 1107

Freire, P. C. C., Ridolfi, A., Kramer, M., et al. 2017, MNRAS, 471, 857

Gautam, T., Ridolfi, A., Freire, P. C. C., et al. 2022, A&A, 664, A54

Gendre, B., Barret, D., & Webb, N. 2003, A&A, 403, L11

Harris, W. E. 1996, AJ, 112, 1487

Healey, S. E., Romani, R. W., Cotter, G., et al. 2008, ApJS, 175, 97

Healey, S. E., Romani, R. W., Taylor, G. B., et al. 2007, ApJS, 171, 61

Heggie, D. C. 1975, MNRAS, 173, 729

H. E. S. S. Collaboration 2011, A&A, 531, L18

H. E. S. S. Collaboration 2013, A&A, 551, A26

Hessels, J. W. T., Ransom, S. M., Stairs, I. H., Kaspi, V. M., & Freire, P. C. C. 2007, ApJ, 670, 363

Hills, J. G. 1975, AJ, 80, 809

Hooper, D., & Linden, T. 2016, JCAP, 2016, 018

Hui, C. Y., Cheng, K. S., & Taam, R. E. 2010, ApJ, 714, 1149

Hui, C. Y., Cheng, K. S., Wang, Y., et al. 2011, ApJ, 726, 100

Johnson, T. J., Guillemot, L., Kerr, M., et al. 2013, ApJ, 778, 106

Johnson, T. J., Ray, P. S., Roy, J., et al. 2015, ApJ, 806, 91

Kabuki, S., Enomoto, R., Bicknell, G. V., et al. 2007, ApJ, 668, 968

Kalapotharakos, C., Wadiasingh, Z., Harding, A. K., & Kazanas, D. 2022, ApJ, 934, 65

Katz, J. I. 1975, Natur, 253, 698

Kerr, M. 2011, ApJ, 732, 38

Kong, A. K. H., Hui, C. Y., & Cheng, K. S. 2010, ApJL, 712, L36

Kovačević, M., Chiaro, G., Cutini, S., & Tosti, G. 2019, MNRAS, 490, 4770

Lloyd, S. J., Chadwick, P. M., & Brown, A. M. 2018, MNRAS, 480, 4782

MAGIC Collaboration 2019, MNRAS, 484, 2876

Massaro, E., Maselli, A., Leto, C., et al. 2015, Ap&SS, 357, 75

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

McCutcheon, M. 2009, arXiv:0907.4974

Mott, A. J., & Freire, P. C. 2003, AAS Meeting 203, 53.07

Pan, Z., Qian, L., Ma, X., et al. 2021, ApJL, 915, L28

Pooley, D., Lewin, W. H. G., Anderson, S. F., et al. 2003, ApJL, 591, L131

Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17

Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191

Smith, D. A., Bruel, P., Cognard, I., et al. 2019, ApJ, 871, 78

Sollima, A., & Baumgardt, H. 2017, MNRAS, 471, 3668

Song, D., Macias, O., Horiuchi, S., Crocker, R. M., & Nataf, D. M. 2021, MNRAS, 507, 5161

Tam, P.-H. T., Hui, C. Y., & Kong, A. K. H. 2016, JASS, 33, 1

Tam, P. H. T., Kong, A. K. H., Hui, C. Y., et al. 2011, ApJ, 729, 90

Tan, G. H., Cornwell, T. J., Dewdney, P. E., & Waterson, M. 2015, in 2015 1st URSI Atlantic Radio Science Conf. (URSI AT-RASC) (Piscataway, NJ: IEEE),

Venter, C., & de Jager, O. C. 2008, ApJL, 680, L125

Venter, C., De Jager, O. C., & Clapson, A. C. 2009, ApJL, 696, L52

Verbunt, F., & Freire, P. C. C. 2014, A&A, 561, A11

Wood, M., Caputo, R., Charles, E., et al. 2017, Proc. 35th Int. Cosmic Ray Conference (ICRC2017) (Trieste: SISSA), 824

Wu, E. M. H., Hui, C. Y., Kong, A. K. H., et al. 2014, ApJL, 788, L40

Wu, J. H. K., Hui, C. Y., Wu, E. M. H., et al. 2013, ApJL, 765, L47

Yuan, M., Ren, C., Zhang, P., Jiang, Z., & Zhang, L. 2022a, RAA, 22, 115013

Yuan, M., Zheng, J., Zhang, P., & Zhang, L. 2022b, RAA, 22, 055019

Yuan, W., Zhang, C., Chen, Y., & Ling, Z. 2022c, in Handbook of X-Ray and Gamma-Ray Astrophysics, ed. C. Bambi & A. Santangelo (New York: Springer), 86

Zhang, L., Freire, P. C. C., Ridolfi, A., et al. 2023a, ApJS, 269, 56

Zhang, P., Xing, Y., & Wang, Z. 2022, ApJL, 935, L36

Zhang, P., Xing, Y., Wang, Z., Wu, W., & Chen, Z. 2023b, ApJ, 945, 70

Zhang, P. F., Xin, Y. L., Fu, L., et al. 2016, MNRAS, 459, 99

Zhou, J. N., Zhang, P. F., Huang, X. Y., et al. 2015, MNRAS, 448, 3215
ORCID iDs

X. Hou https://orcid.org/0000-0003-0933-6101

W. Zhang https://orcid.org/0000-0003-2839-1325

P. C. C. Freire https://orcid.org/0000-0003-1307-9435

D. F. Torres https://orcid.org/0000-0002-1522-9065

J. Ballet https://orcid.org/0000-0002-8784-2977

D. A. Smith https://orcid.org/0000-0002-7833-0275

T. J. Johnson https://orcid.org/0000-0002-2771-472X

M. Kerr https://orcid.org/0000-0002-0893-4073

C. C. Cheung https://orcid.org/0000-0002-4377-0174

L. Guillemot https://orcid.org/0000-0002-9049-8716

J. Li https://orcid.org/0000-0003-1720-9727

L. Zhang https://orcid.org/0000-0001-8539-4237

A. Ridolfi https://orcid.org/0000-0001-6762-2638

P. Wang https://orcid.org/0000-0002-3386-7159

D. Li https://orcid.org/0000-0003-3010-7661

J. Yuan https://orcid.org/0000-0002-5381-6498

N. Wang https://orcid.org/0000-0002-9786-8548

References

Abbate, F., Ridolfi, A., Barr, E. D., et al. 2022, MNRAS, 513, 2292

Abbate, F., Ridolfi, A., Freire, P. C. C., et al. 2023, A&A, 680, A47

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Sci, 325, 845

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, A&A, 524, A75

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53

Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23

Ackermann, M., Ajello, M., Ballet, J., et al. 2012, Sci, 335, 189

Actis, M., Agnetta, G., Aharonian, F., et al. 2011, ExA, 32, 193

Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 499, 273

Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Natur, 300, 728

Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2009, ApJ, 702, 266

Anderson, S. B., Wolsczan, A., Kulkarni, S. R., & Prince, T. A. 1997, ApJ, 482, 870

Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Bahramian, A., Heinke, C. O., Sivakoff, G. R., & Gladstone, J. C. 2013, ApJ, 766, 136

Baldini, L., Ballet, J., Bastieri, D., et al. 2021, ApJS, 256, 13

Ballet, J., Bruel, P., Burnett, T. H., Lott, B. & The Fermi-LAT collaboration 2023, arXiv:2307.12546

Bednarek, W., & Sitarek, J. 2007, MNRAS, 377, 920

Bhattacharya, D., & van den Heuvel, E. P. J. 1991, PhR, 203, 1

Bruel, P. 2019, A&A, 622, A108

Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, arXiv:1810.11394

Cao, Z., della Volpe, D., Liu, S., et al. 2019, arXiv:1905.02773

Chen, K. 1991, Natur, 352, 695

Cheng, K. S., Chernyshov, D. O., Dogiel, V. A., Hui, C. Y., & Kong, A. K. H. 2010, ApJ, 723, 1219

Clapson, A. C., Domainko, W., Jamrozy, M., Dyrda, M., & Eger, P. 2011, A&A, 532, A47

Clark, G. W. 1975, ApJL, 199, L143

D'Abrusco, R., Álvarez Crespo, N., Massaro, F., et al. 2019, ApJS, 242, 4

Dai, S., Johnston, S., Kerr, M., et al. 2023, MNRAS, 521, 2616

de Menezes, R., Cafardo, F., & Nemmen, R. 2019, MNRAS, 486, 851

de Menezes, R., Di Pierro, F., & Chiavassa, A. 2023, MNRAS, 523, 4455

Eger, P., Domainko, W., & Clapson, A. C. 2010, A&A, 513, A66

Feng, L., Cheng, Z., Wang, W., Li, Z., & Chen, Y. 2024, RAA, 24, 025001

Freire, P. C. C., Abdo, A. A., Ajello, M., et al. 2011, Sci, 334, 1107

Freire, P. C. C., Ridolfi, A., Kramer, M., et al. 2017, MNRAS, 471, 857

Gautam, T., Ridolfi, A., Freire, P. C. C., et al. 2022, A&A, 664, A54

Gendre, B., Barret, D., & Webb, N. 2003, A&A, 403, L11

Harris, W. E. 1996, AJ, 112, 1487

Healey, S. E., Romani, R. W., Cotter, G., et al. 2008, ApJS, 175, 97

Healey, S. E., Romani, R. W., Taylor, G. B., et al. 2007, ApJS, 171, 61

Heggie, D. C. 1975, MNRAS, 173, 729

H. E. S. S. Collaboration 2011, A&A, 531, L18

H. E. S. S. Collaboration 2013, A&A, 551, A26

Hessels, J. W. T., Ransom, S. M., Stairs, I. H., Kaspi, V. M., & Freire, P. C. C. 2007, ApJ, 670, 363

Hills, J. G. 1975, AJ, 80, 809

Hooper, D., & Linden, T. 2016, JCAP, 2016, 018

Hui, C. Y., Cheng, K. S., & Taam, R. E. 2010, ApJ, 714, 1149

Hui, C. Y., Cheng, K. S., Wang, Y., et al. 2011, ApJ, 726, 100

Johnson, T. J., Guillemot, L., Kerr, M., et al. 2013, ApJ, 778, 106

Johnson, T. J., Ray, P. S., Roy, J., et al. 2015, ApJ, 806, 91

Kabuki, S., Enomoto, R., Bicknell, G. V., et al. 2007, ApJ, 668, 968

Kalapotharakos, C., Wadiasingh, Z., Harding, A. K., & Kazanas, D. 2022, ApJ, 934, 65

Katz, J. I. 1975, Natur, 253, 698

Kerr, M. 2011, ApJ, 732, 38

Kong, A. K. H., Hui, C. Y., & Cheng, K. S. 2010, ApJL, 712, L36

Kovačević, M., Chiaro, G., Cutini, S., & Tosti, G. 2019, MNRAS, 490, 4770

Lloyd, S. J., Chadwick, P. M., & Brown, A. M. 2018, MNRAS, 480, 4782

MAGIC Collaboration 2019, MNRAS, 484, 2876

Massaro, E., Maselli, A., Leto, C., et al. 2015, Ap&SS, 357, 75

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

McCutcheon, M. 2009, arXiv:0907.4974

Mott, A. J., & Freire, P. C. 2003, AAS Meeting 203, 53.07

Pan, Z., Qian, L., Ma, X., et al. 2021, ApJL, 915, L28

Pooley, D., Lewin, W. H. G., Anderson, S. F., et al. 2003, ApJL, 591, L131

Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17

Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191

Smith, D. A., Bruel, P., Cognard, I., et al. 2019, ApJ, 871, 78

Sollima, A., & Baumgardt, H. 2017, MNRAS, 471, 3668

Song, D., Macias, O., Horiuchi, S., Crocker, R. M., & Nataf, D. M. 2021, MNRAS, 507, 5161

Tam, P.-H. T., Hui, C. Y., & Kong, A. K. H. 2016, JASS, 33, 1

Tam, P. H. T., Kong, A. K. H., Hui, C. Y., et al. 2011, ApJ, 729, 90

Tan, G. H., Cornwell, T. J., Dewdney, P. E., & Waterson, M. 2015, in 2015 1st URSI Atlantic Radio Science Conf. (URSI AT-RASC) (Piscataway, NJ: IEEE),

Venter, C., & de Jager, O. C. 2008, ApJL, 680, L125

Venter, C., De Jager, O. C., & Clapson, A. C. 2009, ApJL, 696, L52

Verbunt, F., & Freire, P. C. C. 2014, A&A, 561, A11

Wood, M., Caputo, R., Charles, E., et al. 2017, Proc. 35th Int. Cosmic Ray Conference (ICRC2017) (Trieste: SISSA), 824

Wu, E. M. H., Hui, C. Y., Kong, A. K. H., et al. 2014, ApJL, 788, L40

Wu, J. H. K., Hui, C. Y., Wu, E. M. H., et al. 2013, ApJL, 765, L47

Yuan, M., Ren, C., Zhang, P., Jiang, Z., & Zhang, L. 2022a, RAA, 22, 115013

Yuan, M., Zheng, J., Zhang, P., & Zhang, L. 2022b, RAA, 22, 055019

Yuan, W., Zhang, C., Chen, Y., & Ling, Z. 2022c, in Handbook of X-Ray and Gamma-Ray Astrophysics, ed. C. Bambi & A. Santangelo (New York: Springer), 86

Zhang, L., Freire, P. C. C., Ridolfi, A., et al. 2023a, ApJS, 269, 56

Zhang, P., Xing, Y., & Wang, Z. 2022, ApJL, 935, L36

Zhang, P., Xing, Y., Wang, Z., Wu, W., & Chen, Z. 2023b, ApJ, 945, 70

Zhang, P. F., Xin, Y. L., Fu, L., et al. 2016, MNRAS, 459, 99

Zhou, J. N., Zhang, P. F., Huang, X. Y., et al. 2015, MNRAS, 448, 3215