Strong constraints on thermal relic dark matter from Fermi-LAT observations of the Galactic Center

The extended excess towards the Galactic Center (GC) in gamma rays inferred from Fermi-LAT observations has been interpreted as being due to dark matter (DM) annihilation. Here, we perform new likelihood analyses of the GC and show that when including templates for the stellar galactic and nuclear bulges, the GC shows no significant detection of a DM annihilation template, even after generous variations in the Galactic diffuse emission models and a wide range of DM halo profiles. We include Galactic diffuse emission models with combinations of 3D inverse Compton maps, variations of interstellar gas maps, and a central source of electrons. For the DM profile, we include both spherical and ellipsoidal DM morphologies and a range of radial profiles from steep cusps to kiloparsec-sized cores, motivated in part by hydrodynamical simulations. Our derived upper limits on the dark matter annihilation flux place strong constraints on DM properties. In the case of the pure $b$-quark annihilation channel, our limits on the annihilation cross section are more stringent than those from the Milky Way dwarfs up to DM masses of $\sim$TeV, and rule out the thermal relic cross section up to $\sim$300 GeV. Better understanding of the DM profile, as well as the Fermi-LAT data at its highest energies, would further improve the sensitivity to DM properties.


I. INTRODUCTION
The particle nature of DM remains one of the most important unresolved questions in astrophysics, cosmology and particle physics. Hierarchical structure formation with cold collisionless or self-interacting dark matter predicts that the Milky Way (MW) GC would contain a large concentration of DM [1,2], providing an avenue for stringent tests of DM annihilation [3][4][5]. After the launch of the Fermi Gamma-Ray Space Telescope, an extended source of gamma ray emission was quickly identified toward the GC and shown to be consistent with the annihilation of thermal weak-interaction-scale DM producing gamma rays [6]. This GC excess (GCE) has since been detected by many followup analyses, which also indicated its potential association with unresolved point sources or new diffuse emission processes [7][8][9][10][11].
A major challenge for establishing DM signals in our MW's GC is the abundant astrophysical activity in the GC region. For example, part of the GCE signal could be explained by gamma-ray emission induced by cosmic rays injected by on-going star formation activity in the GC region [12], cosmic-ray bremsstrahlung off of molecular gas [8], or inverse-Compton emission from leptonic cosmic rays [13,14]. However, these studies are not able to completely explain the data, and still leave the need for a spherical GCE.
Recently, Refs. [15,16] showed that the GCE overwhelmingly prefers the spatially asymmetric morphology of the Galactic stellar bulge-a triaxial bar-like structure extending a few kpc in the GC [17]-over the spher-ically symmetric morphology assumed by a DM origin. The bulge has a radially-varying asymmetry that was not captured in earlier elliptical shape tests conducted on the GCE [8,10]. A detailed study of the robustness for the detection of the bulge, including systematic uncertainties arising from background emissions and other gamma-ray sources, was presented in Ref. [18]. Since the bulge contains a broad mix of star-forming and old stellar populations, this motivates a population of astrophysical gamma-ray emitters such as young pulsars and millisecond pulsars (MSP) as the source of the excess gamma rays. Most significantly, the inclusion of the asymmetric bulge model completely eliminates the need for a spherically symmetric DM component of the GCE [15,16]. This provides an opportunity to substantially improve the sensitivity to test DM properties.
In this article, we present stringent DM limits incorporating recent developments in modeling the bulge and other astrophysical gamma-ray sources in the GC region. To ensure that our limits are robust, we use results from galaxy formation simulations to inform our DM templates, which provides a significant point of departure from previous work. Furthermore, we explore generous variations in models of the gamma-ray emission from cosmic-ray interactions. Even with the considerably larger freedom for the astrophysical emission and DM profiles, our results show that the Fermi Large Area Telescope (Fermi-LAT) observations of the GC provide very stringent constraints on DM annihilation. For twobody final states with hadronic components, we are able to rule out thermal DM up to approximately 300 GeV arXiv:2003.10416v1 [hep-ph] 23 Mar 2020 in mass, surpassing the reach from dwarf satellites of the MW for DM particles with masses less than a TeV.

