EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 708 (April 2026) A&A, 708 (2026) A89 Full HTML
Open Access
Issue
A&A
Volume 708, April 2026
Article Number A89
Number of page(s) 10
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202557963
Published online 31 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

The Galactic Centre (GalC) is a region of intense astrophysical activity, hosting the supermassive black hole (SMBH) Sagittarius A* (Ghez et al. 2008; Gillessen et al. 2017), a dense nuclear stellar cluster (Neumayer et al. 2020), complex interstellar gas and radiation fields, and intense cosmic-ray activity, as reviewed in Genzel et al. (2010). It is expected to be the brightest source of diffuse γ-ray emission in the sky, primarily arising from interactions of high-energy cosmic rays with the interstellar medium (ISM) and radiation fields. However, observations have revealed an anomalous surplus of γ-ray emission in the GeV energy range that exceeds predictions from conventional astrophysical models. This signal, known as the GalC GeV Excess (GCE), was first identified in 2009 using the initial year of data from the Fermi γ-ray Space Telescope’s Large Area Telescope (Fermi-LAT), an instrument with wide field-of-view, high angular resolution, and energy coverage from ≈100 MeV to >300 GeV (Goodenough & Hooper 2009; Hooper & Goodenough 2011). The excess has since been robustly confirmed by independent analyses, including those from the Fermi-LAT Collaboration, using over a decade of accumulated data, establishing it as a spatially extended, spectrally peaked feature accounting for 2% of the total γ-ray flux within a 30◦ radius of the GalC (Ackermann & et al. 2017; Di Mauro 2021). The signal peaks at a few GeV, is approximately spherically symmetric about Sgr A* and is detectable out to ≳10–15◦ (a few kpc) from the centre (Goodenough & Hooper 2009; Hooper & Goodenough 2011; Hooper & Linden 2018; Di Mauro 2021; Cholis et al. 2022).

The origin of the GCE remains a subject of intense debate, with interpretations being broadly divided between exotic physics and astrophysical sources. One leading hypothesis attributes the excess to the annihilation of weakly interacting massive particles (WIMPs), a candidate for dark matter (DM), concentrated in the GalC due to the cuspy density profile of the Milky Way’s (MW) DM halo. This model naturally explains the GCE’s bump-like spectrum (peaking at 1–4 GeV), spherical morphology (consistent with a Navarro-Frenk-White profile with inner slope γ ≈ 1.2), and intensity, which aligns with the thermal relic annihilation cross-section of ⟨σv⟩ ≈ 10−26 cm3 s−1 (Hooper & Goodenough 2011; Daylan et al. 2016; Di Mauro 2021; Cholis et al. 2022). Alternatively, a more conservative explanation posits that the GCE arises from a large, unresolved population of millisecond pulsars (MSPs) in the Galactic bulge, whose collective γ-ray emission mimics the observed spectrum due to a power-law with exponential cut-off at a few GeV, matching the average spectral properties of known Fermi-LAT detected MSPs (Bartels et al. 2016; Lee et al. 2016; Gautam et al. 2022; Holst & Hooper 2025). Eckner et al. (2018) showed that the GCE is not unique to the MW: a similar extended γ-ray emission is also observed in M31. The authors demonstrated that a population of MSPs forming in situ in bulges, with their luminosity scaled by stellar mass, can reproduce both the energetics and the morphology of the observed signal.

A prominent astrophysical route to a bulge population of MSPs invokes the long-term dynamical evolution and disruption of globular clusters (GCs). Early works have shown that a large fraction of the MW’s GC system is destroyed over Hubble time, depositing stars and compact remnants into the bulge (Tremaine et al. 1975; Gnedin & Ostriker 1997; Gnedin et al. 2014). Building on this, Brandt & Kocsis (2015) argued that MSPs delivered by disrupted clusters can quantitatively account for the GCE’s amplitude, spectrum, and approximately spherical morphology.

Direct N-body models that follow infalling clusters with internal mass segregation and tidal stripping support this picture, linking the high-energy emission to compact objects brought in by cluster inspiral and dissolution (Arca-Sedda & Capuzzo-Dolcetta 2014; Arca-Sedda et al. 2018). Independent constraints on the MSP yield of GCs, from population analyses of cluster pulsars (Turk & Lorimer 2013) to recent FAST discoveries of multiple new MSPs in a single GC, underscore that dense clusters can host substantial MSP reservoirs (Turk & Lorimer 2013; Yin et al. 2024). Fragione et al. (2018) performed direct N-body simulations (NBODY7 + SSE + BSE) of GCs evolution with varying initial masses, metallicities, and primordial binary fractions, tracing the formation and spatial distribution of their neutron stars (NSs) populations over 10 Gyr. They showed that most NSs escape from their parent clusters due to the natal kicks received during Type II supernova explosions, as the escape velocity of the clusters is much lower than the typical kick velocity. Assuming 10% of the NSs are recycled into MSPs, the authors demonstrated that the resulting γ-ray luminosity of the GCs matches Fermi observations, supporting the MSP origin of the GCE. Taken together, these results motivate a fresh, fully dynamical assessment of whether GC-borne NSs can indeed reproduce the detailed properties of the GCE.

In our previous work (Kuvatova et al. 2024), we analysed the dynamical evolution of six GCs (Palomar 6, HP 1, NGC 6401, NGC 6642, NGC 6681, and NGC 6981) and estimated the potential contribution of their putative MSPs to the γ-ray excess in the GalC. We found that even when considering a small number of clusters, their MSPs can provide a noticeable contribution to the excess, while not fully accounting for it, which is the motivation behind the present study.

In this work, we revisit the MSP interpretation by running high-resolution, direct N-body simulations of GCs evolving in a time-varying MW potential, explicitly tracing the orbits of all NSs formed in (and stripped from) the clusters. We consider two complementary channels: (i) present-day Galactic GCs on their observed orbits, to test whether repeated pericentre passages and tidal stripping can inject NSs into the central few kiloparsecs where they might be recycled into MSPs and contribute to the GeV signal; and (ii) a population of now-destroyed clusters whose debris was accreted into the central regions, depositing NSs that subsequently power the GCE. Our approach self-consistently captures internal cluster dynamics (e.g. mass segregation) and external tides, and it produces the full phase-space distribution of deposited NSs. We then map these NS populations to the γ-ray emission using empirically calibrated MSP efficiencies and luminosity functions, and confront the resulting surface-brightness profiles and spectra with Fermi-LAT measurements.

