A Fermi-LAT Study of Globular Cluster Dynamical Evolution in the Milky Way: Millisecond Pulsars as the Probe

Li Feng, Zhongqun Cheng, Wei Wang, Zhiyuan Li,4, and Yang Chen,4

1 SchoolofAstronomyand Space Science,NanjingUniversity,Nanjing 210023, China

2 Department of Astronomy,School of Physics andTechnology,WuhanUniversity,Wuhan 430072,China;chengzq@whu.edu.cn

3 WHU-NAOC Joint Center for Astronomy, Wuhan University, Wuhan 430072, China

4 Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University, Nanjing 210023, China

Abstract Using archival Fermi-LAT data with a time span of ∼12 yr, we study the population of Millisecond Pulsars(MSPs)in Globular Clusters(GlCs)and investigate their dependence on cluster dynamical evolution in the Milky Way.We show that the γ-ray luminosity (Lγ) and emissivity (i.e., ∊γ=Lγ/M, with M the cluster mass) are good indicators of the population and abundance of MSPs in GlCs, and they are highly dependent on the dynamical evolution history of the host clusters.Specifically speaking, the dynamically older GlCs with more compact structures are more likely to have larger Lγ and ∊γ,and these trends can be summarized as strong correlations with cluster stellar encounter rate Γ and the specific encounter rate(Λ=Γ/M),with Lγ ∝Γ0.70±0.11 and ∊γ ∝Λ0.73±0.13 for dynamically normal GlCs.However, as GlCs evolve into deep core collapse, these trends are found to be reversed, implying that strong encounters may have lead to the disruption of Low-Mass X-ray Binaries and ejection of MSPs from core-collapsed systems.Besides, the GlCs are found to exhibit larger ∊γ with increasing stellar mass function slope (∊γ ∝10(0.57±0.1)α), decreasing tidal radius (∊ γ ∝Rt-1.0±0.22) and distances from the Galactic Center(GC,∊ γ ∝Rg-c1.13±0.21).These correlations indicate that,as GlCs losing kinetic energy and spiral in toward the GC, tidal stripping and mass segregation have a preference in leading to the loss of normal stars from GlCs, while MSPs are more likely to concentrate to cluster center and be deposited into the GC.Moreover, we gauge ∊γ of GlCs is ∼10–1000 times larger than the Galactic bulge, the latter is thought to reside thousands of unresolved MSPs and may be responsible for the GC γ-ray excess, which supports that GlCs are generous contributors to the population of MSPs in the GC.

Key words: (Galaxy:) globular clusters: general – (stars:) pulsars: general – Galaxy: kinematics and dynamics –Galaxy: center – Galaxy: bulge – gamma-rays: stars – gamma-rays: diffuse background – gamma-rays: galaxies

1.Introduction

Globular Clusters(GlCs)are self-gravitating systems that evolve with dynamical timescale (i.e., two-body relaxation timescale)much smaller than their age(Heggie&Hut 2003),and the binary burning process (i.e., extraction of binary binding energy via dynamical encounters) is the “heating mechanism” that supports the cluster against gravother-thermal collapse(Fregeau et al.2003),which also makes GlCs promising breeding grounds for exotic objects, such as blue straggler stars (BSS, Fregeau et al.2004;Chatterjee et al.2013), coronal active binaries(ABs),cataclysmic variables(CVs;Ivanova et al.2006;Shara&Hurley 2006;Belloni et al.2016,2017,2019),low-mass X-ray binaries(LMXBs;Rasio et al.2000; Ivanova et al.2008; Kremer et al.2018), millisecond pulsars (MSPs; the offspring of LMXBs; Ye et al.2019; Kremer et al.2020b;Ye&Fragione 2022),and gravitational-wave sources made up of compact objects(Clausen et al.2013;Rodriguez et al.2015, 2016; Askar et al.2017; Arca Sedda 2020; Ye et al.2020;Kremer et al.2021).

Depend on their distances from the Galactic Center(GC),the dynamical evolution of GlCs is also subjected to the gravitational tidal field of the Milky Way(Baumgardt&Makino 2003).The expanding cluster halo could be truncated by the external tidal field,and tidal stripping may lead to the loss of stars from GlCs, enhance the energy outflow of GlCs and thereby accelerating the dynamical evolution of the cluster (Gnedin &Ostriker 1997; Gnedin et al.1999).Besides, dynamical friction between GlCs and Galactic background stars may lead to the spiral in of GlCs toward the deep potential well of the Galaxy(Tremaine et al.1975;Moreno et al.2022),where tidal stripping would eventually lead to the fully dissolution of the clusters.GlCs therefore can take their binary burning products into the Inner Galaxy (Fragione et al.2018a; Arca-Sedda et al.2018b),and contribute the cluster mass to the growing Galactic nuclei(Antonini et al.2012; Antonini 2013; Gnedin et al.2014).

Thanks to their exclusive formation process and emission properties,the binary burning products are proved to be superb tracing particles of GlC stellar dynamical interactions and cluster evolution.For example,compared with normal stars,the emission of LMXBs, MSPs, CVs, ABs and BSS are found to be much more luminous either in optical,X-ray,γ-ray or radio band and can be detected and picked out from the dense core of GlCs,making them ideal tracers for studying stellar dynamical encounters (Verbunt & Hut 1987; Pooley et al.2003;Verbunt 2003; Pooley & Hut 2006), two-body relaxation and mass segregation effect of GlCs (Ferraro et al.2012; Cheng et al.2019b, 2019a).The cumulative cluster X-ray luminosity LXand emissivity ∊X(LXper unit stellar mass), proxies of population and abundance of weak X-ray sources(mainly CVs and ABs) in GlCs (Cheng et al.2018; Heinke et al.2020), are also shown to be effective indicators for diagnosing cluster binary burning process (Cheng et al.2018, 2020).

Among all the binary burning products in GlCs, MSPs are outstanding for many reasons.First, compared with BSS, CVs and ABs that can be formed through both primordial binary evolution and dynamical formation channels, it is widely recognized that LMXBs, the progenitors of MSPs, are dynamically formed in GlCs (Verbunt & Hut 1987; Pooley et al.2003; Verbunt 2003), and the abundances (number per unit stellar mass) of LMXBs and MSPs are more than∼100 times higher in GlCs than in the Galactic field (Clark 1975; Katz 1975; Camilo & Rasio 2005; Ransom 2008).Besides, unlike other tracers that could be contaminated by fore- and background sources, MSPs are mainly detected via either radio or γ-ray observations, in which contamination is negligible.Moreover, the lifetime of MSPs is very long(∼1010yr)and their luminosities are observed to be very stable.Taking these aspects together,it is safe to say that MSPs are the best probe of cluster dynamical evolution,especially in tracing the tidal dissolution of GlCs in the Milky Way.

In γ-ray astronomy, both MSPs and GlCs are recognized as important γ-ray emitters in the Galaxy (Abdo et al.2010; Tam et al.2011),and ∼30 GlCs have been identified as point sources in the fourth Fermi Large Area Telescope catalog (4FGL,Abdollahi et al.2020).The γ-ray emission of GlCs is widely assumed to originate from MSPs that reside in the clusters(Cheng et al.2010; Bednarek & Sobczak 2013), and the detection of pulsed gamma-ray emission from two GlCs has further strengthened the case for this connection (Freire et al.2011; Wu et al.2013).In concert with LMXBs, a strong correlation between γ-ray luminosity Lγand the stellar encounter rate Γ has been established among GlCs(Abdo et al.2010;Hui et al.2011; Hooper & Linden 2016; Tam et al.2016; Zhang et al.2016; Lloyd et al.2018; de Menezes et al.2019), which lend a support to the dynamical origin of MSPs in GlCs.

