Strong Evidence that the Galactic Bulge is Shining in Gamma Rays

TL;DR

This paper presents strong evidence that the gamma-ray excess in the Galactic Center is correlated with the stellar distribution of the Galactic bulge, supporting the hypothesis that faint astrophysical sources like millisecond pulsars are responsible.

Contribution

It demonstrates a significant morphological correlation between gamma-ray excess and the Galactic bulge, and constrains MSP formation scenarios, favoring an admixture of primordial and dynamically formed binaries.

Findings

Gamma-ray excess morphology traces the Galactic bulge stellar over-density.

Admixture formation scenario for MSPs is favored over primordial-only at 7.6σ.

Evidence of a high energy tail in the nuclear bulge spectrum suggesting inverse Compton emission.

Abstract

There is growing evidence that the Galactic Center Excess identified in the $\textit{Fermi}$-LAT gamma-ray data arises from a population of faint astrophysical sources. We provide compelling supporting evidence by showing that the morphology of the excess traces the stellar over-density of the Galactic bulge. By adopting a template of the bulge stars obtained from a triaxial 3D fit to the diffuse near-infrared emission, we show that it is detected at high significance. The significance deteriorates when either the position or the orientation of the template is artificially shifted, supporting the correlation of the gamma-ray data with the Galactic bulge. In deriving these results, we have used more sophisticated templates at low-latitudes for the $\textit{Fermi}$ bubbles compared to previous work and the three-dimensional Inverse Compton (IC) maps recently released by the ${\tt GALPROP}$ team. Our results provide strong constraints on Millisecond Pulsar (MSP) formation scenarios proposed to explain the excess. We find that an $\textit{admixture formation}$ scenario, in which some of the relevant binaries are $\textit{primordial}$ and the rest are formed $\textit{dynamically}$, is preferred over a primordial-only formation scenario at $7.6\sigma$ confidence level. Our detailed morphological analysis also disfavors models of the disrupted globular clusters scenario that predict a spherically symmetric distribution of MSPs in the Galactic bulge. For the first time, we report evidence of a high energy tail in the nuclear bulge spectrum that could be the result of IC emission from electrons and positrons injected by a population of MSPs and star formation activity from the same site.

Equations

$${\rho }_{\mathrm{bar}}(R,\varphi ,z)$$

$$\displaystyle\quad\times\begin{cases}1,\,\,\,\mbox{$R\leq R_{\rm end}$}\\ e^{-[(R-R_{\rm end})/h_{\rm end}]^{2}},\,\,\,\mbox{$R>R_{\rm end}$},\end{cases}$$

$$R_{\perp }^{C_{\perp }}$$

$$R_s^{C_{||}}$$

$$\begin{gathered} \text{mbox}bulge\propto \exp (-\{{\left(\frac{(b-b_0)\cos (\theta )+(l-l_0)\sin (\theta )}{\beta }\right)}^2 \\ {+[(l-l_0)\cos (\theta )-(b-b_0)\sin (\theta ){]}^2\}}^{1/2}/\alpha ), \end{gathered}$$

$$p(\mathrm{TS})=2^{-n}\left(\delta (\mathrm{TS})+\sum _{i=1}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right){\chi }_i^2(\mathrm{TS})\right)$$

Full text

1 Introduction

