Disrupted Globular Clusters Can Explain the Galactic Center Gamma Ray Excess

The Fermi satellite has recently detected gamma ray emission from the central regions of our Galaxy. This may be evidence for dark matter particles, a major component of the standard cosmological model, annihilating to produce high-energy photons. We show that the observed signal may instead be generated by millisecond pulsars that formed in dense star clusters in the Galactic halo. Most of these clusters were ultimately disrupted by evaporation and gravitational tides, contributing to a spherical bulge of stars and stellar remnants. The gamma ray amplitude, angular distribution, and spectral signatures of this source may be predicted without free parameters, and are in remarkable agreement with the observations. These gamma rays are from fossil remains of dispersed clusters, telling the history of the Galactic bulge.


INTRODUCTION
While there are strong indications for the existence of cold dark matter from its gravitational effects (e.g. Planck Collaboration et al. 2014), there has not yet been any conclusive direct or indirect detection of the corresponding dark matter particles. One promising avenue to look for these particles is through annihilation in which two dark matter particles (a particle and its antiparticle) convert into high energy photons that we can observe. The dark matter annihilation signal is expected to be strongest where the density of dark matter is highest, i.e., in the centers of galaxies.
Detailed analyses of the Fermi satellite's map of the gamma-ray sky have revealed an excess around the Galactic center peaking at energies of ∼2 GeV (e.g. Hooper & Goodenough 2011;Gordon & Macías 2013;Daylan et al. 2014). This excess appears to be roughly spherical and extends at least ∼10-20 • (1.5-3 kpc) from Sgr A*, the Galaxy's central supermassive black hole. Remarkably, this signal can be interpreted as photons from annihilating ∼30 GeV dark matter particles (Hooper & Goodenough 2011;Daylan et al. 2014). In order to confirm this extraordinary interpretation, one must carefully rule out all other astrophysical sources. Possible alternatives include millisecond pulsars (MSPs), rapidly spinning neutron stars that are observed in other regions of the Galaxy with very similar gamma ray spectra to that of the observed excess (Gordon & Macías 2013;Abazajian & Kaplinghat 2012;Yuan & Zhang 2014;Abazajian 2011;Mirabal 2013;Yuan & Ioka 2015;Petrović et al. 2015); highly magnetized young pulsars created in the innermost nuclear star cluster (O'Leary et al. 2015); injection of cosmic-ray protons (Carlson & Profumo 2014); or cosmic ray outbursts (Petrović et al. 2014). However, it remains to be shown that any of these sources is sufficiently abundant and spatially extended to explain the gamma-ray excess.
Energetic photons have also been observed from within 1 Institute for Advanced Study, Einstein Dr., Princeton, NJ 2 Eötvös University, Pázmány P. s. 1/A, Budapest, Hungary 3 NASA Sagan Fellow the central few pc around Sgr A* itself, extending from soft X-rays to ∼100 TeV gamma rays (Baganoff et al. 2001;Aharonian et al. 2004;Bélanger et al. 2006;Perez et al. 2015). The origin of this emission is subject to debate; see van Eldik (2015) for a review. The region near the event horizon of Sgr A* is likely responsible for bright outbursts in soft X-rays (Baganoff et al. 2001), but this scenario struggles to explain the steady emission at much higher energies. Alternative explanations for the GeV and TeV flux include the supernova remnant Sgr A East (Crocker et al. 2005), though this is strongly disfavored based on its observed offset from the very high energy emission centered on Sgr A* (Acero et al. 2010). Secondary emission from particles accelerated by Sgr A* is another candidate, either in a steady state or from a past burst of accretion (e.g. Atoyan & Dermer 2004;Aharonian & Neronov 2005;Chernyakova et al. 2011). Most of these scenarios cannot account for both the GeV and TeV emission. In contrast, a population of ∼1000 MSPs in the inner few pc could account for the emission from GeV through 100 TeV (Bednarek & Sobczak 2013). None of these scenarios seek to explain the GeV excess extending several kpc from Sgr A*. The pulsar population in the Galactic center has long been sought to test the theory of gravity (Pfahl & Loeb 2004;Liu et al. 2012, and references therein) and the existence of intermediate mass black holes and gravitational waves (Kocsis et al. 2012). Extended multiwavelength observations were conducted which should have detected a significant fraction of the most common second-period pulsars, but only four were seen. This missing pulsar problem indicates that the formation and/or retention of ordinary pulsars may be inefficient in this region Macquart & Kanekar 2015). However, these searches did not significantly constrain the number of MSPs, especially at the relatively large galactocentric distances of 0.1-1 kpc where the gamma ray excess is observed.
In this paper, we argue that the MSPs needed to produce the gamma ray excess were not made under the present conditions of the Galactic bulge, but were produced in dense globular clusters that have since dis-solved. The population of globular clusters constitutes a key component in the theory of galaxy evolution, the formation of galactic bulges, and nuclear star clusters (Tremaine et al. 1975;Arca-Sedda & Capuzzo-Dolcetta 2014;Capuzzo-Dolcetta 1993;Lotz et al. 2001;Bekki et al. 2004;Capuzzo-Dolcetta & Miocchi 2008;Capuzzo-Dolcetta & Mastrobuono-Battisti 2009;Agarwal & Milosavljević 2011;Hartmann et al. 2011;Leigh et al. 2012;Antonini et al. 2012;Antonini 2013Antonini , 2014den Brok et al. 2014;Kruijssen 2014;Perets & Mastrobuono-Battisti 2014;Antonini et al. 2015b,a). The clusters we see today are the ones that have survived throughout the evolution of the Galaxy, and may be a small fraction of the initial cluster population.
We organize the paper as follows. In Section 2, we discuss the relationship of MSPs to globular clusters and the evolution of the Galaxy's population of globular clusters. Section 3 discusses the scaling of the gamma ray luminosity to the predicted population of disrupted clusters, while Section 4 presents the predictions of this model for the Fermi excess. Sections 5 and 6 discuss two objections to a MSP explanation of the excess, the high end of the luminosity function and the average spectrum. We discuss prospects for radio detections of our predicted MSPs in Section 7 and conclude with Section 8.