The paper is organised as follows. In Sect. 2, we describe integration procedure and initial conditions for GCs. Section 3 presents the NS distribution in Galaxy at the present day based on GCs N-body simulations. In Sect. 4, we provide our main results about γ-ray fluxes from the distribution of MSPs based on NSs from the simulations. Section 5 contains our conclusions.

2 Integration procedure, initial mass function, and mass loss

To investigate the spatial distribution of GCs’ NSs around the GalC, we performed detailed full N-body simulations with stellar evolution using the high-order parallel code φ-GPU1 (Berczik et al. 2011, 2013). This code implements a fourth-order Her-mite integration scheme with hierarchical individual block time steps and supports parallel computation via CUDA library. In its current version, the code incorporates an updated stellar evolution library (Banerjee et al. 2020; Kamlah et al. 2022b), which accounts for key evolutionary processes, including stellar mass loss via winds, supernova mechanisms, natal kicks, and more.

Each of our present day clusters was initialised in a state of dynamical equilibrium using a King model (King 1966), characterised by the total initial mass, M, the dimensionless central potential, W0, and the half-mass radius, rhm. These parameters were numerically adjusted to fit the current mass and size of the clusters. The stellar mass distribution followed the Kroupa initial mass function (Kroupa 2001) in the range 0.08–100 M⊙. A more detailed description of the integration procedure and the mass evolution of the investigated GCs can be found in Ishchenko et al. (2024).

For the simulations, we selected 12 real GCs. These clusters were chosen primary because of their close passages near the GalC. The initial physical parameters of GCs used in the simulation are listed in Table 1. The initial kinematic conditions for the cluster centres were derived via individual backward orbital integration of GCs as point mass objects to a lookback time of 8 Gyr from the present state, as presented in Ishchenko et al. (2023a). As external potential in both types simulations we used a MW-like time-variable galactic potential from the IllustrisTNG-100 simulation (Nelson et al. 2019; Ishchenko et al. 2023a).

Figure 1 shows the fraction of lost mass for the real GCs over the 8 Gyr of orbital integration. As seen in the plot, the clusters experience the most rapid mass loss during the first billion years of their evolution. The most prominent mass loss is observed for NGC 1904, NGC 6401, HP 1, and NGC 6981.

In addition to the real clusters, we also simulated a set of theoretical GCs. This sample included 50 clusters with an initial mass of 60 000 M⊙ and 4 more massive clusters with a mass of 180 000 M⊙ (see other parameters at the end of Table 1). All of these theoretical clusters were fully disrupted during the 5 Gyr simulation, contributing their stellar content to the MW’s central region. In one theoretical GC with a mass of 60 000 M⊙, a total of more than 700 NS are formed, and in a more massive one, three times more. More information about initial conditions and modelling process can be found in Kuvatova et al. (2026).

Here, we briefly describe the initial positions and velocities of each cluster’s centre. For generating the mesh of random initial conditions, we use the global specific binding energy and the two components of the specific angular momentum, Ltot and Lz, for each disrupted GC. For this purpose, we integrate back 10 Gyr in our well-tested 411321 Illustris TNG-TVP potential. Using this external potential, we compute the phase-space distribution of the GC centres. We adopt specific binding energies in the range (−19 to −14) × 104 km2 s−2. For the specific angular-momentum space, we chose Ltot = 0–3 × 102 kpc km s−1 and Lz = (−0.5 to +0.5) × 102 kpc km s−1, spanning both prograde and retrograde orbits and selected to cover the range of currently observed (and modelled) Galactic GCs (Ishchenko et al. 2023b, 2024).

Table 1

Initial physical characteristics at –8 Gyr for MW and early theoretical GCs.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Evolution of the mass loss in the GCs in percent due to stellar evolution and orbital tidal stripping for our 12 MW GCs.

3 NS distribution in Galaxy at the present day based on GC simulations

The number of NSs in a GC depends both on its initial mass, range of the IMF, and the individual metallicity. In our simulations, the fraction of NSs formed over the course of the cluster’s evolution ranges from 0.5–0.8% relative to the initial number of stars. At the beginning of GC evolution, the largest number of NSs from the 12 MW GCs with their evolution modelled in Terzan 5 is over 22 k. The lowest number of the NSs evolved in NGC 362 is around 5 k.

However, due to the high natal kick velocity prescription, which is implemented in our N-body code, fewer than one-fifth of them remain bound to their parent clusters (e.g. Lyne & Lorimer 1994; Hansen & Phinney 1997; Podsiadlowski et al. 2004; Banerjee et al. 2020; Kamlah et al. 2022b). Recent studies continue to explore and support diverse natal kick prescriptions (Roberts et al. 2025; Popov et al. 2025). These recent prescriptions have included various non-radial asymmetries in supernova explosions, along with very rapid rotation as well as the influence of magnetic field. As a result (taking in account all factors), the retained NS fraction is only about 0.01–0.1% of the initial stellar population (see more details in Table 2), which is also consistent with previous estimates (Pfahl et al. 2002; Kuranov & Postnov 2006).

We classified the NS as a cluster member if it has a negative binding energy with the nucleus of the cluster (i.e. NSs are bound with the GC). We computed the amount of the NS within a tidal radius. The tidal radius (or ‘Jacobi’ radius) was calculated based on the numerical iteration of the Mtid and rtid values using the equation rtid=[G⋅MtidMGal]1/3,Mathematical equation: ${r_{{\rm{tid}}}} = {\left[ {{{G \cdot {M_{{\rm{tid}}}}} \over {{M_{{\rm{Gal}}}}}}} \right]^{1/3}},$(1)

where G, Mtid, and MGal are respectively the gravitational constant, the cluster tidal mass and the Galaxy enclosed mass at the GC current position. For more details, we refer to King (1962); Just et al. (2009); Ernst et al. (2011).

In Fig. 2, we show the distribution of NSs in the Galaxy and the corresponding velocity values. The top panel shows peaks corresponding to bound NSs within the GC’s rtid boundaries at the present day. As expected, the majority of NSs in the GC are found in Terzan 5, NGC 7078, and NGC 6642. The mean relative distance between the GalC and our set of GCs varies from 0.1– 10 kpc, with a velocity distribution of 10–500 km s−1. We note that the NGC 6642 (red colour) has the lowest velocity values and Terzan 5 has the highest.