II. DARK MATTER LIMITS
To calculate the limits on DM annihilation cross section, we must first generate a likelihood profile for the DM annihilation intensity for a given DM halo model. We consider four classes of MW DM profiles, described in the next Section. The likelihood for each value of the DM annihilation intensity is computed by varying the fluxes of all the background templates such that the loglikelihood is maximized. We use the Fermi UpperLimits tool 1 to perform this maximization and generate the likelihood profile for the DM annihilation intensity.
Our background model contains templates for the following: hadronic emission traced by HI and H2 gas maps divided in four cylindrical concentric rings and two total dust maps, 3D inverse Compton (IC) divided into four or six rings and a 2D IC map with a central source of electrons, an isotropic background, the 4FGL [19] point sources (PS), Fermi Bubbles, Loop I, the Sun and the Moon. Details of the templates, methods employed, likelihood profiles, resulting spectra, as well as our comprehensive checks and analyses of the systematic effects, are all presented in the appendices. Additional tests of the robustness of the preference for the bulge template are discussed in Ref. [18].
The likelihood profile is generated in 15 independent logarithmic-spaced energy bins between 0.667 and 158 GeV, and no broadband spectral shape is assumed for any of the templates. Following this methodology, we are able to marginalize over the uncertainties in the astrophysical backgrounds in a manner that is independent of the uncertainties in the particle physics models. An indicator of the success of our method is that we recover physically-consistent, continuous spectra for all the background templates (see Fig. 6 in appendix C).
With the likelihood profiles in hand, we use Bayes' theorem to calculate a posterior in the annihilation cross section and DM mass parameter space. The flux signal from DM annihilation scales as where dΦ/dE is the differential number flux, σv is the velocity averaged DM cross section times relative velocity, m χ is the DM mass, dN/dE is the gamma-ray energy spectrum, and the J-factor (J) is the integral through the line of sight over the region of interest of the DM density squared, J = dΩ ds ρ 2 (r(s, Ω)). We need to marginalize over this J-factor in order to calculate the posterior for the DM mass and annihilation cross section.
We assume that the dark matter is single-component when calculating the J-factor. If this is not the case, our constraints on σv should be recast as constraints on f 2 DM σv , where f DM is the fraction of cosmological dark matter density in the model being constrained. This is important for thermal relics because f DM scales inversely with the total annihilation cross section in the early Universe and hence the flux decreases for s-wave cross sections larger than the thermal relic cross section.
For comparison purposes, we also calculate the posterior distribution of the DM mass and annihilation cross section for the eight classical MW dwarf spheroidals, with well-determined J-factors. We use the likelihood profiles for the classical dwarfs from Ref. [20] and the uncertainties in the J-factors of the dwarfs are taken from Ref. [21], which are inferred from fits to the stellar kinematic data using generalized NFW profiles. Unlike the GC region, the J-factors for the classical dwarfs are well-constrained by stellar kinematic data because they are dark matter dominated and the region of interest (RoI) ∼ 0.5 • is wellmatched to their stellar half-light radii [22].