On the other hand, the Fermi Large Area Telescope (Fermi-LAT)has discovered an unexpected γ-ray excess around the GC,peaking at a few GeV, with an approximately spherical density profile ∝r-2.4and extending to 10°–20°(1.5–3 kpc)from the GC(Goodenough & Hooper 2009; Hooper & Goodenough 2011a;Hooper & Linden 2011b; Di Mauro 2021).The nature of this“Galactic Center Excess” (GCE) is still not clear, and many promising models have been proposed over the past decade,including dark matter annihilation(Abazajian&Kaplinghat 2012;Daylan et al.2016), emission of thousands of unresolved MSPs(Abazajian 2011; Yuan & Zhang 2014; Bartels et al.2016), or emission from cosmic rays injected at the GC (Carlson &Profumo 2014;Petrović et al.2014;Cholis et al.2015;Gaggero et al.2015).In the MSP scenario, Galactic GlCs fallen into the GC are expected to deposit their population of MSPs,which are then inherited by the Galactic nuclear star cluster and nuclear bulge, and contribute to the observed γ-ray excess (Bednarek &Sobczak 2013; Brandt & Kocsis 2015; Fragione et al.2018a;Arca-Sedda et al.2018b; Abbate et al.2018).

In the present work,we perform a γ-ray study of the Galactic GlCs based on the archival Fermi data.Unlike previous works that only focus on the cumulative GlC γ-ray luminosity Lγ(Abdo et al.2010; Hooper & Linden 2016; de Menezes et al.2019), we also measure the γ-ray emissivity ∊γ(Lγper unit stellar mass) of the GlCs, and explore their dependence on various cluster parameters.The cumulative cluster γ-ray luminosity offers us a chance to study the population of MSPs at the Inner Galaxy, where single MSP is hard to be detected via radio or γ-ray observation.On the other hand, as demonstrated by Cheng et al.(2018) and Heinke et al.(2020)in the X-ray band, the emissivity has been proved to be a reliable indicator of exotic objects abundance in GlCs,which is insensitive to the luminosity function (LF) and can be applied to a large cluster sample in a highly uniform fashion.More importantly, the derived GlC γ-ray emissivity can be directly compared to the stellar γ-ray emissivity of the Galactic nuclear star cluster and nuclear bulge, which are crucial to evaluating the relative abundance of putative MSPs in these environments,thus estimating the possible stellar origins and contribution to the GCE by dissolved GlCs.

The limitation of our approach is that we assume the measured cluster γ-ray luminosity is a good proxy of the population of hosted MSPs,that may not be true for all clusters,especially for GlCs contains few MSPs.However, as we will show below,although the uncertainty of small counts of MSPs in individual GlCs may introduce large scatter to our cluster sample,5In fact, the scatter of derived correlations in this paper is about an order of magnitude larger.but are unlikely to create the correlations and trends observed in this work.

The paper is organized as follows.Section 2 describes the data reduction and analysis that lead to the detection and the measurement of the γ-ray luminosity Lγand emissivity ∊γof individual GlCs.Section 3 explores the dependence of Lγand∊γon various cluster physical properties.A discussion and a summary of our results are presented in Sections 4 and 5,respectively.Throughout this work we quote 1σ errors, unless otherwise stated.

2.γ-Ray Data Analysis

2.1.Data Reduction and Analysis

We analyzed the archival Fermi-LAT data of the 157 GlCs presented in the catalog of Harris (1996).The Fermi data was observed from 2008 August 8 to 2020 November 8 (MET:239 846 401–626486405), with a time span of ∼12 yr.We used the Fermi tools release 2.0 to analyze the data, with the energy band was restricted to 100MeV–300 GeV and divided into 15 logarithmically spaced energy bins.We selected the events with source class (evclass=128, evtype=3) and filter the data with DATAQUAL>0,LATCONFIG==1.A zenith angle cut of<90°and a satellite rocking angle cut of<52°was applied to avoid contamination from the Earth limb.

The region of interest (ROI) of each target was restricted to a 14°×14° rectangular box centered on the optical center of the GlCs.We used the Make4FGLxml.py tool and the 4FGL DR2 catalog to generate the background source list within the ROI.For diffuse background modeling,we adopted the most recent Galactic interstellar emission model gll_iem_v07.fits and the isotropic spectral template iso_P8R3_SOURCE_V3_v1.txt.The instrument response function (IRF) was set as P8R3_SOURCE_V3.

All targets were investigated by means of binned likelihood analysis (gtlike tool—DRMNFB, NEWMENUIT algorithm).To search for the γ-ray emission from the GlCs, we added a putative point source at the optical center of the cluster, and used the gttsmap tool to derive the TS maps.From the TS maps,one can visually inspect the detections of unkonwn sources located within the ROI.We then removed the new sources by adding them to the background source lists and refit the data again.The spectral models of GCs were set as either a powerlaw (PL), a PL plus exponential cut-off (PL+Expcut), or a logparabola (LP) model.During the fitting, all the GlC parameters (i.e., coordinates, photon index, cut off energy and normlization) were set as free, whereas for diffuse backgrounds and point sources located within 5° of the ROI,only the normalization parameter was left free to vary.To quantify the significance of the detection, a test statistic (TS)was calculated, defined asTS= 2 (l ogL1- logL0), where L1(L0) is the maximum likelihood value of the model with(without) the putative source.The chosen criteria for detection was TS>25,corresponding to a significance slightly above 4σ.

2.2.Data Analysis Results

For the 157 GlCs tabulated in the catalog of Harris (1996),about 25% (39/157) of them are found to be γ-ray bright(TS>25)in this work,with two clusters,HP 1(TS=44.3)and Terzan 9(TS=42.1),are first identified as γ-ray emitters.Crosschecking the 39 clusters with the 4FGL-DR3 catalog(Abdollahi et al.2022), five clusters, NGC 362 (TS=16) and NGC 6304(TS = 21.7), are found to have a smaller significance than our detection threshold.Besides, Liller 1 was found to be γ-ray bright by Tam et al.(2011), while our data fitting suggests that the detection is marginal (TS = 21.6).All the 39 GlCs are located within a distance of D=15 kpc from the Earth(Figure 1(a)), and there is no evident observation bias for GlCs near the GC(Figure 1(b))and the Galactic disk(Figure 1(c)),or significant dependence on cluster metallicity (Figure 1(d)).However, it seems that GlCs are more likely to be identified as γ-ray emitters,provided that they are more massive(Figure 1(e))or have larger stellar encounter rate (Figure 1(f)).Interestingly,the γ-ray detection rate of core-collapsed GlCs(12/29)is found to be ∼2 times higher than the dynamically normal GlCs (27/128),which suggest that the dynamical evolution history is also an important factor influence the γ-ray luminosity of GlCs.

The γ-ray data analysis results are summarized in Table 1,which is segmented into two panels according to the core collapse flag in the catalog of Harris (1996).The likelihood fit parameters arranged from left to right are: cluster name,spectral models, test statistic values, PL photon index,exponential cut-off energy Ecand measured GlC energy flux fγ.Adopting the cluster distance(D)presented in the catalog of Baumgardt&Hilker(2018),we also calculated the cumulative GlC γ-ray luminosity with function Lγ=4πD2fγ.We define the GlC γ-ray emissivity as ∊γ=Lγ/M, with M being the cluster mass given in the catalog of Baumgardt & Hilker (2018).Lγand ∊γare listed in Columns (7)and (8) of Table 1.The errors are given in 1σ standard deviation level.

To check the spatial consistency between the γ-ray emission and the optical centers of GlCs,we plot in Figure 2 the TS maps of the 39 clusters.For each GlC, the image was restricted to a 5°×5° box, with the central green circle indicates the tidal radius of the cluster.All maps have evidence for good coincidence between the γ-ray emission and the GlC centroids,except for NGC 1904, where the offset between the optical center and the peak of the γ-ray emission (cayan ellipse) is∼0°.3.Besides the coincidence of spatial coordinates, it is possible that the γ-ray emission is associated with a fore/background source rather than the cluster.We estimated this probability with functionPrandom=πR t2N SROI, where Rtis the GlC tidal radius, SROIis the area of the ROI, and N is the number of γ-ray sources detected within SROI.The values of Prandomare presented in the last column of Table 1.In most clusters,the Prandomvalue is less than ∼10%,the exceptions are NGC 104 (16.9%), NGC 5139 (34.5%), NGC 6624 (18.1%),NGC 6752(32.1%)and NGC 6656(40.1%),where the random detection probability is over ten percent.However, as illustrated in Figure 2, Rtof these clusters are also found to be very large, thus the larger Prandomis more likely to result from the larger tidal radius rather than being associated with a fore/background source.