In bottom panel of Fig. 2, we can see that most of all NS from early destroyed clusters are located in the central region of the Galaxy, less than 1 kpc from the centre. As mentioned above, Sect. 2, we modelled 50 GCs that were initially located in the central region of the Galaxy (at less than 5 kpc) and most of them were located even less than 1 kpc from the centre.

The dependence of velocity as a function of distance is demonstrated in Fig. 3. A small fraction of NSs (1–2% of all NSs) were generated in 12 MW GCs and became unbound at an early stage of the GCs evolution. And after 8 Gyr of evolution, their relative distance exceeded 1000 kpc and even more. These stars had an initial kick velocity of up to 1000 km s−1 at the beginning, as shown by the black line in the bottom panel of Fig. 2. As a result, such stars have the potential to completely escape from our Galaxy.

In Table 2, we summarise all statistical information about the investigated NS. Here, we present the ratio of the initial NS mass to the mass of bound NSs in GCs after 8 Gyr of orbital evolution (i.e. today). In Column 8, we also present the number of observed MSPs in GCs from the web-based online catalogue2.

We also present the theoretical numbers of MSPs for each MW GC based on the study of Turk & Lorimer (2013) and Yin et al. (2024). In the last column, we give the numbers of NSs in the central area of the GalC (see also Fig. 4).

Table 2

Overview of the NS distribution in MW GCs and early theoretical GCs along with the theoretical number of MSPs.

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Position and velocity distribution of the NSs in the Galaxy from the GC’s N-body simulations. Top: real MW GCs. Bottom: early GCs that were destroyed over a 5 Gyr evolution time.
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Dependence of the velocity distribution on distance. Top: real MW GCs. Bottom: destroyed GCs. The colour coding corresponds to the legend in Fig. 1.
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Distribution of NSs in the central region of the Galaxy is less than 1 kpc. This distribution is shown in a three-coordinate projection. From left to right: X vs. Y, Y vs. Z and R vs. Z, where R is the galactic plate radius. Top: NSs from MW GC simulations (different colours correspond to different models, as in Fig. 1). Bottom: NSs from GCs that were destroyed early on.

4 γ-ray fluxes from the distribution of MSPs

4.1 Fitting the MSP-to-NS ratio

From our cluster simulations, we were able to directly obtain the number of NSs located in the central region of the Galaxy within a maximum radius up to 1 kpc (see Table 2, Column 9).

The central NS population is a mixture of stars stripped from different GCs at different epochs. Some of these NSs may end up as MSPs (Bhattacharya & van den Heuvel 1991; Lorimer 2008) with the capacity to emit γ radiation. Our simulations provide the full 6D phase space for all these NSs; however, they do not include binary stellar evolution, which is required to recycle NSs into pulsars. Nevertheless, we can statistically estimate the number of pulsars per NS (based on the observed MSPs in real GCs) to assess the contribution of these MSPs to the observed γ-ray excess.

A cluster’s initial mass largely sets the total number of NSs formed, but this proportionality strictly applies at formation. Figure 5 (top) shows that the total number of NSs correlated with the initial cluster mass. Because all our modelled clusters share the same Kroupa initial mass function and stellar-mass range, the total number of NSs produced in a GC is nearly proportional to its initial mass; any differences between clusters of the same initial mass arise mainly from metallicity. By contrast, the correlation weakens when we examine the number of bound NSs as a function of a cluster’s current mass within the tidal radius (Fig. 5, bottom panel), reflecting the effects of mass loss on different Galactic orbits, stellar stripping, tides, and retention physics.

To infer the fraction of NSs in the GalC that are observable today as MSPs, we adopted an observationally motivated statistical calibration. Specifically, we mapped the number of bound NSs in each modelled cluster to the number of observed MSPs in the corresponding GC (Turk & Lorimer 2013; Yin et al. 2024). Because the observed counts can be incomplete, we also used the predicted total MSP numbers for those clusters (Turk & Lorimer 2013; Yin et al. 2024), which both mitigated any selection biases and allowed for the inclusion of clusters without confirmed MSP detections.

In Fig. 6 (top), observed MSPs in real clusters are shown as coloured points, while theoretical MSP predictions (NMSP) based on a statistical model by Turk & Lorimer (2013) and completed with FAST data (Yin et al. 2024) are shown as open triangles. The corresponding numbers are listed in parentheses in Table 2 (column 8). All displayed MSPs are gravitationally bound to their host clusters. For each cluster, we computed the ratio NMSP/Nbound NS. These values are summarised in Table 3. Here, multiple entries indicate distinct observational or predicted estimates. These values are directly correspond to those shown in Fig. 6 (top). With these ratios, we obtained a median MSP-per-NS coefficient of k≡median(NMSP/Nbound NS)=0.0114−0.0103+0.0173.Mathematical equation: $k \equiv {\rm{median}}\left( {{N_{{\rm{MSP}}}}/{N_{{\rm{bound NS}}}}} \right) = 0.0114_{ - 0.0103}^{ + 0.0173}.$(2)

Figure 6 (lower panel) shows the resulting MSP counts predicted for each cluster (with asterisks). To analyse the ratio of MSPs to bound NSs in selected GCs, we computed the MSP/NS ratio for both observed and theoretical data. The median and median absolute deviation (MAD) were employed instead of the mean and standard deviation, since these statistics are less sensitive to outliers and provide a more robust description of non-Gaussian datasets.

The median represents the central value that divides the dataset into two equal halves, while MAD quantifies the typical deviation of data points from the median, defined as MAD = median (|xi − median(x)|).

  • overall dataset (observed + theoretical): median = 0.011449, MAD = 0.008142;

  • upper level of the predicted values for each cluster (theor2: the greater of the two theoretical estimates per cluster, corresponding to predictions by Yin et al. (2024); theor1 denotes predictions by Turk & Lorimer (2013) and corresponding to the lower of the two theoretically predicted values): mediantheor2 = 0.016154, MADtheor2 = 0.012592;

  • observed values only: medianobs = 0.006608, MADobs = 0.005467.

Based on these results, the lower and upper bounds were defined as lower = medianobs − MADobs = 0.001141, upper = mediantheor2 + MADtheor2 = 0.028746. These values define the grey shaded region in the Fig. 6, representing the statistically expected range of MSP/NS ratios. The red dashed line marks the overall median, while the grey dashed lines correspond to the upper and lower MAD-based limits. Observed MSP ratios are shown as filled circles, while theoretical predictions are shown as open triangles connected by dashed lines for each cluster.