III. DARK MATTER PROFILES
We consider four classes of MW DM profiles: a generalized NFW (gNFW) profile, a cored profile that matches smoothly on to a NFW profile while conserving mass [25], and ellipsoidal versions of both. The gNFW density profile is where ρ is the local DM density, R is the solar radius, r s is the scale radius, and γ is the inner slope. Parameterizing the MW DM profile this way is useful since it allows us to use independent data sets to characterize the uncertainty in each parameter. We allow γ to vary between 0.5 and 1.5. Note that if the GCE is to be explained by a gNFW squared (gNFW 2 ) template, then we need γ ≈ 1.2 [8][9][10]. We adopt a log-normal prior on r s with a mean value of 26 kpc and a width of 0.14 dex, consistent with the ΛCDM concentration-mass relation [26] for halo mass of 10 12 M . We neglect the factor of two uncertainty in MW's halo mass since it is subdominant to the adopted spread. For the local density, we take ρ = 0.28 ± 0.08 GeV/cm 3 as a prior from Zhang et al. (2012) [27], who constrain the local DM density from the vertical motions of K-dwarfs close to the plane of the MW, independent of the DM density at other radii. We note that the local density constraint agrees very well in the ρ s -r s plane with the mass constraint at about 20 kpc obtained using Globular cluster proper motions from Gaia [28,29]. Previous hydrodynamical simulations of the MW with cold DM have typically found the profile to be adiabatically contracted [30]. Our profiles with inner slopes of γ > 1 capture this possibility. However, recent hydrodynamical simulations of MW-like galaxies also show the  ) on the average DM cross section times relative velocity σv for annihilation to bb, among a large number of GDE models and DM distributions considered. The GDE models allow for changes in the interstellar gas, dust, and IC distributions. For both the gNFW (left) and cored (right) DM profiles, we considered spherical and ellipsoidal shapes. For gNFW, the inner slope was also varied. See text and Fig. 2 for details. The dashed black line is the thermal cross section [23]. The H.E.S.S. [24] and stacked dwarfs limits [20] are shown for comparison and do not reflect the different GDE models and DM profiles. All the constraints shown assume that the DM is entirely made up of one kind of particle. If this assumption is relaxed, then the constraints on σv should be weakened by the square of the fraction of DM in the component being constrained. The data files and code necessary to reproduce this figure are available at https://github.com/oscar-macias/Fermi_GC_limits. presence of a core in the DM density profile with a size of roughly a kpc. The origin of this core is not fully understood. Using the Eris simulation [30], Ref. [31] argued that the core is formed in response to the bar, along the lines of ideas proposed earlier [32,33], and not due to feedback. They also noted the supporting fact that a roughly same-size core is present in another simulation identical to Eris but with a lower star formation threshold, which reduces feedback effects dramatically. Further indirect evidence supporting the view that the presence of the bulge, and not the direct impact of the feedback, is causing these kpc-sized cores comes from simulations with a fixed disk and bulge potential that lead to similar cores [34].
We use the cored "Read" profile [25] to investigate the effects of a cored dark matter density. It has a core radius r c that describes the removal of mass from the center to the outer parts due to core formation and the mass asymptotically tends to the NFW profile mass at large radii. The enclosed mass for the cored profile is described by, where we take M NFW (r) to be the NFW profile with γ = 1. We fix the core radius to be 1 kpc in keeping with the discussion of the simulations above and, in order to make a straight-forward one-to-one comparison, we assume the same prior distribution for r s (a mean of 26 kpc and a scatter of 0.14 dex). Note that this neglects the impact of adiabatic contraction, which would increase the inner core density. A better characterization of the the inner density profile of MW halos is likely to lead to stronger results than those presented here. We then use Monte Carlo sampling to calculate the prior uncertainty on the J-factor from the prior uncertainty on these parameters of the MW's DM profile.
The presence of the bulge and bar should also have an impact on the axis ratio of the DM template. The expectation is that the DM density profile is an ellipsoid with the short axis perpendicular to the stellar disk [35]. This flattening of the halo should be due, in part, to the formation of the stellar disk. Moreover, there is likely also a perturbative effect of the bar formation on the halo that induces further flattening [35]. The Eris simulation discussed previously finds a minor-to-major axes ratio of about 0.8 at 1 kpc and intermediate-to-major axes ratio of unity [36].
Given the arguments above, a flattened ellipsoid with a mild radial variation in the density is a reasonable description of the inner kpc of the MW halo. This is very different from the spherical gNFW γ = 1.2 profiles that were used by the bulk of the explorations of the GCE and considered to be representative of the expectations for cold DM. To test for the impact of non-spherical DM distribution, we use two different density ellipsoids with axis ratios of 0.7 (somewhat more flattened than the results in Ref. [36]): one in which the radial profile is the same as the gNFW profile with γ = 1.2 and the other in which the density profile is the same as the cored profile with r c = 1 kpc.   Contours of the posterior that contain 95% of the probability showing the impact of variations in the GDE models and DM spatial morphology. The left panel shows variations around a gNFW profile; the inner-slope γ ranges between [0.5,1.5] and the spatial distribution can be spherical or ellipsoidal (blue lines). The purple lines correspond to the systematic uncertainty arising from different GDE templates. Similarly, the right panel shows variations around the cored profile (spherical or ellipsoidal; blue) and the GDE models (purple). The data files and code are available at https://github.com/oscar-macias/ Fermi_GC_limits