Figure 1.Histogram distributions of GlCs as a function of cluster parameters: distance from the Earth (a), distance from the Galactic center (b), distance from the Galactic disk(c),cluster metallicity(d),cluster mass(e)and the stellar encounter rate(f).The blue lines represent the γ-ray bright GlCs,while clusters without γ-ray detection are shown as green lines.

3.Statistic Relations

As discussed in Section 1, the measured GlC γ-ray luminosity (Lγ) arises from a collection of MSPs reside in the cluster.Generally speaking, there are two dominant factors may influence the number of MSPs of a given GlC: (i) the population (thus the cluster mass M) of neutron stars (NS)hosted by the cluster,and(ii)the stellar dynamical interactions that transfrom NS into LMXBs and MSPs.The second factor also refers as the binary burning process,which is thought to be the essential internal energy source of cluster dynamical evolution.In this regard, the γ-ray emissivity (∊γ=Lγ/M) is better than Lγin tracing the formation efficiency of MSPs in GlCs,and both Lγand ∊γmay serve as sensitive probe to study the dynamical evolution of GlCs.With this in mind, we examine the dependence of Lγand ∊γon various cluster physical parameters in this section.If tidal stripping have played an important effect in the dynamical evolution of GlCs,significant correlations between GlC γ-ray emission and the Galactic environment parameters are also expected.The GlC parameters are mainly taken from Baumgardt & Hilker (2018)or Harris (1996), unless the specific references are quoted.Throughout the paper, we also classified the GlCs into two subgroups (i.e., the dynamically normal and core-collapsed GlCs) according to their labels in Harris (1996).The corecollapsed GlCs are generally considered to be in the late phase of cluster dynamical evolution, and thus are dynamically more older than the normal GlCs.

3.1.Correlations with Cluster Mass

In Figure 3(a),we first explore Lγversus M for all GlCs.The GlC mass is adopted from the online catalog6https://people.smp.uq.edu.au/HolgerBaumgardt/globular/of Baumgardt &Hilker (2018), which is derived by comparing the observed cluster velocity dispersion, surface density, and stellar mass function profiles against N-body simulations.Although the scatter is substantial,a moderate positive correlation between Lγand M is observable, from which we find the Spearman’s rank coefficient r=0.454 and r=0.464 for the total and dynamically normal GlCs, with p = 0.004 and p = 0.017 for random correlation.We fit a power-law function to the correlations and obtained Lγ∝M0.53±0.17and Lγ∝M0.51±0.20for total and the dynamically normal GlCs, respectively.The best fitting functions are shown as the blue and the olive lines in Figure 3 (a),and the dotted curves mark the 95% confidence range.Compared with the dynamically normal GlCs,the dependence of Lγon M is marginal for core-collapsed GlCs, with Spearman’s rank coefficient r = 0.259 and random correlation p-value of p = 0.417.The minimum and median γ-ray luminosity of dynamically normal GlCs is measured to beLγ,min= 5.41×1033erg s-1and Lγ,med=4.1×1034erg s-1, which is comparable to that of the core-collapsed GlCs (i.e.,Lγ,min= 3.45×1033erg s-1and Lγ,med=3.22×1034erg s-1).The exception is the most luminous GlCs (i.e., Lγ≿1035erg s-1), they are dominated by dynamically normal GlCs and thus are responsible for the much larger scatter of Lγin dynamically normal GlCs than the core-collapsed ones.

The dependence of Lγon M is consistent with the work of Hooper & Linden (2016), where massive GlCs are also found to have larger γ-ray luminosities.This trend is also in agreement with Figure 1 (e), supporting a γ-ray detection preference for massive GlCs.However, the negative ∊γ-M correlation suggests that the cluster mass is not the decisive factor influencing the abundance (or formation efficiency) of MSPs in GlCs.We note that a similar ∊γ-M relationship has also been reported by Fragione et al.(2018a)in their Figure 3.

3.2.Correlations with Stellar Encounter Rates

Traditionally, there are two parameters in literature to describe the stellar encounter rate of GlCs.The first is the total stellar encounter rate, Γ ∝∫ρ2/σ, an integral of stellar encounters over the cluster volume, with ρ the stellar density and σ the cluster velocity dispersion (Verbunt & Hut 1987;Verbunt 2003).The second is the specific encounter rate,defined as Λ ≡Γ/M, and measure the chance that a particular star (or binary) undergoes dynamical encounters in GlCs(Pooley & Hut 2006).For this reason, Γ is thought to be a resonable indicator of the total amount of exotic objects that can be dynamically produced in GlCs (Pooley et al.2003;Cheng et al.2018).While Λ is considered to be highly correlated with the abundance (or the formation efficiency) of exotic objects in GlCs(Pooley&Hut 2006;Cheng et al.2018).

In Figure 4(a),we plot Lγversus Γ for each GlC.The value of Γ are adopted from Bahramian et al.(2013), which is integrated over the cluster and normalized to a value of 1000 for NGC 104 (Table 2).It is clear that Lγis highly correlated with Γ for dynamically normal GlCs.The Spearman’s correlation coefficient is r = 0.697, with a random correlation p-value of p=1.1×10-4.If core-collapsed GlCs are taken into account,this correlation is still significant,with r=0.583 and p=1.5×10-4for total GlCs.The best-fitting functions can be written as Lγ∝Γ0.71±0.11(olive curves) and Lγ∝Γ0.73±0.11(blue curves) for the dynamically normal and total GlCs, respectively.These correlations are consistent with the finding of de Menezes et al.(2019), where the possible number of MSPs(N)of GlCs is found to be proportional to Γ,with N ∝Γ0.64±0.15in a sample of 23 GlCs.

According to the Spearman’s rank correlation coefficients,it is clear that Lγdepends on Γ better than M, which suggest that stellar dynamical interaction is the fundamental factor that influence the population of MSPs in GlCs.To minimize the dependence on cluster mass, we also examine ∊γversus Λ in Figure 4 (b).The value of Λ was calculated with equation Λ=Γ/M6,with M6be the cluster mass in units of 106M⊙.Since Γ was estimated based on the V-band luminosity (Bahramian et al.2013), here we calculated M6using the V-band-based magnitude for consistency, following an empirical mass-tomagnitude relation described in Cheng et al.(2018).From Figure 4(b),it is obvious that the dynamically normal GlCs have a larger ∊γwith increasing Λ.The Spearman’s rank correlation coefficient is r = 0.616, with p=1.0×10-3for random correlation.This correlation becomes less significant when taking the core-collapsed GlCs into acount, with r = 0.368 and p=0.025 for the total GlCs.A power-law fitting yields a function of ∊γ∝Λ0.78±0.15(olive curves)for the dynamically normal GlCs.

Compared with dynamically normal GlCs, the corecollapsed GlCs are found to have much large scatter inFigure 4.The difference may arises from the rapid change of cluster structure in these systems.Once GlCs are running out of binary systems and undergo deep core collapse, their structure are not stable any more and the cluster core may suffer from gravothermal oscillations (Fregeau et al.2003), core-collapsed GlCs therefore may create extremely frequent stellar encounter environments.Indeed, some core-collapsed GlCs (e.g., NGC 6624 and NGC 7078)are proved to be clusters with the highest stellar encounter rates in Table 2,supporting a currently deep core collapse in these systems.While some clusters(e.g.,HP 1 and Terzan 1) are found to have the smallest values of Γ and Λ,which may indicate a core bounce phase after the deep core collapse.On the other hand,with extremely large Λ,it is possible for LMXBs to be dynamically disrupted during deep core collapse(Verbunt&Freire 2014),and strong encounters may lead to the ejection of MSPs from the clusters(Section 4.2).These processes may introduce uncertainties to the Lγ-Γ and ∊γ-Λ correlations for core-collapsed GlCs.