We then applied the same coefficient to the NSs stripped from each system, assuming that the MSP-to-NS ratio is approximately conserved during the clusters’ dynamical evolution and subsequent deposition of NSs into the Galactic central bulge.

We extended the calibration to disrupted clusters by using the same coefficient, k, rather than cluster-specific values. Because the now-destroyed GCs lack reliable present-day parameters (e.g. mass, metallicity, encounter rate), a single empirically anchored k provides a consistent mapping from each hypothetical cluster’s neutron-star population to its expected MSP yield. Our simulations supply the full 6D phase-space distribution of deposited NSs, so the MSP-per-NS approach converts these kinematic populations to NMSP without invoking additional, poorly constrained global factors, yielding a robust and uniform estimate for the disrupted population.

Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Relationship between GC mass and NS populations. Colours denote individual clusters as listed in legend of Fig. 2. The full dataset is given in Table 2. Top: initial mass of GCs vs the total number of NSs formed up to the present day. The dashed line shows the best-fit linear relation (slope = 9428.42 NSs per 106M⊙, R2 = 0.87), with the shaded region representing the ±1436 NSs scatter around the fit. Bottom: current mass of GCs vs bound number of NSs. The dashed line shows the best-fit linear relation (slope = 95.29 NSs per 105 M⊙, R2 = 0.11), with the shaded region representing the ±842 NSs scatter around the fit.
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

MSP populations in selected GCs. Top: MSP-to-bound NS ratio. Colours correspond to individual clusters as shown in the legend of Fig. 2, and colour dashed lines connect multiple estimates for the same cluster. Filled circles represent observed MSPs and open triangles show theoretical predictions. The red dashed line indicates the median value, while the shaded region represents the median absolute deviation (MAD) range. Lower panel: number of MSPs in the selected GCs. Star symbols denote MSP numbers predicted in this work, computed as NMSP = 0.0114 × NBound NS based on our simulation results.
Table 3

NMSP/Nbound NS ratios for each cluster used in the calibration for obs/bound NS, theor1/bound NS, and theor2/bound NS.

4.2 Mock observations: NSs from real and destroyed GCs

In Fig. 4, we show distribution of NSs in Cartesian Galactocentric coordinates (X⋆, Y⋆, Z⋆). To convert the individual NSs coordinates to Galactic coordinates (Heliocentric distance – d, Galactic longitude (LON) – ℓ, and latitude (LAT) – b), we used Galactocentric coordinates of the Sun as X⊙ = −8178 pc (Gravity Collaboration 2019), Y⊙ = 0 pc, Z⊙ = 20 pc (Bennett & Bovy 2019). Next, we used Eq. (3) to find the Cartesian Heliocentric coordinates of NSs (x⋆, y⋆, z⋆). After using Eq. (4), we transformed the NSs Heliocentric coordinates into Galactic coordinates. More details on a similar method are given in Kalambay et al. (2022), while Bissekenov et al. (2024) used this method for simulated clusters from Shukirgaliyev et al. (2017, 2018, 2019, 2021): x⋆=X⋆−X⊙,y⋆=Y⋆−Y⊙,z⋆=Z⋆−Z⊙,Mathematical equation: ${{x_ \star } = {X_ \star } - {X_ \odot },\quad {y_ \star } = {Y_ \star } - {Y_ \odot },\quad {z_ \star } = {Z_ \star } - {Z_ \odot },}$(3) d=x⋆2+y⋆2+z⋆2,l=arctan(y⋆x⋆),b=arcsin(z⋆d).Mathematical equation: ${d = \sqrt {x_ \star ^2 + y_ \star ^2 + z_ \star ^2} ,\ell = \arctan \left( {{{{y_ \star }} \over {{x_ \star }}}} \right),b = \arcsin \left( {{{{z_ \star }} \over d}} \right).}$(4)

We combined the NS data from all our simulated clusters to obtain a distribution of NSs around the GalC in the Galactic coordinates. Then, the previously obtained coefficient (k = 0.0114) was multiplied by the number of NSs in the centre of the Galaxy to estimate the number of MSPs within the 1 kpc sphere around the centre. MSPs in GCs have mean γ-ray luminosities roughly ⟨L⟩ ∼ (1 − 8) × 1033 erg/s (Amerio et al. 2025). For the purposes of demonstration, we used three values of luminosity ⟨L⟩ = 1 × 1033 erg/s, 4 × 1033 erg/s, and 8 × 1033 erg/s to estimate the total flux from the GalC. This result is shown in Fig. 7 (top). The same method was applied to NSs from destroyed clusters (Fig. 7, middle), and the total flux (combination of NSs from both simulated real and destroyed theoretical GCs) is shown in the lower panel of Fig. 7.

We want to highlight the interesting feature of the NS distribution from our population of the early destroyed theoretical GCs. The NS distribution from these clusters show a strong axisymmetric signature with the notable over-density in the Galactic disk plane, see Figs. 4 and 7. A similar spatial feature in the γ-ray excess was recently reported by Lipari & Vernetto (2025) and Muru et al. (2025). Here, we also note another less prominent feature perpendicular to the Galactic plane; although this over-density is not as notable as the similar feature in Galactic disk, we like to call attention to it, as it might potentially be detected in future γ-ray observational programs and surveys.

In Fig. 8, we present the cumulative integrated γ-ray flux from the summed NSs in GalC. We summed up all the stripped NSs from the present-day observed GCs with the contribution of NSs from now-destroyed early GCs. We assumed a fixed luminosity of each MSP on three different levels: ⟨L⟩ = 1 × 1033 erg/s, 4 × 1033 erg/s, and 8 × 1033 erg/s. In Fig. 8 (left), we plot the different contributions with different line-styles. Using different colours, we present the various luminosity levels. As we see from the result of our modelling, the summed contribution for the maximum luminosity level is generally quite consistent with the observationally available data Abazajian et al. (2014); Brandt & Kocsis (2015). In this plot, we also present the artificially enhanced by a contribution (of a factor of 2) from the early GC population (grey line: with MSP maximum luminosity level 8 × 1033 erg/s). In this case, we have nearly achieved an ideal correspondence with the observations after subtracting the contribution from the Galactic SMBH (Abazajian et al. 2014), represented by the filled black pentagons.