IV. RESULTS AND DISCUSSION
We find that an emission template that traces stellar mass in the Galactic bulge is preferred in all (independent) energy bins over each of the DM templates considered in this analysis. In none of our maximum likelihood runs-that included a variety of alternative GDE models-was a DM template detected. This allows us to impose strong constraints on σv using the flux likelihood profile for each DM template and GDE combination as described in Section II.
Our constraints on σv are presented in Figs. 1 and 2. The curves correspond to the contours of the posterior in cross section and mass that enclose 95% of the probability. In Fig. 1 we show the maximum cross sections from the set of the 95% limits derived from the gNFW (left panel) and cored(right panel) profiles, while in Fig. 2 we show in more detail how the limits are affected by variations of the GDE model and DM morphology.
First, we explore the impact of the GDE model on the DM limits. The purple lines in both panels of Fig. 2 illustrate the systematic uncertainty that arises from different GDE model assumptions. We explore alternative dust, interstellar gas, 3D IC map composed of four or six independent rings, and a 2D IC map containing an additional central source population of e − (Model B in Ref. [37]). In particular, interstellar gas maps constructed using hydrodynamic simulations [15] or the standard interpolation method [38] and dust maps with E(B − V ) magnitude cuts of either 2 or 5 mag were considered. Details about these models are provided in the appendix and in Ref. [18]. Next, we explore the impacts of the DM profile. The blue lines in Fig. 2 correspond to the systematic uncertainty arising from different MW DM profiles. For the left panel, the blue lines represent 95% C.L. upper limits on the cross sections derived from the gNFW DM templates with a generous γ = [0.5, 1.5] range and an ellipsoidal-NFW case, while the same lines in the right panel are determined by the limits from the cored profile with a core size r c = 1 kpc for spherical and ellipsoidal shapes.
For the gNFW profile we find that γ = 1.2 is the value that attributes the most flux to the squared-NFW profile. This is why we used the γ = 1.2 profile as our baseline when varying the GDE models. Note that while γ = 1.2 is the value with the largest flux for that template, smaller values of γ have smaller J-factors. These two effects compete and in the end the γ = 1.0 value corresponds to the weaker limit for all masses. For all panels shown in Fig. 1 and 2, when DM profile variations are studied we assumed the benchmark GDE model described in the appendix.
To address potential issues of mismodeling and over-or under-subtraction affecting our limits, we have performed a series of injection tests. We found that in the vast majority of cases, our analysis successfully recovers the correct statistical coverage of constraints. However, we found a systematic bias in the last energy bin (34.49 − 158.11 GeV) of our analysis for the cored profile cases considered. Therefore, we removed the last energy bin from our upper limits for all the cored profile cases (see appendix). This explains why the limits shown in the right panels of Figs. 1 and 2 are weaker than those on the left, and become comparable to those from dwarfs for DM masses larger than a TeV.
The results for other channels including W + W − , ZZ, τ + τ − and HH are presented in the appedinx. The qualitative features with respect to variations in the diffuse models and the density profiles are the same as bb channel for these other channels. We do not consider annihilation to e + e − and µ + µ − since the dominant gamma-ray contribution in these cases will arise from the IC process, causing the spatial profile to change from the DM density-squared morphology [39].
A noteworthy aspect of our results is that despite allowing for extensive systematic uncertainty, they provide strong constraints on thermal relic models with DM particle masses smaller than about 300 GeV and they are comparable to the H.E.S.S. constraints for masses around a TeV [24]. It is interesting to compare our results with those of Refs. [4,5], which showed that the inner galaxy constraints could be better than those arising from MW dwarfs, even with a kpc-sized cored DM template. More data, better models for the point sources and the diffuse emission, and the inclusion of the bulge templates have all contributed to making our constraints stronger and more robust. We note that our limits are in reasonable agreement with expected Fermi-LAT limits for the inner Galaxy [40].
Three advances in the future will make our results even more powerful. The first is a deeper understanding of the central density profile of DM in the Milky Way and its correlation with the stellar bulge and disk, which could remove the uncertainty arising from the radial distribution and shape of the DM template. This could allow the properties of the MW dark matter profile in the inner kpc, such as the core size and ellipsoidal shape, to be constrained based on the bulge and disk. In addition, such a study will allow us to include the effect of adiabatic contraction, which has been neglected in our study for the cored profile and could increase the J-factor by up to a factor of 2 (estimated by varying r s by a factor of 2 to account adiabatic contraction). We note that if the core radius r c were larger by a factor of 2, the J-factor for our ROI would only decrease by about 30%. This is because the J-factor for our ROI with a cored profile is dominated by contributions from r > 1 kpc.
The second important advance would be a clear determination of the point source nature of the GCE. While this won't quantitatively change our constraints, it will provide corroborating evidence for the Bulge-GCE connection that our analysis clearly prefers. This may be possible through the non-Poissonian template fitting procedure [41][42][43][44][45][46] and wavelet techniques [47][48][49] to detect clustering of photons or radio detection of point sources responsible for the bulge emission [50], or detection of a significant number of millisecond pulsars (putative sources for the bulge gamma-ray emission) with radio telescopes [51][52][53][54].
The third is further improvements in 3D models of the gas and ISRF maps, which directly feed into the diffuse emission models and determine the residuals from fitting to Fermi-LAT data. For the cored profile, the upper limits for the six-ring 3D ISRF model are evidently more stringent than those for the other GDE models, which exhibit more subtle differences. Since we have chosen our limits to be the weakest among the GDE models, and not the best fitting, a study that includes more extensive set of background models may be able to improve upon our σv limits at high masses by a factor of few.