Table 1 γ-Ray Data Analysis Results of the 27 Dynamically Normal GlCs (Upper Panel) and 12 Core-collapsed GlCs (Lower Panel) Detected in this Work

Figure 2.(Continued.)

Figure 3.GlC γ-ray luminosity (a) and γ-ray emissivity (b) as a function of cluster mass.The olive and red pluses represent the dynamically normal and corecollapsed GlCs separately.The red,olive and blue texts donate the Spearman’s rank correlation coefficient of the core-collapsed,dynamically normal and the total GlC samples,respectively.The solid lines and associated dotted curves are the best fitting functions and the 95%confidence range:red for core-collapsed GlCs,olive for dynamically normal GlCs, and blue for the total sample.

3.3.Correlations with Cluster Structure Parameters

To find out the dependence of GlC γ-ray esmission on GlC dynamical evolution history, we sort the cluster structure parameters and investigate their correlations to Lγand ∊γin Figures 5 and 6.Generally speaking, MSPs are prone to be dynamically formed in the dense core of GlCs.Therefore, we first plot Lγversus GlC core radius Rc, core density ρc, core relaxation time trcand central velocity dispersion σcin the upper panels of Figure 5.No clear correlation exists for the total GlC sample, except that Lγis highly correlated with σc(r = 0.574,p=2.0×10-4) in Figure 5(d).When only the dynamically normal GlCs are considered,this correlation is still evident,with Spearman’s r = 0.672 and random correlation p-value of p=2.4×10-4.We fit this correlation with power-law function,which gives(blue lines) and(olive lines) for the total and dynamically normal GlCs,respectively.However, it seems that the cluster mass is the more fundamental parameter underlying the Lγ–σcrelation,since the dependence of ∊γon σcis less significant in Figure 6(d),and σcis highly correlated with M in GlCs (Table 3).

Unlike the linear relations reported Sections 3.1 and 3.2,the Spearman’s rank correlation coefficients suggest that the dynamically normal GlCs and core-collapsed GlCs are statistically different in Figures 5(a)–(c).As GlCs evolve to older dynamical age(i.e.,with decreasing trc)and become more compact (i.e., with smaller Rcand larger ρc), the dynamically normal GlCs are more likely to exhibit larger Lγ.Whereas for core-collapsed GlCs,this trend is reversed.When the influence of cluster mass is eliminated, this tendency still holds in Figures 6(a)–(c), where the dynamically normal GlCs are measured to have larger ∊γwith increasing dynamical age,which then reverses to smaller values as clusters evolve into core-collapsed GlCs.We suggest that these features may indicate an episodic ejection of MSPs from GlCs.Namely,during the deep core collapse and subsequent core bounce,significant number of MSPs may have been dynamically ejected from the host clusters.The implication of this effect will be addressed in Section 4.2.Here, we find the significant correlations of dynamically normal GlCs can be expressed asand.They are plotted as olive lines in the figures.

Figure 4.GlC γ-ray luminosity as a function of the stellar encounter rate Γ (a) and emissivity as a function of the specific encounter rate Λ (b).The color-coded symbols denate different types of GlCs,as in Figure 3,and the same color coded text gives the Spearman’s rank correlation coefficients:red for core-collapsed,olive for dynamically normal, and blue for the total.The olive and blue solid lines mark the best-fitting funtion of the dynamically normal and the total GlCs, and dotted curves represent the 95% fitting confidence range.

By and large, the dependence of γ-ray esmission on GlC half-mass paramters is in agreement with that on the core paramters.That is, as GlCs evolve to older dynamical ageand become more compact, GlCs are more likely to exhibit larger Lγand ∊γ.The discrepancy may arise from the core collapse,in which the cluster core structure may change dramatically during the gravothermal oscillations.Whereas for the half-mass parameters, their response to core collapse will be more moderate(Fregeau et al.2003),thus both dynamically normal and core-collapsed GlCs are found to exhibit similar correlations in Figures 5(e)–(h) and 6(e)–(h).

In the bottom panels of Figures 5 and 6, we examine the dependence of Lγand ∊γon the cluster tidal radius Rt, stellar mass function (MF) slope α, concentration parameter c and metallicity[Fe/H].Lγis found to be independent of Rt,α and c in Figures 5(i), (j) and (k).However, It is clear that GlCs with smaller Rtand larger α are more likely to have larger ∊γ(Figures 6(i) and (j)).The power-law fitting yieldsandfor the dynamically normal, core-collapsed and total GlCs,respectively.While for the ∊γ–α relation,which can be expressed as log∊γ=aα+bin Figure 6(j), with the best-ftiting parameters a=(0.57±0.10) and b=(29.24±0.07) for total GlCs,a=(0.67±0.14)and b=(29.26±0.09)for dynamically normal GlCs, and a=(0.41±0.15) and b=(29.22±0.11) for core-collapsed GlCs, respectively.

Observationally,many authors argued that LMXBs are more likely to be detected in old, metal-rich GlCs, rather than in young, metal-poor ones (Sivakoff et al.2007; Li et al.2010;Kim et al.2013).As the offspring of LMXB, Figure 5(l) also suggests that MSPs are more likely to be found in metal-rich GlCs.We fit the Lγ–[Fe/H] relation with the function logLγ=a[F e H] +b, the best-fitting parameters are a=(0.53±0.18) and b=(35.16±0.21) for dynamically normal GlCs, a=(0.47±0.14) and b=(35.09±0.17) for total GlCs, respectively.When the influence of cluster mass was eliminated,we find this tendency are still evident in Figure 6(l).The best-ftiting function log∊γ=a[F e H] +bgives parameters of a=(0.58±0.17)and b=(29.71±0.20) for dynamically normal GlCs,a=(0.54±0.18)and b=(29.88±0.26)for core-collapsed GlCs, and a=(0.53±0.13) and b=(29.73±0.13) for total GlCs, respectively.We note that a similarLγ–[Fe/H] relation was also obtained by de Menezes et al.(2019) with 23 GlCs.

Table 2 Integrated Stellar Ecounter Rate (Γ) and Specific Encounter Rate (Λ) of GlCs

3.4.Correlations with the Galactic Environment

Contrary to the contraction of the central regions, the outskirts of GlCs are driven to expand to be truncated by the tidal force of the host galaxy.The tidal radius, defined as the boundary where a GlC may lose its stars to the Galaxy, is subject to its coordinates in the Milky Way (Baumgardt &Makino 2003):.Here G is the gravitation constant,VGthe circular velocity of the Galaxy and Rgcthe distance of the cluster from the GC.The independence of Lγon Rtin Figure 5(i), and the negative ∊γ–Rtrelation in Figure 6(j), suggests that tidal stripping effect is another important factor influence the abundance of MSPs in GlCs.Indeed, this tendency also can be confirmed by the positive∊γ–α correlation, since mass segregation and tidal stripping effects have a preference for the depletion of low-mass stars from GlCs, thus tidally stripped clusters are more likely to exhibit larger stellar MF slope α (Vesperini & Heggie 1997;Baumgardt & Makino 2003).