In Fig. 8 (right panel), we present the same cumulative distribution, calculated using the minimum and maximum coefficients (theor1 and theor2) from Turk & Lorimer (2013) and Yin et al. (2024), respectively. As can be seen, when adopting the maximum average flux 8 × 1033 erg/s from a single MSP according to Amerio et al. (2025) and using the maximum coefficient k (red line), even with the current number of early (now already disrupted) GCs, we are still able to account for the required overall level of the GCE.

For the lower luminosity levels, we see the significant deficit of our theoretical cumulative flux compared with the observations. Such a flux deficit can mainly be explained by our applied assumptions: (i) on the total numbers of early GC; (ii) on the initial mass of these clusters; (iii) on the disruption rate of such systems (possibly accounting for intensive dynamical friction in initially more massive clusters). Additionally, in this work, we focussed on a scenario where MSPs were deposited into the GalC by GCs. However, pulsars could also form in situ, given the complex star formation history of the central region (Schödel et al. 2020; Nogueras-Lara et al. 2021). Moreover, Panamarev et al. (2019) showed that binaries in the inner few parsecs can harden efficiently (see their Fig. 10), potentially boosting MSP formation.

5 Conclusions and discussion

We estimated the contribution to the γ-ray flux in the Galactic Centre from MSPs that originated in MW GCs. We used a suite of direct N-body simulations of observed GCs evolving in a time-varying MW potential to track the individual NSs tidally stripped from their host clusters and deposited in the central kpc region of the Galaxy. We also generated a set of hypothetical GCs that underwent complete tidal disruption and did not survive to the present day. We then estimated the number of MSPs per NS for both present-day and hypothetical GCs and used the mean γ-ray luminosity per pulsar to compute their net flux. In our numerical investigation, we found that:

  • The total NS production scales with the initial GC mass (Fig. 5, top), as expected for an applied common IMF. In contrast, the correlation between the bound NS and the current GC mass is weak (Fig. 5, bottom), reflecting the fact that there is a significant kick velocity for the remnants during the formation of NS (Banerjee et al. 2020; Kamlah et al. 2022b). As a consequence, the slight dependence of the number of retained NSs on the cluster with the mass of the cluster;

  • Calibrating the bound-NS counts in real clusters to their observed versus predicted MSP content yields the sample median coefficient k = median(NMSP/Nbound NS) ≃ 0.0114 (Table 3, Fig. 6). Applying the same k to all stripped NSs in the Galactic center provides a uniform mapping to NMSP for both the sources of surviving and disrupted GC systems. The cumulative flux profile (Fig. 8, left panel) matches the observed trend when adopting higher MSP luminosities (⟨L⟩ ∼ 8 × 1033 erg s−1) and is brought into near-ideal agreement after subtracting the SMBH component (Abazajian et al. 2014; Brandt & Kocsis 2015) if the early-GC contribution is modestly enhanced (by factor of ∼2). For a more detailed analysis of the cumulative distribution of γ-ray flux around the GalC, we also employed two alternative ratios of the number of MSPs to the total number of NSs in observed Galactic stellar systems, namely, the coefficient k. One of them, theor1 from Turk & Lorimer (2013), corresponds to the minimum value of this coefficient. The other, theor2 from Yin et al. (2024), corresponds to the maximum value of k. The results of these calculations are shown in Fig. 8 (right panel). According to our estimates, even with the current number of early (i.e., now already disrupted) GCs, assuming the maximum theor2 individual k values for each observed GC, we were able to reproduce the overall level of γ–ray flux around the GalC at the observed level (Abazajian et al. 2014; Brandt & Kocsis 2015);

  • For lower ⟨L⟩ (e.g. (1–4) × 1033 erg s−1), the model under-produces the cumulative flux, pointing to uncertainties in the early GC populations, initial mass spectrum, phase-space distribution, disruption efficiencies, and dynamical friction. These are the parameters that directly modulate the stripped NS reservoir;

  • Mock-observation maps within 1 kpc (Fig. 7) show that the MSP population inferred from present-day GCs alone already contributes substantially to the GalC γ-ray signal. Adding the disrupted population boosts both the amplitude and central concentration;

  • NSs deposited by early destroyed GCs exhibit an axisymmetric signature with a pronounced over-density along the disk plane (Figs. 4 and 7) are qualitatively consistent with recent reports of a similar structure in the γ-ray excess morphology Lipari & Vernetto (2025); Muru et al. (2025). We also note another less prominent feature perpendicular to the Galactic plane, which can potentially be detected in future γ-ray observational surveys.

Thus, our modelling favours the MSP origin of the GCE over the DM-annihilation scenario. To estimate the net flux from pulsars, we assumed a mean γ-ray luminosity of ⟨L⟩ = 8 × 1033 erg s−1. This assumption is critical for the total flux calculation. For example, Fragione et al. (2018) found that the observed GCE flux can be reproduced using a more moderate mean luminosity of ⟨L⟩ = 2 × 1033 erg s−1. The main difference between our calculations and those performed by Fragione et al. (2018) is the more self-consistent treatment we adopted with respect to the number of NSs retained in the Galactic Centre region (<1 kpc).

In that work, the authors assumed a relatively large and fixed retained fraction of ∼10% of the total number of produced NSs. In our investigation, we modelled each NS orbit individually and dynamically compute the retained fraction for each cluster. In this way, we found that only ∼6% of NSs remain within the central 1 kpc region of the Galaxy. Additionally, including a larger sample of potentially destroyed clusters could increase the number of deposited NSs in the central region, thereby relaxing the requirement for a relatively high average per-pulsar luminosity to match the observed flux. Finally, the MSP population in the Galactic Centre might also be supplemented by in situ formation, given the complex star formation history of the region (Schödel et al. 2020; Nogueras-Lara et al. 2021) and efficient dynamical hardening of binaries in the inner few parsecs could alsio promote MSP production (e.g., Panamarev et al. 2019). A future study that jointly accounts for these effects could help mitigate the associated uncertainties.

The following extensions will make our predictions more precise. First, incorporating a fully time-varying Galactic potential directly in a code which features accurate treatment of binary stars, such as NBODY6++GPU (Wang et al. 2015, 2016; Panamarev et al. 2019; Kamlah et al. 2022a), together with binary stellar evolution (to follow LMXB stages, recycling, and MSP birth+aging) would remove our reliance on a global k and enable cluster-by-cluster predictions of NMSP and their γ-ray luminosities. Second, implementing dynamical friction self-consistently (or coupling it to an evolving live background) should tighten the disrupted-GC contribution and central concentration. Third, we assumed a coefficient k and also fixed the per-pulsar average luminosity ⟨L⟩. A hierarchical treatment that fits a cluster-dependent k, assuming individual natal-kick retention, and a full MSP luminosity function would reduce systematics tied to selection effects in current MSP/LMXB catalogues. Finally, forward-modelling MSP beaming geometry, duty cycles, and magnetospheric cut-offs would further reduce biases in the flux maps and in the cumulative profile.