V. CONCLUSION
The detection in the Fermi-LAT data of a spatiallyconcentrated excess of gamma-ray emission in the MW potentially consistent with DM annihilation [6,8,10] has sparked great interest in the sources of high-energy emission in the GC. At the same time, the Fermi-LAT data has spurred steady progress in our understanding of the gamma-ray emission from our galaxy over the past decade. With detected sources that are consistent with the Fermi bubbles, 4FGL point sources, detailed IC emission maps, disk gas, and, most importantly, the emission from the stellar Galactic bulge and nuclear bulge, there is no significant excess in the GC that may be attributed to DM annihilation.
Our results strongly favor the idea that the excess emission in the GC at GeV energies is dominantly of astrophysical origin, consistent with other analyses that show the statistical distribution of photon counts is more consistent with point sources than spatially-uniform DM annihilation [41,43,46,47], studies of the asymmetry of the GCE [15,16,18], and detailed models of MSP stellar remnants that are consistent with the flux and spectrum in the region (e.g., [55,56]).
While gamma-ray emission from DM annihilation in the GC is still possible, the flux would have to be below that of the GCE, and with parameters consistent with the exclusion regions of Fig. 1. In arriving at this conclusion, we allowed for a variety of DM templates. These include ellipsoidal profiles with a kpc-sized core that we suggest, based on existing simulations of the MW, are closest to the true prediction for the density profile of cold dark matter. We explored in detail the robustness of our results to variations in the GDE models arising from new sources of relativistic e ± , new 3D IC templates, and changes to the standard gas maps. Our results provide stringent constraints on models of thermal relic dark matter with masses up to a few hundred GeV and prompt annihilation to Standard Model particles.