Therefore, we also investigate the dependence of GlC γ-ray emission on the Galactic environment parameters in Figures 7 and 8.The coordinates of GlCs are adopted from the online catalog of Baumgardt et al.(2019), with X point from the Galactic center toward the Sun, Y point in the direction of Galactic rotation at the Solar position, and Z point toward the North Galactic pole,respectively.In Figures 7(a)and 8(a),both Lγand ∊γare found to decrease moderately with increasing X.However, these features may reflect an observational selection effect of GlC γ-ray detection by Fermi-LAT at the position of Sun(i.e.,with X ≈8 kpc,Y=0 kpc and Z=0 kpc),rather than a physical dependence of GlC γ-ray emission on X.Nevertheless,it seems that Lγtend to reach the maximum values near the GC in Y and Z direction (Figures 7(b) and (c)), and this tendency is still evident when the influence of cluster mass was considered in Figures 8(b)and(c).We further plot in Figures 7(d)–(f)Lγas a function of cluster distance to the GC (Rgc), the perigalactic distance (Rgc,per) and apogalactic distance (Rgc,apo) of GlC orbit in the Milky Way.Although with large scatter, the mild anticorrelations between Lγand the GlC distances to GC is in agreement with the finding of Figures 7(b) and (c), suggesting that more MSPs tend to be dynamically formed in GlCs closer to the GC.This tendency is even more significant in Figures 8(d)–(f), where the best-fitting functions can be expressed as∊γ∝Rg-c1.13±0.21,∊γ∝Rg-c0,p.7e4r±0.13and∊γ∝Rg-c1,a.2p1o±0.20for total GlCs,∊γ∝Rg-c1.43±0.28,∊γ∝Rg-c0,p.7e3r±0.16and∊γ∝Rg-c1,a.5p2o±0.27for dynamically normal GlCs, and∊γ∝Rg-c0.78±0.34,∊γ∝Rg-c1,p.0e6r±0.15and∊γ∝Rg-c0,a.8p4o±0.35for core-collapsed GlCs,respectively.

For comparison,it will be interesting to mark in Figures 8(d)–(f)the γ-ray emissivity of stars in the Inner Galaxy,since tidally disrupted GlCs are thought to have a contribution to at least part of the stars therein.With a total mass of (1.4±0.6)×109M⊙for nuclear bulge stars and (0.91±0.7)×1010M⊙for boxy bulge stars,Bartels et al.(2018)showed that the emssion profile of GCE can be better fitted by stellar mass in the boxy and nuclear bulge rather than the dark matter profiles,and the γ-ray emissivity of the Galactic boxy bulge and nuclear bulge is measured to be∊γ= (1.9 ± 0.2) ×1027erg s-1M⊙-1and∊γ= (1.1 ± 0.6) ×1027erg s-1M⊙-1separately.We mark the boxy bulge and nuclear bulge with blue and olive rectangles in Figures 8(d)–(f), with vertical ranges indicate γ-ray emissivity errors and the horizontal ranges roughly denote the extent of the boxy bulge (∼10°) and nuclear bulge (∼200 pc).We also displayed the best fit of GCE emission per unit stellar mass as red dashed horizontal line (∊γ= 1.8 ×1027erg s-1M⊙-1, Bartels et al.2018) in the figures.

Figure 5.Dependence of the γ-ray luminosity on various physical properties of GlCs.Top panels,from left to right:cluster core radius Rc,core densitylog(ρ c ),core relaxation timelog( trc ),and central velocity dispersion σc.Middle panels,from left to right:cluster half-mass–radius Rh,density inside half-mass–radiuslog(ρ h ),halfmass relaxation timelog( trh ), and central escape velocity vesc.Bottom panels, from left to right: cluster tidal radius Rt, global mass function slope α, central concentration c, and metallicity [Fe/H].All these parameters are adopted from the on line catalog of Baumgardt & Hilker (2018), exceptlog( trc ), c and [Fe/H] are adopted from Harris (1996).The color-coded symbols, solid lines and texts have the same meanings as in Figure 3.For simplicity, we have omitted the confidence curves of the best-fitting functions.

Figure 6.Dependence of the γ-ray emissivity on various physical properties of GlCs.The layout of the panels,color-coded symbols,lines and comment texts in each panel are same as in Figure 5.

From Figures 8(d)–(f), it can be seen that GlCs have a γ-ray emissivity ∼10–1000 times of higher than the Galactic boxy bulge and nuclear bulge,which indicates that tidally disrupted GlCs may enhance the γ-ray emissivity of the Inner Galaxy greatly.More importantly, as GlCs inspiral in toward the GC, it seems that normal stars are more likely to be stripped off the clusters by Galactic tidal force,while MSPs are preferentially to be deposited to the Inner Galaxy,thus resulting in a negative dependence of ∊γon Rgc, Rgc,perand Rgc,apoin GlCs.These features suggest that, in addition to the tidal stripping,mass seggregation also may plays an important role in allocating the objects that can be delivered to the GC.Therefore,an estimation of GlC contribution to the GCE must take these effects into account.We leave the implication of this problem to be addressed in Section 4.4.

We end this section by emphasizing that the GlC parameters are not independent of each other.For example, there is a significant correlation between Γ and M in GlCs (Cheng et al.2018), and massive GlCs are usually found to have larger central velocity dispersion σc,escape velocity vesc,core density ρc,and smaller half-mass(core)relaxation time trh(trc).On the other hand, significant dependencce of cluster paramerters on Rcand Rgchas been reported in literature (Djorgovski &Meylan 1994).These mutual relations suggest that the GlC parameters are highly coupled with each other, and they are closely connected with the dynamical evolution of GlCs in the Milky Way tidal field(Meylan&Heggie 1997).In Table 3,wesummarized the Spearman’s rank coefficients of Lγ, ∊γ, M and Rgcon various GlC parameters, the corresponding random correlation p-values are also listed in the table.

Table 3 Tested Correlations and Coefficients

Table 3(Continued)

4.Discussion

By exploring the dependence of Lγand ∊γon various GlC parameters in Section 3,we have established many correlations and trends between GlC dynamical evolution and their binary burning products.In particular, the meaured γ-ray emission is mainly contributed by the MSPs, which facilitates us to trace the dynamical formation of MSPs within GlCs,their migration with host clusters in the Milky Way, and the final settlement and spatial distribution surrounding the GC.As introduced in Section 1, these processes are particularly important for the MSP interpretation of the GCE.Here we discuss our findings and conclusions as follows.

4.1.GlC Dynamical Evolution and Formation of MSPs

As self-gravitating systems,the dynamical evolution of GlCs is balanced by the production of energy in the core and the outflow of energy from the cluster, and two-body relaxation is the fundamental process that dominates the transportation of energy and mass in GlCs.Stars are driven to reach a state of energy equipartition by two-body relaxation, massive stars (or binaries) therefore tend to lose energy and drop to the cluster center (thus influences the types of binary burning products),whereas lower-mass stars tend to gain energy and move faster,and they will migrate outward (thus leading to the energy outflow in GlCs) and drives the cluster envelope to expand(Heggie & Hut 2003).Binary burning is thought to be the internal energy source that supports the cluster against gravothermal collapse, and the energy production rate of GlCs is sensitive to the encounter rate of the hard binaries (Cheng et al.2018):

where fbis the fraction of stars in binary, Ab∝a/σ2is the binary encounter cross section, with a being the orbital separation of the binary.In practice, it is hard to determine the distribution of a by observation, while fbis found to vary not too much among GlCs (i.e., fb∼1%–20%;Milone et al.2012).Therefore, Γbis usually simplified as Γ ∝∫ρ2/σ in literature.

As binaries are disrupted or been driven to evolve to harder systems by dynamical ecounters, both fband Abwill become smaller and the energy production rate of GlCs may derease according to Equation (1).However, this may not happen in a real cluster.In fact, GlCs contract their cores to increase the stellar density, which enhances the stellar encounter rates (i.e.,Γ and Λ)and thus the energy production rate.Accordingly,the cluster two-body relaxation time was adjusted to smaller values, to enhance the energy outflow rate and maintain the balance between the energy production in the cores and the outflow of energy from the systems.As a result,the dynamical evolution of GlCs will become more and more rapid, untill dynamical ejection of close binaries become important (see Section 4.2) and the cluster evolve into a core-collapsed GlC.