MILLISECOND PULSARS FROM DISRUPTED
GLOBULAR CLUSTERS MSPs are thought to be "recycled" pulsars, spun up by the accretion of material from a close binary companion (Bhattacharya & van den Heuvel 1991). These close binaries are formed and driven to smaller separations in dense stellar environments where the rate of stellar dynamical encounters is high. Once these interactions have sufficiently decreased the binary separation, the neutron star's companion transfers material and angular momentum, and reduces the neutron star's magnetic field. This phase lasts ∼10 7 -10 9 years and is visible, for lowmass companions, as a low mass X-ray binary (LMXB) (Ivanova et al. 2008). A long-lived MSP remains after the mass transfer stops. While the strong magnetic braking in ordinary pulsars leads to a rapid spindown, MSPs persist for ∼10 10 years (Bhattacharya & van den Heuvel 1991); these pulsars can long outlive their birth clusters.
The highest abundances of MSPs are found in globular clusters, the old dense stellar islands orbiting in the galactic halo. There are several indications that a fraction of the stellar mass of galactic bulges may have been formed by dissolving globular clusters (Tremaine et al. 1975;Arca-Sedda & Capuzzo-Dolcetta 2014;Gnedin et al. 2014). The distribution of old globular clusters within galaxies and the geometry of galactic bulges are both approximately spherical. The number density of globular clusters increases inwards on kpc scales, but shows a relative decrease within the galactic bulge. A tight correlation is observed between the mass of the galactic bulge and the number of globular clusters (Harris et al. 2014).
Massive globular clusters spiral in towards the Galactic center due to dynamical friction. In the central kpc, the tidal gravitational field of the Galaxy may exceed the attractive field of the cluster stars, stripping the cluster from its outskirts and eventually down to its core.
The cluster then spills its entire contents, including the MSPs in its core, into a spherical shell about the Galactic center. Since MSPs are long-lived, they remain bright in gamma rays after the cluster is disrupted. The high dynamical encounter rates needed to form new LMXBs and new MSPs are, however, strongly suppressed after the disruption of the globular cluster. Therefore, the MSP population will be frozen at the time and orbit of the cluster's disruption while LMXBs (precursors to MSPs) will burn out within ∼10 8 years. As a result, the ratio of LMXBs to MSPs from disrupted globular clusters becomes much lower in the galactic bulge than in the surviving globular clusters we observe today. Indeed, LMXBs are observed to be rare in the bulge (Revnivtsev et al. 2008;Cholis et al. 2015).
We model the distribution of globular clusters and the Galactic bulge following Gnedin et al. (2014), who account for mass loss from passive stellar evolution, cluster evaporation, infall due to dynamical friction, and tidal disruption; we adopt all of their fiducial parameters. This simple model was set up to reproduce the radial and mass distribution of extant globular clusters in the halo with no eye towards reproducing the Fermi excess. The distribution of mass from dissolved globular clusters is shown in Figure 3 of Gnedin et al. (2014); we use this result directly. It is a cored radial profile with ρ(r) ∼ r −2.2 and an enclosed mass ∼10 8 M ⊙ at 1 kpc. We assume that most gamma-ray sources in globular clusters formed sufficiently quickly that the luminosity per unit mass in a population of disrupted clusters may be approximated by the observed value in surviving globular clusters. Their radial distribution is set by their orbits at their points of disruption in the model. We do not tune the fiducial model of Gnedin et al. (2014); our approach has no free parameters.
The Appendix gives a brief overview of the Gnedin et al. (2014) model and calculations, with equations giving the characteristic disruption timescales. We refer the reader to that paper for a more thorough discussion.
3. SCALING THE GAMMA RAY LUMINOSITY We compute the gamma ray luminosity per unit stellar mass for the globular clusters studied by Abdo et al. (2010) and listed in Table 1. Of these eleven clusters, Abdo et al. (2010) reported gamma ray detections for eight. We use the more recent 2 GeV fluxes measured by Cholis et al. (2015); this data set includes fluxes for two of the three globular clusters that were previously undetected. We adopt the cluster distances compiled by Abdo et al. (2010) and convert the absolute V magnitudes of Harris 1996Harris (2010 to mass by assuming a mass-to-light ratio of 3 M ⊙ /L ⊙ , appropriate for an old, slightly metal-poor population (Maraston 2005). The mass of Terzan 5 is uncertain due to its very large extinction (Lanzoni et al. 2010), with those authors favoring a value 2 × 10 6 M ⊙ compared to the 3 × 10 5 M ⊙ implied by the absolute V magnitude given by Harris (1996), but this has little effect on our results.
We neglect systematic variations in the number of MSPs per unit globular cluster mass as a function of the cluster properties and evolutionary stage. The formation rate and the total number of MSPs are observed to correlate with the rate of encounters, not simply the cluster mass (Hui et al. 2010;Bahramian et al. 2013), while the encounter rates are strongly affected by core collapse and the primordial binary fraction which are not well understood (Binney & Tremaine 2008). Indeed, the clusters listed in Table 1 have only a weak correlation between stellar mass and gamma ray luminosity. Future studies are needed to examine more detailed models of the MSPs within a population of globular clusters.

THE PREDICTED FERMI EXCESS
With an independent model of the population of disrupted globular clusters (Gnedin et al. 2014) and a scaling of total globular cluster mass to gamma ray luminosity (Section 3), we may compute the predicted Fermi GeV signal. Our results for the integrated flux within a circular region around the Galactic center, and the approximate number of enclosed MSPs, are shown in Figure 1 as a function of the angular distance to the center. Figure 2 shows the differential flux within circular annuli. Each figure shows the recent measurements of the Fermi excess to be in excellent agreement with our predictions.
The integral flux ( Figure 1) contains two particularly noteworthy results. First, the red line in Figure 1 is a discontinous set of dark matter profiles fit at different annuli: the best annihilation fit to the data is not a physically meaningful dark matter profile. At separations >0.   (2013), Daylan et al. (2014) and Abazajian et al. (2014), compared to the prediction (solid blue curve) from disrupted globular clusters, assuming the same gamma ray luminosity per unit mass as for intact clusters. The yellow hatching shows a factor of two uncertainty; the right axis shows the approximate number of enclosed MSPs. The gamma ray signal from scaling the disrupted globular clusters of Gnedin et al. (2014) correctly predicts all measurements, including an unresolved source around Sgr A* seen by Abazajian et al. (2014), with no free parameters. The black open stars include Sgr A* and the filled pentagons exclude it; we have also shown the disrupted globular cluster prediction with the inner 0. • 1 masked for comparison (blue dotted-dashed curve). We interpret this unresolved source as emission from MSPs in a nuclear star cluster. Daylan et al. (2014) and Gordon & Macías (2013) include this unresolved source in their diffuse fits.  Calore et al. (2015), compared to the prediction from scaling the Gnedin et al. (2014) disrupted globular clusters to the same gamma ray luminosity per unit mass as for intact clusters. We have included a factor of two uncertainty (yellow hatching) on the globular cluster prediction (blue curve). The blue curve and yellow hatching are not fitted to the data; they include no free parameters.
ters provide a physical motivation for the form of this signal. Second, Abazajian et al. (2014) fit for an unresolved source at the Galactic center, Sgr A*, which can be explained as MSP emission in a nuclear star cluster. We have added this central source to the diffuse emission fitted by Abazajian et al. (2014) to obtain the open black stars in Figure 1. We have placed the left-most black star at 0. • 1, roughly the angular resolution of the Fermi telescope. The actual extent of the nuclear star cluster is predicted to be just a few pc, or ∼0. • 02. To facilitate a direct comparison, we have also removed emission from the inner 0. • 1 of the disrupted cluster prediction, and indicated the resulting diffuse gamma ray signal with a dotted-dashed blue line. Our prediction matches the gamma-ray flux at all separations.
The existence of strong gamma ray emission from a nuclear star cluster reported by Abazajian et al. (2014) further supports the disrupted globular cluster hypothesis. Faucher-Giguère & Loeb (2011) found that the total encounter rate in the inner parsec of the Galaxy is similar to that of the globular cluster Terzan 5, so the formation rate of MSPs (and resulting gamma ray luminosity) may be similar. The 2 GeV flux of Terzan 5 would be just 5 × 10 −9 GeV cm −2 s −1 at 8.3 kpc. Unless the gravitational potential of Sgr A* retains a much higher number of neutron stars than in globular clusters, this scenario can only explain ∼10% of the flux seen by Abazajian et al. (2014) and shown in Figure 1. MSPs deposited in the nuclear star cluster by massive globular clusters were also suggested by Bednarek & Sobczak (2013) as an explanation of the observed TeV photons from around Sgr A*. Those authors required a population of ∼1000-3000 MSPs, fully consistent with our results both in the Galactic center and at larger separations. Recent 20-40 keV X-ray observations by the NuS-TAR satellite can also be explained by a large population of MSPs in the inner few pc (Perez et al. 2015).

THE MAXIMUM LUMINOSITY OF MSPS
One objection to a population of bulge MSPs as the source of the Fermi excess is the paucity of individually identified high luminosity gamma ray pulsars detected as point sources within ∼10 • of the Galactic center . Lee et al. (2015) found evidence that the GeV excess is from unresolved point sources with a 1.9-12 GeV flux cutoff at ∼1.5-2 × 10 −10 photons cm −2 s −1 , for a 0.1-100 GeV luminosity of ∼1.5-2.5 × 10 34 erg s −1 at 8.3 kpc. Cholis et al. (2014) found a very hard luminosity function for MSPs, with most of the luminosity contributed by objects above a few 10 34 erg s −1 , which should have been detected as Fermi point sources. We reanalyze the data of Cholis et al. (2014) and reexamine the cutoff at high luminosities. Cholis et al. (2014) use a sample of 59 field MSPs to determine the luminosity function (listed in their Table  IV). The distances they adopt are taken from the ATNF pulsar database (Manchester et al. 2005), available at http://www.atnf.csiro.au/people/pulsar/psrcat/, which uses the model Galaxy of Taylor & Cordes (1993), hereafter TC93, to convert observed dispersion measures into distances.
The Fermi Second Pulsar Catalog (2PC) distances (Abdo et al. 2013)   |b| > 10 • to determine the luminosity function; the revised distances thus have a large effect on the results. Table 2 lists the NE2001 distances for all pulsars without 2PC distances listed in Cholis et al. (2014). In the case of J1843−1113, Hobbs et al. (2004) have provided a distance using the TC93 model Galaxy. We have taken their dispersion measure and converted it into the tabulated NE2001 distance estimate. One MSP, J1137+7528, does not appear in any database. The only reference to this pulsar that we were able to find (apart from the Fermi Third Point Source catalog, Acero et al. 2015) was in a table of pulsar data available at http://astro.phys.wvu.edu/GalacticMSPs/GalacticMSPs.txt. This site lists a dispersion measure of 29.2 pc cm −3 which, according to the NE2001 free electron model, implies a distance of 1.5 kpc.
We also note that these new distances have significant uncertainties. The brightest MSP in our sample, for example, is J0614−3329, which was identified with Fermi by Ransom et al. (2011). This is a generally unremarkable MSP apart from the very high gamma ray efficiency (greater than unity) implied by its dispersion measure distance of 1.9 kpc. As Ransom et al. (2011) note, J0614 is nearly tangent to the Gum nebula where the NE2001 model has a very steep gradient in dispersion measure (and hence, in derived distance); those authors suggest that the true distance is likely to be closer by a factor of 2 or more. Decreasing its distance by such a factor would reduce its integrated gamma ray luminosity to ∼10 34 erg s −1 and remove the need to invoke a very large gamma ray efficiency and/or beaming factor.
Other pulsars, including those just below the cutoff, could also have distance errors. Random errors, however, will tend to smear out a distribution; a deconvolution should sharpen the cutoff. Also, an anomalously  Cholis et al. (2015) sample with |b| > 10 • adopting the Taylor & Cordes (1993) (TC93, blue line) and an updated model (Cordes & Lazio 2002 of the Galaxy's dispersion measure. The latter leaves just two MSPs sufficiently luminous to detect as point sources near the Galactic center. One, J0614−3329 is likely to be at least a factor of 4 less luminous than shown (Ransom et al. 2011); the other, J0218+4232, is just marginally (∼1σ) more luminous than the Galactic center detection threshold (Abdo et al. 2013). The vertical shading shows the luminosity cutoff needed to reproduce the gamma ray excess with an unresolved source population (Lee et al. 2015); Bartels et al. (2015) favor a somewhat higher cutoff.
large distance (by a factor of 2, say) requires simply an unmodeled clump of ionized gas along the line-of-sight. An anomalously small distance requires a void or bubble; the required void for the same fractional distance error becomes larger with increasing distance. A factor of 2 error in a 1 kpc distance would require a 1 kpc completely empty void (to go along with the 1 kpc of properly modeled free electron density). Figure 3 shows the cutoff at high luminosities in the Cholis et al. (2015) MSP sample using both the TC93 and the NE2001 distances. For reference, we have also shown the approximate luminosity cutoff favored by Lee et al. (2015) of ∼1.5-2 photons cm −2 s −1 between 1.9 and 12 GeV, transformed into 0.1-100 GeV flux assuming the well-measured spectrum of J0614−3329 and placed at a distance of 8.3 kpc. There are only two MSPs above this cutoff. One, J0614−3329 itself, was discussed above; Ransom et al. (2011) argue that its true distance is likely to be at least a factor of ∼2 smaller than that implied by its dispersion measure. The other is J0218+4232, for which Abdo et al. (2013) quote a factor of 2 error in the luminosity as a result of dispersion measure uncertainties in the NE2001 model, making it more luminous than 2 × 10 34 erg s −1 by just ∼1σ. The updated pulsar distances support the luminosity cutoff needed by Lee et al. (2015) and Bartels et al. (2015) to account for the GeV excess with unresolved point sources. Bartels et al. (2015) favor a cutoff at higher luminosities than Lee et al. (2015) (though over a different bandpass).

THE AVERAGE SPECTRUM OF MSPS
We now turn to the spectrum of an unresolved distribution of Fermi MSPs, estimating the integrated light using the field MSPs listed in Cholis et al. (2014). The spectrum of MSPs has been suggested to differ modestly from that of the observed GeV excess , arguing against their ability to produce the gamma rays around the Galactic center. We use the field MSPs under the assumption that many of them have similar origins to the one we suggest for the central bulge population, in the cores of stellar clusters disrupted long ago. We also account for the fact that Fermi's sensitivity is a strong function of photon energy (Acero et al. 2015). Following Cholis et al. (2014), we do not apply a cut in Galactic latitude b (as we and they did for the high luminosity cutoff, Section 5). Applying such a cut would make the spectra we show in this section slightly harder, and lessen the tension with the Daylan et al. (2014) spectrum of the GeV excess.
The MSP sample of Cholis et al. (2014) is not selected at a single frequency: for a fixed 2 GeV flux density with little emission from higher energy photons, very soft sources are easier to detect than harder sources. Only the soft sources supply enough photons at energies <1 GeV to contribute significantly to their detectability. We therefore create a 2 GeV-selected sample from the Cholis et al. (2014) MSPs by including only those with a 1-3 GeV test statistic of 8 2 , corresponding to a signal-to-noise ratio of at least 8 in this band. Of the 59 MSPs, 45 pass this cut. This sample of 45 MSPs should be complete independently of spectral shape.
Under the (unlikely) assumption that our field MSPs form a flux-limited sample of a spatially uniform population, we can derive a weighted average of the individual MSP spectra that matches the spectrum expected for a population at fixed distance. Uniformly weighting each MSP's spectrum is equivalent to assuming a volume-limited survey. These two scenarios, fluxlimited/spatially uniform and volume-limited, almost certainly bracket the truth. The rest of our calculations use an average or composite of these two limiting cases.
Assuming a flux-limited survey, the total gamma-ray flux from Fermi-resolved MSPs in a luminosity range where dN/dL is the differential number density and r max ∝ L/F min is the maximum radius out to which the source could be detected. For a source population around the Galactic center, we wish to know the total luminosity per unit volume, which is simply given by We combine Equations (1) and (2) to obtain The integrated spectrum from an unresolved MSP population near the Galactic center may therefore be estimated by scaling the observed flux density of each object by that object's luminosity to the −1/2 power, and simply adding all of the scaled flux densities together. Figure 4 shows the results, with all spectra normalized to their 1.7 GeV values. The blue dotted-dashed line shows the result for simply adding together all MSP spectra without selecting them by 1-3 GeV flux and without scaling by luminosity. The red line is the average of the spectra with and without weighting by L −1/2 , i.e, assuming volume-limited and flux-limited samples, respectively. The blue and orange hatching show the 1σ and 2σ uncertainties in the red spectrum as estimated from bootstrap resampling of the 45 MSPs. For this exercise, we have adopted the fitted spectra in Table I of Cholis et al. (2014) and have neglected measurement errors, fitting errors, and distance errors.
The difference in Figure 4 between the scaled and unscaled spectra results from a correlation between luminosity and spectral index. Distance errors will tend to blur this correlation; the MSP spectrum of a population at a single distance is likely to be slightly harder than the red line in Figure 4. Including this effect and adding measurement errors would not bring the MSP spectrum into perfect agreement with the Galactic center excess, but it could bring the 1σ discrepancy to as little as ∼20-30% at 500 MeV. Selecting only those MSPs with |b| > 10 • (38 of the 45 that pass our 1-3 GeV signal-to-noise cut) would also marginally improve the agreement with the spectrum of the GeV excess.
The discrepancy between our estimated average MSP spectrum and the GeV excess is only significant at the lowest energies (<800 MeV) where Fermi's sensitivity is rapidly falling. Uncertainties in Galactic diffuse emission are largest here (Calore et al. 2015). As a result, there are spectrally correlated systematic errors in the spectrum of the GeV excess not shown in the black stars of Figure 4. Systematic errors can be quite large, and can also arise from the method of masking point sources and from the assumed morphology of the excess, among other aspects of the fitting (Daylan et al. 2014;Calore et al. 2015). Figure 4 also shows the systematic errors from varying the diffuse backgrounds as estimated by Calore et al. (2015). These gray and gold hatched regions neglect statistical errors.

PROSPECTS FOR RADIO DETECTIONS
Our results show that a population of disrupted globular clusters, which must exist to explain the current clusters, naturally predicts a field population of MSPs in the Galaxy's inner few kpc. These MSPs satisfy the spatial, spectral, and luminosity requirements imposed by the Fermi observations. A large population of MSPs in a nuclear star cluster is another necessary consequence of a population of disrupted massive globular clusters. Such a population explains the 20-40 keV X-ray emission seen by NuSTAR (Perez et al. 2015) and implies that many of the unidentified Chandra point sources may be MSPs (Muno et al. 2004;Perez et al. 2015). Astro-H (Takahashi et al. 2010) will also be sensitive to highenergy X-rays, and could confirm the NuSTAR results. A population of ∼1000 MSPs around Sgr A* can also explain the observed TeV emission by inverse Compton scattering of the dense interstellar radiation field (Bednarek & Sobczak 2013).
Radio observations could individually detect our predicted MSPs and confirm their identities. However, the bulk of the radio observations to date have focused not on scales of tens to thousands of pc, where most of our predicted MSPs lie, but in the innermost pc. This was motivated by theoretical estimates predicting ∼100-1000 pulsars formed in situ within , and is the average of the spectra expected for a population at uniform distance assuming the Cholis et al. (2014) to be volume-limited and flux-limited. These scenarios almost certainly bracket the truth. The blue and orange hatching show 1σ and 2σ sample variances as estimated using bootstrap resampling. We have neglected errors in the MSP distances and in the spectral measurements; both would tend to alleviate the discrepancy with the observed Galactic center excess (Daylan et al. 2014). The error bars on the Daylan et al. (2014) fits are only statistical; systematic errors (which are spectrally correlated) are neglected. The gold and gray hatching show 1σ and 2σ systematic uncertainties (neglecting statistical errors) as estimated by Calore et al. (2015). 0.02 pc of Sgr A* (Pfahl & Loeb 2004). More recently, Faucher-Giguère & Loeb (2011) noted that the encounter rate in the inner 1 pc of the central star cluster is comparable to that of the globular cluster Terzan 5 (which has many MSPs), and estimated that up to ∼1200 MSPs may be present in this region due to the deeper gravitational potential well of Sgr A*. The disrupted globular cluster scenario instead predicts these MSPs to be found over a larger region: we predict ∼1,000 MSPs within 3 pc of Sgr A*, and a further ∼1,000 MSPs within 300 pc (2 • , see Figure 1).
MSP observations towards the Galactic center are extremely challenging because of the large dispersion measures. Radio pulses at a frequency ν are broadened by an amount τ = (1.3 ± 0.2)(ν/GHz) −3.8±0.2 (with τ in seconds, Spitler et al. 2014), implying that MSPs may not be observed below ∼8 GHz. The radio intensity of pulsars scales steeply with frequency (I ∝ ν −1.6 to ν −1.8 , Kramer et al. 1998), so high-frequency detections require extended integration times.
While discovering and timing MSPs 0.001 pc from the central supermassive black hole would offer tantalizing measurements of general relativity and tests of alternative theories of gravity (Wex & Kopeikin 1999;Kramer et al. 2004;Cordes et al. 2004;Pfahl & Loeb 2004;Liu et al. 2012), discovering MSPs further out within 10 pc would also be invaluable. Such MSPs could be used to measure the properties of the nuclear star cluster, find intermediate mass black holes, and measure the gravitational waves of the Galactic center (Kocsis et al. 2012). If the nuclear star cluster indeed formed from disrupted globular clusters, we predict ∼1000 MSPs within ∼10 pc of the Galactic center; such a population of MSPs can account for the unresolved Fermi flux seen by Abazajian et al. (2014). Future high frequency radio surveys will have the frequency coverage and sensitivity needed to detect this MSP population (Chennamangalam & Lorimer 2014;Macquart & Kanekar 2015). Even larger radio surveys such as the square kilometer array (SKA) on 100 pc to 2 kpc scales are required to confirm or disprove the disrupted globular cluster origin of MSPs in the Galactic bulge.

CONCLUSIONS
The Fermi Galactic center excess is in excellent agreement with independent predictions of the population of MSPs produced in disrupted globular clusters. This astrophysical model appears to fit the observations as well as dark matter annihilation, but without any free parameters. MSPs from disrupted clusters also provide an excellent match to the observed emission near Sgr A* from hard X-rays through very hard gamma rays. If the bulge indeed contains a large population of stars from long-dead clusters, such MSPs form a background that must necessarily be present in the Fermi data.
The observed emission extends at least ∼2 kpc from the Galactic center (Hooper & Slatyer 2013), far from the nuclear star cluster around Sgr A* where dynamical formation of MSPs is plausible. LMXBs burn out after the disruption of globular clusters, reducing their relative numbers in the galactic bulge, consistent with the lack of LMXB observations (c.f. Cholis et al. 2015). We conclude that the dominant MSP population is not likely to have formed under the current conditions in the bulge, but was deposited by dissolving globular clusters. If the Fermi excess is indeed the relic of a previous large population of globular clusters, it provides the first direct evidence for their existence, and strongly supports the theory for the globular cluster origin of the nuclear star cluster. Future radio observations may be directly sensitive to these MSPs and could offer decisive evidence of a broad distribution of MSPs deposited by globular clusters.
While our results disfavor a dark matter interpretation of the GeV excess, they show that Fermi can offer a new probe of the formation history of the bulge, and of the evolution of the Galaxy's globular cluster system. Our reevaluation of field MSP luminosities, combined with the results of Lee et al. (2015) and Bartels et al. (2015), suggest that we will soon begin to resolve the brightest of these fossils.
The authors thank Scott Tremaine, Doron Kushnir, Mariangela Lisanti, Neal Dalal, and Oleg Gnedin for very useful discussions, and an anonymous referee for helpful suggestions. B.K. gratefully acknowledges support from the W.M. Keck Foundation Fund of the Institute for Advanced Study, NASA grants NNX11AF29G and 13-ATP13-0056, and NSF grant AST-1406166. This work was performed in part under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. This work was supported in part by the European Research Council under the European Union's Horizon 2020 Programme, ERC-2014-STG grant GalNUC 638435. relaxation times (Gnedin et al. 1999), which is between 0.1 and several Gyr, longer than the orbital time within 2 kpc (see Eq. A1). The tidal debris from globular clusters overlapping with the gamma-ray excess may remain spherical in a wider class of models than Gnedin et al. (2014).