VI. ACKNOWLEDGMENTS
We are grateful to Troy Porter for his very helpful comments about GALPROP v56 and for making the new 3D ISRF data publicly available. We thank Shin'ichiro Ando, Francesca Calore, Roland Crocker, Douglas Finkbeiner, Chris Gordon We used 8 years (August 4, 2008−August 2, 2016) of P8R3 ULTRACLEANVETO data recorded by Fermi -LAT. We chose these particular time cuts because these are exactly the same used in the construction of the 4FGL catalog [19], thus making our point source modeling completely self-consistent. Note that if a bigger amount of data had been chosen, then a dedicated point source search would have been necessary for this work. Events with measured energies between 667 MeV and 158 GeV were divided into 15 logarithmic energy bins. In order to minimize contamination from the Earth atmosphere we only considered photons detected at zenith angles larger than 90 • . Moreover, we employed the recommended data quality filters (DATA − QUAL>0)&&(LAT − CONFIG==1) and restricted the analysis to a square region of 40 • × 40 • around the GC. Our study was carried out using the standard Fermitools v1.0.1 2 analysis framework, and instrument response functions (IRFs) P8R3 − ULTRACLEANVETO − V2. The gamma-ray background and foreground model used in this work is similar to that developed in Ref. [15]. However, this has been further improved by including new three-dimensional (3D) Inverse Compton (IC) maps [57] and a more robust low-latitude Fermi Bubbles (FB) template [18,58]. In particular, the 3D IC maps were modeled using the 3D Interstellar Radiation Field (ISRF) data available with the recent GALPROP v56 [57, 59], though conventional 2D IC maps were also tested in our analysis of the systematic uncertainties in the GDE model. Furthermore, the 3D IC maps have been divided in several rings (see Tab. I) and their corresponding normalization floated during the fits to account for the impact of CR density uncertainties. As for the Fermi bubbles component, we have included the map recently developed in Ref. [18]. In that study, the structured FB template of Ref.
[58] was further modified by an inpainting algorithm to help restore image processing artifacts due to point source masks used in its derivation.  We also included templates for the Sun and the Moon that match our photon event class and cuts (available in the Fermi fourth catalog of point sources 4FGL [19]), an isotropic component (iso − P8R3 − ULTRACLEANVETO − V2 − v1.txt) and an emission model map for Loop I [15,60]. The diffuse gammaray emission resulting from the interaction of energetic cosmic ray (CR) particles with the interstellar medium (ISM) was modeled as a linear combination of atomic and molecular hydrogen gas templates divided into four concentric rings (0 − 3.5, 3. The gamma-ray point sources present in our RoI were modeled using the 4FGL catalog [19]. There is a total of 487 point-like and extended sources inside our RoI.  Due to limitations in the maximum number of parameters that can be reliably fitted in a given run within Fermitools, we have employed the hybrid fitting procedure implemented in Ref. [18]. Specifically, we varied the normalization of each of the 120 brightest point sources in our RoI, while for the remaining 367 sources we constructed a point source population template whose normalization was allowed to vary at each energy bin. This is a good approximation given that the amount of data utilized in the present study is the same used in the construction of the 4FGL catalog. The point source population template was constructed by using the best-fit spectra in the 4FGL and convolving it with the The systematic uncertainties in the gas-correlated emission were studied using alternative model templates. In particular, the interstellar gas maps included in our benchmark model were obtained from a suite of hydrodynamic simulations of interstellar gas flow [65]. However, we also considered interpolated gas templates that reproduce those used in the construction of the official Fermi diffuse emission model [38]. Reference [15] 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 GC region. Note that this result has been independently confirmed with the non-Poissonian template fitting pipeline [46]. Some previous GCE analyses estimated the systematic uncertainties associated to the IC component by using the results of a GALPROP propagation parameter scan in Ref. [66]. However, that study was restricted to a selected set of CR injection and propagation scenarios that assumed 2D Galactocentric cylindrically symmetric geometry for the Galaxy. Although this assumption is physically sensible and has allowed to gain deep insights into the gamma-ray sky, it is expected to introduce a bias to GCE studies since the 2D IC models fail to incorporate the non-axisymmetric characteristics of the stellar distribution in the MW, such as the spiral arms and bar [57,67]. Indeed, the most recent release of the GAL-PROP code [57] has introduced more realistic 3D spatial models for the CR source and ISRF densities. These include sophisticated templates for the spiral arms, the bulge/bar complex and warped stellar/dust disk [67]. In the present study we have reproduced the results in Ref. [57] and included in our analysis one of their main 3D IC models named F98-SA50 (see Table 3 of Ref. [57]). The choice of this particular model had no impact in our results since we have divided the 3D IC map in four or six rings and allowed their normalization to float in the fits in order to account for uncertainties in the CR densities. To allow for a greater range of systematics, we also included a 2D IC map containing an additional central source population of e − (Model B in Ref. [37]). A summary of the foreground/background models considered in this study is shown in Tab. I.
In Refs. [15,18], a subset of us showed that the GCE spatial morphology was better explained by the stellar nuclear bulge [61] and galactic bulge [63] templates than by a spherically symmetric excess map given by a gNFW 2 profile, as would be consistent with annihilating DM [1] (e.g., Table I of [15]). The preference for the bulge model was typically ∼ 10σ or higher. Similar results have been quantitatively obtained also by Ref. [16,18]. Given the high significance of these results, we add the nuclear bulge and galactic bulge templates to our benchmark model for the GC. Note that these templates are detected in addition to the 3D IC templates that already contain the galactic bulge as a source of photons and CRs. This can be interpreted as gamma-ray sources distributed according to the galactic bulge. For example, in the MSP scenario, the prompt gamma-ray emission would still be required to be accounted for even while their secondary emission is modeled by the 3D IC maps. It is also worth noting that a nuclear bulge component has not yet been included in the GALPROP Galaxy model.
The analysis procedure used here is similar to that of Ref. [15,18]. We employed a bin-by-bin fitting method in which a separate maximum-likelihood was run at each energy bin. In each of these maximum-likelihood runs we simultaneously varied the normalization of all the GDE templates, the 120 brightest 4FGL point sources, and the point source population template containing the remaining 367 point sources. We used the pyLikelihood tool to vary a total of 146 parameters in the fits and ensure they converged.
Appendix B: Spatial maps for the GCE Detailed specifications of the templates for the GCE are given in Sec. II of Ref. [18]. Here we provide a brief description of the templates considered with an emphasis on those that are new in the present work.
We used two types of spatial models in our analysis of the morphology of the GCE: Stellar density and DM density (squared) maps. For the bulge stars we included the "boxy bulge" (Fig. 3-a) model -obtained in Ref.
[63] from a fit to diffuse infrared COBE/DIRBE data -as well as the "nuclear bulge" [61] (Fig. 3-b) which is a stellar density map of the inner 400 pc. For DM we used the density distributions given by gNFW (Fig. 3-c) or cored (Fig. 3-e) profiles, already described in the main text. However, in our analysis we also included a DM halo shape that departs from the commonly assumed spherical symmetry; we considered an oblate halo shape with its longer axis aligned with the Galactic disk. Indeed, collisionless N-body simulations predict ellipsoidal halos with density profile minor-to-major axis ratios ∼ 0.4−0.6 (e.g. [68, 69]). Moreover, hydrodynamical simulations have shown that baryonic dissipation can mitigate this halo shape contraction by making the DM halos more spherical (see e.g. [70] and references therein). In practice the halo shape contraction is implemented in our analysis by making a transformation in the Galactic distance (introduced in Eq. 2 and 3) of the form r → r , where r is given by x, y, and z are Galactocentric Cartesian coordinates and a, b, and c are the major, intermediate, and minor axis scale lengths. We have opted for assuming a minor-tomajor axis ratio c/a ∼ 0.7 and intermediate-to-minor b/a ∼ 1, which are the best values found in a recent study [36] (based on the results of the Eris simulations). The actual DM templates included in our maximum likelihood runs were constructed by performing a line-ofsight integral of the density squared profiles. Figure 3

-(d)
and -(f) show the NFW and cored profiles after implementation of the above halo shape contraction.

Appendix C: Evaluation of the Systematic Uncertainties in the DM limits
The systematic uncertainties were evaluated by repeating our DM limits procedure with variants of the background/foreground emission model. In particular, the log-likelihood scans for the DM source were performed with a bin-by-bin fitting method in which the different templates were fitted independently in small energy bins. This helps to mitigate the impact on the results of the assumed spectrum of the several templates. At each energy bin the differential DM flux was assumed to be described by a simple power-law of the form N 0 (E/2002.3 MeV) −2 . In our procedure we consecutively fixed the DM flux normalization N 0 to a given value in the range [10 −17 , 10 −7 ] MeV −1 cm −2 s −1 and maximized the likelihood to find the best-fit normalization parameters for all the remaining extended and pointlike sources included in the model. This was done with the UpperLimits tool within Fermitools. We started the scans with the benchmark model described in the main text, but also applied the procedure to variants of the foreground/background model. Specifically, we ran independent log-likelihood scans in which we replaced the benchmark 3D IC map (divided in four rings) by an alternative 3D IC template (divided in six rings). In addition, we considered a 2D IC model that contains an extra electrons-only [37] source population in the GC. The spatial distribution of the additional source of electrons used in the construction of this 2D IC model can be seen in Fig.  13 of Ref. [37]. The uncertainties introduced by some of the assumptions in the creation of the hydrodynamic gas and dust templates were investigated in the same manner. Since the amount of dust traced by the E(B − V ) extinction map is not accurate in regions of high extinction, we utilized dust map templates constructed with two different magnitude cuts; 5 mag (benchmark model) and 2 mag. To encompass a greater range of systematic uncertainties we also included the interpolated gas maps that reproduce the ones in Ref. [38].
The results of our scans for the benchmark and alternative background/foreground models are shown in Fig. 4. Regardless of the background/foreground model assumed, we find that a putative DM source starts to significantly worsen the fits for DM fluxes in the range ∼ 5 × 10 −11 ph cm −2 s −1 and 2 × 10 −8 ph cm −2 s −1 , depending on the energy bin. To have a better understanding of the constraining power of each of our energy bins, we have also displayed the 95% C.L flux upper limits in Fig. 5. These were computed by requiring a change in each profile log-likelihood of 2.71/2 from their maximum. We remind the reader that our background/foreground model includes templates for the spatial distribution of the bulge stars. As thoroughly discussed in previous studies [15,16,18], once the stellar bulge models are included to the fits, DM-like spatial models are strongly disfavored. Our current analysis leveraged on this fact to impose some of the strongest constraints on self-annihilating DM models.
We note that the Sun and the Moon contribute extended gamma-ray emission in our RoI, not accounting for this emission can bias the spectra of other sources included in our analysis. Templates describing gamma-rays originating from the Sun and the Moon need to be independently constructed to match the specific data cuts adopted in the analysis (photon event type, maximum zenith angle cut, energy and time range). However, constructing newer Sun and Moon templates is bounded by computational costs. As a compromise between computational requirements and photon statistics we used the same data cuts as in the 4FGL [19] for which there are appropriate Sun and Moon templates readily available. We note that this is another important factor justifying the amount of Fermi data included in this analysis.
The best-fit spectra of the benchmark background and foreground models are shown in Fig. 6. As can be seen, our bin-by-bin method produces stable and physically sensible spectra for the model components considered in this work. For display purposes, we have combined the spectra of different sources in groups. It can be observed that the GCE is replaced by the stellar bulge templates. And importantly, the bulge is found to be spectrally distinct to our Fermi bubbles map. Figure 7 shows the latitudinal and longitudinal flux profiles of the various components included in the fits in comparison with the Fermi-LAT data. There are noticeable differences in the shape of the "galactic bulge" component between the latitudinal and longitudinal flux profiles. This is due to the oblateness of both the boxy bulge and the nuclear bulge templates. It is of importance for this study that the background model components are spatially and spectrally different to the expected galactic bulge emission as this helps preventing possible degeneracies that could impact our log-likelihood scans.  [2.8, 11.8] and [11.8, 158.1] GeV, respectively. All the panels were smoothed with a Gaussian filter of radius 0.5 • . Although the spatial resolution of the LAT is higher than this for energies > ∼ 10 GeV, this choice is motivated by limitations in some of our background/foreground model templates. For example, the distribution of the atomic hydrogen column density was derived from the Leiden-Argentine-Bonn 21 cm galactic HI composite survey [71] which itself has a spatial resolution 0.5 • .
The panels in the third column of of Fig. 8 show the fractional residuals, (Data − Model)/Model, for the benchmark model (see Fig. 6 for the spectrum). It can be seen in the first three energy windows ([0.6, 1.1], [1.1, 2.8] and [2.8, 11.8] GeV) that the model mostly underpredicts the data at the < ∼ 10% level, with the exception of some more localized negative residuals that reach up to ∼ 20%. However, in the last energy window ([11.8, 158.1] GeV) the regions of under/over-prediction can reach to the ∼ 30% level. Interestingly, these fractional residual images (specially the last energy window) bear some resemblance to the low latitude Fermi bubbles counterpart (e.g., Fig.8 bottom-right of Ref. [58]). It should be noticed that the Fermi bubbles template used in this analysis is an inpainted version of the original Fermi bubbles template obtained in Ref. [58]. In that study, a Spectral Component Analysis [72] was applied to data in the [1,10] GeV energy range in order to reconstruct a morphological template with photons having the same spectrum as that of the Fermi bubbles in the high-latitude region. It is possible that if the same image reconstruction technique is applied to data that includes the [11.8, 158.1] GeV energy range, regions of under/over-prediction in our last energy window will be ameliorated. A more thorough investigation of the Fermi bubbles template in our last  . Profiles of the bin-by-bin log-likelihood function used to test for a putative DM source in the 40 • × 40 • region of the GC. Each profile shows the log-likelihood ratio between "background/foreground"+"DM source" model and the background/foreground-only-model. The bin-by-bin log-likelihood was calculated by scanning the flux normalization of the DM source within each energy bin in the range 10 −17 and 10 −7 ph cm −2 s −1 . When running the scan, the flux normalization of the background sources were varied in the fit while their spectral slope fixed to −2. The best-fit spectrum of the background/foreground model components included in the fit can be seen in Fig. 6. Within each energy bin, the line colors denote an alternative spatial template for the DM source (cored, ellipsoidal-cored, NFW and ellipsoidal-NFW, see Tab. I for descriptions) or an alternative gamma-ray background/foreground model. We changed the benchmark 3D IC map (divided in four rings) by an alternative 3D IC map divided in six rings, and 2D IC map containing a central source of e − [37]. We also varied the magnitude cut used in the construction of the gas maps and explored the results obtained with the Interpolated gas maps. Unless otherwise stated, we conservatively assume an NFW-squared (γ = 1.2) density profile for the DM map since this model has the largest log-likelihood.
energy window is beyond the scope of our current work and we leave this interesting possibility for future a analysis.
We note that when a DM-like template is included as a model for the positive residual, this is unable to account for all of the residual emission. In this sense our DM limits should be seen as conservative. Even though the residual emission does not appear spherically symmetric distributed, our fitting procedure allows sufficient freedom to the DM template to try account for most of the residual photons.
Our main concern in this section was to investigate the extent at which the computed DM constraints depend on the specific fore-/background model assumed. It was not   our aim to perform an exhaustive search for an alternative foreground model that matches the LAT data best in the GC region. Indeed, in Ref. [18] we have shown that GDE models that assume the hydrodynamical gas and the new 3D IC maps are better fits to the Fermi data. However, here we used the different variations in   the fore-/background model for the purpose of testing the impact they had in our limits and estimating their expected variance.

Appendix D: Dark Matter Injection and Recovery tests
Given that our upper limit procedure allows for all the sources to vary in the fits 3 , it is crucial to verify that our foreground/background model would not absorb a DM signal if one were present in the data. For this, we have artificially injected DM signals of different characteristics into the real data and consecutively applied our upper limits procedure to each augmented data set.
Our tests are similar to those carried out in Ref. [42,73]; we have simulated DM injections by taking a random Poisson draw of DM maps generated for a range of DM masses and annihilation cross sections. In particular, we considered self-annihilating DM in thebb channel, DM masses of 10, 25, 100, and 500 GeV, annihilation cross sections in the range [10 −27 , 3 × 10 −25 ] cm 3 /s, and two different DM spatial morphologies (gNFW γ = 1.2 and cored profiles). For a given realization, we obtained the 95% CL flux upper limits by requiring a change in the log-likelihood of 2.71/2 from the best-fitting point.
The results of our DM injection tests are presented in Fig. 9 (gNFW γ = 1.2) and Fig. 10 (cored profile). In each panel, the black line shows the DM signal that was injected into the data and the red arrows display the 95% CL flux upper limits recovered with our log-likelihood profile scan method. As can be seen from these figures, for a large majority of our realizations, the recovered binby-bin flux upper limits have the correct statistical coverage. There are a few cases in which our upper limits are below the injected DM signal; for most of those we nonetheless obtain that the upper limits weaken in a way that is consistent with the strength of the injected signal. The only exception to this pattern was observed in the highest energy bin (37.49 − 158.11 GeV) for the DM injections corresponding to the cored profile. In this case it was found that the flux upper limits did not have the correct statistical coverage for all our high DM mass injection trials. It is possible that this is due to a combination of complicating factors: First, the cored profile is much flatter than the gNFW profile. Second, in the highest energy bin the statics are low. Degeneracies between the injected DM signal and the GDE model components appear to be difficult to resolve under these conditions.