The center of the Milky Way Galaxy is a complex environment displaying an extremely rich phenomenology of energetic cosmic-ray production and interactions. Many types of transient and persistent cosmic-ray source candidates populate the region, including a supermassive black hole [1], active star formation [2], supernova remnants [3], pulsars and black holes [4], as well as new physics involving particle dark matter annihilations or decays (see, e.g., [5] for a review). The bulk of the gamma-ray signal from the Galactic Center arises from the interactions of energetic cosmic rays with the interstellar gas and radiation fields. However, there is growing evidence that there exists an excess of extended gamma rays observed on top of models of cosmic-ray interactions in the Galactic Center. Dubbed the Galactic Center Excess (GCE), this signal has been detected by multiple analyses of the Fermi-LAT data [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

The origin(s) of the GCE remain elusive. The main properties of the GCE—being centrally peaked with a roughly spherically symmetric spatial morphology, having an energy spectrum peaking at a few GeV, and total luminosity of $\sim 10^{37}$ erg/s—can be accommodated within a scenario where weakly interacting massive particle (WIMP) dark matter of $\mathcal{O}(10\,{\rm GeV})$ mass distributed in a slightly steep Navarro-Frenk-White (NFW)-like profile self-annihilate. However, the GCE might also arise from a population of gamma-ray emitting pulsars, in the form of either old “recycled” millisecond pulsars [18, 13], or young pulsar remnants from Galactic Center supernovae [19, 20]. The gamma-ray spectra of known pulsars are similar to the GCE, and the capture of pulsars onto the Galactic Center region by Globular Cluster disruptions could explain the GCE’s quasi-spherical spatial distribution [21, 22]. Other explanations for the GCE discussed in the literature include outbursts of cosmic-ray production by the central supermassive black hole (e.g., [23]).

While pulsars are a natural astrophysical candidate, there is ongoing debate regarding the consistency of the luminosity function required to explain the GCE with that measured for the pulsar population observed elsewhere as point sources [24, 25, 26, 27]. Nevertheless, there are strengthening indications that the GCE may have an astrophysical origin. Cosmic-ray interactions in the Galactic Center region are notoriously complex to model. New analyses have cast doubt on some of the main properties claimed previously for the GCE. First, the energy spectrum of the GCE is subject to large systematic uncertainties arising from the incomplete understanding of cosmic-ray interactions in the Galactic Center region (e.g., [13, 14, 15, 17]). A variety of scenarios might, therefore, describe the GCE based on the energy spectrum alone. Second, it has been argued that the photon count distribution of the GCE is more consistent with arising from a population of faint point sources rather than being truly diffuse [28, 29]. This observation supports an astrophysical origin in the form of faint point sources over a microphysics origin in the form of dark matter particle annihilation. Third, there are growing reports of a spatial variation of the GCE energy spectrum with latitude [30] showing a high-energy tail away from the plane, similar to that recently reported by Ref. [31]. Again, this disfavors a dark matter origin because particle properties would be expected to be spatially invariant (although secondary gamma-ray mechanisms might remain spatially-dependent).

Finally, it has recently been argued that the spatial morphology of the GCE is better matched to the asymmetric (stellar) bulge of the Milky Way than the spherically symmetric distribution expected for dark matter [32, 33]. The Galactic bulge is a triaxial bar-like structure in the central region of our Galaxy extending in size to a few kpc [34, 35, 36]. While the exact morphological properties of the bulge still remain under investigation, studies agree that it is asymmetric on our sky plane given that its major axis approaches us at positive Galactic longitudes (see, e.g., [37]). Reference [32] adopted two bulge maps to analyze the GCE: (i) The X-shaped Galactic bulge template of Ness & Lang [38], and (ii) the boxy shaped bulge of Freudenreich [35]. They found that the data overwhelmingly preferred either of the bulge maps over the spherically symmetric distribution expected by a dark matter origin. The bulge houses a broad mix of stellar populations from old to star forming [39, 40], and should contain ample candidates of astrophysical gamma-ray emitters.

In this paper, we put the reported bulge-correlation of the Galactic Center Excess to further tests. First, we include the recent GALPROP three-dimensional (3D) Inverse Compton (IC) maps constructed in Ref. [41] in our pipeline. The new IC maps do not assume the Galactocentric cylindrical symmetry adopted in all two-dimensional (2D) IC maps constructed with previous GALPROP versions (v54 or older) and contain sophisticated model templates for the bulge/bar, spiral arms and stellar disk. Second, compared to previous work in Ref. [32], our region of interest (RoI) is about a factor of $7$ larger. Third, we translate and rotate the bulge template to investigate how robustly it is detected in the Fermi-LAT data. We find that the Galactic bulge is detected strongly in gamma rays only when consistently positioned and aligned within one or two degrees of the known stellar bulge position [35].

Fourth, we perform a morphological analysis of the stellar bulge distribution in order to test the different MSP formation scenarios that have been discussed in the literature. We consider an * admixture formation* scenario in which some of relevant binaries are primordial and the remaining ones are the result of stellar interactions. We find that the admixture formation scenario is preferred to primordial-only formation with a confidence of $7.6\sigma$. Finally, we present the gamma-ray spectra of the boxy bulge and the nuclear bulge. For the first time, we report evidence of a high-energy tail in the nuclear bulge spectrum that could be the result of IC emission from electrons and positrons injected by a population of MSPs and/or star formation activity from the same site. The dark matter implications of our new results will be presented in a separate article.

The structure of this article is as follows. In Sec. 2 we provide details about the data, model templates, and pipeline used in this study. In Sec. 3 we show our results for the empirical tests of rotation and translation of the stellar templates as well as our morphological analysis. In Sec. 4 we present the gamma-ray properties of the Galactic bulge. Finally, a discussion and our conclusions are given in Sec. 5.

2 Data and Methods

We first describe our modeling procedure, starting with a description of data selection, the templates used in our analysis, and finally the analysis method.

2.1 Data selection

We used $\sim 7$ years (August 4, 2008$-$September 4, 2015) of Pass 8 ULTRACLEANVETO class photons with reconstructed energy in the 667 MeV$-$158 GeV range. Photons detected at zenith angles larger than 90*∘* were excised to limit the contamination from $\gamma$-rays generated by cosmic-ray interactions in the upper layers of the atmosphere. Moreover, we filtered the data using the recommended specifications (DATA*-QUAL$>$0)&&(LAT-CONFIG==1). Fermi Science Tools v10r0p5 and instrument response functions (IRFs) P8R2-ULTRACLEANVETO-*V6 were used for this analysis. In addition, the analysis was restricted to a square region of $40^{\circ}\times 40^{\circ}$ centered at Galactic coordinates $(l,b)=(0,0)$ with a spatial binning of $0.2^{\circ}$.

2.2 Baseline model templates

We started with a similar baseline model to that developed in Ref. [32]. In summary, this consists of the point-like sources listed in the 2FIG catalog [45], specialized templates for the Sun and the Moon, an isotropic component (iso*-P8R2-ULTRACLEANVETO-V6-*v06.txt), a standard 2D IC template generated with GALPROP v56 [41], and an emission map for Loop I. In addition, the interstellar gas-correlated photons were modeled with a linear combination of atomic and molecular hydrogen gas templates divided into concentric rings (see Section 2.7). Furthermore, these gas maps were obtained from a suite of hydrodynamical simulations of interstellar gas flow. This contrasts with gas maps in GALPROP and the official Fermi diffuse emission model that are constructed with an interpolation technique. Reference [32] showed in detail that there are important morphological differences between the interpolated and hydrodynamic gas maps and that the latter provides a significantly better fit to the gamma-ray data in the Galactic Center region. See Table 1 for a summary of the model templates considered in this work.

2.3 Bulge templates

Analyses [48, 49, 35] of the near-infrared emission measured by the DIRBE instrument on board of the COBE satellite were the first to uncover the non-spherical nature of the Galactic bulge. More recent studies of the Galactic bulge stars [50, 51] have firmly established that the bulk of bulge stellar mass forms a Boxy/Peanut (B/P) spatial morphology (this geometry is also sometimes called an “X-shape” in the literature). Using $N$-body simulations of disk galaxy formation [52], it has been found that the bulge hosts a rapid initial starburst component, as well as a second component which forms after disk build-up through dynamical instabilities [53]. The former is thought to be associated with a classical bulge structure111Classical bulges are known to have roughly spherically symmetric morphologies. while the latter with the B/P structure. Studies [54] of the internal kinematics of the stars of the Galactic bulge have determined that these are consistent with at least 90% of the bulge mass being contained in the B/P structure. Interestingly, the detection of a classical bulge component in the Milky Way has been claimed in [55]. However, as noted above, this is expected to harbor a sub-dominant fraction of the Galactic bulge stellar mass.

Further evidence for an X-shaped or B/P bulge morphology was deduced from the bifurcation in red clump giant counts in the 2MASS and OGLE-III surveys [56, 57, 58]. The analysis in Ref. [38] found, by constructing an image of the Milky Way bulge from an independent co-adding of publicly available WISE data, that the Milky Way bulge shows a distinct X-shaped structure in the projected stellar distribution that is different from previous works in that it has a more pronounced pinch-in effect at longitude ($l\sim 0^{\circ}$).

It has been posited in previous studies (e.g. [18]) that the GCE could be due to a population of unresolved sources such as MSPs, which are known to have GeV-peaked gamma-ray spectra. Since MSPs can be generated from old stellar populations, here we test several different stellar templates for the Galactic bulge in the Fermi data.

2.3.1 Boxy/Peanut bulge

Reference [35] derived a parametric model for the spatial morphology of the Galactic bulge by fitting to COBE/DIRBE near-infrared ($1.25$–$4.9$ $\mu$m) data. The best-fitting model obtained in Ref. [35] was Model S, which is a ${\rm~{}sech~{}}$-squared function on the bar radial spatial profile. The density distribution of the boxy bulge for Model S is written as222Note that there was a typo in the argument of the $\exp$ function in Eq. (14) of Ref. [35] which has been corrected in our Eq. (2.1). We have confirmed this in private communication with H. Freudenreich.

[TABLE]

where $R_{s}$ is the effective radius:

[TABLE]

and ($X^{\prime}$, $Y^{\prime}$, $Z^{\prime}$) is the Boxy bulge frame of reference. The bulge axis scale-lenghts are given by $a_{x}=1.696$ kpc, $a_{y}=0.6426$ kpc and $a_{z}=0.4425$ kpc, the bulge cut-off radius is $R_{\rm end}=3.128$ kpc and the scale height $h_{\rm end}=0.461$ kpc. In addition, the bulge face-on and edge-on shape parameters are $C_{||}=3.501$ and $C_{\perp}=1.574$ respectively.

In order to transform from the Galactocentric frame333The Galactocentric $y$ axis is assumed to be positive in the Sun-Galactic Center direction and the $xy$-plane is parallel to the plane of the Galaxy. $(x,y,z)$, to the bulge frame ($X^{\prime}$, $Y^{\prime}$, $Z^{\prime}$), we need to perform two successive passive rotations of coordinates: the first is a clockwise rotation $\theta_{0}=13.79^{\circ}$ around the Galactocentric $z$ axis, and the second is a clockwise rotation of $\beta_{\rm bulge}=0.023^{\circ}$ around the new axis $y^{\prime}$. The Sun is assumed to be at $z_{0}=16.46$ pc and $R_{0}=8.5$ kpc in the Galactocentric frame. Lastly, we performed a line-of-sight integration of the density function in Eq. (2.1) and normalized the map appropriately for analyses within the Fermi Science Tools.

2.3.2 X-shaped bulge

We reproduced the method described in Ref. [38] and applied it to the raw infrared WISE image of the Milky Way made available in Ref. [38]. We employed the 3.4 (W1) and 4.6 (W2) micron WISE data and masked out pixels with negative flux, as well as the top and bottom 5% of pixels based on W1$-$W2 color. We fitted an exponential disk model to the W1 and W2 bands and then subtracted the best fit models from the W1 and W2 data, respectively. Namely, we assumed

[TABLE]

where $\alpha=1.2$  kpc, $\beta=0.35$, $l_{0}=-0.45^{\circ}$, $b_{0}=0.11^{\circ}$, and $\theta=0.06^{\circ}$. We then applied a median filter of radius $1.7^{\circ}$ to the masked versions of both residual templates. In order to include this template in our gamma-ray pipeline, all pixels values with negative fluxes were set to zero in each median filtered exponential subtracted map. Finally, we computed the average of the two resultant maps and normalized it to flux unity for analyses with the Fermi Science Tools.

2.3.3 Nuclear Bulge

Another stellar structure that is evident from infrared observations corresponds to the so called Nuclear Bulge (NB); this is thought to consist of a spherical Nuclear Stellar Cluster centered at the position of the supermassive black hole and a Nuclear Stellar Disk with radius $\sim 230$ pc and scale height $\sim 45$ pc. The gamma-ray emission from this stellar component was also studied in Ref. [32]. This template was constructed from a near-infrared stellar density measurement [46] of the central region of our Galaxy ($|l|\geq 3^{\circ}$ and $|b|\geq 1^{\circ}$). More details about this template can be found in Ref. [32].

2.4 Spherically symmetric template

In addition to the stellar bulge templates, we also consider a map describing the potential gamma-ray emission from dark matter annihilation or a spherically symmetric distribution of MSPs. We model this by the square of a NFW density profile with an inner slope of 1.2, which has been shown to describe the GCE well in previous works (see Ref. [10] for details on the precise parameters choices used in this work).

2.5 Fermi Bubbles

Reference [42] carried out a reanalysis of the low latitude counterpart of the Fermi bubbles (FB). The residual map found in that study has a spatial morphology whose boundaries are well described by two catenary curves of the form $10.5^{\circ}\times(\cosh((l-1^{\circ})/10.5^{\circ})-1^{\circ})$ and $8.7^{\circ}\times(\cosh((l+1.7^{\circ})/8.7^{\circ})-1^{\circ})$ for the Northern and the Southern bubbles, respectively. This flat FB template (hereafter simply called “Catenary”) is displayed in the left panel of Fig. 1.

Recently, Ref. [43] derived an all-sky template for the Fermi bubbles using a spectral component analysis technique [60]. This method assumes that the FB spectrum at high latitudes is the same as the spectrum at low latitudes in the energy range $[1,10]$ GeV. We show their “Structured FB” 444The Strucutured FB map used in this work corresponds to Fig.12 of Ref. [43]. The authors of that reference constructed the map with their 3 component spectral decomposition technique. template for our RoI in the central panel of Fig. 1. it is evident that this FB template contains missing data (see green contours in Fig. 1) due to the need of a point source mask template.

In order to ameliorate the artifacts present in the Structured FB map, we have applied an inpainting method to the FB image. The inpainting algorithm utilizes Laplace Interpolation [59] which is a specialized interpolation technique for restoring missing data on a grid. Specifically, this algorithm solves the Poisson equation in the masked pixels of an image using as boundary conditions the pixel values in the unmasked portion of the image. The inpainted FB template is shown in the right panel of Fig. 1.

2.6 Inverse Compton Models

Reference [32] investigated the GCE systematics associated to the IC component with the use of the 2D Galactocentric cylindrically symmetric interstellar radiation field (ISRF) model attainable in GALPROP v54 [61]. The authors of Ref. [62] compared 128 different GALPROP models with all sky Fermi-LAT data finding a range of possible propagation parameter setups that are consistent with local CRs and gamma-ray measurements. Reference [32] reconstructed a subset of the IC maps in Ref. [62] that were found to have the largest morphological and spectral variations in order to estimate the impact of the IC component in the best-fit GCE emission.

In this article we improve upon the pipeline of Ref. [32] by using the two recently developed 3D ISRF models made publicly available in the most recent release of GALPROP v56 [41]. These are based in two different Galaxy wide dust and stellar distribution models labeled F98 and R12: the F98 model is based on the work by Freudenreich (1998) [35] and its bulge component is described in more detail in Sec. 2.3, while the R12 is based on the study by Robitaille et al. (2012) [47]. The corresponding spectral intensities for the ISRF were calculated by Ref. [41] with the Numerical Kalculation for Interstellar Emission (FRaNKIE) code. Although both models assume different stellar luminosities and dust densities, their predicted local intensities are consistent with near-infrared to far-infrared observations.

The IC maps investigated in [41] utilized three different CR propagation setups (see Table 3 of that publication). These were labeled SA0, SA50 and SA100 according to the proportion of CR luminosity injected by the spiral arms. We have reproduced the main IC skymaps presented in that study and included them in our step-wise statistical procedure described in the following section. In particular, for each Galaxy model (F98 and R12) we have utilized the three proposed CRs propagation setups. These are named: F98-SA0, F98-SA50, F98-SA100, R12-SA0, R12-SA50 and R12-SA100. Figure 2 displays the fractional residual differences between a standard 2D IC map extracted from Ref. [62] (Std-SA0) and the 3D model F98-SA0.

2.7 Analysis pipeline

In Ref. [32] some of us showed that the GCE spatial morphology was better explained by the stellar nuclear bulge and Boxy/Peanut (or X-shaped) bulge templates than by a spherically symmetric excess map given by an NFW2 profile. This was done by performing a set of maximum-likelihood fits summarized in Table I of that work [32]. In the present analysis we use a similar step-wise statistical procedure to establish whether a certain extended template for the GCE is statistically preferred over another.

The analysis pipeline used here to study the Fermi-LAT data utilizes the inner $40^{\circ}\times 40^{\circ}$ of the Galactic Center. Similar to Ref. [32], we employed a bin-by-bin analysis method in which the full Fermi-LAT data were divided into 15 logarithmically spaced energy bins where, for each energy bin, we performed a separate maximum-likelihood fit with pyLikelihood within the Fermi Science Tools. In each of these maximum-likelihood runs we simultaneously fitted the normalization of all included diffuse emission components while keeping their spectral slope as a fixed parameter. Due to the small size of the bins, our results were not sensitive to the assumed spectral shapes of the extended and point-like sources. For simplicity we therefore used power-law spectra with fixed index of negative two.

The point source modeling was done with the 2FIG catalog [45]. We considered only the point sources that were confidently found with the two different interstellar emission models considered in that work. There is a total of 288 such point sources inside our RoI. Due to limitations in the maximum number of parameters that can be reliably fitted in a given run with Fermi Science Tools, we have opted to follow a hybrid fitting approach analogous to that employed in Ref. [43]. Namely, we independently varied the normalizations of the 100 brightest 2FIG point sources in the RoI, while for the remaining point sources we constructed a single point source template in each energy bin. This is a reasonable approximation given that the amount of data utilized in the present study is very similar to the one used by the authors of the 2FIG catalog. To construct the combined point source template, we used the best-fit spectra in the 2FIG catalog and convolved it with the Fermi point spread function at each energy bin with the gtmodel tool. The resulting point source template was appropriately normalized and included in the fit along with the other extended and point-like sources. Other extended sources (W28, W30, RXJ1713, HESSJ1825-137, HESSJ1632-478 and HESSJ1837-069) that exist in our RoI were taken from the 3FGL catalog [63] and varied independently in the fits.

The significance of each source was evaluated using the test statistic $\mbox{TS}=2(\ln\mathcal{L}_{1}-\ln\mathcal{L}_{0})$, where $\mathcal{L}_{0}$ is the likelihoods of the background (null hypothesis) and $\mathcal{L}_{1}$ is the likelihood of the source being tested plus background. The total TS of a source is given by $\mbox{TS}=\sum_{i=1}^{n}\mbox{TS}_{i}$ where $\mbox{TS}_{i}$ is the TS for bin $i$ and there are $n=15$ bins. Since a source has now 15 degrees of freedom (non-negative normalization in each energy bin) instead of two in a standard broad-band analysis (non-negative normalization and index), the detection threshold condition has to be adjusted accordingly. As discussed in Ref. [32], this can be done by utilizing the formula given in case 9 of [64]:

[TABLE]

where $\delta$ is the Dirac delta function, $\binom{n}{i}$ is a binomial coefficient with $n=15$, and $\chi_{i}^{2}$ is a $\chi^{2}$ distribution with $i$ degrees of freedom. It follows (see Ref. [32]) that a threshold of $\mbox{TS}\geq 34.8$ corresponds to a 4$\sigma$ detection of a new extended source.

3 Results

Using the initial baseline model described in Sec. 2.2, the first step of our procedure is to find the best-fit normalization parameters for all free extended and point-like sources555The fitting of the 2FIG point sources in carried out with a hybrid approach described in Sec. 2.7. in the RoI. The results are shown in Table 2. In this first iteration we seek to find out which of the three alternative Fermi bubbles maps (see Fig. 1) is preferred by the data. Three fits were done, each time adding a different FB template: Catenary [42], Structured FB [43] and Structured FB corrected with the inpainting method. As summarized in this table, the Structured FB (Inpainted) map was the one found to improve the fit the most. Therefore, for the next iterations we add this FB template and re-optimize the baseline model.

In the next stage of our procedure we perform six fits in which we replace the standard 2D IC map by one of the new 3D IC templates: F98-SA0, F98-SA50, F98-SA100, R12-SA0, R12-SA50 or R12-SA100 (which were recently proposed in Ref. [41]; see Sec. 2.6). A comparison of the log-likelihood values obtained for each 3D IC template versus the standard 2D IC is displayed in Fig. 3 (larger is better). Interestingly, gamma-ray data from the inner $40^{\circ}\times 40^{\circ}$ prefer some of the 3D IC models to the Standard 2D one, with F98-SA0 being particularly favored. However, we note that although the new 3D IC maps are more physically motivated and much more realistic than the Standard 2D IC maps, the assumed CR propagation models and CR luminosities have not been tuned to match the gamma-ray data. Rather, they were obtained from a maximum-likelihood fit to local CR data [41]. So, in the remainder of this study we use the 3D IC maps in two different ways: first, we evaluate the impact of the new 3D IC maps on the results obtained in Ref. [32]. And second, we replace the Standard 2D IC model by F98-SA0 in our baseline when performing empirical tests of the bulge-correlated gamma rays (described following sections).

Table 3 shows the log-likelihood values obtained by adding to the fits the new extended templates (bulge models and preferred Fermi bubbles map). There are seven row sets in which for each different set an alternative 3D IC map is included in the baseline instead of the Standard 2D IC map: the first row set considers the Standard 2D IC model, the second F98-SA0, the third F98-SA50, and so forth, until the last row set that uses R12-SA100. Regardless of the IC map assumed, we find that a bulge model is still needed by the data. As can be seen, this holds even after inclusion of the inpainted structured FB map. Figure 4 shows a comparison of the TS values obtained for the boxy bulge + NB and NFW2 models respectively.

In each row set shown in Table 3, the methodology is to keep the component with highest log-likelihood value and move to the next iteration. In all the cases considered in this study we find that the stellar models have larger TS values than the NFW2 map. Remarkably, when the boxy bulge + NB is included in the base model, there is no statistically significant need for an NFW2 template (the significance of the NFW2 map is always $\lesssim 4\sigma$). The stellar templates evidently absorb much of the residual photons and the data no longer need a map describing a dark matter distribution. Finally, even when the NFW2 template is added to the base model, we find very strong evidence for the nuclear bulge + boxy bulge template. We therefore conclude that the data most prefer the nuclear bulge + boxy bulge. The low statistical significances obtained for the NFW2 model when different IC maps were assumed can be seen in Fig. 5. We note that this is consistent with our previous work [32] where a smaller RoI and only the standard 2D IC maps were used.

3.1 Rotation

Given the statistical preference for the stellar bulge templates under consideration, we continue to investigate whether the likelihoods can be improved by allowing more flexibility to the templates. In particular, we consider rotations and translations of the boxy bulge map. We even allow changes that are inconsistent with stellar data, which serve to test the interpretation of the gamma-ray correlations. Although the nuclear bulge is also robustly detected, we only perform these tests on the boxy bulge as this template is crucial to explain the GCE at higher latitudes. In order to test the rotational symmetry of the boxy bulge morphology we rotate the template about a line defined by $l=0^{\circ}$ (see inset image of Fig. 6) in increments of $9^{\circ}$ and compute the maximum log-likelihood at each turn. The results from this test are shown in Fig. 6 and indicate that the template in its original orientation is preferred. Note that the boxy bulge is East-West asymmetric (this is elegantly illustrated in Fig. 9 of Ref. [52]), and thus the $\Delta\log(\mathcal{L})$ profile looks different for negative and positive rotations. This indicates that the data has constraining power to statistically distinguish the East-West asymmetries present in the boxy bulge template. The maximum likelihood rotation was found to be $\sim 9^{\circ}\pm 2^{\circ}$ with significance $3.8\sigma$ ($\Delta\log(\mathcal{L})=7.2$ for one extra parameter). Therefore the data is consistent with the original orientation of the Boxy bulge image.

It is instructive to compare the maximum likelihoods of the original Boxy bulge image for $\pm 90^{\circ}$ rotations. In the case of a clockwise rotation of $\theta_{0}=90^{\circ}$, we find that the original orientation is a better model than the rotated one with a confidence of $15.8\sigma$ ($\Delta\log(\mathcal{L})=124.6$ for one additional parameter). For the counterclockwise rotation of the same angle ( $\theta_{0}=-90^{\circ}$), we find that the original orientation is preferred with a statistical significance of $7.5\sigma$ ($\Delta\log(\mathcal{L})=27.9$ for one extra parameter). Additionally, the profile displayed in Fig. 6 shows that a rotation around $180^{\circ}$ is disfavored at $13.0\sigma$ ($\Delta\log(\mathcal{L})=84.5$ for one extra parameter). This demonstrates that the data can also statistically differentiate the North-South bulge asymmetries.

3.2 Translation

Next, we translate the boxy bulge template along the longitudinal and latitudinal directions in increments of $0.2^{\circ}$ and compute the log-likelihood value at each translated position respectively. The results from these tests are shown in Figs. 7 and 8.

In the longitudinal direction case (Fig. 7), we find evidence for a translation of $l\sim-1.7^{\circ}$, with a broad $\Delta\log(\mathcal{L})$ curve encompassing the range $l\sim[0^{\circ},-3.5^{\circ}]$, while in the latitudinal case (Fig. 8) we find a translation at $b\sim 0.7^{\circ}$ with a $\Delta\log(\mathcal{L})$ curve that is comparatively smaller ($b\sim[0^{\circ},1.7^{\circ}]$). The boxy bulge images translated by $l\sim-1.7^{\circ}$ and $b\sim 0.7^{\circ}$ are preferred with $12.5\sigma$ ($\Delta\log(\mathcal{L})=77.7$) and $9.2\sigma$ ($\Delta\log(\mathcal{L})=42.0$ for one additional parameter) respectively. The preference for a shifted central position arises due to the presence of residuals on the Galactic west and north sides of the Galactic Center that are able to be absorbed by the boxy bulge template when it is translated along the latitude and longitude directions. Therefore, the translations appear real in the sense it is warranted by the gamma-ray data. However, in the construction of the 2FIG catalog [45], a group of several new point sources separated by less than $0.6^{\circ}$ were found around the region $l\sim 0^{\circ}$ and $b\sim[0^{\circ},2^{\circ}]$. Due to those new point sources not being easily individually resolved, the procedure of the authors of Ref.  [45] excluded that group of new point source candidates from the 2FIG catalog. It is likely that this exclusion drives the bulge translation shift and that with their inclusion in the form of, for example, an empirical patch of positive residual emission at $l\sim 0^{\circ}$ and $b\sim[0^{\circ},2^{\circ}]$, the boxy bulge latitudinal shift could be ameliorated. This is an interesting possibility that should be pursued in future work since it requires a dedicated focus along the Galactic plane, rather than the larger-scale Bulge emission of this study.

With regards to the longitudinal shift, recent gamma-ray analyses of the inner Galaxy [42, 43, 65] have also found excesses of gamma radiation off the Galactic Center at negative longitudes. Specifically, Refs. [42, 43] have both found that the base of the FB is a few degrees off the Galactic Center toward negative longitudes (see also central panel of Fig. 1). In addition, Ref. [65] has pointed out that the GCE is potentially offset by $-4^{\circ}$. However, when the authors run their pipeline with the Galactic plane masked, the offset becomes $-1^{\circ}$. The fact that this offset is sensitive to the applied masks indicates that the excess emission is mainly due to excess gamma-ray emission at very low latitudes ($b\sim 0^{\circ}$) but that there is also some evidence for a $\sim\mathcal{O}(1^{\circ})$ shift of the gamma-ray bulge from high latitude data.

In this context, we note that the distribution of molecular clouds and young stars in the so-called “central molecular zone” (CMZ) is highly asymmetric [66, 67, 68, 69]. Close to 2/3 of the molecular gas is at positive longitude—Galactic East near Sgr A⋆—while about 2/3 of the young stars are at negative longitudes—Galactic West near Sgr A⋆ 666See for example the bright 24 $\mu$m emission shown in Fig.1 of Ref. [69] Represented by blue color in that figure, this channel shows many bright point sources at $l\sim[0^{\circ},-2^{\circ}]$ which correspond to a population of massive young stellar objects and evolved high mass stars.. The asymmetric nature of the disposition of the molecular gas and young stars in the CMZ implies that the ISRF density distribution and the locations of CRs accelerators like supernovae and young massive stars, should be modeled with templates capable of accounting for these asymmetries. However, the triaxial function used in the derivation of the boxy bulge template adopted in this work [35] is not able to account for deviations of the triaxial symmetry in the stellar bulge data and, by construction, traces mostly the oldest bulge stars [35]. Nonetheless, Ref. [50] performed a translation analysis in their fits to red giant clump stars and found that the data always preferred an offset toward negative longitudes (see Table 1 in Ref. [50]). While there are still quantitative differences, with the best-fit offsets of Ref. [50] being smaller than those we find in the gamma rays, it is very likely that these are mostly due to gamma-ray correlations with star formation sites along the disk or in the far Galactic East region of the nuclear bulge. As a consequence, high level emission from these sites might be biasing our fits to larger shift values than those found by Ref. [50]. The asymmetry of the gamma-ray bulge is a very interesting topic that should be more carefully investigated in future work, using bulge templates that go beyond the functional forms adopted in this study 777 We observe that the 3D IC templates [41] used in this work do not yet account for the stellar populations of the nuclear bulge and assume triaxial geometry (F98 or R12) for the bulge/bar structure. .

4 Gamma-ray properties of the bulge

4.1 Gamma-ray bulge morphology

Globular clusters are old and extremely dense stellar agglomerations that might provide important clues about the formation mechanisms of MSPs in the Galactic bulge [9]. There is strong observational evidence that some of the most massive globular clusters in the Milky Way contain large numbers of MSPs [70]. However, given the relatively low escape velocities in globular clusters ($\bar{v}_{\rm escape}\lesssim 60$ km/s) [71], it is not straightforward to reconcile the presence of MSPs in these objects with empirical estimates of typical neutron star (NS) kick velocities ($\bar{v}_{\rm kick}\sim 90$ km/s) [72].

In order to solve this “NS retention problem” [73] in globular clusters, two scenarios have been proposed: First is the primordial formation scenario [74], which posits that NSs form in massive binaries making them more easy to retain. This is because the kick-momentum transferred to the NS is shared with a massive companion, which leads to a reduction in the velocity of the post-supernova compact object. The number of MSPs is expected to scale linearly with the total stellar mass in this scenario. Second is the dynamical formation scenario, which rests on the observed positive correlations between the stellar encounter rate $\Gamma_{c}$ and the number of low-mass X-ray binaries (thought to be the progenitors of MSPs) in globular clusters [75]. Such observations indicate that after NSs are born, they can be captured by massive stars to form binaries which in turn can evolve to become MSPs. Since $\Gamma_{c}\propto\rho^{2}$, where $\rho$ is the stellar density, it follows that in this picture the number of MSPs scales with the square of the stellar density.

In spite of the very different stellar densities and dynamical histories of globular clusters and the Galactic bulge, similar MSPs formation mechanisms might be operating in the Galactic Center [9, 13, 76]. The spatial distribution of low-mass X-ray binaries in M31 provides evidence for the view that both dynamical and primordial components exist [77], with the dynamical component dominating in the center. The close match between M31 X-ray binaries and the GCE morphology was previously used by Ref. [13] to argue for unresolved MSPs being the dominant source of the GCE.

Here, we perform a morphological analysis of the GCE data to find out whether the spatial distribution of photons is better explained by the primordial or dynamical formation scenario. To this end, we modify the boxy bulge density function in Eq. (2.1) to include an additional slope parameter $s$ in the form $\rho^{s}_{\rm bar}(R,\phi,z)$. By varying $s$ in the range $[0.5,2.5]$ in steps of 0.1 and integrating the resulting density function along the line-of-sight, we create additional stellar maps which we individually pass through our pipeline routine. We note that in doing so the nuclear bulge map is left unmodified as we are focusing on the larger scale bulge component.

The results of the morphological analysis are presented in Fig. 9. Panels (a) and (b) show the longitudinal and latitudinal dependencies of the boxy bulge + NB models for various slope values: the primordial formation ($s=1.0$), dynamical formation ($s=2.0$) and the best-fit slope ($s=1.4$). There are appreciable differences between the modified boxy bulge + NB profiles and the NFW2, which helps to explain why the maximum-likelihood analysis assigns largely different log-likelihood values to bulge models compared to spherically symmetric ones (e.g., Table 3). Interestingly, we find that the density slope that best fits the data is $s=1.38^{+0.06}_{-0.05}$, which is preferred over the primordial scenario with a statistical significance of $7.6\sigma$ ($\Delta\log(\mathcal{L})=28.8$ for one additional parameter), see panel (c) of the same figure. Therefore, we find significant evidence for an admixture formation scenario in which a fraction of the MSPs are formed through the primordial channel and the remainder are formed dynamically. Figure 9-(c) also shows that the fit worsens when stellar density models are assumed with slopes $s<1.0$ or $s>2.0$. For instance, a slope $s=0.9$ is statistically disfavored with $6.8\sigma$ ($\Delta\log(\mathcal{L})=23.1$ for one additional parameter) significance. Likewise, slopes of $s=2.1$ and $s=2.2$ are disfavored with $5.0\sigma$ ($\Delta\log(\mathcal{L})=12.3$) and $6.6\sigma$ ($\Delta\log(\mathcal{L})=21.5$) significance, respectively. This is an important consistency check that gives further support to the MSPs hypothesis for the GCE.

It has been argued that the putative population of Galactic bulge MSPs could have been the result of depositions from tidally disrupted globular clusters [21, 22, 78]. However, a fundamental prediction of this model seems to be a spherically symmetric distribution of Galactic bulge MSPs [78], which we have demonstrated to be highly disfavored by the data. There is a growing consensus that the origin of the Boxy/Peanut bulge morphology is the result of the dynamical evolution of a disk stars via buckling instabilities [79, 52]. In this picture, the Galactic bulge is predominantly composed of disk stars now on bar orbits. It is possible that by changing the initial conditions of their $N$-body simulations, the disrupted globular cluster scenario could reproduce the Boxy/Peanut bulge morphology that the gamma-ray data requires. Nevertheless, such simulations would still need to explain the nearly isotropic distribution of the surviving population of Milky Way globular clusters. We note that accounting for these two seemingly irreconcilable observations in the disrupted globular cluster scenario [21, 22, 78] appears to be a difficult task, thus disfavoring current versions of this model for the GCE.

4.2 Gamma-ray bulge spectrum

The gamma-ray luminosity888In the computation of the luminosities we assumed 8.25 kpc for the distance to the Galactic Center. of the boxy bulge model ($s=1$) was found to be $(2.2\pm 0.4)\times 10^{37}$ erg s*-1* for $E\geq 100$ MeV, as estimated for the whole $40^{\circ}\times 40^{\circ}$ region of interest. A power-law with an exponential cut-off (PLE) model was preferred to a simple power law (PL) with a confidence of $2.6\sigma$. As the amplitude is restricted to be non-negative, a $\chi^{2}/2$ distribution rather than the $\chi^{2}$ distribution is needed. We computed the best-fit spectral parameters with a $\chi^{2}$ fitting method to the inferred fluxes at each energy bin. We found $E_{\rm cut}=15\pm 14$ GeV and $\Gamma=2.1\pm 0.2$, which as can be seen from Fig. 5 of Ref. [11], are consistent with the resolved MSPs and also the GCE obtained using an NFW2 template. In Fig. 10, we show the spectra of the boxy bulge component in the inner $40^{\circ}\times 40^{\circ}$ around the Galactic Center. The yellow (green) error bars show the statistical (systematic) uncertainties. We have estimated the systematic uncertainties by computing the dispersion of the best-fit bin fluxes obtained when all the alternative 3D and standard-2D IC models are included in the maximum-likelihood analyses respectively. Our approach is consistent with earlier analyses by the Fermi-LAT collaboration [17] that argued for a large impact in the GCE spectrum due to IC modeling. After rescaling the GCE data from Ref. [15] to our RoI and comparing to our boxy bulge spectrum (red data points), we find very good agreement with their results. We note that in that work the GCE was obtained by masking the inner $|b|<2^{\circ}$ region of the Galactic plane and assuming an NFW2 spatial model.

The quality of the fit to the inner $40^{\circ}\times 40^{\circ}$ of the GC is illustrated in Fig. 11. The different columns show the data, best-fitting model and fractional residuals $(Data-Model)/Model$, respectively. The rows show different energy ranges, as labeled above the panels. All maps were smoothed using a gaussian filter of $0.5^{\circ}$ radius for display purposes. The first two columns demonstrate that there is in general a good agreement between the model and the data. As can be seen from the third column, our model underpredicts the data at the $\lesssim 20\%$ level for energy bins $[0.6,1.1]$, $[1.1,2.8]$ and $[2.8,11.8]$ GeV. The level of agreement of the model with the data is downgraded to the $\sim 30\%$ level for energies $[11.8,158.1]$ GeV. It is possible that this enhanced residual emission is due to an imperfect FB template for this particular energy range. We observe that the authors of Ref. [43] derived their FB template using gamma-ray data with energies between $[1,10]$ GeV. Application of their technique may be capable of better accommodating the small-scale features of the data in the $[11.8,158.1]$ GeV energy range. Similarly to what is shown in the second column of Fig. 11, we also show in Fig. 12 a comparison of the data and the best fit model components in a way that the relative importance of each of the components included in the fit can immediately be seen. Both these figures show that our model provides a good fit to the data.

A fit to the nuclear bulge template preferred a two component spectral model over a single one. Initially, we attempted to fit a PL to only the Fermi-LAT data taken from the nuclear bulge. We found a best-fit slope of $\alpha=2.40\pm 0.03$ and normalization $N_{0}=(9.8\pm 2.4)\times 10^{-11}$ MeV*-1* cm*-2* s*-1* with $\mbox{p-value}=0.045$ ($\chi^{2}=22.7$ for $15-2$ degrees of freedom). Although this p-value passes the usual acceptance condition $\mbox{p-value}>0.001$, when the resulting PL formula is extrapolated to higher energies, it clearly overshoots the H.E.S.S. measurements from the same region (left-hand side panel of Fig. 13). We also tested a PLE spectrum and obtained $\mbox{p-value}=0.0003$ ($\chi^{2}_{\rm PLE}=36.7$ for $15-3$ degrees of freedom ), which is less than the usual acceptance value.

Next, we included the H.E.S.S. data points [80] to our $\chi^{2}$ test so as to make the model consistent with higher energy data from the same region. We therefore opted to add a PLE spectral model to account for data at around $\sim 1$ GeV and a PL model for $E\geq 11$ GeV. This composite model consisting of a PLE+PL was preferred to a single PL with $3.4\sigma$ statistical significance ($\Delta\chi^{2}=\chi^{2}_{\rm PL}-\chi^{2}_{\rm PLE}$=17.7 for four extra parameters). Best fitting spectral parameters corresponding to the PLE component were found to be $\Gamma=1.6\pm 0.2$, $E_{\rm cut}=5.3\pm 1.4$ GeV and $L=(3.9\pm 0.5)\times 10^{36}$ erg s*-1* calculated in the energy range $E=[100,500000]$ MeV. As for the PL component we found $N_{0}=(5.6\pm 1.5)\times 10^{-11}$ MeV*-1* cm*-2* s*-1*, $\alpha=2.4\pm 0.1$ for $E\geq 11$ GeV 999This value is motivated by the known typical $E_{\rm cut}$ of the resolved MSPs. . While the H.E.S.S ROI ($|l|<0.8^{\circ}$ and $|b|<0.3^{\circ}$) is indeed close to that of the nuclear bulge template, the H.E.S.S ROI could be slightly smaller. A dedicated observation of exactly the nuclear bulge region with H.E.S.S or the forthcoming Cherenkov Telescope Array (CTA) could make the consistency across the full gamma-ray energy range much better. Similarly to the boxy bulge component, the best-fit PLE parameters are consistent with spectrum of resolved MSPs. These results are summarized in the right-hand side panel of Fig. 13.

A simple astrophysical model that can match the characteristics of the preferred PLE+PL model is an unresolved population of MSPs in the nuclear bulge in addition to a population of energetic electrons from star formation activity and the same MSPs. Prompt gamma-rays from MSPs could account for the PLE spectrum seen at Fermi-LAT energies and a population of electrons injected from the same region could explain both the high-energy tail observed at energies $>10$ GeV and the H.E.S.S. Galactic ridge spectrum (see for example Ref. [81] for a variation of this picture). Hadronic gamma-rays are also expected to contribute to the Galactic ridge at some level. A combination of hadronic and leptonic scenarios for this region have been explored in Ref. [82]. However, our new results for the nuclear bulge region motivates a much more detailed modeling that self-consistently accounts for hadronic gamma-rays from the CMZ region as well as a new population of electrons from MSPs and star formation activity. This particular topic will be addressed in future work.

It is useful to compare the mass-to-light ratios of our hypothesized boxy bulge and nuclear bulge MSP populations with others residing in different systems101010We caution that this comparison assumes that the ratio of MSP mass to stellar mass is constant across the different regions considered, which can only be approximately correct.. The stellar mass of the Galactic bulge is estimated to be around $1.5\times 10^{10}$ solar masses ($M_{\odot}$) and for the nuclear bulge $1.4\times 10^{9}$ $M_{\odot}$. Bringing together the bulge mass estimates with our best-fit luminosities of the boxy bulge and nuclear bulge components, the luminosity per stellar mass for the boxy bulge is found to be $(1.5\pm 0.3)\times 10^{27}$ erg s*-1* $M_{\odot}^{-1}$ while for the nuclear bulge it is found to be $(2.8\pm 0.4)\times 10^{27}$ erg s*-1* $M_{\odot}^{-1}$. Thus, we find a luminosity per stellar mass for the nuclear bulge that is $\sim 2$ times higher than for the boxy bulge. This is qualitatively consistent with a picture in which the MSPs of the boxy bulge are systematically older than those in the nuclear bulge, and hence less luminous per stellar mass. The higher luminosity per stellar mass of the nuclear bulge might also be due to additional energy being injected from star formation activity on top of an old MSP population.

5 Discussion and conclusions

In this work we have performed a detailed reevaluation of the systematic uncertainties affecting the properties of the Galactic Center Excess (GCE). We have utilized hydrodynamical interstellar gas and dust maps [32] in order to account for uncertainties associated with gas-correlated gamma-ray emission. In addition, we have employed, for the first time, the new 3D inverse-Compton (IC) maps111111Some of the new 3D IC maps dramatically improve the total log-likelihood of baseline model. This demonstrates the great constraining potential of the Fermi-LAT data from the inner Galaxy for constructing improved CRs and ISRF models. that abandon the assumption of Galactocentric cylindrical symmetry adopted in most previous analyses of the same region. Improved low-latitude Fermi bubbles maps have also been included. By considering this wide range of maps and uncertainties, we have confirmed recent studies [32, 33] that reported a positive correlation between the stellar distribution of the Galactic bulge and the GCE.

A recent study by the Fermi collaboration [17] has raised some concerns with regard to whether the GCE is the result of incomplete modeling of the IC component in the inner Galaxy. However, our results reveal that such concerns are unwarranted. Even though the new 3D IC maps [41] include the bulge/bar for the CR source and ISRF densities, an additional stellar bulge component in the form of the nuclear bulge + boxy bulge is still required by the Fermi-LAT data. This is consistent with the hypothesis of an unresolved population of MSPs in the Galactic Center region. Since CRs interacting with interstellar gas and the ISRF produce IC morphologies at $\sim 1$ GeV that are different to both those of the CR source distribution and the target material, statistical analyses of the GCE would still be required to add template model maps to account for the prompt gamma-ray emission from MSPs that closely match their source distribution.

Although the X-shape [38] structure of the Galactic bulge is estimated to contain some $\sim 40$% of the mass of the bulge stars, these stars spend only a small fraction of their orbital time within the X arms [52, 79]. Therefore, at any given time, about 5–10% of bulge stars are observed to be within the X arms [50, 79, 84]. Under the assumption that gamma-ray sources are kinematically similar to stars, we would therefore expect that the X arms would contribute $\sim 10\%$ of the Galactic bulge gamma-ray luminosity. In the maximum-likelihood analyses presented in Table 3, we found that there is statistical evidence for the X-shape bulge in gamma-rays, but that this was weaker than that of the whole boxy bulge structure. In our step-wise statistical procedure, only new sources that increased the log-likelihood the most were allowed to be included in the subsequent baseline model.

In Ref. [32], using a RoI of size $15^{\circ}\times 15^{\circ}$ and only the standard 2D IC maps, some of us showed that the GCE was better fit by either the nuclear bulge + X-bulge or nuclear bulge + boxy bulge. In that work both stellar templates were found to have similar statistical significances. Now with our larger RoI and improved IC maps used in this analysis, we are finding that the significance of the nuclear bulge + X-bulge is smaller than that for the nuclear bulge + boxy bulge template. However, we note that it is not at present possible to nest the boxy bulge template with the X-bulge template in the same maximum-likelihood fit, which is necessary for a robust comparison. This is because the boxy bulge model already contains a substantial fraction—if not all—of the X-bulge stars. Our new results strongly encourage a careful revaluation of the luminosities associated to the X-bulge and boxy bulge stellar components using stellar maps that are physically self-consistent.

In order to test whether the stellar models can be improved by allowing more flexibility to the templates, we rotated the boxy bulge template around the central position in increments of $9^{\circ}$ and computed the log-likelihood at each orientation. The results from this test indicate that the boxy bulge in its original orientation is consistent and that the data is sufficiently constraining to differentiate the East-West and North-South asymmetries present in the stellar bulge model (see Fig. 6).

Similarly, with the aim of testing the GCE for any possible translational shifts with respect to the location of the stellar bulge, we have moved the central position of the boxy bulge along latitude and longitude in steps of $0.4^{\circ}$. Statistically significant shifts were found toward positive latitudes at $b=0.72^{\circ}\pm 0.07^{\circ}$ and negative longitudes at $l=-1.69^{\circ}\pm 0.12^{\circ}$ (see Figs. 7 and 8). The former is relatively small and is likely due to a positive excess resulting from the suppression of a cluster of point sources in the 2FIG [45] catalog around the region $l=0^{\circ}$ and $b\sim[0^{\circ},2^{\circ}]$. However, the later prefers a shifted center at larger angular distance and increased statistical confidence. This result is in line with a similar such analysis performed with the red clump giant stellar data as seen in Table 1 of Ref. [50], where a Galactic westward translation of around 0.3∘ was found. Although the preferred shifted positions in the two studies differ by $\Delta l\sim 1.4^{\circ}$, it is worth noting that in the region around $b\sim 0^{\circ}$ and $l\sim[-0^{\circ},-2^{\circ}]$ there is a population of young massive stellar objects that might be emitting significant gamma-ray flux which, by construction, the analysis of Ref. [50] does not account for. This is because the red clump giants are systematically older stellar objects. In our view the results of these empirical tests emphasize the validity of the proposal for a stellar origin for the GCE, within the systematic limitations of the current maps. In summary, we find that the gamma-ray data highly prefers a Galactic bulge model that is consistently positioned and aligned within one or two degrees of the known Galactic bulge [35]. All other translational and rotational variations result in significantly poorer fits to the data.

Perhaps the most striking result to emerge from the GCE data is that any model predicting a spherically symmetric morphology is robustly ruled out (see Fig. 5). This, of course, includes the annihilation of weakly-interacting dark matter models but also current versions of the disrupted globular cluster scenario [21, 22, 78]. Any successful disrupted globular cluster model must explain both the spherical distribution of surviving globular clusters in the Galaxy as well as the boxy/peanut distribution of MSPs in the bulge. Taken together, these two observations seem to disfavor this scenario (at least in its current form [21, 22, 78]) as the main formation mechanism for MSPs in the Galactic bulge.

We have also tested the primordial and dynamical formation scenarios for the Galactic bulge MSPs. In our detailed morphological analysis of the GCE, we have found a significant statistical preference ($7.6\sigma$) for an admixture formation scenario in which a fraction of the MSPs are formed through the primordial channel ($52\%\pm 23\%$) and the remaining one through the dynamical channel. Due to systematic uncertainties in the Galactic diffuse emission model, our maximum-likelihood method assigns almost identical log-likelihood values to a pure primordial formation as well as to a pure dynamical formation scenario (see Fig. 9). However, models whose spatial morphology start to depart from those predicted by these two scenarios are statistically disfavored by the data. This is an important check that lends further support to the MSP theory for the GCE because we expect both components to exist in the MSP population [77].

Reference [85] used the results of a correlation between globular cluster gamma-ray luminosity and the stellar encounter rate [75] to infer the energetics and morphology of the dynamical formation scenario for MSPs in the Galactic bulge. That study predicts a contribution to the bulge luminosity relative to the primordial gamma-ray luminosity at the percent level ($\lesssim 5\%$). In contrast, using our maximum-likelihood procedure we find that the best-fit contribution for this channel is $\sim 48\%\pm 21\%$ (see Fig. 15). This was obtained by simultaneously including in the fit a $\rho^{1}_{\rm bar}(R,\phi,z)$ and $\rho^{2}_{\rm bar}(R,\phi,z)$ map and then estimating their corresponding best-fit gamma-ray luminosities121212The fractional luminosity error was obtained by adding in quadrature the inferred energy fluxes at each energy bin. We also note that in this case there are 15 new degrees of freedom (one for each band of the $\rho_{\rm bar}^{2}$ template) while in estimating $s$ there was only one new degree of freedom.. Figure 15 shows how the flux fractions computed this way adds up to approximately our best-fit stellar density slope $s\sim 1.4$ obtained with the slope $s$ scan of Fig. 9. Although our detailed morphological analysis does not seem to corroborate the luminosity fraction predicted by Ref. [85], we caution that a more direct comparison with that study is not at present possible given that the authors have assumed a spherically symmetric model (see Sec. 3 in [85]) for the distribution of the bulge stars. Such stellar distributions are thought to be not realistic [52, 79] and are robustly ruled out by the gamma-ray data [32, 33].

In agreement with most previous analyses of the GCE [9, 10, 11, 12, 13, 14, 15, 16, 17], we have found that the spectrum of the boxy bulge is consistent with that expected from an unresolved population of MSPs in the Galactic bulge (see Fig. 10). In particular, we confirm earlier results [31] that claimed the detection for a high energy component above $\sim 10$ GeV associated to the GCE. As can be seen in Fig. 14, this high-energy tail has a different spectrum from that of the Fermi Bubbles spectrum. Such a component could be the result of IC emission of electrons and positrons injected by the same Galactic bulge MSP population [31].

Furthermore, we identified a corresponding high-energy tail also in the nuclear bulge template (Fig. 13). As far as we are aware, this is the first time that a high-energy tail has been detected in this smaller angular region. A fit to the Fermi and H.E.S.S. observations from this site preferred a two component spectral model (although with a low confidence of $3.5\sigma$) comprising a power-law with an exponential cutoff model plus a simple power-law for energies $\gtrsim 10$ GeV. An important implication of this result is that energetic electrons emitted by MSPs, and possibly star forming activity, could potentially contribute a significant fraction of the TeV gamma-rays detected by H.E.S.S.

Our results show that the Galactic Center Excess in gamma rays is correlated with the bulge of the Milky Way. The most likely explanation for this correlation is that the unresolved bulge millisecond pulsars are responsible for the Galactic Center Excess.

Acknowledgments

We are grateful to Kevork Abazajian, Shin’ichiro Ando, Anthony M. Brown, Henry Freudenreich, Ryan Keeley, Dmitry Malyshev, Harrison Ploeg, Troy Porter, Nicholas Rodd, Tracy Slatyer, and Deheng Song for useful comments and discussions during the project. This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. OM was partially supported by JSPS KAKENHI Grant Numbers JP17H04836, JP18H04340 and JP18H04578. SH thanks support by the U.S. Department of Energy under award number de-sc0018327. MK was supported by NSF Grant No. PHY1620638. DMN was supported by the Allan C. and Dorothy H. Davis Fellowship. OM acknowledges the high performance computer center at Kavli IPMU for providing computational resources and specially to Tatsuhiro Yano for technical support that have contributed to the results reported in this paper.