Despite these modelling choices, in its current state, our work favours an MSP origin of the GalC γ-ray excess over DM annihilation. This is primarily because the combined contribution of MSPs delivered by surviving and disrupted GCs naturally reproduces both the amplitude and concentration of the observed signal under reasonable assumptions.

Further observational leverage will come from deep bulge radio pulsar searches (e.g. FAST, MeerKAT, SKA) targeted along our predicted over-densities (Han et al. 2021, 2024; Frail et al. 2024; Keane et al. 2015); TeV-halo mapping with CTA to test the integrated MSP population (Hooper & Linden 2018); improved MeV–GeV spectroscopy and diffuse-background control with next-generation γ-ray missions (Caputo et al. 2022; De Angelis et al. 2017); updated LMXB inventories (e.g. eROSITA/Chandra/NuSTAR) to anchor progenitors (Jonker et al. 2011; Wevers et al. 2016; Hong et al. 2016); continuous-GW searches constraining the bulge MSP population (Miller & Zhao 2023); and refined bulge and bar stellar maps to compare with our predicted morphology (Wegg & Gerhard 2013; Valenti et al. 2016). These data will allow for tighter, physically anchored priors to be set on MSP demographics, thereby improving our understanding of the Galactic γ-ray excess content.

Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

γ-ray flux maps from simulated NSs within 1 kpc in Galactic coordinates. The flux is computed assuming a MSP luminosity of ⟨L⟩ ∼ 1 × 1033 erg/s (left), ⟨L⟩ ∼ 4 × 1033 erg/s (middle), and ⟨L⟩ ∼ 8 × 1033 erg/s (right), scaled by empirical factors for observed GCs (top), disrupted GCs (middle), and the combination of both populations (bottom). The colour scale indicates the total γ-ray flux in each bin, shown on a logarithmic scale between 10−12 and 10−16 erg s−1 cm−2. To obtain a smoother flux map, we applied a Gaussian filter with a smoothing parameter of σ = 1, which corresponds to a smoothing scale of approximately 0.5°.
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Cumulative γ-ray flux as a function of projected radius from the GalC. The bottom axis shows the angular distance Ψ in degrees, while the top axis indicates the corresponding projected distance in kiloparsecs, assuming a Sun–GalC distance of 8178 pc. Left panel: results for three different luminosities: L = 1 × 1033 erg s−1 (blue), L = 4 × 1033 erg s−1 (orange), and L = 8 × 1033 erg s−1 (green). Right: case with fixed luminosity of L = 8 × 1033 erg s−1, but different values of the scaling parameter k: the median value k = 0.0114 (same green as on the left), the maximum individual value, corresponding to theor2 from Yin et al. (2024, shown in red), and the minimum individual value, corresponding to theor1 from Turk & Lorimer (2013, shown in sky blue,) as listed in Table 3 for each cluster. The legends in the figure: MW is the cumulative flux calculated only based on the selected real MW clusters (see Table 1); Destr. refers to fluxes calculated only on the basis of the theoretical, already destroyed clusters; MW+Destr. is the sum of the flux from real MW plus the Destr. theoretical clusters; MW + 2 × Destr. is the same as the above sum, but with the contribution from the ‘Destr.’ clusters artificially enhanced by a factor of 2.

Acknowledgements

The authors thank the anonymous referee for a very constructive report and suggestions that helped significantly improve the quality of the manuscript. This research has been funded by the Science Committee of the Ministry of Science and Higher Education, Republic of Kazakhstan (Grant No. AP26102895). PB and MI appreciate the support from the Polish Academy of Sciences Special Program and the U.S. National Academy of Sciences Special Program under the Long-term program to support Ukrainian research teams grant No. PAN.BFB.S.BWZ.329.022.2023. We also gratefully acknowledge the Polish high-performance computing infrastructure PLGrid (HPC Center: ACK Cyfronet AGH) for providing computer facilities and support within computational grant No. PLG/2026/019243.