With larger mass than normal stars, NS are expected to concentrate to cluster center and take part in the binary burning process,thus leading to the formation of LMXBs and MSPs in GlCs.However, recent N-body simulations suggest that many GlCs may host significant population of primordial stellar mass black holes (BHs), even evolve to an age larger than ∼12 Gyr(Breen&Heggie 2013;Morscher et al.2013,2015;Wang et al.2016).Compared with NS, the BHs are expected to drop to cluster center first and form a high-density BH subsystem(BHS), where the BH burning process (i.e., the dynamical hardening of BH binaries, Kremer et al.2020b) may dominate the energy production of GlCs and lead to the formation of gravitational wave sources (Rodriguez et al.2016;Askar et al.2017; Anagnostou et al.2020; Antonini & Gieles 2020).The thermal coupling of BHS with GlCs is sufficient to reshape the structure of the cluster (Merritt et al.2004; Mackey et al.2007, 2008; Giersz et al.2019) and quench or slow down the mass segregation of less-massive objects (such as NS, binaries and BSS) in GlCs (Fragione et al.2018b; Weatherford et al.2018).Besides, through exchange encounters involving a BH,binaries are preferentially been transformed into BH binaries,BHS therefore may suppress the formation of other binary burning products in GlCs, and clusters with a large number of BHs are more likely to have a lower formation efficiency of LMXBs,MSPs,CVs,ABs,and BSS(Fragione et al.2018b;Ye et al.2019; Kremer et al.2020a).

The mass segregation delay of heavy objects in GlCs has been confirmed by the radial distribution study of X-ray sources,MSPs and BBS in 47 Tuc and Terzan 5(Ferraro et al.2012; Cheng et al.2019b, 2019a), where massive objects (i.e.,MSPs, Bright X-ray sources or BSS) are found to drop to cluster center earlier and be more centrally concentrated than the less massive ones (i.e., Faint X-ray sources and normal stars).The exception is M28, Cheng et al.(2020) reported an abnormal deficiency of X-ray sources in the cluster central region with respect to its outskirts, and argued that an early phase of primordial binary disruption by BHS may have happened in M28 (Cheng et al.2020).On the other hand,significant evidence of BH burning has been reported in GlC ω Cenaturi,in which the dynamical formation of X-ray sources is found to be highly suppressed,and BHS is the essential internal energy source that supports the cluster against collapse(Cheng et al.2020).Indeed, with ∼4.6% of the cluster mass being invisible, ω Cenaturi is thought to be the cluster with the heaviest BHS in the Galaxy (Baumgardt et al.2019).In Table 1,the γ-ray emissivity of ω Cenaturi is also the smallest among the 39 GlCs, suggesting a low formation efficiency of MSPs in this cluster.

Figure 7.Dependence of GlC γ-ray luminosity on Galacitic environment parameters.Top panels for cluster coordonates:with X point from GC toward the Sun(a),Y in direction of the Solar motion (b) and Z point toward the North Galactic pole (c).Bottom panels for the cluster orbital parameters in the Milky Way:Rg c = X 2 + Y 2 +Z2 represents the current GlC distance to the GC(d),while Rgc,per and Rgc,apo mark the perigalactic(e)and apogalactic(f)distances of the GlC orbit to the GC.The color-coded symbols and texts have the same meanings as in Figure 3.

The evolution fate of the BHS is beyond the scope of this paper.It is very likely that the majority of the BHs will be ejected from the cluster gradually, and host GlCs tend to contract their cores and increase the stellar density(Arca Sedda et al.2018a).As a consequence, NS starts to concentrate to cluster center and gradually take over the binary burning process, which then results in the formation a large number of LMXBs and MSPs in GlCs.Indeed,we have demonstrated that Lγ, the proxy of the population of MSPs in GlCs, is highly correlated with the cluter stellar encounter rate,7According to Equation (1), the dynamical encounters between invisible compact objects are not included in Γ,thus the energy input of BH burning are also not counted here.Γ, with Lγ∝Γ0.71±0.11in dynamically normal GlCs (Figure 4(a)).Besides, the GlC γ-ray emissivity ∊γ, a proxy of the MSP abundance of GlCs, is highly correlated with the cluster specific encounter rate Λ, with ∊γ∝Λ0.78±0.15in Figure 4(b).These correlations provide strong evidence for the dynamical formation of MSPs in GlCs.

On the other hand,by examing the dependence of Lγand ∊γon the cluster structure paramters, we have tested the relationship between MSP formation and the dynamcial evolution of GlCs.For example,both Lγand ∊γare anti-correlated with GlC core radius Rcand half-mass–radius Rhin dynamically normal GlCs (Figures 5 and 6).In particular, the ∊γ–Rcand ∊γ–Rhcorrelations can be expressed asand∊γ∝in Figures 6(a) and (e), respectively.Furthermore, Γ is strongly dependent on the cluster stellar density ρ,with Γ ∝∫ρ2/σ.We have also found strong correlations between Lγand ρcand ρhin dynamically normal GlCs, withandin Figures 5(b) and (f).These trends are in agreement with the dynamical evolution of GlCs, that is, as GlCs contract their central regions to smaller radius and become more dense,more MSPs are expected to be dynamically formed in the clusters.

Figure 8.Dependence of GlC γ-ray esmissivity on Galacitic environment parameters.The layout of the panels are same as in Figure 7.The color-coded symbols and texts have the same meanings as in Figure 3.The measured γ-ray emissivity of the Galactic boxy bulge and nuclear bulge(Bartels et al.2018)are shown as blue and olive rectangles in the botom panels, while the red dashed horizontal line represents the best fitting γ-ray esmissivity of the Galactic boxy bulge and nuclear bulge.

Finally, it is necessary to emphasize that the dynamical formation of MSPs is far from finished in many dynamically normal GlCs.As confirmed by the radial distribution studies of X-ray sources (Cheng et al.2019b, 2019a, 2020) and BSS(Ferraro et al.2012) in GlCs, it is possible that there are considerable numbers of primordial binaries and NS exist in the outskirts of the dynamically young GlCs.Under the effect of mass segregation, these objects tend to drop to the dense core of the cluster, where they could be transformed into LMXBs and MSPs by the binary burning process, thereby increasing the abundance of MSPs in dynamically normal GlCs.As a consequence, dynamically normal GlCs with an older dynamical age(i.e.,td∝ortd∝)are more likely to exhibit a larger abundance of MSPs, with the ∊γ-trccorrelation can be written asin Figure 6(c), and the ∊γ–trhcorrelation can be expressed asin Figure 6(g),respectively.

4.2.Dynamical Disruption of LMXBs and Ejection of MSPs in Core-Collapsed GlCs

From Equation(1),it can be seen that the evolution of GlCs to core collapse is sensitively depends on the encounter cross section (i.e., Ab∝a/σ2) of the binaries.With larger Ab, the primordial binaries are expected to take part in the binary burning process first and are been transformed into close binaries, which may lead to the formation of many exotic objects, such as LMXBs (e.g., through exchange encounters involving NS), MSPs, CVs, ABs and BSS in GCs.However,when primordial binaries were exhuasted,GlCs tend to contract their cores and increase the stellar density, thereby increasing the specific encounter rate (Λ) of stars in the cluster core.The enhancement of Λ is helpful for the extraction energy from harder binaries, and even very close binary systems, such as LMXBs, MSPs, CVs and ABs, may suffer strong dynamical encounters in core-collapsed GlCs.Nevertheless, it is argured that the “burning” of very hard binaries is not sufficient to terminate the deep core collapse of GlCs.The cluster cores are expected to contract further, untill dynamical interactions between single stars take effect and lead to the formation of new binaries in GlCs(i.e.,binary formation either via the“twobody binaries”or the“three-body binaries”channels,Heggie&Hut 2003).The sudden introduction of a large number of binaries in extremely dense core is sufficient to reverse the cluster core collapse to core bounce, and trigger gravothermal oscillations in GlCs.

