Evidence for inverse Compton emission from globular clusters
Deheng Song, Oscar Macias, Shunsaku Horiuchi, Roland M. Crocker, David M. Nataf
MMNRAS , 1–19 (2021) Preprint 2 February 2021 Compiled using MNRAS L A TEX style file v3.0
Evidence for inverse Compton emission from globular clusters
Deheng Song, ★ Oscar Macias, , , † Shunsaku Horiuchi, ‡ Roland M. Crocker andDavid M. Nataf Center for Neutrino Physics, Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA Kavli Institute for the Physics and Mathematics of the Universe (WPI),University of Tokyo, Kashiwa, Chiba 277-8583, Japan GRAPPA Institute, Institute of Physics, University of Amsterdam, 1098 XH Amsterdam, The Netherlands Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia Center for Astrophysical Sciences and Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Millisecond pulsars are very likely the main source of gamma ray emission from globularclusters. However, the relative contributions of two separate emission processes–curvatureradiation from millisecond pulsar magnetospheres vs. inverse Compton emission from rela-tivistic pairs launched into the globular cluster environment by millisecond pulsars–has longbeen unclear. To address this, we search for evidence of inverse Compton emission in 8-year
Fermi -LAT data from the directions of 157 Milky Way globular clusters. We find a mildlystatistically significant (3.8 𝜎 ) correlation between the measured globular cluster gamma-rayluminosities and their photon field energy densities. However, this may also be explained by ahidden correlation between the photon field densities and the stellar encounter rates of globularclusters. Analysed in toto , we demonstrate that the 𝛾 -ray emission of globular clusters can beresolved spectrally into two components: i) an exponentially cut-off power law and ii) a purepower law. The latter component–which we uncover at a significance of 8.2 𝜎 –is most naturallyinterpreted as inverse Compton emission by cosmic-ray electrons and positrons injected bymillisecond pulsars. We find the luminosity of this inverse Compton component is comparableto, or slightly smaller than, the luminosity of the curved component, suggesting the fractionof millisecond pulsar spin-down luminosity into relativistic leptons is similar to the fractionof the spin-down luminosity into prompt magnetospheric radiation. Key words: gamma-rays:general – globular clusters: general – pulsars:general
Over two dozen globular clusters (GCs) have been detected in 𝛾 raysin Fermi
Large Area Telescope (LAT) data (Abdo et al. 2009a, 2010;Kong et al. 2010; Tam et al. 2011; Zhou et al. 2015; Zhang et al.2016). The millisecond pulsar (MSP) populations of those GCs arebelieved to be the main source of this 𝛾 -ray emission. In particular,MSPs have been firmly established as 𝛾 -ray sources (Verbunt et al.1996; Abdo et al. 2009b; Espinoza et al. 2013; Abdo et al. 2013)and a large fraction of them have been discovered in GCs (Camilo& Rasio 2005). Recently, Fermi -LAT detected 𝛾 -ray pulsations intwo GCs (Freire et al. 2011; Johnson et al. 2013), providing furthersupport for this scenario.The high-energy emission from MSPs emerges from the pri-mary electrons accelerated by them and subsequent radiation by ★ [email protected] † [email protected] ‡ [email protected] secondary, relativistic electrons and positrons ( 𝑒 ± ) pair created intheir magnetospheres. In particular, Harding et al. (2005) studied thecurvature radiation (CR) of primary electrons within MSP magne-tospheres with a focus on GeV-scale emission. Bednarek & Sitarek(2007) then considered a scenario in which 𝑒 ± , injected by MSPs,gradually diffuse through a GC, up-scattering ambient photons,and thus producing inverse Compton (IC) 𝛾 -ray emission in theGeV − TeV energy range. Venter et al. (2009) calculated the CR andIC spectra for an ensemble of MSPs in the GCs 47 Tucanae andTerzan 5. Cheng et al. (2010) found that the spectra of 47 Tucanaeand seven other GCs can be explained by IC alone, invoking back-ground photons from the cosmic microwave background (CMB)or Galactic infrared/stellar radiation. In general, the GeV emissionmechanism of MSPs remains in contention with CR, IC, and syn-chrotron radiation all proposed (Harding 2021).Here, motivated by the increasing number of GCs detectedin 𝛾 rays, we perform a collective statistical study of their prop-erties in order to gain insight into the nature of their high-energy © a r X i v : . [ a s t r o - ph . H E ] J a n Song et al. emission. Our particular aim is to investigate the importance of thecontribution of IC emission to the overall 𝛾 -ray emission of GCs.Relations between the detected GC gamma-ray luminositiesand properties of GC can be used to probe the origins of 𝛾 -rayemission and their underlying sources. For example, correlationswith the photon field energy density of GCs could unveil the poten-tial contribution from IC, and correlations with the stellar encounterrate and metallicity could provide insight into the dynamical for-mation of MPSs in GCs. Previous work here includes a study byAbdo et al. (2010) that reported a correlation between the 𝛾 -rayluminosity 𝐿 𝛾 and the stellar encounter rate of eight GCs. Hui et al.(2011) studied a group of 15 𝛾 -ray emitting GCs with 2 years ofFermi data and found a positive correlation between 𝐿 𝛾 and, re-spectively, encounter rate, metallicity, and Galactic photon energydensity. Hooper & Linden (2016) studied 25 GCs using 85 monthsof Fermi data, and found that the 𝛾 -ray luminosity function of MSPsin GCs is consistent with that applying to MSPs detected in the field.Lloyd et al. (2018) studied 𝛾 -ray emission from high-latitude GCsand its connection to their X-ray emission. de Menezes et al. (2019)reanalysed 9 years of Fermi data and found 23 𝛾 -ray emitting GCs;they found that the metallicity only mildly contributes to 𝐿 𝛾 whilea very high encounter rate seemed to reduce the 𝐿 𝛾 from GCs.In parallel, modeling of GCs’ observed broadband spectralenergy distributions provides a handle on their CR and IC emissions.Recently, Kopp et al. (2013) and Ndiyavala et al. (2019) modelledthe multiwavelength emission from MSPs considering a potentialCR origin for GeV and IC emissions for TeV 𝛾 rays, as well assynchrotron radiation for the radio and X-ray wavebands. Thesemodels are successful in explaining the multiwavelength spectraof Terzan 5. However, Terzan 5 is the only GC (perhaps) detectedabove TeV energies (H. E. S. S. Collaboration et al. 2011). Detailedspectral modelling similar to that presented by Kopp et al. (2013)and Ndiyavala et al. (2019) is difficult for other GCs at present dueto a lack of TeV 𝛾 -ray data. Although Fermi -LAT is sensitive to 𝛾 -rays of up to ≈ Fermi -LAT fourth source cata-log (Abdollahi et al. 2020) (4FGL) , 30 GCs have been detected inGeV 𝛾 rays. With such a number, we can begin to carefully study thenature of the 𝛾 -ray emission from GCs through a population study.In this paper, we repeat the bin-by-bin analysis of the 4FGL data forthe 157 known Milky Way GCs in the Harris (1996) catalog . Wesearch for correlations between the 𝛾 -ray luminosity of the GCs andother parameters of the GCs to probe which are good proxies for the 𝛾 -ray luminosity and study potential IC contributions; to this end,we consider the photon field energy densities, the stellar encounterrate, and the metallicity of the GCs. Unlike previous studies of cor-relations of the GC 𝛾 -ray emissions, we consider also the upperlimits placed by null detections, which we implement via a Kendall 𝜏 coefficient test statistic. Furthermore, we also look for evidence ofIC from the spectra of GCs. For the first time, we implement a uni-versal two-component model to study the spectra of 𝛾 -ray-detectedGCs. The two-component model comprises a CR component, which The Fermi collaboration has recently released an incremental update ofthe fourth source catalog (Ballet et al. 2020) (4FGL-DR2, for Data Release2). The new catalog uses 10 years of data, a 25% increase with respect tothe 4FGL. However, only 1 new GC (NGC 362) has been detected. Giventhis marginal change, we retain use of the 4FGL catalog constructed withthe 8-year data set. is spectrally curved, plus an IC component modeled as a power lawin the energy range of interest.The remainder of the paper is as follows: In Section 2 we dis-cuss the choice of GC samples and the data analysis procedure.Section 3 presents the methodology and results of our correlationanalysis. Section 4 describes the spectral analysis method and re-ports the 𝑒 ± injection efficiency in the GCs. We discuss the impli-cations of our results in Section 5 and conclude in Section 6. 𝛾 -RAY EMISSION FROM GLOBULAR CLUSTERS In this section, we describe our choice of GC sample and the GCs’ 𝛾 -ray-related parameters. The Fermi data analysis process is reportedas well. For GCs with a 𝛾 -ray counterpart in the 4FGL, we updatetheir spectral parameters through a maximum likelihood procedure.For those not detected in the 4FGL, we estimate their 95% C.L. 𝛾 -ray upper limits. We consider the Harris (1996) catalog (2010 edition), which con-tains identifications and basic parameters for 157 GCs in the MilkyWay. Here, we reanalyse publicly available
Fermi -LAT data fromthe direction of all GCs in the Harris (1996) catalog. Figure 1 showsthe spatial distribution of the GCs. The top panel shows the pro-jected direction of the GCs on the celestial plane while the bottomtwo panels display their 3D coordinates. The GCs which are de-tected in the 4FGL are marked by red stars while null detections areindicated by green circles. Most 𝛾 -ray GCs are near the Sun (yellowcircle) or located in the Galactic bulge (assumed to be sphere of 3kpc radius, grey circular area).The origin of the 𝛾 -ray emission from GCs can be studied bycomparing its dependency on GC properties. IC emission is sensi-tive to the ambient photon field on which the 𝑒 ± scatter. Hui et al.(2011) reported a positive correlation between the 𝛾 -ray luminosity 𝐿 𝛾 and the photon field energy density at the cluster location, indi-cating an IC contribution. In the present work, we improve upon theGalactic radiation field model used by Hui et al. (2011) by extract-ing the energy density of the interstellar radiation at the locations ofthe GCs from the three-dimensional interstellar radiation model in GALPROP v56 (Porter et al. 2017; Jóhannesson et al. 2018). Thisis a fully 3D model that combines the CMB, infrared, and opticalphotons of the Milky Way, denoted as 𝑢 MW . In addition, photonsfrom stars in the GCs are expected to make a dominant contributionto the total, ambient radiation field. We estimate this component by 𝑢 GC = 𝐿 ∗ 𝜋𝑐𝑟 ℎ , (1)where 𝐿 ∗ and 𝑟 ℎ are the stellar luminosity and the half-light radius ofthe GC. The total photon field energy density is 𝑢 Total = 𝑢 MW + 𝑢 GC .A potential correlation between the 𝛾 -ray luminosity 𝐿 𝛾 andthe stellar encounter rate has been studied as a way to probe thedynamic formation of MSPs in GCs (Abdo et al. 2010; Hui et al.2011; de Menezes et al. 2019). In the present work, we adopt thestellar encounter rate estimated by Bahramian et al. (2013), whichis defined as Γ 𝑐 = 𝜋𝜎 𝑐 ∫ 𝜌 ( 𝑟 ) 𝑟 𝑑𝑟, (2) http://galprop.stanford.edu/ MNRAS000
Fermi -LAT data fromthe direction of all GCs in the Harris (1996) catalog. Figure 1 showsthe spatial distribution of the GCs. The top panel shows the pro-jected direction of the GCs on the celestial plane while the bottomtwo panels display their 3D coordinates. The GCs which are de-tected in the 4FGL are marked by red stars while null detections areindicated by green circles. Most 𝛾 -ray GCs are near the Sun (yellowcircle) or located in the Galactic bulge (assumed to be sphere of 3kpc radius, grey circular area).The origin of the 𝛾 -ray emission from GCs can be studied bycomparing its dependency on GC properties. IC emission is sensi-tive to the ambient photon field on which the 𝑒 ± scatter. Hui et al.(2011) reported a positive correlation between the 𝛾 -ray luminosity 𝐿 𝛾 and the photon field energy density at the cluster location, indi-cating an IC contribution. In the present work, we improve upon theGalactic radiation field model used by Hui et al. (2011) by extract-ing the energy density of the interstellar radiation at the locations ofthe GCs from the three-dimensional interstellar radiation model in GALPROP v56 (Porter et al. 2017; Jóhannesson et al. 2018). Thisis a fully 3D model that combines the CMB, infrared, and opticalphotons of the Milky Way, denoted as 𝑢 MW . In addition, photonsfrom stars in the GCs are expected to make a dominant contributionto the total, ambient radiation field. We estimate this component by 𝑢 GC = 𝐿 ∗ 𝜋𝑐𝑟 ℎ , (1)where 𝐿 ∗ and 𝑟 ℎ are the stellar luminosity and the half-light radius ofthe GC. The total photon field energy density is 𝑢 Total = 𝑢 MW + 𝑢 GC .A potential correlation between the 𝛾 -ray luminosity 𝐿 𝛾 andthe stellar encounter rate has been studied as a way to probe thedynamic formation of MSPs in GCs (Abdo et al. 2010; Hui et al.2011; de Menezes et al. 2019). In the present work, we adopt thestellar encounter rate estimated by Bahramian et al. (2013), whichis defined as Γ 𝑐 = 𝜋𝜎 𝑐 ∫ 𝜌 ( 𝑟 ) 𝑟 𝑑𝑟, (2) http://galprop.stanford.edu/ MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters -60°-30°0°+30° +60° -120°-60°0°+60°+120°
20 15 10 5 0 5 10 15 20505 z [ k p c ] Sun20 15 10 5 0 5 10 15 20x [kpc]15105051015 y [ k p c ] SunGCs not detected in raysGCs detected in the 4FGL catalog
Figure 1.
Spatial distribution of the 157 Milky Way GCs in the Harris (1996)catalog. The top panel shows the all sky spatial distribution in Galacticcoordinate. The middle and bottom panel display the three-dimensional(3D) Cartesian coordinates of the GCs, in which the Sun (yellow circle) islocated in the negative x-axis. The 𝛾 -ray-detected GCs in the 4FGL catalogare shown as red stars, and the GCs not detected in 𝛾 rays are shown asgreen circles. Most 𝛾 -ray-detected GCs are located near the Sun or in theGalactic bulge (grey shaded area in the middle and bottom panel). where 𝜎 𝑐 is the velocity dispersion at the core radius, 𝜌 ( 𝑟 ) is thestellar density profile of the cluster, and the line-of-sight integrationis performed along the half-light radius.Additionally, it has been argued (Hui et al. 2011; de Menezeset al. 2019) that high metallicity [Fe/H] in the GCs could enhance thedynamical formation of MSP. The outer convective zone of metal-rich stars enables magnetic braking, which assists orbital shrinkingduring binary formation. In this analysis, we use the GCs metal-licities reported in the Harris (1996) catalog, which summarizesspectroscopic or photometric estimates in the literature.In summary, we consider the empirical dependence of the inferred 𝛾 -ray luminosity of GCs on four parameters, namely, 𝑢 MW , 𝑢 Total , Γ 𝑐 , and [Fe/H]. We summarize the values of these parametersfor the 30 𝛾 -ray-detected GCs in Table 1. Also included are the stellarmasses M ∗ of the GCs adopted from Baumgardt (2017), Sollima& Baumgardt (2017), and Baumgardt & Hilker (2018). They areestimated from N-body modelling of the velocity dispersion andsurface density profiles. In Section 3, we study the correlationsbetween the 𝛾 -ray emission of GCs and these parameters. We use 8 years of
Fermi -LAT data, from 2008 August 4 to2016 August 2. This constitutes the same data as the 4FGL.The newest Pass 8 data release is applied. As recommended bythe
Fermi -LAT data analysis documentation , the event class forthe analysis is "P8 Source" class (evclass=128) and the eventtype is "FRONT+BACK" (evtype=3). We use a 90 ◦ zenith an-gle cut to remove Earth limb events and filter the data by(DATA_QUAL>0)&&(LAT_CONFIG==1). The corresponding in-strument response function is P8R3_SOURCE_V2 . For our analysis,the
Fermipy software version is used, together with the
FermiScience Tools version .For 30 GCs detected by the 4FGL, we simply reanalyse the4FGL 𝛾 -ray source. We use a 10 ◦ by 10 ◦ Region-of-Interest aroundthe source with a 0.1 ◦ by 0.1 ◦ bin size. Photons from 300 MeVto 500 GeV are analysed using 9 logarithmic bins. Given that weuse a different Region-of-Interest size and photon class comparedto those adopted in the construction of the 4FGL, additional pointsources might emerge in our Region-of-Interest. However, sincewe use the same observation time as in the 4FGL, the impact ofthose potential new sources is expected to be minimal. Therefore,we only include known 4FGL sources in our analysis. As rec-ommended by the Fermi team, we re-run a maximum likelihoodanalysis that starts from the best-fit parameter values found in the4FGL and updating accordingly. The most recent Galactic inter-stellar emission model gll_iem_v07 and the isotropic component iso_P8R3_SOURCE_V2_v1 are employed as fore/backgrounds withfree-floating normalization. We have followed the Fermipy recom-mended procedure and fixed the spectral parameters of the sourceswith TS < 10 and 10 < Npred < 100 to their 4FGL values. However,the spectral parameters of the 4FGL sources lying within 5 ◦ of theRegion-of-Interest center are allowed to float freely. The MINUIT algorithm is used to determine the best-fit parameters of the sourcesfor each energy bin independently.For the 127 additional GCs in the Harris (1996) catalog without4FGL detections, we estimate the 95% C.L. 𝛾 -ray upper limits fromtheir locations. More specifically, we place a point source at thecoordinates of those GCs. The point source is assumed to havea power-law spectrum 𝑑𝑁 / 𝑑𝐸 ∼ 𝐸 − Γ with fixed index Γ = 𝐿 𝛾 for 30 𝛾 -emitting GCs. For eachGC, the total photon flux is summed over the bin-by-bin fluxes fromthe Fermi analysis. The statistical error of the total flux is addedquadratically from the bin-by-bin flux errors. The energy flux isestimated similarly, then 𝐿 𝛾 = 𝜋𝑅 (cid:12) × ( energy flux ) . We ignore http://fermi.gsfc.nasa.gov/ssc/data/analysis/ http://fermipy.readthedocs.io/en/latest/quickstart.html MNRAS , 1–19 (2021)
Song et al.
Table 1.
Parameters and data analysis results for 30 𝛾 -ray-detected GCs.Name Γ 𝑐 a [Fe/H] b M ∗ c 𝑢 MWd 𝑢 Totale 𝑅 (cid:12) f Flux g 𝐿 𝛾 g TS h (10 𝑀 (cid:12) ) (eV cm − ) (eV cm − ) (kpc) (10 − ph cm − s − ) (10 erg s − )2MS-GC01 ... ... ... 1.79 7.14 3.60 1 . ± .
26 3 . ± .
81 153.82GLIMPSE01 ... ... ... 1.55 30.23 4.20 2 . ± .
28 8 . ± .
94 535.61GLIMPSE02 ... -0.33 ... 2.61 >2.61 5.50 2 . ± .
25 11 . ± .
57 318.41M 62 2470.00 -1.18 6.76 2.14 293.14 6.80 0 . ± .
09 9 . ± .
89 1012.19M 80 937.00 -1.75 2.82 0.92 276.86 10.00 0 . ± .
07 4 . ± .
39 94.83NGC 104 1000.00 -0.72 8.13 0.55 31.12 4.50 1 . ± .
07 5 . ± .
34 4853.63NGC 1904 126.00 -1.60 1.66 0.29 173.13 12.90 0 . ± .
04 2 . ± .
98 23.84NGC 2808 1210.00 -1.14 8.13 0.38 467.35 9.60 0 . ± .
06 3 . ± .
03 90.30NGC 5904 120.00 -1.29 3.63 0.55 56.46 7.50 0 . ± .
04 1 . ± .
47 39.07NGC 6139 407.00 -1.65 3.47 1.15 161.34 10.10 0 . ± .
09 5 . ± .
19 59.29NGC 6218 18.10 -1.37 0.83 0.90 14.94 4.80 0 . ± .
05 0 . ± .
20 33.92NGC 6304 150.00 -0.45 1.62 2.33 23.95 5.90 0 . ± .
03 1 . ± .
42 21.71NGC 6316 131.00 -0.45 3.63 1.88 270.82 10.40 0 . ± .
11 10 . ± .
13 207.99NGC 6341 265.00 -2.31 3.09 0.42 97.31 8.30 0 . ± .
04 0 . ± .
37 15.84NGC 6388 1770.00 -0.55 10.47 1.30 1127.11 9.90 0 . ± .
09 18 . ± .
63 970.86NGC 6397 146.00 -2.02 0.89 0.92 3.75 2.30 0 . ± .
06 0 . ± .
05 17.21NGC 6402 106.00 -1.28 7.41 0.87 136.26 9.30 0 . ± .
08 3 . ± .
18 51.16NGC 6440 1750.00 -0.36 5.01 2.50 721.93 8.50 0 . ± .
13 10 . ± .
97 259.55NGC 6441 3150.00 -0.46 11.75 1.36 1148.78 11.60 0 . ± .
10 18 . ± .
62 363.50NGC 6528 233.00 -0.11 0.59 4.00 158.13 7.90 0 . ± .
03 2 . ± .
75 31.27NGC 6541 567.00 -1.81 2.51 1.48 120.84 7.50 0 . ± .
06 2 . ± .
63 77.12NGC 6652 805.00 -0.81 0.47 1.22 106.17 10.00 0 . ± .
06 4 . ± .
03 120.53NGC 6717 46.10 -1.26 0.36 1.56 22.38 7.10 0 . ± .
07 2 . ± .
63 70.85NGC 6752 374.00 -1.54 2.29 0.83 18.59 4.00 0 . ± .
05 0 . ± .
12 157.19NGC 6838 2.05 -0.78 0.54 0.85 4.14 4.00 0 . ± .
08 0 . ± .
18 40.13NGC 7078 6460.00 -2.37 4.90 0.39 248.98 10.40 0 . ± .
05 2 . ± .
78 46.55Omega Cen 144.00 -1.53 33.11 0.71 27.35 5.20 0 . ± .
07 3 . ± .
42 747.94Terzan 1 0.63 -1.03 2.95 4.06 4.27 6.70 0 . ± .
02 2 . ± .
73 62.48Terzan 2 19.60 -0.69 0.33 4.23 9.33 7.50 0 . ± .
05 3 . ± .
00 42.65Terzan 5 1400.00 -0.23 6.17 4.80 98.73 6.90 3 . ± .
20 38 . ± .
51 3740.32 a Stellar encounter rate computed using equation (2). The numerical values are normalized by the encounter rate of NGC 104, which is set to1000. b Metallicity. c Stellar mass. d Galactic photon field energy density. e Total photon field energy density, defined as the sum of the Galactic photon field and the photons from stars in the GC. f Distance from the Sun. g 𝛾 -ray flux and luminosity between 300 MeV to 500 GeV. h Test statistic. the uncertainties on 𝑅 (cid:12) for the GCs since they are either unavail-able in the Harris (1996) catalog or estimated only at percentagelevel (Baumgardt 2017) and so make a negligible contribution tothe overall error of 𝐿 𝛾 . For the parameters and flux upper limits of127 additional GCs, see Appendix A. In this section, we investigate the correlation between 𝐿 𝛾 ’s and otherGC observables. However, GCs not yet detected in 𝛾 rays and sampleselection effects must be taken into account to properly determinethe significance of any apparent correlations. We use the Kendall 𝜏 coefficient as the test statistic for estimating the significance of thecorrelations, and the expectation-maximization (EM) algorithm forthe linear regression of the correlations. Both methods allow us toproperly incorporate the luminosity upper limits − implied by GCsnot detected in the 4FGL − into our statistical analysis. To study the correlations between the 𝐿 𝛾 ’s and the other GC ob-servables, we assume a linear relation in logarithmic space of theform:log ( 𝐿 𝛾 ) = 𝑎 log ( 𝑋 ) + 𝑏, (3)where 𝐿 𝛾 is the gamma-ray luminosity of the GC, 𝑋 is the inde-pendent observable considered, and 𝑎 and 𝑏 are parameters to bedetermined.We use an EM algorithm (Feigelson & Nelson 1985; Isobeet al. 1986; Lavalley et al. 1992) to find the maximum likelihoodestimates of the parameters 𝑎 and 𝑏 . In contrast to the standardmaximum likelihood method, the EM algorithm is designed to beused with censored data, i.e., data consisting of both measurementsand limits. Upper limits must be properly incorporated in the cor-relation analyses so as to obtain statistically robust results. Briefly,the implementation of the EM algorithm is done as follows: first,the expected values of the censored data are estimated based on MNRAS000
51 3740.32 a Stellar encounter rate computed using equation (2). The numerical values are normalized by the encounter rate of NGC 104, which is set to1000. b Metallicity. c Stellar mass. d Galactic photon field energy density. e Total photon field energy density, defined as the sum of the Galactic photon field and the photons from stars in the GC. f Distance from the Sun. g 𝛾 -ray flux and luminosity between 300 MeV to 500 GeV. h Test statistic. the uncertainties on 𝑅 (cid:12) for the GCs since they are either unavail-able in the Harris (1996) catalog or estimated only at percentagelevel (Baumgardt 2017) and so make a negligible contribution tothe overall error of 𝐿 𝛾 . For the parameters and flux upper limits of127 additional GCs, see Appendix A. In this section, we investigate the correlation between 𝐿 𝛾 ’s and otherGC observables. However, GCs not yet detected in 𝛾 rays and sampleselection effects must be taken into account to properly determinethe significance of any apparent correlations. We use the Kendall 𝜏 coefficient as the test statistic for estimating the significance of thecorrelations, and the expectation-maximization (EM) algorithm forthe linear regression of the correlations. Both methods allow us toproperly incorporate the luminosity upper limits − implied by GCsnot detected in the 4FGL − into our statistical analysis. To study the correlations between the 𝐿 𝛾 ’s and the other GC ob-servables, we assume a linear relation in logarithmic space of theform:log ( 𝐿 𝛾 ) = 𝑎 log ( 𝑋 ) + 𝑏, (3)where 𝐿 𝛾 is the gamma-ray luminosity of the GC, 𝑋 is the inde-pendent observable considered, and 𝑎 and 𝑏 are parameters to bedetermined.We use an EM algorithm (Feigelson & Nelson 1985; Isobeet al. 1986; Lavalley et al. 1992) to find the maximum likelihoodestimates of the parameters 𝑎 and 𝑏 . In contrast to the standardmaximum likelihood method, the EM algorithm is designed to beused with censored data, i.e., data consisting of both measurementsand limits. Upper limits must be properly incorporated in the cor-relation analyses so as to obtain statistically robust results. Briefly,the implementation of the EM algorithm is done as follows: first,the expected values of the censored data are estimated based on MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters the regression parameters and the variance of the uncensored data.Second, a least-squares fit is performed and the variance is updated.Lastly, the procedure is repeated until convergence is achieved on 𝑎 , 𝑏 , and the variance. Using the EM algorithm, we are able toutilize the complete data set (including upper limits) in estimatingrelations between the 𝐿 𝛾 and the other observables. 𝜏 coefficient and significance While the EM algorithm allows us to estimate the linear relationsbetween the 𝐿 𝛾 ’s and other GC observables, we are also interestedin determining the statistical significance of those relations. To thatend, we apply the generalized Kendall 𝜏 rank correlation test andperform Monte Carlo (MC) simulations to determine the signifi-cance of each correlation studied with the EM algorithm.The Kendall 𝜏 rank correlation coefficient (also referred to asthe Kendall 𝜏 coefficient) is a non-parametric statistical test that hasbeen used to study multi-wavelength correlations of star-forminggalaxies (Ackermann et al. 2012; Ajello et al. 2020), and misalignedactive galactic nuclei (Di Mauro et al. 2014). It has been generalizedto include upper limits in the statistical procedure (Ackermann et al.2012). Therefore, we can calculate the Kendall 𝜏 coefficient usingall available information concerning GCs (measurements and upperlimits).To estimate the significance of the correlations, we adopt a sim-ilar procedure as advanced previously in the literature (Ackermannet al. 2012). Namely, the null hypothesis assumes no correlationbetween 𝐿 𝛾 and 𝑋 . A set of null hypothesis samples is generatedby repeating the following steps: (1) randomly exchange 𝐿 𝛾 of twoGCs while preserving their locations; (2) if the energy fluxes ofthe GCs after exchanging the 𝐿 𝛾 are above the detection thresholdof Fermi -LAT, the exchange is kept ; and (3) we perform a largenumber of exchanges, until obtaining a nearly uniform 𝐿 𝛾 sample(including corrections from applying the detection threshold) over 𝑋 , as required by the null hypothesis. In Appendix B, we discussthe number of exchanges needed to generate the null hypothesissample.For each correlation, we generate 10 null hypothesis samplesand calculate their Kendall 𝜏 coefficients. For a large number ofsamples, the coefficients can be fitted to a normal distribution (Efron& Petrosian 1999), { ˆ 𝜏 𝑖 } ∼ 𝑁 ( 𝜈 , 𝜎 ) , (4)where { ˆ 𝜏 𝑖 } represents the distribution of the 𝜏 coefficients from thenull hypothesis sample, and 𝑁 ( 𝜈 , 𝜎 ) is a normal distribution withmean 𝜈 and standard deviation 𝜎 . For each correlation, we cancompare the observed value of 𝜏 with the corresponding normaldistribution from the MC results and compute the significance, 𝜎 = 𝜏 − 𝜈 𝜎 . (5)Figure 2 shows an example of the 𝐿 𝛾 – 𝑢 Total data set. The bluehistogram shows the probability density of Kendall 𝜏 coefficients This step guarantees the detectability of the null hypothesis samples. Itis crucial to apply realistic estimates of the detection threshold so that thenull hypothesis samples are valid. Ackermann et al. (2012) and Ajello et al.(2020) have used the minimum fluxes in their data. Since we are using thesame amount of data as the 4FGL, we take advantage of the spatial map ofthe 8-year LAT detection threshold published with the 4FGL. We expect thisto be a more rigorous way of generating the samples since the map includesthe spatial dependence of the LAT threshold. .
05 0 .
06 0 .
07 0 .
08 0 . τ coefficient10 − P r o b a b ili t y d e n s i t y vs u Total
Figure 2.
Probability density distribution of the Kendall 𝜏 coefficients forthe 𝐿 𝛾 – 𝑢 Total data set. The blue histogram corresponds to the density of theKendall 𝜏 coefficients of the null hypothesis samples. The dash-dotted lineshows the best-fit normally-distributed probability density function for thenull hypothesis. The red vertical line indicates the Kendall 𝜏 coefficient forthe actual data. Table 2.
Summary of correlations between 𝐿 𝛾 and four astrophysical pa-rameters of the GCs. The best-fit parameters 𝑎, b, and the correspondingvariance of 𝐿 𝛾 are found using the EM algorithm. The significance of thecorrelations is found by MC simulations with Kendall 𝜏 coefficients.Correlation 𝑎 𝑏 √ Variance Significancevs Γ 𝑐 ± ± ± 𝜎 vs 𝑢 Total ± ± ± 𝜎 vs [Fe/H] 0.35 ± ± ± 𝜎 vs 𝑢 MW ± ± ± 𝜎 of the null hypothesis samples. The dash-dotted line is the best fitnormal distribution of the probability density, which has 𝜈 = . 𝜎 = . 𝜏 coefficient of real data is 0.093,shown by the red vertical line. The real data is about 3.8 𝜎 awayfrom the center of the null hypothesis distribution. The top (bottom) panel of Figure 3 shows the correlations between 𝐿 𝛾 and 𝑢 MW ( 𝑢 Total ). GCs with measured 𝛾 -ray luminosity areshown in red, while GCs with upper limits are shown in blue. We finda very small slope for the 𝐿 𝛾 - 𝑢 MW correlation, with 𝑎 = . ± . 𝐿 𝛾 - 𝑢 MW correlation is found to be1.5 𝜎 . When the total photon field is considered, we find a 𝐿 𝛾 – 𝑢 Total correlation with 𝑎 = . ± .
09. In this case, the significanceincreases to 3.8 𝜎 . The 𝐿 𝛾 – 𝑢 Total correlation is mostly driven by 𝑢 GC , the photon field from the starlight in the GCs (see equation (1)).As shown by Table 1, 𝑢 Total is much greater than 𝑢 MW due to thedominant contribution from 𝑢 GC .We also investigate the correlation of the 𝐿 𝛾 ’s with the stellar MNRAS , 1–19 (2021)
Song et al. u MW [eV cm − ]10 L γ [ e r g s − ] a = 0 . ± . u Total [eV cm − ]10 L γ [ e r g s − ] a = 0 . ± . Figure 3.
Correlations between 𝐿 𝛾 and the photon field energy densities.The top panel shows the 𝐿 𝛾 - 𝑢 MW correlation and the bottom panel showsthe 𝐿 𝛾 - 𝑢 Total correlation. GCs with measured 𝛾 rays are shown in red,while GCs with upper limits are shown in blue. The best-fit correlations(black solid lines) are calculated using the EM algorithm discussed in Sec-tion 3.1, with 1 𝜎 uncertainties included as the gray shaded bands. We finda shallow correlation between 𝐿 𝛾 and 𝑢 MW with 𝑎 = . ± .
26. Thecorrelation between 𝐿 𝛾 and 𝑢 Total is more significant, with 𝑎 = . ± . encounter rate ( Γ 𝑐 ) and GC metallicities ([Fe/H]). These observ-ables are argued to berelated to the formation of MSPs and mayprovide a proxy for the total number of MSPs in GCs. Figure 4shows the 𝐿 𝛾 – Γ 𝑐 correlation (top panel) and the 𝐿 𝛾 –[Fe/H] corre-lation (bottom pannel) obtained with the EM algorithm. We find apositive correlation between the 𝐿 𝛾 and Γ 𝑐 , with 𝑎 = . ± . 𝜏 test yields a 6.4 𝜎 statistical significance.Similarly, we find a correlation between 𝐿 𝛾 and [Fe/H] with the − Γ c L γ [ e r g s − ] a = 0 . ± . − − [Fe/H]10 L γ [ e r g s − ] a = 0 . ± . Figure 4.
Same as Figure 3, but correlated with the encounter rate (toppanel) and metallicity (bottom panel). best-fit value 𝑎 = . ± .
14. However, the statistical significanceof the correlation is only 1.8 𝜎 .We summarize the best-fit correlation results and their respec-tive statistical significance in Table 2. Positive and statistically significant correlations are obtained in boththe 𝐿 𝛾 – 𝑢 Total and the 𝐿 𝛾 – Γ 𝑐 space. The positive 𝐿 𝛾 – 𝑢 Total corre-lation could indicate a significant contribution from IC emission. Ifthe 𝑒 ± injected by MSPs lose energy through multiple comparableprocesses, e.g., IC and synchrotron radiaton, the 𝐿 𝛾 is proportionalto the IC energy loss rates, which is linearly proportional to the 𝑢 Total . In the extreme limit where all the 𝑒 ± injected by MSPs losetheir energy through IC, the 𝐿 𝛾 is constrained by the energy injec-tion rate of 𝑒 ± by MSPs and the 𝑢 Total would have less impact. Since
MNRAS000
MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters u Total [eV cm − ] − − l og Γ c a = 1 . ± . Figure 5.
Hidden correlation between 𝑢 Total and Γ 𝑐 . GCs with higher en-counter rates tend to have higher total photon field energy densities. Thered line shows the relation between 𝑢 Total and Γ 𝑐 based on a least-squiresmethod in logarithmic space. we find a preference for a non-linear correlation ( 𝑎 = . ± . 𝛾 rays are unlikely all originated from IC radiation.However, the positive correlation between 𝐿 𝛾 and 𝑢 Total couldalternatively be driven by the 𝐿 𝛾 – Γ 𝑐 correlation. Here, we investi-gate a potential hidden correlation between 𝑢 Total and Γ 𝑐 in orderto better understand the nature of our detections. Figure 5 shows the 𝑢 Total and Γ 𝑐 values for our sample of GCs. It is apparent that the 𝑢 Total tends to be higher for GCs with higher encounter rates. Sincethese data are uncensored, we simply estimate the correlation usingthe Spearman coefficient: we find 0.72, confirming a strong corre-lation. This result is not surprising because a higher photon densityimplies higher stellar density which implies higher encounter rates(see also equation (2)).An important implication of this result is that the 𝐿 𝛾 – 𝑢 Total and the 𝐿 𝛾 – Γ 𝑐 correlations are not necessarily independent. Usinga simple least squares method in the logarithmic space, we find therelation between Γ 𝑐 and 𝑢 Total to be Γ 𝑐 ∝ 𝑢 . ± . . (6)As reported in Table 2, the correlation between 𝐿 𝛾 and Γ 𝑐 has apower index 𝑎 = . ± .
10. Based on the hidden relation between Γ 𝑐 and 𝑢 Total , the projected correlation between 𝐿 𝛾 and 𝑢 Total wouldhave an index 𝑎 = . ± .
15. Within the uncertainty, this projectedresult is consistent with the directly measured correlation between 𝐿 𝛾 and 𝑢 Total found in real data, 𝑎 = . ± .
09. Therefore, thepositive correlation between 𝐿 𝛾 and 𝑢 Total could be evidence for IC,or alternatively, an indirect effect of the 𝐿 𝛾 – Γ 𝑐 correlation which isconnected to the dynamic formation of MSPs. The correlation foundbetween 𝐿 𝛾 and 𝑢 Total cannot be considered as concrete evidencefor IC due to this ambiguity implicated by the hidden correlation.However, as we discuss next, evidence for IC emission in GCsmay still be revealed from the detailed spectral properties of theseobjects.
Motivated by the correlations detected in the previous section, weperform a spectral analysis of the 30 GCs detected in the 4FGLcatalog, with the aim of finding further evidence for IC emission.First, we model the spectra of the GCs individually and comparetheir spectral parameters with those describing the field MSPs. Sec-ond, we fit the GCs spectra with universal spectral models whichphenomenologically describe possible IC emission. Lastly, we usethe detected IC component to constrain the 𝑒 ± injection efficiencyin the GCs. We consider two possible mechanisms of 𝛾 -ray emission, not mu-tually exclusive: CR, and IC up-scattering of starlight.Detailed theoretical models predict that the maximum energyof the 𝑒 ± accelerated by the MSPs is limited by the CR in the pulsarmagnetosphere. The predicted CR spectrum exhibits an energy cut-off which is related to the 𝑒 ± Lorentz factor (Harding et al. 2005).For this reason, we model the GC 𝛾 -ray spectrum–as predicted byCR models–using a power law with an exponential cut-off (PLE) ofthe form: (cid:20) 𝑑𝑁𝑑𝐸 (cid:21) CR = 𝑁 (cid:18) 𝐸𝐸 (cid:19) − Γ exp (cid:18) − 𝐸𝐸 cut (cid:19) , (7)where 𝑁 is the normalization factor, Γ is the spectral index, 𝐸 isthe scaling energy, and 𝐸 cut is the energy cutoff.The 𝑒 ± may also leave the MSPs through open magnetic fieldlines and diffuse into the GC medium. Escaping pairs may up-scatter ambient photons and produce IC emission. The spectrum ofthe IC is determined by the 𝑒 ± spectrum and the ambient photonfield. Theoretical studies (Harding & Muslimov 2011) show thatthe MSPs can inject 𝑒 ± with Lorentz factors 𝛾 𝑒 ± > efficiently.Given ambient photons of 𝐸 ∼ 𝛾 𝑒 ± 𝐸 = 𝑒 ± . In the Thomson regime, the IC spectrum resulting from theinteraction of power-law-like 𝑒 ± with ambient photons followinga black-body radiation distribution (Blumenthal & Gould 1970) isstill a power law in 𝛾 -ray energy. We consider this spectral form asa phenomenological description of the IC model. Specifically , (cid:20) 𝑑𝑁𝑑𝐸 (cid:21) IC = 𝑁 (cid:18) 𝐸𝐸 (cid:19) − Γ . (8)We first estimate the GCs’ spectral parameters using a maxi-mum likelihood method. For this, we use the CR model and the ICmodel separately. We perform a 𝜒 test using the bin-by-bin 𝛾 -rayfluxes (9 energy bins from 300 MeV to 500 GeV) of each GC andthe CR and IC emission models. Therefore, we define 𝜒 = ∑︁ 𝑖 ( 𝐹 𝑖 data − 𝐹 𝑖 model ) ( Δ 𝐹 𝑖 data ) + ( 𝑓 𝑖 ref 𝐹 𝑖 data ) , (9)where 𝐹 𝑖 data and Δ 𝐹 𝑖 data are the measured fluxes and flux uncer-tainties obtained at each independent energy bin, 𝐹 𝑖 model are thepredicted fluxes (either the CR or IC models). We allow all model For the maximum 𝛾 -ray energy (hundreds of GeV) and the photon field(starlight) we considered, the IC is in transition from the Thomson regimeto the Klein-Nishina regime with the Thomson regime still an adequateapproximation.MNRAS , 1–19 (2021) Song et al. parameters to be free (i.e., normalization, power-law index, and cut-off energy for CR, and normalization and power-law index for IC).The 𝑓 𝑖 ref values encapsulate the systematic uncertainties on the ef-fective area of the LAT. We follow the values reported in the 4FGLcatalog (Abdollahi et al. 2020), and set 𝑓 ref to 0.05 for the first threeenergy bins, 0.06 for the fourth bin, and 0.1 for the last five bins.The significance of the spectral curvature is estimated by com-puting the difference of the best-fit 𝜒 between the IC and theCR models, TS curve = 𝜒 − 𝜒 . We apply a 2 𝜎 threshold todetermine the type of spectrum: for GCs with TS curve ≥
4, theirPLE spectra are reported. Otherwise, the power-law spectra are re-ported. Note that this is a lower threshold than the 4FGL, whichrequires TS curve ≥ 𝜎 threshold adequate in our analysis since our fits generate finite 𝐸 cut ’s within uncertainties for all GCs with TS curve ≥
4. For thoseGCs with TS curve <
4, the fits only generate lower limits for 𝐸 cut .Table 3 summarizes the best-fit parameters of the spectra for 30 𝛾 -ray-detected GCs, sorted by their TS curve . The majority prefercurved spectra, with only 5 preferring simple power law spectra.Figure 6 shows the spectra for 2 GCs as examples. The top panelshows the spectrum of NGC 6397, which is best fit by a simple powerlaw, while the lower panel shows the spectrum for NGC 6541, whichprefers an exponential cut-off at ∼
350 MeV with TS curve = Fermi -LAT has detected more than 200 pulsars. Most ofthese have been found to have a curved spectrum with best-fit energycutoffs of the order of a few GeV. Therefore, their 𝛾 -ray emission islikely dominated by a CR process. Nevertheless, Hooper & Linden(2018) find that many MSPs could be surrounded by TeV halosof IC. The IC emission may also extend to the GeV energy range.Figure 7 compares the distribution of the spectral parameters of108 field MSPs in the 4FGL (red dots) with the 𝛾 -ray GCs (bluedots), assuming a PLE spectra. The 1 𝜎 uncertainties of the best-fit parameters are also shown. We find that within uncertainties,the spectral distribution of the GCs and the field MSPs are verysimilar. However, given the starlight in GCs typically contributesa much larger photon field energy density than for field MSPs, ICemission may still provide a sizeable contribution to the overall GC 𝛾 -ray emission. The results from individual spectral fit cannot ruleout the presence of IC for the following reasons: (1) There are 5GCs for which the spectra shows no obvious energy cutoffs. This ishard to explain using the CR emission model alone. (2) Many GCshave energy bins above 10 GeV detected even though their spectrahave cutoffs of order a GeV (see Appendix C). These high-energymeasurements may be indicative of an IC component. The curvature of the GCs’ spectra at around a few GeV’s, as well astheir similarity to the field MSPs spectra, support the hypothesis thatthe GeV 𝛾 -ray emission from most GCs is due to mainly local CRemission from MSPs within GCs. However, IC may still contributesub-dominantly, especially at the high-energy end. To probe thispossibility, we perform a reduced 𝜒 analysis in which we fit, bin-by-bin, the GCs’ spectra using a linear combination of the spectral Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 639710 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6541
Figure 6.
Best-fit spectra (blue solid line) and the 1 𝜎 uncertainties (blueband) for two GCs. The bin-by-bin fluxes from the Fermi data analysisare included as black points. The top panel shows the spectrum of NGC6397 as a simple power law because the PLE model is only slightly favored(TS curve = ∼ GeV with TS curve = components introduced in equation 7 and 8. Specifically; 𝑑𝑁𝑑𝐸 = (cid:20) 𝑑𝑁𝑑𝐸 (cid:21) CR + (cid:20) 𝑑𝑁𝑑𝐸 (cid:21) IC (10) = 𝑁 (cid:18) 𝐸𝐸 (cid:19) − Γ exp (cid:18) − 𝐸𝐸 cut (cid:19) + 𝑁 (cid:18) 𝐸𝐸 (cid:19) − Γ . Fitting such a two-component model to each GC’s bin-by-bindata is difficult since the GC spectra only contains 9 energy bins, andmany high energy bins only provide upper limits. On the other hand,typical GCs can host close to ∼
20 MSPs (Ye et al. 2019) each. So,as a simplifying approximation, we hypothesise that the 𝛾 -ray and 𝑒 ± injection from the collection of MSPs in each GC to be similarto one another. Then, we can fit a common or universal spectrum to MNRAS000
20 MSPs (Ye et al. 2019) each. So,as a simplifying approximation, we hypothesise that the 𝛾 -ray and 𝑒 ± injection from the collection of MSPs in each GC to be similarto one another. Then, we can fit a common or universal spectrum to MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters Table 3.
Spectral parameters for 30 𝛾 -ray-detected GCs from the individual fits, ordered from the least curved to the most curved. For GCs with TS curve < 𝜎 ), the best-fit simple power laws (PL) are reported. For the rest GCs, the power laws with an exponential cutoff (PLE) are reported.Name Γ log (cid:16) 𝐸 cut MeV (cid:17) 𝜒 /d.o.f. Type a TS curve ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± a Spectrum type: PL for power law; PLE for power law with an exponentialcut-off. all the 𝛾 -ray detected GCs, i.e., one set of spectral shape parametersin the two component model in equation (10) for all the GCs. Morespecifically, we tie the Γ , Γ and 𝐸 cut across all GCs considered(hereafter referred to as the “universal model”). The normalizationfactors 𝑁 and 𝑁 are allowed to float for each GC as these shoulddepend on the number of MSPs and the photon field energy densityin the GCs.We perform the universal fit by minimizing the total 𝜒 of 30 𝛾 -ray-detected GCs, 𝜒 ( Γ , Γ , 𝐸 cut ) = ∑︁ 𝑖 𝜒 𝑖 ( Γ , Γ , 𝐸 cut , 𝑁 𝑖 , 𝑁 𝑖 ) . (11)In practice, we assign the same Γ , Γ , and 𝐸 cut to all 𝛾 -ray GCs andperform a minimum 𝜒 for each different object. However, duringthe fit, we free the 𝑁 𝑖 and 𝑁 𝑖 parameters. By scanning the parameterspace of Γ , Γ , and 𝐸 cut , we find the values that minimize the total 𝜒 for the two-component model. These are, Γ = . ± . , Γ = . ± . , log (cid:18) 𝐸 cut MeV (cid:19) = . ± . , for which we find a 𝜒 / d . o . f = / = .
99 (we have 30 × + . o . f = × − − − = 𝜎 contours and correlated uncertainties forthe parameters Γ , Γ , and 𝐸 cut as found in this procedure.In order to compute the statistical significance of the IC compo-nent, it is necessary to define the null hypothesis. This correspondsto the universal model containing only the CR component (seeequation 7). Again, we tie Γ and 𝐸 cut across all GCs and allow thenormalization factors to individually vary. We find that the best-fitparameters for the CR-only model are: Γ = . ± . , log (cid:18) 𝐸 cut MeV (cid:19) = . ± . . In this case, we find a 𝜒 / d . o . f = / = .
47 (the nullhypothesis has 30 + 2 free parameters so the d . o . f = × − − − = 𝜎 level ( Δ 𝜒 = −
204 for 31 d.o.f [1 power-law indexplus 30 normalization factors]). It is useful to compare the best-fitspectral results of the CR component for the universal models withthe best-fit spectral parameters of the MSPs in the 4FGL catalog.As seen in Figure 9, although the CR-only model (null hypothesis)has larger Γ and higher 𝐸 cut than the CR component from the two-component model, our results for both models are compatible withthe field MSPs, up to statistical uncertainties.The universal fitting procedure used in this section is similar MNRAS , 1–19 (2021) Song et al. . . . . . E cut [MeV] − Γ MSPsGCs
Figure 7.
Spectral parameters of the 𝛾 -ray emission. The PLE spectrumis assumed, with Γ and 𝐸 cut as free parameters. Both the GCs (blue dots)and the MSPs (red dots) detected in the 4FGL are included. The errorbars represent the 1 𝜎 parameter uncertainties. Within uncertainties, thedistribution of the GCs’ spectra is very similar to that of the MSPs. to a stacking analysis. This method is usually applied to explore thecharacteristics of an astrophysical population, especially one that isundetected. Numerous studies have shown that this technique canincrease the detection sensitivity to such population characteristics.So, even though there is good statistical evidence for the IC compo-nent in the universal fit, this might not be apparent from individualfitting of the two-component model.We show examples of the spectra obtained in the universal fitof the two-component model for NGC 6397 and NGC 6541 in Fig-ure 10. As can be seen, the solutions for the CR and IC componentslook physically plausible. The spectra also include 1 𝜎 bow-tie er-rors, which immediately reveal the level of statistical support for theCR and IC components, respectively. For comparison, the resultsshown in Figure 6 presented a traditional (Abdollahi et al. 2020)spectral curvature analysis applied to NGC 6397 and NGC 6541,individually.We show some additional noteworthy results of the universalfit in Figure 11. Here, we see that in the case of GC 2MS-GC01,the IC model is sufficient to explain the bin-by-bin spectrum overthe full energy range, but we also display the estimated 95% C.L.upper limit for the normalization of the CR model. By contrast, inthe case of GC M 80, we find that the data is best described by theCR model alone, and we show the 95% C.L. upper limit for thenormalization of the IC component. These examples might indicatespecial conditions of the environment of the GC.For 19 GCs (out of the 30 GCs included in the universal fit),we find good statistical support for both the CR and IC models. Forthe remaining 11 GCs we find that only one component is sufficientto explain the spectrum: 7 GCs require only the CR model and theother 4 GCs require only the IC model. The two-component spectralresults for all 30 GCs are shown in Appendix C.To explain the best-fit power law index of the IC component( Γ = . ± . 𝑒 ± spectrum would have anindex of 4 . ± .
50. The minimum 𝑒 ± energy required to maintain a power law IC in the energy range of our analysis (300 MeV) is (cid:46)
10 GeV given that the upscattered photon field has energies ∼ 𝑒 ± pair cascades from pulsar polarcaps. For typical MSP parameters, they show that the injected 𝑒 ± flux decreases by ∼ −
10 orders of magnitude when the 𝑒 ± energyincreases from ∼
10 GeV to ∼ 𝑒 ± spectrum we foundis in line with their results. The relative normalization of the CR and IC components probesan important property of the MSPs: the 𝛾 -ray production and 𝑒 ± injection efficiencies, respectively. Indeed, the spin-down energy ofMSPs can be injected into 𝛾 rays and 𝑒 ± . While prompt 𝛾 rays aremainly produced by CR in the magnetosphere, the 𝑒 ± can propagateinto the interstellar environment. We can write down the followingempirical relations, 𝐿 CR = 𝑓 𝛾 𝐿 sd , (12) 𝐿 𝑒 ± = 𝑓 𝑒 ± 𝐿 sd , (13)where 𝐿 sd is the spin-down luminosity and the 𝑓 ’s are efficiencyparameters.Assuming that the 𝛾 -ray emission is the superposition of CRand IC processes, we have that 𝐿 𝛾 = 𝐿 CR + 𝐿 IC . (14)However, 𝑒 ± can also lose energy via synchrotron radiation. We cancompare the relative strength of the sychrotron radiation vs the ICemission through (cid:164) 𝐸 SR (cid:164) 𝐸 IC = 𝑢 𝑢 , (15)where (cid:164) 𝐸 SR and (cid:164) 𝐸 IC are the synchrotron and IC energy loss rates,respectively, 𝑢 B is the magnetic energy density, and 𝑢 rad is theradiation field energy density. We note that this relation assumes thatthe 𝑒 ± lose all their energy within the GCs. We provide justificationsfor this assumption in Appendix D. For a typical GC the magneticfield is estimated to be 𝐵 (cid:46) 𝜇 G (Bednarek & Sitarek 2007), sowe expect to have a 𝑢 B = ( 𝜇 G ) /( 𝜇 ) = . − , whichis much smaller than the total radiation field of most GCs shown inTable 1. Thus in the usual instance when IC is the leading energyloss process, we have that 𝐿 IC (cid:39) 𝐿 𝑒 ± . (16)Since no GC is detected as an extended source by the Fermi -LAT, the energy carried away to the interstellar medium by 𝑒 ± propagation is expected to be small. Thus, we can use the followingapproximate scaling relation, 𝑓 𝑒 ± 𝑓 𝛾 (cid:39) 𝐿 IC 𝐿 CR . (17)Using this , we estimate the ratios 𝑓 𝑒 ± / 𝑓 𝛾 for all 𝛾 -ray emitting GCsin Table 4. These are found to be in the range ≈ . − .
04. Notethat for some GCs we present only upper or lower limits as onlyone component is detected. The measurement of pulsars by
Fermi -LAT estimated the 𝑓 𝛾 efficiency from observations of pulsars and We discuss caveats to the approximation used in equation (17) in Ap-pendix D. MNRAS000
Fermi -LAT estimated the 𝑓 𝛾 efficiency from observations of pulsars and We discuss caveats to the approximation used in equation (17) in Ap-pendix D. MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters − . . . . . . . . . . l og E c u t [ M e V ] − . . . . . . . . . . . . Γ . . . . E cut [MeV]2 . . . . . . Γ Figure 8.
The projected parameter space of the universal model fit, as illustrated in equation (10). The blue shaded contours show the 1 𝜎 , 2 𝜎 , and 3 𝜎 confidence levels for the two-component model. The crosses indicate the best fit values for Γ , Γ , and 𝐸 cut . On the top-left panel, the red shaded region showsthe best-fit value and 3 𝜎 confidence levels for the null hypothesis, which includes only the CR model. found that on average, 𝑓 𝛾 ∼ 𝑒 ± efficiency 𝑓 𝑒 ± was also estimated to be around 10% from TeV observationsof nearby pulsars (Hooper et al. 2017; Hooper & Linden 2018) andthe Galactic center (Bednarek & Sobczak 2013), although MAGICCollaboration et al. (2019) claims 𝑓 𝑒 ± is at the percentage level forone GC they observed (NGC 7078), and Sudoh et al. (2020) suggest 𝑓 𝑒 ± ∼
90% on the basis of the radio continuum emission detectedfrom galaxies with low specific star formation rates.For the CR and IC luminosities, we integrate the best-fit two-component spectra from 300 MeV to 500 GeV, the same energyrange used in the Fermi data analysis. For the IC emission, theminimum 𝑒 ± injection energy probed by this energy range is (cid:46) 𝑒 ± pair cascades fromMSPs and proposed several theoretical models. Their Figure 10shows that their predicted pair spectra peak at ∼ GeV and extendto (cid:38)
TeV. This roughly corresponds to the Fermi energy range weassume. If the 𝑒 ± injection spectra extend to lower energy, they willlead to higher 𝐿 IC . Therefore, the choice of 𝛾 -ray energy range willcontribute as systematic uncertainties on the estimated 𝑓 𝑒 ± / 𝑓 𝛾 . For example, we verify that the 𝑓 𝑒 ± / 𝑓 𝛾 would be ∼ 𝛾 -ray energy is assumed to be 30 MeV. We have found strong positive correlations between 𝐿 𝛾 , the stellarencounter rate Γ 𝑐 , and the total photon field energy density 𝑢 Total ofGCs. The latter correlation may indicate a significant contributionof IC upscattering of ambient starlight to the total 𝛾 -ray emission ofGCs. However, we showed in Figure 5 that the 𝑢 Total also increaseswith Γ 𝑐 . So, the detection of the 𝐿 𝛾 – 𝑢 Total correlation alone doesnot unambiguously demonstrate the presence of IC emission inGCs . On the other hand, corroborating evidence for IC emissionwas found from the universal two-component fit, wherein we were We also analyzed other potential hidden correlations, but no obviouscorrelations with other parameters such as the interstellar radiation field andthe distance from the Sun were found (see Appendix D).MNRAS , 1–19 (2021) Song et al. . . . . . E cut [MeV] − Γ MSPs
Figure 9.
Same as Figure 7, but with best-fit parameters of GC replaced bythose obtained from the universal models. The 3 𝜎 contours are shown for theCR-only model (red) and the CR component from the two-component model(blue), as in Figure 8. The MSP parameters are included in the background(yellow dots). able to estimate, separately, the luminosities of the CR and ICcomponents of most GCs. The ratios of the luminosities betweenthe CR and IC components were found to be comparable. Thisimplies that MSPs in GCs can potentially inject 𝑒 ± as efficiently asthey inject prompt magnetospheric 𝛾 rays.Overall, our correlation results in the 𝐿 𝛾 - Γ 𝑐 plane are consis-tent with those in Hui et al. (2011) and de Menezes et al. (2019),though it is important to note that our method is more statisticallyrobust since we include GCs with detection limits which were pre-viously neglected. In particular, our high significance (6 . 𝜎 ) detec-tion of a 𝐿 𝛾 – Γ 𝑐 correlation naively supports a dynamic formationscenario for MSPs in GCs. However, as pointed out earlier, thiscorrelation may not be independent due to the hidden correlationof 𝑢 Total and Γ 𝑐 . On the other hand, we have not found an obviouscorrelation between 𝑓 𝑒 ± / 𝑓 𝛾 and 𝑢 Total (see Appendix D). The lackof this latter correlation may indicate that IC is, in fact, the leadingenergy loss process for 𝑒 ± in GCs: in the limit of IC dominance, theIC luminosity of GCs already saturates the power going into freshly-injected 𝑒 ± pairs, so “dialling-up” the light field energy density hasno effect on the IC luminosity. Thus, in this situation of IC dom-inance, we expect, at most, only a weak correlation between 𝐿 𝛾 and 𝑢 Total and we would anticipate that the 𝐿 𝛾 - Γ 𝑐 correlation is thefundamental one (while the 𝐿 𝛾 - 𝑢 Total correlation is caused by thefact that GCs with higher stellar encounter rate naturally have higherstellar density which leads to higher photon field density). With theuncertainties of the data and the number of variables involved, it ischallenging to statistically confirm this scenario. Overall, however,our results are consistent with there being both a significant rolefor dynamical formation of MSPs in GCs and for the presence of asignificant contribution of IC to the overall 𝛾 -ray emission of GCs.Previous studies (Hui et al. 2011) found a positive correlationbetween 𝐿 𝛾 and 𝑢 MW , and 𝐿 𝛾 and [Fe/H]. However, our study doesnot confirm these results. The former discrepancy is possibly dueto the different interstellar radiation field models assumed in these Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6397 CRIC10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6541 CRIC
Figure 10.
The best-fit two-component spectra for NGC 6397 (top panel)and NGC 6541 (bottom panel). The spectra are fit to a universal shape withsame Γ , Γ , and 𝐸 cut for all GCs. Only the normalizations of the twocomponents are allowed to vary between GCs. The best-fit parameters forthe CR component (red line with shaded band) is Γ = . ± .
44 andlog ( 𝐸 cut / MeV ) = . ± .
16. The best fit parameter for the IC component(blue line with shaded band) is Γ = . ± .
25. The black dashed lineindicates the total of the two components. works, or it could be due to the more limited sample data usedin Hui et al. (2011). Specifically, while we have used the most up-to-date interstellar radiation field for the Milky Way–which is the3D radiation field model in GALPROP v56 (Porter et al. 2017) − Hui et al. (2011) used the 2D radiation field model in GALPROPv54. Also, as explained above, our correlation study includes 30 𝛾 -ray-detected GCs , as well as the luminosity upper limits from the127 non-detected ones, thus covering the entire GC Harris (1996)catalog. As for the latter discrepancy, similar results for the 𝐿 𝛾 –[Fe/H] correlation were obtained by de Menezes et al. (2019), whichalso found low statistical evidence for this correlation. MNRAS000
25. The black dashed lineindicates the total of the two components. works, or it could be due to the more limited sample data usedin Hui et al. (2011). Specifically, while we have used the most up-to-date interstellar radiation field for the Milky Way–which is the3D radiation field model in GALPROP v56 (Porter et al. 2017) − Hui et al. (2011) used the 2D radiation field model in GALPROPv54. Also, as explained above, our correlation study includes 30 𝛾 -ray-detected GCs , as well as the luminosity upper limits from the127 non-detected ones, thus covering the entire GC Harris (1996)catalog. As for the latter discrepancy, similar results for the 𝐿 𝛾 –[Fe/H] correlation were obtained by de Menezes et al. (2019), whichalso found low statistical evidence for this correlation. MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters Energy [MeV]10 − − − − − − E d N / d E [ e r g c m − s − ] Energy [MeV]10 E d N / d E [ e r g c m s ] M 80 CRIC
Figure 11.
Spectra of 2MS-GC01 (top panel) and M 80 (bottom panel) fromthe two-component fit. The 2MS-GC01 prefers only the power-law model(blue solid line). The 95% upper limit on the normalization of the CR modelis shown by the red shaded region. In contrast, only the CR model is detectedfor M 80 (red solid curve). The 95% upper limit on the normalization of theIC model is shown by the blue shaded region.
The emission from a putative population of about 40 ,
000 (Ploeget al. 2020) unresolved MSPs in the Galactic Center region is cur-rently the preferred explanation for the Fermi GeV excess (Maciaset al. 2018; Bartels et al. 2018; Macias et al. 2019; Abazajian et al.2020). Since GCs also contain large numbers of unresolved MSPs,it is useful to compare the light-to-mass ratios for these two systemsso as to obtain additional clues for the physical processes causingthe observed high-energy 𝛾 -ray emissions in their directions. InFigure 12, we show the relation between 𝐿 𝛾 and the stellar massfor several different systems. The blue dots show the sample of the 𝛾 -ray detected GCs in this work. The nuclear bulge (orange dot)has a stellar mass around 1 . × M (cid:12) and a 𝛾 -ray luminosity Table 4. 𝛾 -ray luminosity for IC and CR components and the ratios between 𝑓 𝑒 ± and 𝑓 𝛾 . For GCs with only one component detected, the 95% C.L. upperlimits are reported for another component.Name 𝐿 IC 𝐿 CR 𝑓 𝑒 ± / 𝑓 𝛾 (10 erg s − ) (10 erg s − )GLIMPSE02 10.90 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± of ( . ± . ) × erg s − and the boxy bulge (green dot) has1 . × M (cid:12) and ( . ± . ) × erg s − (Macias et al. 2019).The combination of the nuclear bulge and the boxy bulge is re-sponsible for the Galactic center GeV excess. Also included are theGalactic disk (red dot) luminosity predicted by Bartels et al. (2018)and the M31 galaxy (purple star) (Ackermann et al. 2017). Thedot-dashed line shows the 𝛾 -ray luminosity-to-stellar-mass relationimplied for the nuclear bulge and the boxy bulge, which is 2 × erg s − M − (cid:12) .As can be seen in Figure 12, the luminosities of the detectedsample of GCs exceed the luminosities expected based on the bulgecorrelations. In total, the GC samples have a stellar mass of ∼ . × M (cid:12) and a 𝛾 -ray luminosity of ∼ . × erg s − .This means that GCs systematically emit ∼
50 times more 𝛾 raysper stellar mass than other objects such as the nuclear bulge andthe Galactic bulge. The GCs have long been known for producingMSPs efficiently (Camilo & Rasio 2005). On average GCs make up ∼ .
05% of the total number of stars in the Milky Way (Ye et al.2019), but more than one-third of the known MPSs are found in thesesystems (Manchester et al. 2005). Our observations support thisscenario. This is also consistent with the larger stellar densities andlarger stellar encounter rates in GCs than in the Galactic bulge. Wealso note in passing that a large fraction of 𝛾 -ray-detected GCs arelocated in the Galactic bulge region (see Figure 1), which means thatit is possible that at least some of the unresolved MSPs contributing MNRAS , 1–19 (2021) Song et al. Stellar Mass [ M (cid:12) ]10 L γ [ e r g s − ] M31GCsNBBBDisk
Figure 12.
Relations between 𝐿 𝛾 and the stellar mass for several systems.The results for the nuclear bulge (orange dot) and boxy bulge (green dot) areadopted from Macias et al. (2019). The Galactic disk (red dot) luminosityis predicted by Bartels et al. (2018), and the M31 (purple star) are adoptedfrom Ackermann et al. (2017). The blue dots show the data for 30 𝛾 -raydetected GCs. The dash-dotted line is the relation implied by the nuclearbulge and boxy bulge. to the Fermi GeV excess are hosted by GCs in the Galactic bulgeregion. Using the universal two-component fit, we identified a power lawcomponent with a slope of 2 . ± .
25 from the spectra of 𝛾 -ray-detected GCs. The power law component can be plausibly explainedby IC emission from GCs. The fact that this power law componentis rather soft may explain why most GCs are not detected in the TeVenergy range. In order to explore this more closely, we extrapolatethe high energy tail of the GCs spectra to TeV energies in Figure 13.In this figure, the black line shows the extrapolated fluxes for Terzan5, and the gray band shows the range of extrapolated fluxes for theother 23 GCs with a detected IC component. Above 100 GeV, Fermi -LAT only find upper limits (blue arrows) for Terzan 5. The red dotsare the H.E.S.S. measurements from the direction of Terzan 5. Itis interesting to note that the extrapolated spectrum for Terzan 5is about one order of magnitude lower than the H.E.S.S. measure-ments from the same object. This discrepancy might be explainedby the fact that the 𝛾 -ray source reported by H.E.S.S is misalignedwith the center of Terzan 5 so that this association could be a chancecoincidence. However, such a coincidence with known objects hasbeen estimated to be improbable ( ∼ − ) (H. E. S. S. Collaborationet al. 2011). If the H.E.S.S. source is indeed associated to Terzan5, it could be that 𝑒 ± injection spectrum from MSPs has a spectralbreak at approximately 1 TeV. Note that a substantial fraction ofstars in Terzan 5 have been identified as young and centrally con-centrated (Ferraro et al. 2016; Gotta et al. 2020), which could leadto a larger number of younger pulsars. The H.E.S.S. measurementscould be explained if these young pulsars have higher energy 𝑒 ± cutoffs. Therefore, Terzan 5 may not be representative compared to − Energy [TeV]10 − − − − − − − − E d N / d E [ e r g c m − s − ] C T A h L H A A S O y r Terzan 523 GCsH.E.S.S. 2011Fermi
Figure 13.
Extrapolated spectra for Terzan 5 (black line) and 23 GCs with theIC component detected (grey region). The H.E.S.S. measurement of Terzan5 is shown by red dots with error bars. The
Fermi -LAT only find upperlimits (blue arrows) for Terzan 5 in this energy range. Also included are thesensitivities for 100-hour CTA South (orange line) and 1-year LHAASO(green line). other GCs which are dominated by old stellar systems. However,Nataf et al. (2019) also find that the abundance variations amongTerzan 5 is indeed consistent with a regular globular cluster. Al-ternatively, the TeV 𝛾 rays could originate from sources other thanMSPs (e.g., hadronic emission from supernova remnants). Furtherinvestigation of those scenarios, though very interesting, is beyondthe scope of this work.We also include in Figure 13 the sensitivities to point-likesources for the next generation 𝛾 ray observatories. The orangeline shows the sensitivity for the Cherenkov Telescope Array(CTA) − South assuming 100 hours of observation time. The greenline shows the 1-year sensitivity of the Large High Altitude AirShower Observatory (LHAASO). The extrapolated IC fluxes areclose to the 100-hour CTA sensitivity. It is clear that it will bedifficult for the next-generation TeV 𝛾 -ray telescopes to actuallydetect each individual GC considered in our study. This might re-quire a much more ambitious observation strategy that increasesthe sensitivity by factor of a few at the TeV energy range. Efforts tomeasure the diffuse IC emission from the putative MSP populationresponsible for the Fermi GeV excess have been made and are veryencouraging; see Song et al. (2019). Alternatively, Bednarek et al.(2016) studied TeV 𝛾 -ray emission from MSPs taking into accountthe advection of 𝑒 ± with the wind from the GC. They showed thatCTA can constrain models incorporating such effects. We have reanalyzed
Fermi -LAT data in the energy range between300 MeV and 500 GeV from the direction of 157 GCs in the Harris(1996) catalog. Using the same data cuts adopted in the constructionof the 4FGL catalog, we confirmed the detection of 30 GCs in 𝛾 rays, and updated the 𝛾 -ray spectral parameters for the sample of MNRAS000
Fermi -LAT data in the energy range between300 MeV and 500 GeV from the direction of 157 GCs in the Harris(1996) catalog. Using the same data cuts adopted in the constructionof the 4FGL catalog, we confirmed the detection of 30 GCs in 𝛾 rays, and updated the 𝛾 -ray spectral parameters for the sample of MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters detected objects. We also estimated the 95% C.L. luminosity upperlimits for the sample of 127 undetected GCs in the 4FGL catalog.The main objective of our reanalysis was to find evidence for ICemission from 𝑒 ± injected by MSPs in GCs. This was done usingtwo different methodologies. First, we searched for correlationsof the 𝛾 -ray luminosities with other GCs properties. Second, weperformed a spectral analysis of the GCs with a universal fit methodthat enhances the sensitivity to the high energy tail of the spectra.Specifically:1. Using an expectation-maximization algorithm that properlyincorporates null detections ( 𝐿 𝛾 upper limits) in the pipeline, wefound a correlation between 𝐿 𝛾 and the GCs’ total photon fieldenergy density 𝑢 Total of the form,log (cid:18) 𝐿 𝛾 erg s − (cid:19) = ( . ± . ) log (cid:18) 𝑢 Total eV cm − (cid:19) +( . ± . ) . (18)Using the Kendall 𝜏 coefficient as the test statistic we determinedthis correlation to have a 3 . 𝜎 significance. The total photon field isdominated by the stellar light of the GCs ( 𝑢 GC ), and we find a muchweaker correlation (below the 2 𝜎 level) when only the photon fieldat the location of the GC ( 𝑢 MW ) is used. In addition, we obtained astrong correlation (at 6.4 𝜎 significance) between 𝐿 𝛾 and the stellarencounter rate Γ 𝑐 , which is given bylog (cid:18) 𝐿 𝛾 erg s − (cid:19) = ( . ± . ) log ( Γ 𝑐 ) + ( . ± . ) . (19)Finally, we found only weak evidence (below the 2 𝜎 level) for acorrelations between 𝐿 𝛾 and the stellar metallicity [Fe/H].2. We revealed a hidden correlation between 𝑢 Total and Γ 𝑐 , whichimplies that the 𝐿 𝛾 – 𝑢 Total and 𝐿 𝛾 – Γ 𝑐 correlations are not entirelyindependent. However, as described below, we find spectral evi-dence for IC emission. The correlation results are consistent withthere being both a significant role for dynamical formation of MSPsin GCs and the for the presence of a significant contribution of ICto the observed 𝛾 -ray luminosity.3. We applied a universal spectral fit to the sample of 30 GCs inthe 4FGL catalog and searched for evidence of an IC componenton top of a curvature radiation model–accounting for the MSPsprompt emission in the GCs. We found that the extra power-law ICcomponent is preferred at the 8.2 𝜎 significance over the curvatureradiation model only. The best-fit power law index of the IC com-ponent was found to be 2 . ± .
25. This implies a power-law 𝑒 ± spectrum with an index of 4 . ± .
50 and a minimum energy aslow as ∼
10 GeV.4. We estimated the 𝑒 ± injection efficiency 𝑓 𝑒 ± for MSPs residingin GCs. We determined the IC 𝛾 -ray luminosities over 300 MeV to500 GeV, which roughly corresponds to 𝑒 ± energies from 10 GeV to1 TeV. We found the fraction of MSP spin-down energy injected to 𝑒 ± is comparable to or slightly smaller than that injected to 𝛾 rays, 𝑓 𝑒 ± (cid:46) 𝑓 𝛾 and is therefore at (cid:46)
10% level. This parameter has beenestimated in different environments, such as nearby pulsars (Hooperet al. 2017; Hooper & Linden 2018), the Galactic center (Bednarek& Sobczak 2013), individual GCs (MAGIC Collaboration et al.2019), and galaxies with low specific star formation rate (Sudohet al. 2020). Our results provide new insights into the 𝑓 𝑒 ± parameterbased on the universal properties of 𝛾 -ray-detected GCs in the MilkyWay.In summary, our analysis reveals strong evidence for IC emis-sion in Fermi -LAT GCs. This is indicative of 𝑒 ± injected by MSPshosted by such systems. Although the Fermi -LAT sensitivity for en-ergies larger than 10 GeV is not sufficiently high to claim a detection in each individual GC, we employed a universal fit method with thebin-by-bin spectra of the sample of detected objects and were ableto increase the sensitivity to the IC component. Our results also ex-plain why it is difficult to detect GCs with TeV 𝛾 -ray telescopes: wehave obtained a very soft spectra for the high energy tail of the GCpopulation. It is possible that with a more aggressive observationcampaign such objects could be detected by forthcoming TeV tele-scopes [see (Ndiyavala et al. 2018) for a recent sensitivity analysis]such as CTA (Cherenkov Telescope Array Consortium et al. 2019)and LHAASO (Bai et al. 2019). Globular clusters remain some ofthe most important systems within which to search for and studymillisecond pulsars. We have shown the potential of extracting crit-ical knowledge from 𝛾 -ray data of globular clusters with advancedstatistical tools and intensive modelling. ACKNOWLEDGEMENTS
We thank Shin’ichiro Ando and Holger Baumgardt for discus-sions. The work of D.S. is supported by the U.S. Department ofEnergy under the award number DE-SC0020262. O.M. acknowl-edges support by JSPS KAKENHI Grant Numbers JP17H04836,JP18H04340, JP18H04578, and JP20K14463. The work of S.H. issupported by the U.S. Department of Energy under the award num-ber DE-SC0020262 and NSF Grant numbers AST-1908960 andPHY-1914409. This work was supported by World Premier Interna-tional Research Center Initiative (WPI Initiative), MEXT, Japan.R.M.C. acknowledges support from the Australian Governmentthrough the Australian Research Council for grant DP190101258shared with Prof. Mark Krumholz at the ANU. D.M.N. acknowl-edges support from NASA under award Number 80NSSC19K0589.
REFERENCES
Abazajian K. N., Horiuchi S., Kaplinghat M., Keeley R. E., Macias O., 2020,Phys. Rev. D, 102, 043012Abdo A. A., et al., 2009a, Science, 325, 845Abdo A. A., et al., 2009b, ApJ, 699, 1171Abdo A. A., et al., 2010, A&A, 524, A75Abdo A. A., et al., 2013, ApJS, 208, 17Abdollahi S., et al., 2020, ApJS, 247, 33Ackermann M., et al., 2012, ApJ, 755, 164Ackermann M., et al., 2017, ApJ, 836, 208Ajello M., Di Mauro M., Paliya V. S., Garrappa S., 2020, ApJ, 894, 88Bahramian A., Heinke C. O., Sivakoff G. R., Gladstone J. C., 2013, ApJ,766, 136Bai X., et al., 2019, arXiv e-prints, p. arXiv:1905.02773Ballet J., Burnett T. H., Digel S. W., Lott B., 2020, arXiv e-prints, p.arXiv:2005.11208Bartels R., Storm E., Weniger C., Calore F., 2018, Nature Astronomy, 2,819Baumgardt H., 2017, MNRAS, 464, 2174Baumgardt H., Hilker M., 2018, MNRAS, 478, 1520Bednarek W., Sitarek J., 2007, MNRAS, 377, 920Bednarek W., Sobczak T., 2013, MNRAS, 435, L14Bednarek W., Sitarek J., Sobczak T., 2016, MNRAS, 458, 1083Blumenthal G. R., Gould R. J., 1970, Reviews of Modern Physics, 42, 237Camilo F., Rasio F. A., 2005, in Rasio F. A., Stairs I. H., eds, AstronomicalSociety of the Pacific Conference Series Vol. 328, Binary Radio Pulsars.p. 147 ( arXiv:astro-ph/0501226 )Cheng K. S., Chernyshov D. O., Dogiel V. A., Hui C. Y., Kong A. K. H.,2010, ApJ, 723, 1219Cherenkov Telescope Array Consortium et al., 2019, Science with theCherenkov Telescope Array, doi:10.1142/10986.MNRAS , 1–19 (2021) Song et al.
Di Mauro M., Calore F., Donato F., Ajello M., Latronico L., 2014, ApJ, 780,161Efron B., Petrosian V., 1999, Journal of the American Statistical Association,94, 824Espinoza C. M., et al., 2013, MNRAS, 430, 571Feigelson E. D., Nelson P. I., 1985, ApJ, 293, 192Ferraro F. R., Massari D., Dalessandro E., Lanzoni B., Origlia L., RichR. M., Mucciarelli A., 2016, ApJ, 828, 75Freire P. C. C., et al., 2011, Science, 334, 1107Gotta V., Mauro F., Moni Bidin C., Geisler D., Ferraro F., 2020, Boletin dela Asociacion Argentina de Astronomia La Plata Argentina, 61C, 90H. E. S. S. Collaboration et al., 2011, A&A, 531, L18Harding A. K., 2021, arXiv e-prints, p. arXiv:2101.05751Harding A. K., Muslimov A. G., 2011, ApJ, 743, 181Harding A. K., Usov V. V., Muslimov A. G., 2005, ApJ, 622, 531Harris W. E., 1996, AJ, 112, 1487Hooper D., Linden T., 2016, J. Cosmology Astropart. Phys., 2016, 018Hooper D., Linden T., 2018, Phys. Rev. D, 98, 043005Hooper D., Cholis I., Linden T., Fang K., 2017, Phys. Rev. D, 96, 103013Hui C. Y., Cheng K. S., Wang Y., Tam P. H. T., Kong A. K. H., ChernyshovD. O., Dogiel V. A., 2011, ApJ, 726, 100Isobe T., Feigelson E. D., Nelson P. I., 1986, ApJ, 306, 490Jóhannesson G., Porter T. A., Moskalenko I. V., 2018, ApJ, 856, 45Johnson T. J., et al., 2013, ApJ, 778, 106Kong A. K. H., Hui C. Y., Cheng K. S., 2010, ApJ, 712, L36Kopp A., Venter C., Büsching I., de Jager O. C., 2013, ApJ, 779, 126Lavalley M. P., Isobe T., Feigelson E. D., 1992, in Bulletin of the AmericanAstronomical Society. pp 839–840Lloyd S. J., Chadwick P. M., Brown A. M., 2018, MNRAS, 480, 4782MAGIC Collaboration et al., 2019, MNRAS, 484, 2876Macias O., Gordon C., Crocker R. M., Coleman B., Paterson D., HoriuchiS., Pohl M., 2018, Nature Astronomy, 2, 387Macias O., Horiuchi S., Kaplinghat M., Gordon C., Crocker R. M., NatafD. M., 2019, J. Cosmology Astropart. Phys., 2019, 042Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993Nataf D. M., et al., 2019, AJ, 158, 14Ndiyavala H., Krüger P. P., Venter C., 2018, MNRAS, 473, 897Ndiyavala H., Venter C., Johnson T. J., Harding A. K., Smith D. A., Eger P.,Kopp A., van der Walt D. J., 2019, ApJ, 880, 53Ploeg H., Gordon C., Crocker R., Macias O., 2020, J. Cosmology Astropart.Phys., 2020, 035Porter T. A., Jóhannesson G., Moskalenko I. V., 2017, ApJ, 846, 67Sollima A., Baumgardt H., 2017, MNRAS, 471, 3668Song D., Macias O., Horiuchi S., 2019, Phys. Rev. D, 99, 123020Sudoh T., Linden T., Beacom J. F., 2020, arXiv e-prints, p. arXiv:2005.08982Tam P. H. T., Kong A. K. H., Hui C. Y., Cheng K. S., Li C., Lu T. N., 2011,ApJ, 729, 90Venter C., De Jager O. C., Clapson A. C., 2009, ApJ, 696, L52Verbunt F., Kuiper L., Belloni T., Johnston H. M., de Bruyn A. G., HermsenW., van der Klis M., 1996, A&A, 311, L9Ye C. S., Kremer K., Chatterjee S., Rodriguez C. L., Rasio F. A., 2019, ApJ,877, 122Zhang P. F., Xin Y. L., Fu L., Zhou J. N., Yan J. Z., Liu Q. Z., Zhang L.,2016, MNRAS, 459, 99Zhou J. N., Zhang P. F., Huang X. Y., Li X., Liang Y. F., Fu L., Yan J. Z.,Liu Q. Z., 2015, MNRAS, 448, 3215de Menezes R., Cafardo F., Nemmen R., 2019, MNRAS, 486, 851
APPENDIX A: GLOBULAR CLUSTERS NOT DETECTEDIN THE 4FGL CATALOG
Table A1 reports the parameters and 𝛾 -ray analysis results for 127additional GCs in the Harris (1996) catalog with no counterpartdetected in the 4FGL. For fluxes and 𝐿 𝛾 , we report their 95% C.L.upper limits by placing a putative point source at the sky location of the GC and running a maximum-likelihood procedure in whichwe assume a power-law spectrum with a spectral slope of Γ = − 𝐿 𝛾 upper limits in the correlationanalysis shown in Section 3. The EM algorithm uses the 𝐿 𝛾 upperlimits to perform maximum likelihood estimates and find the best-fitparameters for the correlations with the other GC observables. Thesignificance of the correlations is estimated with the generalizedKendall 𝜏 coefficient as the test statistic, which also includes theluminosity upper limits. APPENDIX B: KENDALL 𝜏 TEST
During the MC simulations that determine the significance of cor-relations between 𝐿 𝛾 and other observables, the null hypothesissamples are generated by repetitively exchanging 𝐿 𝛾 of two GCswhile keeping their location fixed. A large number of exchanges isneeded to fully randomize the data and fulfill the requirement ofthe null hypothesis. A previous study (Ackermann et al. 2012) thatused this technique with star forming galaxies found that 1000 ex-changes were enough in their case. In this section, we test differentnumbers of exchanges and their impact on the significance of thecorrelations. Figure B1 shows the calculated significance of eachcorrelation when different numbers of exchanges are used to gen-erate null hypothesis samples. We find that the significance usuallyincreases with the number of exchanges until some large exchangenumber. This is expected as a small number of exchanges will gen-erate samples too close to the real data, which is unwanted for thenull hypothesis. The dependence of significance on the exchangesdisappears from a large number of exchanges when the samplesare fully randomized and represents the null hypothesis. We findthat about 10 exchanges are needed to generate the null hypothesissamples and achieve converged significance for four correlationsestimated. APPENDIX C: SPECTRAL ENERGY DISTRIBUTIONS
Figure C1 shows the best-fit spectra for 30 𝛾 -ray-detected GCs fromthe individual fits. We report the 𝛾 -ray fluxes of the GCs from 300MeV to 500 GeV in 9 logarithmic energy bins (black dots). Thespectra are sorted by their spectral curvature (TS curve ). For five ofthem (2MS-GC01, NGC 1904, NGC 6397, NGC 7078, and NGC5904), we report the best-fit spectra from a power law since theyhave TS curve < 4. For the rest of them, we present the best-fit spectrausing a PLE function. In general, we find the cut-off energy 𝐸 cut tobe around ∼ 𝜒 in the universal fit. The normalization factors ofthe two-component model are allowed to float. We report 19 GCswith both components detected (nonzero normalization factors).Four GCs are only fitted to the IC component. The rest of themare fitted to the CR component. For GCs with only one componentdetected, we include the 95% upper limit on the normalization ofthe non-detected component. MNRAS000
Figure C1 shows the best-fit spectra for 30 𝛾 -ray-detected GCs fromthe individual fits. We report the 𝛾 -ray fluxes of the GCs from 300MeV to 500 GeV in 9 logarithmic energy bins (black dots). Thespectra are sorted by their spectral curvature (TS curve ). For five ofthem (2MS-GC01, NGC 1904, NGC 6397, NGC 7078, and NGC5904), we report the best-fit spectra from a power law since theyhave TS curve < 4. For the rest of them, we present the best-fit spectrausing a PLE function. In general, we find the cut-off energy 𝐸 cut tobe around ∼ 𝜒 in the universal fit. The normalization factors ofthe two-component model are allowed to float. We report 19 GCswith both components detected (nonzero normalization factors).Four GCs are only fitted to the IC component. The rest of themare fitted to the CR component. For GCs with only one componentdetected, we include the 95% upper limit on the normalization ofthe non-detected component. MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters Table A1.
Parameters and data analysis results of 127 GCs not detected in the 4FGL. Notations same as Table 1.Name Γ 𝑐 [Fe/H] 𝑀 ∗ 𝑢 MW 𝑢 Total 𝑅 (cid:12) Flux 95% UL 𝐿 𝛾
95% UL(10 𝑀 (cid:12) ) (eV cm − ) (eV cm − ) (kpc) (10 − ph cm − s − ) (10 erg s − )2MS-GC02 ... -1.08 0.26 2.56 17.79 4.90 < 0.57 < 2.33AM 1 0.01 -1.70 0.32 0.25 24.57 123.30 < 0.04 < 457.85AM 4 0.00 -1.30 0.01 0.27 1.77 32.20 < 0.09 < 38.62Arp 2 0.01 -1.75 0.37 0.28 2.47 28.60 < 0.06 < 27.28BH 176 0.15 0.00 0.48 0.34 3.06 18.90 < 0.16 < 23.19BH 261 ... -1.30 0.31 3.29 11.51 6.50 < 0.07 < 1.54Djorg 1 9.04 -1.51 0.78 0.97 13.81 13.70 < 0.16 < 13.22Djorg 2 46.40 -0.65 0.56 4.01 34.01 6.30 < 0.09 < 1.72E 3 0.08 -0.83 0.03 0.46 0.99 8.10 < 0.16 < 3.92Eridanus 0.03 -1.43 0.11 0.25 28.17 90.10 < 0.07 < 271.32ESO280-SC06 ... -1.80 0.26 0.33 4.54 21.40 < 0.06 < 17.45ESO452-SC11 1.72 -1.50 0.06 1.44 9.94 8.30 < 0.10 < 4.18FSR 1735 ... ... 0.58 1.47 173.86 9.80 < 0.23 < 5.79HP 1 2.75 -1.00 1.62 4.89 6.99 8.20 < 0.17 < 4.27IC 1257 ... -1.70 0.65 0.30 8.01 25.00 < 0.21 < 44.16IC 1276 7.78 -0.75 0.72 1.58 5.88 5.40 < 0.17 < 1.47IC 4499 0.94 -1.53 1.29 0.31 15.50 18.80 < 0.09 < 15.22Ko 1 ... ... ... 0.25 2.94 48.30 < 0.03 < 68.99Ko 2 ... ... ... 0.25 1.89 34.70 < 0.02 < 32.76Liller 1 391.00 -0.33 6.61 4.38 >4.38 8.20 < 0.56 < 9.34Lynga 7 ... -1.01 1.02 1.24 17.13 8.00 < 0.06 < 2.39NGC 1261 17.90 -1.27 1.74 0.30 149.73 16.30 < 0.03 < 7.73NGC 1851 1910.00 -1.18 2.82 0.30 433.15 12.10 < 0.14 < 6.77NGC 2298 5.37 -1.92 0.54 0.30 18.54 10.80 < 0.04 < 8.44NGC 2419 3.37 -2.15 14.45 0.25 388.12 82.60 < 0.08 < 202.00NGC 288 1.23 -1.32 1.20 0.35 5.64 8.90 < 0.02 < 2.96NGC 3201 8.45 -1.59 1.45 0.52 5.73 4.90 < 0.12 < 0.90NGC 362 569.00 -1.26 3.39 0.42 184.01 8.60 < 0.06 < 3.07NGC 4147 14.90 -1.80 0.29 0.28 67.12 19.30 < 0.03 < 11.28NGC 4372 2.28 -2.17 2.19 0.65 5.13 5.80 < 0.14 < 1.62NGC 4590 4.58 -2.23 1.29 0.40 20.79 10.30 < 0.03 < 3.16NGC 4833 25.00 -1.85 2.04 0.67 17.40 6.60 < 0.10 < 1.45NGC 5024 28.50 -2.10 4.27 0.30 93.40 17.90 < 0.06 < 14.33NGC 5053 0.15 -2.27 0.62 0.31 4.20 17.40 < 0.09 < 10.59NGC 5272 167.00 -1.50 3.63 0.35 35.36 10.20 < 0.04 < 3.41NGC 5286 569.00 -1.69 3.80 0.47 308.67 11.70 < 0.13 < 6.82NGC 5466 0.18 -1.98 0.52 0.32 6.45 16.00 < 0.02 < 7.68NGC 5634 24.50 -1.88 2.00 0.29 84.71 25.20 < 0.04 < 18.05NGC 5694 205.00 -1.98 3.89 0.27 444.24 35.00 < 0.10 < 45.81NGC 5824 1220.00 -1.91 8.51 0.27 897.80 32.10 < 0.03 < 28.45NGC 5897 1.16 -1.90 1.55 0.49 10.12 12.50 < 0.05 < 5.04NGC 5927 251.00 -0.49 3.47 1.08 58.71 7.70 < 0.05 < 2.00NGC 5946 122.00 -1.29 1.15 0.85 50.13 10.60 < 0.06 < 3.34NGC 5986 56.10 -1.59 3.31 0.79 130.52 10.40 < 0.03 < 3.01NGC 6101 1.21 -1.98 1.29 0.38 28.76 15.40 < 0.04 < 7.32NGC 6121 34.50 -1.16 0.89 0.91 3.01 2.20 < 0.07 < 0.21NGC 6144 3.61 -1.76 0.48 1.14 11.98 8.90 < 0.06 < 2.38NGC 6171 9.65 -1.02 0.78 1.08 13.42 6.40 < 0.06 < 1.17NGC 6205 89.90 -1.53 4.47 0.46 48.74 7.10 < 0.04 < 1.57NGC 6229 49.90 -1.47 2.95 0.27 677.70 30.50 < 0.01 < 51.41NGC 6235 7.11 -1.28 1.12 0.92 18.12 11.50 < 0.30 < 6.36NGC 6254 42.80 -1.56 1.86 0.98 14.51 4.40 < 0.05 < 0.60NGC 6256 242.00 -1.02 1.07 1.66 53.00 10.30 < 0.11 < 3.63NGC 6273 246.00 -1.74 6.46 1.93 136.93 8.80 < 0.06 < 2.28NGC 6284 797.00 -1.26 2.40 0.56 184.37 15.30 < 0.03 < 6.50NGC 6287 52.30 -2.10 1.32 1.49 85.63 9.40 < 0.21 < 4.99NGC 6293 1220.00 -1.99 1.38 1.99 87.63 9.50 < 0.12 < 4.45NGC 6325 189.00 -1.25 0.72 2.42 82.73 7.80 < 0.13 < 1.89NGC 6333 153.00 -1.77 3.16 1.82 87.91 7.90 < 0.04 < 1.79NGC 6342 83.70 -0.55 0.60 1.82 38.20 8.50 < 0.16 < 3.06NGC 6352 12.50 -0.64 0.55 1.40 6.23 5.60 < 0.12 < 0.96NGC 6355 130.00 -1.37 1.17 2.74 117.16 9.20 < 0.13 < 3.13NGC 6356 110.00 -0.40 3.80 0.56 203.09 15.10 < 0.06 < 6.52NGC 6362 3.57 -0.99 1.07 0.75 8.27 7.60 < 0.03 < 1.69MNRAS , 1–19 (2021) Song et al.
Table A1 – continued Name Γ 𝑐 [Fe/H] 𝑀 ∗ 𝑢 MW 𝑢 Total 𝑅 (cid:12) Flux 𝐿 𝛾 (10 𝑀 (cid:12) ) (eV cm − ) (eV cm − ) (kpc) (10 − ph cm − s − ) (10 erg s − )NGC 6366 4.06 -0.59 0.59 1.04 2.26 3.50 < 0.11 < 0.68NGC 6380 96.20 -0.75 3.02 1.64 97.36 10.90 < 0.12 < 6.66NGC 6401 60.20 -1.02 2.75 1.80 22.57 10.60 < 0.11 < 3.65NGC 6426 2.68 -2.15 0.63 0.32 29.15 20.60 < 0.04 < 11.25NGC 6453 183.00 -1.50 2.34 1.57 210.77 11.60 < 0.12 < 6.25NGC 6496 2.02 -0.46 0.83 0.94 39.16 11.30 < 0.10 < 3.58NGC 6517 661.00 -1.23 3.02 1.00 419.34 10.60 < 0.22 < 7.06NGC 6522 467.00 -1.34 2.14 4.15 64.33 7.70 < 0.08 < 2.51NGC 6535 1.42 -1.79 0.13 1.11 6.87 6.80 < 0.27 < 2.11NGC 6539 271.00 -0.63 2.40 1.34 38.89 7.80 < 0.24 < 2.39NGC 6540 263.00 -1.35 0.43 2.62 >2.62 5.30 < 0.13 < 1.42NGC 6544 462.00 -1.40 1.15 1.43 22.80 3.00 < 0.13 < 0.37NGC 6553 103.00 -0.18 3.02 3.40 66.75 6.00 < 0.06 < 1.38NGC 6558 109.00 -1.32 0.39 3.23 7.50 7.40 < 0.07 < 1.53NGC 6569 72.80 -0.76 2.40 1.49 169.48 10.90 < 0.08 < 4.65NGC 6584 22.00 -1.50 1.17 0.56 117.73 13.50 < 0.13 < 7.88NGC 6624 1080.00 -0.44 0.62 2.37 79.61 7.90 < 0.07 < 2.21NGC 6626 688.00 -1.32 2.82 2.41 27.22 5.50 < 0.08 < 1.17NGC 6637 92.40 -0.64 1.48 1.78 86.29 8.80 < 0.09 < 5.98NGC 6638 103.00 -0.95 1.74 1.60 143.62 9.40 < 0.07 < 3.50NGC 6642 112.00 -1.26 0.25 1.97 47.34 8.10 < 0.04 < 2.00NGC 6656 92.40 -1.70 4.07 1.26 12.92 3.20 < 0.05 < 0.35NGC 6681 964.00 -1.62 1.12 1.33 74.61 9.00 < 0.12 < 3.64NGC 6712 33.40 -1.02 1.20 1.44 31.07 6.90 < 0.19 < 2.60NGC 6715 2030.00 -1.49 15.85 0.29 765.61 26.50 < 0.09 < 27.82NGC 6723 13.90 -1.10 1.74 1.18 31.53 8.70 < 0.11 < 3.51NGC 6749 38.50 -1.60 0.79 1.17 21.91 7.90 < 0.09 < 1.83NGC 6760 77.00 -0.40 2.57 1.11 45.55 7.40 < 0.10 < 2.01NGC 6779 25.90 -1.98 1.66 0.48 40.35 9.40 < 0.15 < 2.52NGC 6809 3.65 -1.94 1.86 1.01 7.99 5.40 < 0.10 < 1.09NGC 6864 370.00 -1.29 3.98 0.32 663.99 20.90 < 0.09 < 14.80NGC 6934 31.80 -1.47 1.41 0.34 105.49 15.60 < 0.05 < 8.08NGC 6981 4.12 -1.42 0.69 0.34 40.01 17.00 < 0.06 < 8.68NGC 7006 6.46 -1.52 1.48 0.25 316.89 41.20 < 0.07 < 47.47NGC 7089 441.00 -1.65 5.01 0.39 191.32 11.50 < 0.02 < 4.30NGC 7099 366.00 -2.27 1.26 0.50 47.69 8.10 < 0.14 < 2.71NGC 7492 0.86 -1.78 0.29 0.25 8.61 26.30 < 0.09 < 34.99Pal 1 0.05 -0.65 0.02 0.30 2.82 11.10 < 0.04 < 2.79Pal 10 75.80 -0.10 0.55 0.90 11.97 5.90 < 0.08 < 1.05Pal 11 12.40 -0.40 0.11 0.48 14.90 13.40 < 0.05 < 7.89Pal 12 0.54 -0.85 0.06 0.32 1.41 19.00 < 0.04 < 11.78Pal 13 0.00 -1.88 0.03 0.27 13.18 26.00 < 0.07 < 20.02Pal 14 0.00 -1.62 0.15 0.25 3.18 76.50 < 0.06 < 263.46Pal 15 0.02 -2.07 0.42 0.25 7.18 45.10 < 0.06 < 57.78Pal 2 457.00 -1.42 2.29 0.25 323.49 27.20 < 0.07 < 21.58Pal 3 0.03 -1.63 0.23 0.25 23.67 92.50 < 0.10 < 319.08Pal 4 0.03 -1.41 0.28 0.25 51.34 108.70 < 0.04 < 352.19Pal 5 0.01 -1.41 0.18 0.30 1.12 23.20 < 0.07 < 15.55Pal 6 40.20 -0.91 1.00 4.07 23.00 5.80 < 0.23 < 2.53Pal 8 3.03 -0.37 0.58 0.79 25.72 12.80 < 0.15 < 5.03Pyxis ... -1.20 0.23 0.25 >0.25 39.40 < 0.17 < 67.61Rup 106 0.36 -1.68 0.35 0.29 16.78 21.20 < 0.08 < 18.51Terzan 10 4430.00 -1.00 2.45 3.42 10.99 5.80 < 0.03 < 0.94Terzan 12 35.80 -0.50 0.01 2.34 6.57 4.80 < 0.12 < 0.85Terzan 3 0.89 -0.74 0.58 1.36 4.20 8.20 < 0.11 < 3.29Terzan 4 ... -1.41 0.76 4.12 5.07 7.20 < 0.23 < 1.94Terzan 6 1300.00 -0.56 0.89 4.23 298.38 6.80 < 0.13 < 2.35Terzan 7 0.62 -0.32 0.19 0.31 9.23 22.80 < 0.04 < 15.89Terzan 8 0.08 -2.16 0.54 0.29 6.49 26.30 < 0.13 < 24.96Terzan 9 2.79 -1.05 0.04 4.88 7.51 7.10 < 0.10 < 2.67Ton 2 6.45 -0.70 0.26 2.51 11.62 8.20 < 0.13 < 3.40UKS 1 100.00 -0.64 0.78 4.47 >4.47 7.80 < 0.45 < 6.28Whiting 1 ... -0.70 0.02 0.25 10.69 30.10 < 0.05 < 28.00MNRAS000
Table A1 – continued Name Γ 𝑐 [Fe/H] 𝑀 ∗ 𝑢 MW 𝑢 Total 𝑅 (cid:12) Flux 𝐿 𝛾 (10 𝑀 (cid:12) ) (eV cm − ) (eV cm − ) (kpc) (10 − ph cm − s − ) (10 erg s − )NGC 6366 4.06 -0.59 0.59 1.04 2.26 3.50 < 0.11 < 0.68NGC 6380 96.20 -0.75 3.02 1.64 97.36 10.90 < 0.12 < 6.66NGC 6401 60.20 -1.02 2.75 1.80 22.57 10.60 < 0.11 < 3.65NGC 6426 2.68 -2.15 0.63 0.32 29.15 20.60 < 0.04 < 11.25NGC 6453 183.00 -1.50 2.34 1.57 210.77 11.60 < 0.12 < 6.25NGC 6496 2.02 -0.46 0.83 0.94 39.16 11.30 < 0.10 < 3.58NGC 6517 661.00 -1.23 3.02 1.00 419.34 10.60 < 0.22 < 7.06NGC 6522 467.00 -1.34 2.14 4.15 64.33 7.70 < 0.08 < 2.51NGC 6535 1.42 -1.79 0.13 1.11 6.87 6.80 < 0.27 < 2.11NGC 6539 271.00 -0.63 2.40 1.34 38.89 7.80 < 0.24 < 2.39NGC 6540 263.00 -1.35 0.43 2.62 >2.62 5.30 < 0.13 < 1.42NGC 6544 462.00 -1.40 1.15 1.43 22.80 3.00 < 0.13 < 0.37NGC 6553 103.00 -0.18 3.02 3.40 66.75 6.00 < 0.06 < 1.38NGC 6558 109.00 -1.32 0.39 3.23 7.50 7.40 < 0.07 < 1.53NGC 6569 72.80 -0.76 2.40 1.49 169.48 10.90 < 0.08 < 4.65NGC 6584 22.00 -1.50 1.17 0.56 117.73 13.50 < 0.13 < 7.88NGC 6624 1080.00 -0.44 0.62 2.37 79.61 7.90 < 0.07 < 2.21NGC 6626 688.00 -1.32 2.82 2.41 27.22 5.50 < 0.08 < 1.17NGC 6637 92.40 -0.64 1.48 1.78 86.29 8.80 < 0.09 < 5.98NGC 6638 103.00 -0.95 1.74 1.60 143.62 9.40 < 0.07 < 3.50NGC 6642 112.00 -1.26 0.25 1.97 47.34 8.10 < 0.04 < 2.00NGC 6656 92.40 -1.70 4.07 1.26 12.92 3.20 < 0.05 < 0.35NGC 6681 964.00 -1.62 1.12 1.33 74.61 9.00 < 0.12 < 3.64NGC 6712 33.40 -1.02 1.20 1.44 31.07 6.90 < 0.19 < 2.60NGC 6715 2030.00 -1.49 15.85 0.29 765.61 26.50 < 0.09 < 27.82NGC 6723 13.90 -1.10 1.74 1.18 31.53 8.70 < 0.11 < 3.51NGC 6749 38.50 -1.60 0.79 1.17 21.91 7.90 < 0.09 < 1.83NGC 6760 77.00 -0.40 2.57 1.11 45.55 7.40 < 0.10 < 2.01NGC 6779 25.90 -1.98 1.66 0.48 40.35 9.40 < 0.15 < 2.52NGC 6809 3.65 -1.94 1.86 1.01 7.99 5.40 < 0.10 < 1.09NGC 6864 370.00 -1.29 3.98 0.32 663.99 20.90 < 0.09 < 14.80NGC 6934 31.80 -1.47 1.41 0.34 105.49 15.60 < 0.05 < 8.08NGC 6981 4.12 -1.42 0.69 0.34 40.01 17.00 < 0.06 < 8.68NGC 7006 6.46 -1.52 1.48 0.25 316.89 41.20 < 0.07 < 47.47NGC 7089 441.00 -1.65 5.01 0.39 191.32 11.50 < 0.02 < 4.30NGC 7099 366.00 -2.27 1.26 0.50 47.69 8.10 < 0.14 < 2.71NGC 7492 0.86 -1.78 0.29 0.25 8.61 26.30 < 0.09 < 34.99Pal 1 0.05 -0.65 0.02 0.30 2.82 11.10 < 0.04 < 2.79Pal 10 75.80 -0.10 0.55 0.90 11.97 5.90 < 0.08 < 1.05Pal 11 12.40 -0.40 0.11 0.48 14.90 13.40 < 0.05 < 7.89Pal 12 0.54 -0.85 0.06 0.32 1.41 19.00 < 0.04 < 11.78Pal 13 0.00 -1.88 0.03 0.27 13.18 26.00 < 0.07 < 20.02Pal 14 0.00 -1.62 0.15 0.25 3.18 76.50 < 0.06 < 263.46Pal 15 0.02 -2.07 0.42 0.25 7.18 45.10 < 0.06 < 57.78Pal 2 457.00 -1.42 2.29 0.25 323.49 27.20 < 0.07 < 21.58Pal 3 0.03 -1.63 0.23 0.25 23.67 92.50 < 0.10 < 319.08Pal 4 0.03 -1.41 0.28 0.25 51.34 108.70 < 0.04 < 352.19Pal 5 0.01 -1.41 0.18 0.30 1.12 23.20 < 0.07 < 15.55Pal 6 40.20 -0.91 1.00 4.07 23.00 5.80 < 0.23 < 2.53Pal 8 3.03 -0.37 0.58 0.79 25.72 12.80 < 0.15 < 5.03Pyxis ... -1.20 0.23 0.25 >0.25 39.40 < 0.17 < 67.61Rup 106 0.36 -1.68 0.35 0.29 16.78 21.20 < 0.08 < 18.51Terzan 10 4430.00 -1.00 2.45 3.42 10.99 5.80 < 0.03 < 0.94Terzan 12 35.80 -0.50 0.01 2.34 6.57 4.80 < 0.12 < 0.85Terzan 3 0.89 -0.74 0.58 1.36 4.20 8.20 < 0.11 < 3.29Terzan 4 ... -1.41 0.76 4.12 5.07 7.20 < 0.23 < 1.94Terzan 6 1300.00 -0.56 0.89 4.23 298.38 6.80 < 0.13 < 2.35Terzan 7 0.62 -0.32 0.19 0.31 9.23 22.80 < 0.04 < 15.89Terzan 8 0.08 -2.16 0.54 0.29 6.49 26.30 < 0.13 < 24.96Terzan 9 2.79 -1.05 0.04 4.88 7.51 7.10 < 0.10 < 2.67Ton 2 6.45 -0.70 0.26 2.51 11.62 8.20 < 0.13 < 3.40UKS 1 100.00 -0.64 0.78 4.47 >4.47 7.80 < 0.45 < 6.28Whiting 1 ... -0.70 0.02 0.25 10.69 30.10 < 0.05 < 28.00MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters . . . . . . S i g n i fi c a n ce [ σ ] L γ vs Γ c . . . . . . . . S i g n i fi c a n ce [ σ ] L γ vs u Total . . . . . . S i g n i fi c a n ce [ σ ] L γ vs [Fe/H] − . . . . . . . . S i g n i fi c a n ce [ σ ] L γ vs u MW Figure B1.
Tests of the number of 𝐿 𝛾 exchanges during the MC simulations. The significance of correlations increases with the number of exchanges untilthe number is large enough to generate fully randomized samples. For the four observables investigated, the typical number of exchanges needed to achieveconverged significance is ∼ . APPENDIX D: SYSTEMATIC UNCERTAINTY ONELECTRON/POSITRON LUMINOSITY
In Section 4.3, we discussed the implication for the 𝑒 ± injectionefficiency from modelling the universal GC spectra. We argued thatthe 𝑒 ± luminosity from the MSPs in the GCs can be approximated bythe luminosity of the IC component in our two-component model.Figure D1 plots the 𝐿 IC / 𝐿 CR against the total radiation field 𝑢 Total of the GCs (left panel) and their distance 𝑅 (cid:12) from the Sun(right panel). It is shown that the ratios of 𝐿 IC / 𝐿 CR are ∼ 𝑂 ( ) .In both cases, the 𝐿 IC / 𝐿 CR are scattered and show no obviouscorrelations with 𝑢 Total or 𝑅 (cid:12) . This supports the idea that IC is theleading energy loss process for 𝑒 ± injected by MSPs diffusing inthe GCs. The propagating 𝑒 ± eventually lose all their energy wheninteracting with the background photon field (IC) or magnetic field(synchrotron radiation). The 𝐿 IC would be suppressed for small 𝑢 Total when the synchrotron energy loss is comparable with the IC,which is not observed. For similar reason, we can rule out large uncertainties due to thepoint spread function effects of the
Fermi -LAT. In the present work,we are only considering the point-source luminosity of the GCs.For GCs closer to the Sun, the GC spatial extensions correspond tolarger angular separations. The 𝑒 ± injected by MSPs could leak andcarry away energy from the GCs. Should that happen, we may finda positive correlation between 𝐿 IC / 𝐿 CR and 𝑅 (cid:12) since the apertureswould only see partial IC from the 𝑒 ± . Within uncertainties, this isnot observed.Given the above, it is a good approximation to assume that 𝑓 𝑒 ± 𝑓 𝛾 (cid:39) 𝐿 IC 𝐿 CR , (D1)as we have done in this study. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS , 1–19 (2021) Song et al. Energy [MeV]10 E d N / d E [ e r g c m s ] Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 1904 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6397 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 707810 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 5904 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6341 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6541 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 652810 Energy [MeV]10 E d N / d E [ e r g c m s ] GLIMPSE02 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6717 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6218 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 640210 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 2808 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6139 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6838 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6304
Figure C1.
Best-fit spectra and 1 𝜎 uncertainty (blue line and shaded area) for 30 𝛾 -ray-detected GCs from the individual fits. The bin-by-bin fluxes fromFermi data analysis are include by the black dots with error bars. The spectra are sorted by their curvature. For GCs with TS curve <
4, we report their best-fitPL spectrum. For the rest, the best-fit PLE spectra are shown. MNRAS000
4, we report their best-fitPL spectrum. For the rest, the best-fit PLE spectra are shown. MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters Energy [MeV]10 E d N / d E [ e r g c m s ] M 80 10 Energy [MeV]10 E d N / d E [ e r g c m s ] Terzan 2 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6440 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 644110 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6652 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6316 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6752 10 Energy [MeV]10 E d N / d E [ e r g c m s ] Terzan 110 Energy [MeV]10 E d N / d E [ e r g c m s ] GLIMPSE01 10 Energy [MeV]10 E d N / d E [ e r g c m s ] M 62 10 Energy [MeV]10 E d N / d E [ e r g c m s ] Omega Cen 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 638810 Energy [MeV]10 E d N / d E [ e r g c m s ] Terzan 5 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 104
Figure C1 – continued MNRAS , 1–19 (2021) Song et al. Energy [MeV]10 − − − − − − E d N / d E [ e r g c m − s − ] Energy [MeV]10 − − − − − − E d N / d E [ e r g c m − s − ] NGC 1904 ICCR Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6397 CRIC Energy [MeV]10 − − − − − − E d N / d E [ e r g c m − s − ] NGC 7078 ICCR Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 5904 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6341 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6541 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6528 CRIC Energy [MeV]10 − − − − − − E d N / d E [ e r g c m − s − ] GLIMPSE02 ICCR Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6717 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6218 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6402 CRIC10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 2808 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6139 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6838 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6304 CRIC
Figure C2.
Best-fit spectra and 1 𝜎 uncertainty for 30 𝛾 -ray-detected GCs from the universal two-component fits (see section 4.2). The best-fit parameters forthe CR component (red line with shaded band) is Γ = . ± .
44 and log ( 𝐸 cut / MeV ) = . ± .
16. The best fit parameter for the IC component (blue linewith shaded band) is Γ = . ± .
25. The black dashed line shows the combined CR and IC components. When only one component detected, the 95% C.L.upper limit for the CR (IC) component is shown by the red (blue) shaded area. MNRAS000
25. The black dashed line shows the combined CR and IC components. When only one component detected, the 95% C.L.upper limit for the CR (IC) component is shown by the red (blue) shaded area. MNRAS000 , 1–19 (2021) vidence for inverse Compton emission from globular clusters Energy [MeV]10 E d N / d E [ e r g c m s ] M 80 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] Terzan 2 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6440 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6441 CRIC10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6652 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6316 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6752 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] Terzan 1 CRIC Energy [MeV]10 E d N / d E [ e r g c m s ] GLIMPSE01 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] M 62 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] Omega Cen CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 6388 CRIC10 Energy [MeV]10 E d N / d E [ e r g c m s ] Terzan 5 CRIC 10 Energy [MeV]10 E d N / d E [ e r g c m s ] NGC 104 CRIC
Figure C2 – continued MNRAS , 1–19 (2021) Song et al. u Total [eV cm − ]10 − L I C / L C R R (cid:12) [kpc]10 − L I C / L C R Figure D1.
Distribution of the ratio 𝐿 IC / 𝐿 CR compared with the total photon field energy density 𝑢 Total (left panel) and the distance from the Sun 𝑅 (cid:12) (rightpanel). The ratios are shown to be ∼ O(1) among GCs, with uncertainties presented. In both cases, the ratios show no obvious correlations with the otherparameters. MNRAS000