References

  1. Abazajian, K. N., Canac, N., Horiuchi, S., & Kaplinghat, M. 2014, Phys. Rev. D, 90, 023526 [Google Scholar]
  2. Ackermann, M., Ajello, M., Albert, A., et al. 2017, ApJ, 840, 43 [NASA ADS] [CrossRef] [Google Scholar]
  3. Amerio, A., Hooper, D., & Linden, T. 2025, J. Cosmology Astropart. Phys., 2025, 106 [Google Scholar]
  4. Arca-Sedda, M., & Capuzzo-Dolcetta, R. 2014, MNRAS, 444, 3738 [NASA ADS] [CrossRef] [Google Scholar]
  5. Arca-Sedda, M., Kocsis, B., & Brandt, T. D. 2018, MNRAS, 479, 900 [NASA ADS] [CrossRef] [Google Scholar]
  6. Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bartels, R., Krishnamurthy, S., & Weniger, C. 2016, Phys. Rev. Lett., 116, 051102 [CrossRef] [PubMed] [Google Scholar]
  8. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
  9. Berczik, P., Nitadori, K., Zhong, S., et al. 2011, in International Conference on High Performance Computing, HPC-UA 2011, 8 [Google Scholar]
  10. Berczik, P., Spurzem, R., & Wang, L. 2013, in Third International Conference on High Performance Computing, HPC-UA 2013, 52 [Google Scholar]
  11. Bhattacharya, D., & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1 [Google Scholar]
  12. Bissekenov, A., Kalambay, M., Abdikamalov, E., et al. 2024, A&A, 689, A282 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Brandt, T. D., & Kocsis, B. 2015, ApJ, 812, 15 [Google Scholar]
  14. Caputo, R., Ajello, M., Kierans, C. A., et al. 2022, J. Astron. Telesc. Instrum. Syst., 8, 044003 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cholis, I., Zhong, Y.-M., McDermott, S. D., & Surdutovich, J. P. 2022, Phys. Rev. D, 105, 103023 [CrossRef] [Google Scholar]
  16. Daylan, T., Finkbeiner, D. P., Hooper, D., et al. 2016, Phys. Dark Universe, 12, 1 [NASA ADS] [CrossRef] [Google Scholar]
  17. De Angelis, A., Tatischeff, V., Tavani, M., et al. 2017, Exp. Astron., 44, 25 [NASA ADS] [CrossRef] [Google Scholar]
  18. Di Mauro, M. 2021, Phys. Rev. D, 103, 063029 [Google Scholar]
  19. Eckner, C., Hou, X., Serpico, P. D., et al. 2018, ApJ, 862, 79 [NASA ADS] [CrossRef] [Google Scholar]
  20. Ernst, A., Just, A., Berczik, P., & Olczak, C. 2011, A&A, 536, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Fragione, G., Pavlík, V., & Banerjee, S. 2018, MNRAS, 480, 4955 [NASA ADS] [Google Scholar]
  22. Frail, D. A., Polisensky, E., Hyman, S. D., et al. 2024, ApJ, 975, 34 [Google Scholar]
  23. Gautam, A., Crocker, R. M., Ferrario, L., et al. 2022, Nat. Astron., 6, 703 [CrossRef] [Google Scholar]
  24. Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, Rev. Mod. Phys., 82, 3121 [Google Scholar]
  25. Ghez, A. M., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 1044 [Google Scholar]
  26. Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30 [Google Scholar]
  27. Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223 [Google Scholar]
  28. Gnedin, O. Y., Ostriker, J. P., & Tremaine, S. 2014, ApJ, 785, 71 [Google Scholar]
  29. Goodenough, L., & Hooper, D. 2009, Fermilab | Technical Publications [arXiv:0910.2998] [Google Scholar]
  30. Gravity Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Han, J. L., Wang, C., Wang, P. F., et al. 2021, Res. Astron. Astrophys., 21, 107 [CrossRef] [Google Scholar]
  32. Han, J. L., Zhou, D. J., Wang, P. F., et al. 2024, Res. Astron. Astrophys., 24, 125008 [Google Scholar]
  33. Hansen, B. M. S., & Phinney, E. S. 1997, MNRAS, 291, 569 [NASA ADS] [Google Scholar]
  34. Holst, I., & Hooper, D. 2025, Phys. Rev. D, 111, 023048 [Google Scholar]
  35. Hong, J., Mori, K., Hailey, C. J., et al. 2016, ApJ, 825, 132 [Google Scholar]
  36. Hooper, D., & Goodenough, L. 2011, Phys. Lett. B, 697, 412 [NASA ADS] [CrossRef] [Google Scholar]
  37. Hooper, D., & Linden, T. 2018, Phys. Rev. D, 98, 043005 [NASA ADS] [CrossRef] [Google Scholar]
  38. Ishchenko, M., Sobolenko, M., Berczik, P., et al. 2023a, A&A, 673, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Ishchenko, M., Sobolenko, M., Kuvatova, D., Panamarev, T., & Berczik, P. 2023b, A&A, 674, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Ishchenko, M., Berczik, P., Panamarev, T., et al. 2024, A&A, 689, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Jonker, P. G., Bassa, C. G., Nelemans, G., et al. 2011, ApJS, 194, 18 [NASA ADS] [CrossRef] [Google Scholar]
  42. Just, A., Berczik, P., Petrov, M. I., & Ernst, A. 2009, MNRAS, 392, 969 [Google Scholar]
  43. Kalambay, M. T., Naurzbayeva, A. Z., Otebay, A. B., et al. 2022, Rec. Contrib. Phys., 83, 4 [Google Scholar]
  44. Kamlah, A. W. H., Leveque, A., Spurzem, R., et al. 2022a, MNRAS, 511, 4060 [NASA ADS] [CrossRef] [Google Scholar]
  45. Kamlah, A. W. H., Spurzem, R., Berczik, P., et al. 2022b, MNRAS, 516, 3266 [CrossRef] [Google Scholar]
  46. Keane, E., Bhattacharyya, B., Kramer, M., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 40 [Google Scholar]
  47. King, I. 1962, AJ, 67, 471 [Google Scholar]
  48. King, I. R. 1966, AJ, 71, 64 [Google Scholar]
  49. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  50. Kuranov, A. G., & Postnov, K. A. 2006, Astron. Lett., 32, 393 [Google Scholar]
  51. Kuvatova, D., Panamarev, T., Ishchenko, M., Gluchshenko, A., & Berczik, P. 2024, Rec. Contrib. Phys., 90, 14 [Google Scholar]
  52. Kuvatova, D., Ishchenko, M., Berczik, P., et al. 2026, A&A, 706, A220 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Lee, S. K., Lisanti, M., Safdi, B. R., Slatyer, T. R., & Xue, W. 2016, Phys. Rev. Lett., 116, 051103 [NASA ADS] [CrossRef] [Google Scholar]
  54. Lipari, P., & Vernetto, S. 2025, Phys. Rev. D, 111, 123040 [Google Scholar]
  55. Lorimer, D. R. 2008, Liv. Rev. Relativ., 11, 8 [NASA ADS] [CrossRef] [Google Scholar]
  56. Lyne, A. G., & Lorimer, D. R. 1994, Nature, 369, 127 [NASA ADS] [CrossRef] [Google Scholar]
  57. Miller, A. L., & Zhao, Y. 2023, Phys. Rev. Lett., 131, 081401 [NASA ADS] [CrossRef] [Google Scholar]
  58. Muru, M. M., Silk, J., Libeskind, N. I., Gottlöber, S., & Hoffman, Y. 2025, Phys. Rev. Lett., 135, 161005 [Google Scholar]
  59. Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computat. Astrophys. Cosmol., 6, 2 [NASA ADS] [CrossRef] [Google Scholar]
  60. Neumayer, N., Seth, A., & Böker, T. 2020, A&A Rev., 28, 4 [Google Scholar]
  61. Nogueras-Lara, F., Schödel, R., & Neumayer, N. 2021, ApJ, 920, 97 [NASA ADS] [CrossRef] [Google Scholar]
  62. Panamarev, T., Just, A., Spurzem, R., et al. 2019, MNRAS, 484, 3279 [NASA ADS] [CrossRef] [Google Scholar]
  63. Pfahl, E., Rappaport, S., & Podsiadlowski, P. 2002, ApJ, 573, 283 [NASA ADS] [CrossRef] [Google Scholar]
  64. Podsiadlowski, P., Langer, N., Poelarends, A. J. T., et al. 2004, ApJ, 612, 1044 [NASA ADS] [CrossRef] [Google Scholar]
  65. Popov, S., Müller, B., & Mandel, I. 2025, New A Rev., 101, 101734 [Google Scholar]
  66. Roberts, D., Gieles, M., Erkal, D., & Sanders, J. L. 2025, MNRAS, 538, 454 [Google Scholar]
  67. Schödel, R., Nogueras-Lara, F., Gallego-Cano, E., et al. 2020, A&A, 641, A102 [Google Scholar]
  68. Shukirgaliyev, B., Parmentier, G., Berczik, P., & Just, A. 2017, A&A, 605, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Shukirgaliyev, B., Parmentier, G., Just, A., & Berczik, P. 2018, ApJ, 863, 171 [Google Scholar]
  70. Shukirgaliyev, B., Parmentier, G., Berczik, P., & Just, A. 2019, MNRAS, 486, 1045 [NASA ADS] [CrossRef] [Google Scholar]
  71. Shukirgaliyev, B., Otebay, A., Sobolenko, M., et al. 2021, A&A, 654, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Tremaine, S. D., Ostriker, J. P., & Spitzer, L., J. 1975, ApJ, 196, 407 [Google Scholar]
  73. Turk, P. J., & Lorimer, D. R. 2013, MNRAS, 436, 3720 [Google Scholar]
  74. Valenti, E., Zoccali, M., Gonzalez, O. A., et al. 2016, A&A, 587, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070 [NASA ADS] [CrossRef] [Google Scholar]
  76. Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450 [Google Scholar]
  77. Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874 [Google Scholar]
  78. Wevers, T., Hodgkin, S. T., Jonker, P. G., et al. 2016, MNRAS, 458, 4530 [Google Scholar]
  79. Yin, D., Zhang, L.-y., Qian, L., et al. 2024, ApJ, 969, L7 [NASA ADS] [CrossRef] [Google Scholar]