Before the discussion of the dynamical evolution fate of LMXBs and MSPs in core-collapsed GlCs,it is necessary to have an intuitive understanding on the dynamic effects of binary burning process of GlCs.Accoding to the Hills–Heggie law,through binary–single encounters, hard binaries tend to increase their binding energy (Eb=Gm2/2a) and become harder.The average increasement of Ebper encounter is roughly proportional to its inital value, with δEb≃0.4–0.5 Ebfor equal-mass encounters(Heggie 1975;Hills 1975).The release of the binding energy is transformed into the kinetic energy of the ejected star,and the remaining binay will receive a recoil velocity as well,with.Therefore, hard binaries tend to receive larger recoil velocity after strong encounters, which may lead to the recoil of the binaries out of the core into the halo,and sometimes out of the host cluster (Hut et al.1992).

With typical central escape velocity of 10 km s-1≾vesc≾100 km s-1, LMXBs and MSPs may receive a large recoil velocity and escape the core-collapsed GlCs after suffering strong encounters.On the other hand, the specific encounter rate Λ of deep core collapse GlCs is extremely high,it is possible for LMXBs to be dynamically disrupted (i.e., via the binary-binary encounters) or be greatly modified (i.e.,exchange encounters, etc.) before they otherwise could evolve into MSPs (Verbunt & Freire 2014).The net effect is that the population of MSPs tend to be gradually exhuasted in corecollapsed GlCs.Indeed, compared with dynamically normal GlCs, the core-collapsed GlCs are found to exhibit smaller Lγand ∊γwith decreasing Rc,trc,and increasing ρcin Figures 5(a)–(c) and 6(a)–(c).Two well-known core-collapsed GCs, NGC 6397 and NGC 6752, are proved to be the clusters with the smallest Lγand ∊γin Table 1.

The decreases of MSP populations in core-collapsed GlCs have also been reported by de Menezes et al.(2019),which was interpreted as the disruption of LMXBs by the extremely large specific encounter rates in the dense cores.However, these authors have neglected the importance of MSP ejection in deep core collapse clusters.The lifetime of MSPs (∼1010yr) is at least one order of magnitudes larger than the cluster core relaxation time (trc∼105–9yr), and once MSPs are formed via the recycling process of LMXBs, it is hard to change their properties through direct stellar dynamical interactions.Therefore,the decline of Lγand ∊γfrom dynamically normal GlCs to core-collapsed GlCs in Figures 5(a)–(c) and 6(a)–(c) may suggests an episodic ejection of MSPs during the deep core collapse phase, and the injection of MSPs into the Milky Way may plays an important role in shaping the spatial extent of the GCE.Nevertheless, it is also possible for some MSPs to be retained by the host cluster (Hut et al.1992), and the retention of these MSPs may be responsible for the ∼2 times larger detection rate of γ-ray emission in core-collapsed GlCs than in dynamically normal GlCs by Fermi-LAT (Section 2.2).

4.3.Tidal Stripping and Deposition of MSPs into the Galactic Center

In addition to the internal effects,the structure of GlCs is also subjected to the external influence of the host galaxy,which may affect the dynamical evolution of GlCs in two ways.First, the expanding envelopes of GlCs are expected to be truncated by the gravitational tidal field of the host galaxy, and compared with clusters in isolation, GlCs in the tidal field may suffer an enhanced rate of energy outflow and mass loss.The tidal stripping therefore has a net effect in accelerating the dynamical evolution of the GlCs.Second, through dynamical friction with background stars,GlCs tend to lose energy and spiral in toward the center of the galaxy, where tidal stripping and interactions with the dense nucleus will eventually lead to the dissolution of the clusters (Tremaine et al.1975).In particular, during the passages close to the Galactic bulge or passing through the dense disk,GlCs may suffer gravitational shocks and these two effects could be enhanced greatly.The bulge shocking and disk shocking tend to heat up the GlC outskirts and increase evaporation of the cluster (Gnedin & Ostriker 1997), which in turn accelerate the core collapse and shorten the destruction time of GlCs(Gnedin et al.1999).Taking these aspects together,it is suggested that GlCs are the ancient building blocks of our Galaxy,and the inner galactic structures,such as the nuclear star cluster and the nuclear bulge, are assembled at least in part, by the merger of tidally disrupted GlCs that spiraled into the deep gravitational well of the GC(Antonini et al.2012;Antonini 2013;Gnedin et al.2014).

Observationally, the dynamical evolution of GlCs and its coupling with host galaxy has been widely confirmed in literature, either via the stellar streams of tidal stripping (Leon et al.2000; Sollima 2020), or the dependence of cluster structure parameters on GlC positions within the Galaxy(Djorgovski & Meylan 1994).As shown in Figure 9, with decreasing distance from the GC, GlCs are found to be more concentrated and have smaller and denser cores.These trends are consistent with the dynamical evolution of GlCs in the Galactic environment, namley, as GlCs udergo orbital decay and spiral in toward the GC,the clusters tend to suffer stronger tidal stripping and will be more dynamically evolved.Besides,it is proposed that mass segregation and tidal stripping effects have a preference in evaporating the low-mass stars from GlCs,which is strong enough to turn cluster initially increasing MF into MF that decrease toward the low-mass end (Vesperini &Heggie 1997;Baumgardt&Makino 2003;Lamers et al.2013),and the global MF slope α is strongly anti-correlated with the half-mass relaxation time trhin GlCs (Baumgardt & Sollima 2017; Sollima & Baumgardt 2017).

The relationship of GlC dynamical evolution and its coupling with the Galactic environment can also be confirmed by our γ-ray data.For example,we showed that there is a mild anti-correlation between GlC γ-ray luminosity Lγand their current distances Rgc, perigalactic distances Rgc,perand apogalactic distances Rgc,apofrom the GC in Figure 7.Compared with the strong dependence of Lγon Γ in Figure 4(a), these correlations suggest that tidal stripping is a secondary factor affecting the binary burning process of GlCs and the populations of MSPs formed therein.On the other hand, tidal stripping plays a crucial role in diminishing the mass of the GlCs and leading to the final dissolution of the clusters, which can be further confirmed by the much stronger anti-correlations between GlC γ-ray emissivity ∊γand Rgc,Rgc,perand Rgc,apoin Figure 8.In fact,we also showed that there is no evident dependence of Lγon GlC tidal radius Rtin Figure 5(i),but ∊γis found to strongly anti-correlated with Rtin Figure 6(i).Therefore, the loss of cluster mass is the most plausible reason for these anti-correlations.

Finally, it is necessary to emphasize the importance of mass segregation in affecting the tidal stripping and dissolution of GlCs.Unlike the low-mass normal stars that tend to migrate outward in GlCs and will be stripped off the clusters preferentially, LMXBs and MSPs are more likely to segregate to GlC center and could be retained by the clusters for much longer time.As a result, GlCs with increasing MS slope α are found to exhibit larger ∊γin Figure 6(j),and as GlCs udergo orbital decay and spiral in toward the Inner Galaxy, LMXBs and MSPs are more likely to be delivered to the GC, leading to a strong anti-correlations between∊γand Rgc, Rgc,perand Rgc,apoin Figure 8.

4.4.Contribution of GlC MSPs to the GCE

As introduced in Section 1, it is proposed that the GCE may arising from a large number (∼103–104) of unresolved MSPs residing in the GC (Abazajian 2011; Yuan & Zhang 2014;Bartels et al.2016).The generation of MSPs could be arise from the in situ star formation and evolution at GC (Yuan &Zhang 2014; Eckner et al.2018; Gautam et al.2022), or inherited from GlCs that were brought to the Inner Galaxy by dynamical friction and assembled the Galactic nuclear star cluster and bulge (Bednarek & Sobczak 2013; Brandt &Kocsis 2015; Fragione et al.2018a; Arca-Sedda et al.2018b;Abbate et al.2018).Compared with the dynamical channel within GlCs,the in situ formation of MSPs at GC is disfavoured because of LMXBs,the progenitors of MSPs,are observed to be quite rare in the bulge of our Galaxy than the prediction of the MSP interpretation (Cholis et al.2015; Haggard et al.2017).Besides, Boodram & Heinke (2022) argued that the natal velocity kicks received by newly formed NS during supernova may lead to a lower abundance of MSPs(thus a lower synthetic γ-ray surface brightness) in the central degrees (R ≾150 pc) of the Galatic bulge with respect to its outskirts, and which is inconsistent with the measured γ-ray surface brightness profile of the GCE (Boodram & Heinke 2022).