1

N-body code φ−GPU : https://github.com/berczik/phi-GPU-mole

All Tables

Table 1

Initial physical characteristics at –8 Gyr for MW and early theoretical GCs.

In the text
Table 2

Overview of the NS distribution in MW GCs and early theoretical GCs along with the theoretical number of MSPs.

In the text
Table 3

NMSP/Nbound NS ratios for each cluster used in the calibration for obs/bound NS, theor1/bound NS, and theor2/bound NS.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Evolution of the mass loss in the GCs in percent due to stellar evolution and orbital tidal stripping for our 12 MW GCs.
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Position and velocity distribution of the NSs in the Galaxy from the GC’s N-body simulations. Top: real MW GCs. Bottom: early GCs that were destroyed over a 5 Gyr evolution time.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Dependence of the velocity distribution on distance. Top: real MW GCs. Bottom: destroyed GCs. The colour coding corresponds to the legend in Fig. 1.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Distribution of NSs in the central region of the Galaxy is less than 1 kpc. This distribution is shown in a three-coordinate projection. From left to right: X vs. Y, Y vs. Z and R vs. Z, where R is the galactic plate radius. Top: NSs from MW GC simulations (different colours correspond to different models, as in Fig. 1). Bottom: NSs from GCs that were destroyed early on.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Relationship between GC mass and NS populations. Colours denote individual clusters as listed in legend of Fig. 2. The full dataset is given in Table 2. Top: initial mass of GCs vs the total number of NSs formed up to the present day. The dashed line shows the best-fit linear relation (slope = 9428.42 NSs per 106M⊙, R2 = 0.87), with the shaded region representing the ±1436 NSs scatter around the fit. Bottom: current mass of GCs vs bound number of NSs. The dashed line shows the best-fit linear relation (slope = 95.29 NSs per 105 M⊙, R2 = 0.11), with the shaded region representing the ±842 NSs scatter around the fit.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

MSP populations in selected GCs. Top: MSP-to-bound NS ratio. Colours correspond to individual clusters as shown in the legend of Fig. 2, and colour dashed lines connect multiple estimates for the same cluster. Filled circles represent observed MSPs and open triangles show theoretical predictions. The red dashed line indicates the median value, while the shaded region represents the median absolute deviation (MAD) range. Lower panel: number of MSPs in the selected GCs. Star symbols denote MSP numbers predicted in this work, computed as NMSP = 0.0114 × NBound NS based on our simulation results.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

γ-ray flux maps from simulated NSs within 1 kpc in Galactic coordinates. The flux is computed assuming a MSP luminosity of ⟨L⟩ ∼ 1 × 1033 erg/s (left), ⟨L⟩ ∼ 4 × 1033 erg/s (middle), and ⟨L⟩ ∼ 8 × 1033 erg/s (right), scaled by empirical factors for observed GCs (top), disrupted GCs (middle), and the combination of both populations (bottom). The colour scale indicates the total γ-ray flux in each bin, shown on a logarithmic scale between 10−12 and 10−16 erg s−1 cm−2. To obtain a smoother flux map, we applied a Gaussian filter with a smoothing parameter of σ = 1, which corresponds to a smoothing scale of approximately 0.5°.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Cumulative γ-ray flux as a function of projected radius from the GalC. The bottom axis shows the angular distance Ψ in degrees, while the top axis indicates the corresponding projected distance in kiloparsecs, assuming a Sun–GalC distance of 8178 pc. Left panel: results for three different luminosities: L = 1 × 1033 erg s−1 (blue), L = 4 × 1033 erg s−1 (orange), and L = 8 × 1033 erg s−1 (green). Right: case with fixed luminosity of L = 8 × 1033 erg s−1, but different values of the scaling parameter k: the median value k = 0.0114 (same green as on the left), the maximum individual value, corresponding to theor2 from Yin et al. (2024, shown in red), and the minimum individual value, corresponding to theor1 from Turk & Lorimer (2013, shown in sky blue,) as listed in Table 3 for each cluster. The legends in the figure: MW is the cumulative flux calculated only based on the selected real MW clusters (see Table 1); Destr. refers to fluxes calculated only on the basis of the theoretical, already destroyed clusters; MW+Destr. is the sum of the flux from real MW plus the Destr. theoretical clusters; MW + 2 × Destr. is the same as the above sum, but with the contribution from the ‘Destr.’ clusters artificially enhanced by a factor of 2.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.