On the other hand,we have demonstrated that the dynamical channel of MSPs are very efficient during cluster dynamical evolution, and as GlCs losing kinetic energy and sipiral in toward the GC,MSPs are preferentially to be deposited into the the deep gravitational potential well of the Galaxy,8In fact, it can be seen from Figure 8(e) that the GlC MSPs could be effectively delivered to a distance of R ∼100 pc from the Galactic center.which may avoid the problems reported in the work of Boodram&Heinke(2022).The spatial distribution of GlCs are found to be spherically distributed around the GC (Baumgardt &Hilker 2018), and their velocity dispersion is highly isotropic within the central Galaxy (i.e., R ≾10 kpc;Vasiliev 2019).Therefore, the tidal stripping and ejection of MSPs from GlCs to the Milky Way are also expected to create a spherically distributed MSP population in the Inner Galaxy,and the radial distribution profile resambles that of the GCE (Brandt &Kocsis 2015; Fragione et al.2018a; Arca-Sedda et al.2018b;Abbate et al.2018).Moreover,considering the lifetime(∼1010yr)of MSPs is much longer than that(∼108yr)of LMXBs,the very small LMXBs to MSPs ratio observed in the Galatic bulge can be interpreted as the dynamical disruption of LMXBs during the deep core collapse phase of cluster evolution(Section 4.2), and the sudden interruption of the dynamical formation of LMXBs since host GlCs were tidally disrupted(Brandt & Kocsis 2015).

Figure 9.GlC parameters as a function of cluster distance from GC.All the 157 GlCs are shown in this figure,with olive and red pluses representing the dynamically normal and core-collapsed GlCs separately.The red,olive and blue texts donate the Spearman’s rank correlation coeficients of the core-collapsed,dynamically normal and the total GlC samples, respectively.

To estimate of the contribution of Galactic GlC MSPs to the GCE, the knowledge of the luminosity function (LF) of MSPs is particularly important.By utilizing the cluster stellar encounter rate Γ to estimate the formation rate of MSPs in GlCs,Hooper&Linden(2016)constrain the γ-ray LF of MSPs in GlCs and found that they are as luminous as the Galactic field MSPs,with log-normal LF of L0≃8.8×1033erg s-1and σ ≃0.62.Since MSPs would spin down quickly and fade away in γ-ray luminosity, Hooper &Linden (2016)therefore argued that the GlC MSPs can account for only a few percent or less of the observed GCE.Nevertheless, the LF of Hooper & Linden(2016)also predict that most of the GlCs will be dominated by only one or two MSPs in γ-ray luminosity, which is inconsistent with the radio survey result of MSPs in GlCs,where many clusters detected by Fermi are found to host far more than one MSPs.9http://www.naic.edu/～pfreire/GCpsr.htmlOn the other hand, by using the cluster∊γto scale the stellar mass deposited by GlCs at GC,Brandt&Kocsis (2015) showed that a reasonable disrupted GlC mass,such as that calculated in Gnedin et al.(2014),may account for the total γ-ray emission of the GCE.

However, both the works of Hooper & Linden (2016) and Brandt & Kocsis (2015) have neglected the influence of cluster dynamical evolution on the formation rate of MSPs in GlCs.As dicussed in Sections 4.1 and 4.2, the dynamical formation of MSPs is far from finished in dynamically young GlCs,while the population of MSPs in core-collapsed GlCs is underestimated,since strong encounters may have lead to the ejection of many MSPs from these systems.These effects can also be confirmed by the large variance of measured ∊γin GlCs.From Figures 8(e)–(f), it can be seen that ∊γof GlCs is 1–3 orders of magnitude larger than that of the Galactic bulge stars,the large variance of∊γsuggests that the final population of MSPs created by individual GlCs could be more than tens times higher than that currently hosted in the cluster.Indeed, by assuming a cluster ∊γof log∊γ= 32.66 ± 0.06 - (0.63 ±0.11) logMfor GlCs,Fragione et al.(2018a) simulating the spirial in and tidal stripping of GlCs in the Milky Way.Although their model has neglected the essential features of cluster dynamical evolution and stellar dynamics (i.e., internal energy sources, stellar mass segregation, core collapse, etc.) within GlCs, the empirical∊γ-M relation of Fragione et al.(2018a) is consistent with our results in Figure 3(b), and as GlCs spiral in toward GC and get less massive, the ∊γof remaining cluster debris will become larger and larger, as displayed in Figures 8(e)–(f).More importantly, Fragione et al.(2018a) showed that the cluster ∊γderived γ-ray luminosity is about one order of magnitude larger than the observed GCE,and spin down of MSPs will reproduce a γ-ray luminosity consistent with the GCE.

5.Summary

We present a γ-ray study of the 157 Milky Way GlCs based on the archival Fermi-LAT data with a time span of ∼12 yr.By examing the dependence of GlC γ-ray luminosity (Lγ) and emissivity(∊γ=Lγ/M)on various cluster parameters,our main findings are as follows.

1.39 GlCs are found to be γ-ray bright(i.e.,with TS>25)in our work, which corresponds to about 25% (39/157) of the total Milky Way GlCs.The detection rate of core-collapsed GlCs (12/29) is ∼2 times higher than the dynamically normal GlCs (27/128).

2.Compared with cluster mass (M), the GlC γ-ray emission is highly correlated with the stellar encounter rate (Γ) and the specific encounter rate (Λ=Γ/M), with Lγ∝Γ0.71±0.11and∊γ∝Λ0.78±0.15in dynamically normal GlCs.These correlations provide strong evidence for the dynamical formation of MSPs in GlCs.

3.The formation of MSPs is also highly dependent on the dynamical evolution history of the host GlCs.As GlCs evolve to older dynamcial age (i.e.,td∝tr-c1ortd∝tr-h1) and become more compact (i.e., with smaller Rcor Rh,and larger ρcor ρh),the dynamically normal GlCs are expected to produce more MSPs, thereby exhibit larger Lγand ∊γin Figures 5 and 6.

4.However, compared with dynamically normal GlCs, the core-collapsed GlCs are found to exhibit decreasing Lγand ∊γas their cores undergo deep core collapse.This feature may imply that even LMXBs could be dynamically disrupted or be greatly modified in extremely dense cluster cores, and strong encounters may lead to the ejection of MSPs from corecollapsed GlCs.

5.With decreasing tidal radius(Rt)and distance(Rgc)from the GC, GlCs are found to have increasing γ-ray emissivity, with∊γ∝Rt-1.0±0.22and∊γ∝Rg-c1.13±0.21.Besides, ∊γis positively correlated with GlC stellar mass function slope α, with∊γ∝10(0.57±0.1)α.These correlations suggest that both tidal stripping and mass segregation effects are improtant factors influencing the abundance of MSPs in GlCs, and as GlCs undergo orbital decay and spiral in toward the GC, MSPs are more likely to be deposited into the GC rather than normal stars.

6.We gauge the cluster ∊γis about 1–3 orders of magnitude larger than that of the Galactic bulge stars, which implies that GlCs may enhance the γ-ray emissivity of the Galactic bulge greatly.The large variance of cluster ∊γmay arise from the ongoing dynamical formation (or ejection) of MSPs in GlCs,and the different degree of tidal stripping among the clusters.More importantly, the ∊γrelations derived in our paper are in agreement with the empirical relation adopted in the simulation of Fragione et al.(2018a), which states that tidally disrupted GlCs may provide a natural astrophysical explanation to the observed GCE.

Acknowledgments

This work is supported by the Youth Program of the National Natural Science Foundation of China No.12 003 017.